Development of a three dimension multi-physics code for molten salt fast reactor

Development of a three dimension multi-physics code for molten salt fast reactor 1 2 first-author 1 corresp Corresponding author, daizhimin@sinap.ac.cn Corresponding author, daizhimin@sinap.ac.cn daizhimin@sinap.ac.cn Center for Thorium Molten Salt Reactor System, Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai201800, China Center for Thorium Molten Salt Reactor System, Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai201800, China University of Chinese Academy of Sciences, Beijing100049, China University of Chinese Academy of Sciences, Beijing100049, China Molten Salt Reactor (MSR) was selected as one of the six innovative nuclear reactors by the Generation IV International Forum (GIF). The circulating-fuel in the can-type molten salt fast reactor makes the neutronics and thermo-hydraulics of the reactor strongly coupled and different from that of traditional solid-fuel reactors. In the present paper, a new coupling model is presented that physically describes the inherent relations between the neutron flux, the delayed neutron precursor, the heat transfer and the turbulent flow. Based on the model, integrating nuclear data processing, CAD modeling, structured and unstructured mesh technology, data analysis and visualization application, a three dimension steady state simulation code system (MSR3DS) for the can-type molten salt fast reactor is developed and validated. In order to demonstrate the ability of the code, the three dimension distributions of the velocity, the neutron flux, the delayed neutron precursor and the temperature were obtained for the simplified MOlten Salt Advanced Reactor Transmuter (MOSART) using this code. The results indicate that the MSR3DS code can provide a feasible description of multi-physical coupling phenomena in can-type molten salt fast reactor. Furthermore, the code can well predict the flow effect of fuel salt and the transport effect of the turbulent diffusion. Molten salt fast reactor Turbulent model Delayed neutron precursor Neutronics Thermo-hydraulics Turbulent diffusion I. INTRODUCTION 1 The first molten salt reactor (MSR) was developed in late 1940s as the aircraft propulsion at Oak Ridge National Laboratory (ORNL). The Aaircraft Reactor Experiment (ARE) [ 1 1 ] in 1954 operated successfully. In 1965, the Molten Salt Reactor Experiment (MSRE) [ 2 2 ] goes critical, and after six months of successful operation the enriched U-235 was removed and replaced by denatured U-233 as fissile fuel. In October 1966, the MSRE was the first reactor that reached criticality with U-233. A detailed 1000 MWe engineering conceptual design of a Molten Salt Breeder Reactor (MSBR) [ 3 3 ] was developed in 1970s. Even though the concept looked promising, the studies were stopped in 1976. MSR has the potential of meeting the goals of Generation IV reactors and high-level waste transmutation program, being better than solid fuel reactors. Consequently, many countries are interested in developing new concept MSRs, such as the FUJI series [ 4 4 ], Actinides Molten Salt TransmutER (AMSTER) [ 5 5 ], MOSART [ 6 6 ], Fast Spectrum Molten Salt Reactor (FS-MSR) [ 7 7 ], European Molten Salt Fast Reactor (MSFR) [ 8 8 ], and Small Mobile Molten Salt Reactor (SM-MSR) [ 9 9 ]. In 2010, Chinese Academy of Sciences restarted the thorium molten salt reactor project and established the Center for Thorium Molten Salt Reactor System (CTMSRS) in Shanghai Institute of Applied Physics, and CTMSRS finished the conceptual design of a Liquid Fuel Thorium Molten Salt Reactor (LF-TMSR). The Can-type Molten Salt Fast Reactor (CMSFR) can be employed to consume actinides from light water reactor (LWR) fuel or, alternatively, to extend fissile resource availability through U/Pu and Th/U breeding. CMSFRs are highly flexible and can be configured into modified open or full-recycle configuration. Since 2005 the GIF has focused on CMSFR from thermal spectrum molten salt reactor. Many works have been conducted to investigate the complex behavior of the CMSFR. Wang et al . [ 10 10 ] extended the SIMMER-III code for simulating the MOSART with extra thermo-hydraulic and neutronic modules at two-dimension (2D) axial-symmetric geometry. Nicolino et al . [ 11 11 ] developed a new approach to describe the strong co upling between neutronics and thermo-fluid dynamics with particular focus to the MOSART molten salt fast reactor at 2D axial-symmetric geometry. Zhang et al . [ 12 12 ] provided a theoretical model of flow coupling the heat transfer and neutronics models at 2D axial-symmetric geometry and calculated of the steady characteristics of a molten salt fast reactor without graphite moderator in the core. In all the studies, the transport effect of turbulent diffusion wasn’t taken into account for the delayed neutron precursor (DNP) concentration. In this paper, we present a new model to describe the inherent relations between the neutronics and the thermohydraulics. Based on this model, a three dimension steady state simulation code (MSR3DS) for the molten salt fast reactor is developed and validated. Using the multi-physical coupling code, the three dimension distributions of the velocity, neutron fluxes, DNPs and the temperature are obtained for a simplified MOSART. II. MATHEMATICAL MODEL 1 More recently, several fluoride molten salt fast reactors have been proposed as part of the Gen IV program. These include the European MOSART reactor using LiF/NaF/BeF 2 /(TRU)F 3 as a fuel salt, and an European MSFR concept using LiF/NaF/(TRU)F 3 or LiF/NaF/(U+Th)F 4 . In this paper, in order to establish the theoretical model and develop the simulation code, a core configuration shown in Fig. 1 Fig. 1 is adopted based on the core configuration of the MOSART. The main parameters of the core are list in Table 1 Table 1 . The mass proportions at equilibrium in finite critical core with 20 cm graphite reflector are given in Table 2 Table 2 . Table 3 Table 3 lists six group delayed neutron fractions and the precursor decay constants. 1 1 Parameters 1 1 Values 1 1 Thermal power(MW) 1 1 2400 1 1 Reflector thickness(m) 1 1 0.2 1 1 Inlet temperature(K) 1 1 873.15 1 1 Inlet velocity(m⋅s -1 ) 1 1 0.5 1 1 Resident time out of the core(s) 1 1 3.94 1 1 Compositions 1 1 Atomic number density/10 -24 cm -3 1 1 Li-6 1 1 3.856E-06 1 1 Li-7 1 1 4.301E-03 1 1 Na-23 1 1 1.774E-02 1 1 Be-9 1 1 8.030E-03 1 1 F-19 1 1 3.856E-02 1 1 Np-237 1 1 5.713E-06 1 1 Np-239 1 1 4.407E-11 1 1 U-234 1 1 3.993E-07 1 1 U-235 1 1 3.378E-07 1 1 U-236 1 1 2.522E-07 1 1 U-237 1 1 1.617E-09 1 1 U-238 1 1 2.172E-09 1 1 Pu-236 1 1 4.843E-12 1 1 Pu-238 1 1 2.088E-05 1 1 Pu-239 1 1 4.111E-05 1 1 Pu-240 1 1 7.585E-05 1 1 Pu-241 1 1 3.322E-05 1 1 Pu-242 1 1 4.954E-05 1 1 Pu-244 1 1 1.049E-11 1 1 Am-241 1 1 3.728E-06 1 1 Am-242m 1 1 1.266E-07 1 1 Am-243 1 1 1.789E-05 1 1 Cm-242 1 1 7.802E-07 1 1 Cm-243 1 1 1.359E-07 1 1 Cm-244 1 1 2.849E-05 1 1 Cm-245 1 1 1.286E-05 1 1 Cm-246 1 1 9.889E-06 1 1 Cm-247 1 1 3.473E-06 1 1 Cm-248 1 1 1.202E-06 1 1 Bk-249 1 1 2.336E-07 1 1 Cf-249 1 1 2.878E-07 1 1 Cf-250 1 1 3.215E-07 1 1 Cf-251 1 1 2.169E-07 1 1 Zr-93 1 1 1.121E-06 1 1 Nd-143 1 1 1.252E-06 1 1 Nd-145 1 1 8.960E-07 1 1 Nd-147 1 1 3.608E-08 1 1 Pm-147 1 1 3.011E-07 1 1 Sm-149 1 1 2.055E-07 1 1 Sm-150 1 1 2.402E-07 1 1 Sm-151 1 1 1.554E-07 1 1 Sm-152 1 1 1.736E-07 1 1 Eu-153 1 1 1.966E-07 1 1 Eu-154 1 1 9.447E-08 1 1 Eu-155 1 1 8.558E-08 1 1 Gd-157 1 1 2.378E-08 1 1 B-10 1 1 5.282E-07 1 1 Groups 1 1 λ i (s -1 ) 1 1 β i (10 -5 ) 1 1 C1 1 1 0.0124 1 1 22.3 1 1 C2 1 1 0.0305 1 1 145.7 1 1 C3 1 1 0.111 1 1 130.7 1 1 C4 1 1 0.301 1 1 262.8 1 1 C5 1 1 1.14 1 1 76.6 1 1 C6 1 1 3.01 1 1 28.0 A. Thermohydraulics model 2 The fuel salt flow in the molten salt fast reactor without moderator in the core like MOSART and MSFR is considered as turbulent flow. Their Reynolds numbers in the core region can be higher than 10 5 . Many models can be used to describe turbulent flow. The most accurate method is the direct numerical simulation (DNS), which computes the mean flow and all turbulent velocity fluctuations. But the DNS model is highly costly in terms of computing resources, so the method is not used for industrial flow computations. The large eddy simulation is an intermediate form of turbulence calculations which tracks the behavior of larger eddies. The effects on the resolved flow due to the smallest, unresolved eddies are included by means of a so-called sub-grid scale model in the large eddy method. So, the demands on computing resources in terms of storage and volume of calculations are large. The Reynolds-averaged Navier-Stokes (RANS) equations are focused on the mean flow and the effects of turbulence on mean flow properties. The computing resources required for reasonably accurate flow computations are modest, so this approach has been the mainstay of engineering flow calculations in engineering practice. In this study, the flow, heat transfer and turbulent characteristics in the core were obtained by solving the following three-dimensional, incompressible, steady state Navier-Stokes equations and turbulence model. The RANS equations with Boussiesq’s closure hypothesis and the standard k-ε turbulence model [ 13 13 ] were adopted. The continuity equation is: <math display="block" id="M1"><mrow><mo>∇</mo><mo>⋅</mo><mi>U</mi><mo>=</mo><mn>0</mn><mn>.</mn></mrow></math> The momentum conservation equation is: <math display="block" id="M2"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mi>U</mi><mo>⋅</mo><mo>∇</mo><mo stretchy="false">)</mo><mi>U</mi><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>−</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>p</mi><mtext>I</mtext><mo>+</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><mi>k</mi><mtext>I</mtext></mrow><mo>)</mo></mrow></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mo stretchy="false">(</mo><mi>μ</mi><mo>+</mo><msub><mi>μ</mi><mi>t</mi></msub><mo stretchy="false">)</mo><mrow><mo>(</mo><mrow><mo>∇</mo><mi>U</mi><mo>+</mo><msup><mrow><mo stretchy="false">(</mo><mo>∇</mo><mi>U</mi><mo stretchy="false">)</mo></mrow><mi>T</mi></msup></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The transport equations for k and ε in the standard k-ε turbulence model are: <math display="block" id="M3"><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><mi>U</mi><mo>⋅</mo><mo>∇</mo><mi>k</mi><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>k</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>k</mi><mo>+</mo><mi>G</mi><mo>−</mo><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><mi>ε</mi><mn>,</mn></mrow></math> and <math display="block" id="M4"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><mi>U</mi><mo>⋅</mo><mo>∇</mo><mi>ε</mi><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>ε</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>ε</mi><mo>+</mo><msub><mi>C</mi><mrow><mi>ε</mi><mn>1</mn></mrow></msub><mfrac><mi>ε</mi><mi>k</mi></mfrac><mi>G</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>–</mo><msub><mi>C</mi><mrow><mi>ε</mi><mn>2</mn></mrow></msub><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><mfrac><mrow><msup><mi>ε</mi><mn>2</mn></msup></mrow><mi>k</mi></mfrac><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> <math id="M5"><mrow><msub><mi>μ</mi><mi>t</mi></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="M6"><mrow><mi>G</mi><mo>:</mo><msub><mi>ν</mi><mi>t</mi></msub><msup><mi>S</mi><mn>2</mn></msup></mrow></math> <math id="M7"><mrow><mi>S</mi><mo>:</mo><mo>∇</mo><mi>U</mi><mo>+</mo><msup><mrow><mo stretchy="false">(</mo><mo>∇</mo><mi>U</mi><mo stretchy="false">)</mo></mrow><mi>T</mi></msup></mrow></math> <math id="M8"><mrow><msub><mi>C</mi><mrow><mi>ε</mi><mn>1</mn></mrow></msub><mo>=</mo><mn>1.44</mn></mrow></math> , <math id="M9"><mrow><msub><mi>C</mi><mrow><mi>ε</mi><mn>2</mn></mrow></msub><mo>=</mo><mn>1.92</mn></mrow></math> , <math id="M10"><mrow><msub><mi>C</mi><mi>μ</mi></msub><mo>=</mo><mn>0.09</mn></mrow></math> <math id="M11"><mrow><msub><mi>σ</mi><mi>k</mi></msub><mo>=</mo><mn>1.0</mn></mrow></math> , <math id="M12"><mrow><msub><mi>σ</mi><mi>ε</mi></msub><mo>=</mo><mn>1.3</mn></mrow></math> . Where, U represents the mean velocity vector, and k , γ , ρ fs , μ , μ t , and I are turbulent kinetic energy, turbulent dissipation rate, fuel salt density, dynamic viscosity, turbulent dynamic viscosity and Identity matrix, respectively. The equations contain five constants: C ε 1 , C ε 2 , C μ , σ k and σ ε . The standard k - ε turbulence model employs values for the constants that are arrived at by comprehensive data fitting for a wide range of turbulent flows. The energy conservation equation of the fuel salt expressed by temperature: <math display="block" id="M13"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><msub><mi>c</mi><mrow><mtext>p</mtext><mn>,</mn><mtext>fs</mtext></mrow></msub><mo>∇</mo><mo>⋅</mo><mo stretchy="false">(</mo><mi>U</mi><msub><mi>T</mi><mrow><mtext>fs</mtext></mrow></msub><mo stretchy="false">)</mo><mo>=</mo><mo>∇</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><msub><mi>λ</mi><mrow><mi>T</mi><mn>,</mn><mtext>fs</mtext></mrow></msub><mo>+</mo><msub><mi>c</mi><mrow><mtext>p</mtext><mn>,</mn><mtext>fs</mtext></mrow></msub><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><mi>P</mi><msub><mi>r</mi><mi>t</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><msub><mi>T</mi><mrow><mtext>fs</mtext></mrow></msub></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>s</mi><mrow><mtext>fs</mtext></mrow></msub><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The energy conservation equation of the reflector expressed by temperature: <math display="block" id="M14"><mrow><mo>∇</mo><mo>⋅</mo><msub><mi>λ</mi><mrow><mi>T</mi><mn>,</mn><mtext>ref</mtext></mrow></msub><mo>∇</mo><msub><mi>T</mi><mrow><mtext>ref</mtext></mrow></msub><mo>+</mo><msub><mi>s</mi><mrow><mtext>ref</mtext></mrow></msub><mo>=</mo><mn>0</mn><mn>.</mn></mrow></math> The inner heat sources of fuel salt ( s fs ) and the reflector ( s ref ) in Eqs. (7) and (8) are calculated using the neutron fission reactions. <math display="block" id="M15"><mrow><msub><mi>s</mi><mrow><mtext>fs</mtext></mrow></msub><mo>=</mo><mi>γ</mi><msub><mi>E</mi><mi>f</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mo stretchy="false">(</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mo>⋅</mo><msub><mi>Σ</mi><mrow><mi>f</mi><mn>.</mn><mi>g</mi></mrow></msub><mo stretchy="false">)</mo><mn>,</mn></mrow></math> <math display="block" id="M16"><mrow><msub><mi>s</mi><mrow><mtext>ref</mtext></mrow></msub><mo>=</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>γ</mi><mo stretchy="false">)</mo><msub><mi>E</mi><mtext>f</mtext></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mtext>G</mtext></munderover></mstyle><mo stretchy="false">(</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mo>⋅</mo><msub><mi>Σ</mi><mrow><mtext>f</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">)</mo><mn>,</mn></mrow></math> where T fs , T ref , λ T ,fs , λ T , ref , c p,fs and Pr t respectively represent the temperature of the fuel salt and the reflector, thermal conductivity of the fuel salt and the reflector, specific heat capacity of the fuel salt and turbulent Prandtl number; E f , ϕ g , ∑ f, g and γ are energy released from each fission reaction, neutron flux for group g , fission cross-section for group g and the fraction of power released into fuel molten salt, respectively. B. Neutronincs model 2 According to the basic conservation of the neutron number and DNP concentration in a control volume with the multi-group diffusion theory, the neutronics model of CMSFR can be derived. The diffusion equation of the neutron for group g [ 11 11 ]: <math display="block" id="M17"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><msub><mi>u</mi><mi>g</mi></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>ϕ</mi><mi>g</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac></mrow></mtd><mtd columnalign="left"><mrow><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>u</mi><mi>g</mi></msub></mrow></mfrac><mo>∇</mo><mo>⋅</mo><mo stretchy="false">[</mo><mi>U</mi><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>=</mo><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mrow><mtext>n</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>∇</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>r</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Λ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>χ</mi><mrow><mtext>p</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mtext>G</mtext></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>Σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mn>,</mn><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mtext>I</mtext></munderover></mstyle><msub><mi>χ</mi><mrow><mi>d</mi><mn>,</mn><mi>i</mi><mn>,</mn><mi>g</mi></mrow></msub><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><menclose notation="updiagonalstrike"><mo>=</mo></menclose><mi>g</mi></mrow><mtext>G</mtext></munderover></mstyle><msub><mi>Σ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The delayed neutron precursors are classified into six groups by half-life periods. The balance equation of the DNP concentration for group i : <math display="block" id="M18"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac></mrow></mtd><mtd columnalign="left"><mrow><mo>+</mo><mo>∇</mo><mo>⋅</mo><mo stretchy="false">[</mo><mi>U</mi><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><msub><mi>D</mi><mrow><mtext>c</mtext><mn>,</mn><mi>i</mi></mrow></msub><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><msub><mrow><mtext>Sc</mtext></mrow><mtext>t</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>β</mi><mi>i</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mtext>I</mtext></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>Σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mn>,</mn><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The second items in Eqs. (9) and (10) indicate the flow effect of the fuel molten salt on the neutron fluxes and DNPs. The third item in Eq. (10) include two parts, the <math id="M19"><mrow><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mrow><mtext>c</mtext><mn>,</mn><mi>i</mi></mrow></msub><mo>∇</mo><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></math> part represents the molecular diffusion of the DNP group i and the <math id="M20"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><msub><mrow><mtext>Sc</mtext></mrow><mtext>t</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></math> part represents turbulent diffusion of the DNP group i . In liquids and for turbulent flow, compared to turbulent diffusion, the molecular diffusion can be neglected. In above equations, u g , C i , D n, g , D c, i , β i , β , ν , <math id="M21"><mrow><msub><mo>∑</mo><mrow><mi>g</mi><mo>’</mo><mo>→</mo><mi>g</mi></mrow></msub></mrow></math> , χ p, g , χ p,i, g and Sc t respectively represent neutron mean velocity for group g , the DNP for group i , neutron diffusion coefficient for energy group g , DNP diffusion coefficient for group i , the fraction of delayed neutron for group i , the total fraction of delayed neutron, the average number of neutrons produced in energy group g , scatter cross section from energy group g ’ to group g , the fission spectrums of prompt neutron for group g , the fission spectrums of delayed neutron for energy group g and delayed neutron group i , and turbulent Schmidt number. In steady state conditions, the time-dependent item can be removed from Eqs. (9) and (10) and the effective multiplication factor k eff is introduced. Therefore, the 3D steady state model of neutronics for molten salt fast reactor can be obtained. <math display="block" id="M22"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><msub><mi>u</mi><mi>g</mi></msub></mrow></mfrac><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mi>U</mi><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mrow><mtext>n</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>∇</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>r</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mfrac><mrow><msub><mi>χ</mi><mrow><mtext>p</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mtext>G</mtext></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>Σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mn>,</mn><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mtext>I</mtext></munderover></mstyle><msub><mi>χ</mi><mrow><mtext>d</mtext><mn>,</mn><mi>i</mi><mn>,</mn><mi>g</mi></mrow></msub><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><menclose notation="updiagonalstrike"><mo>=</mo></menclose><mi>g</mi></mrow><mtext>G</mtext></munderover></mstyle><msub><mi>Σ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M23"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><mo stretchy="false">[</mo><mi>U</mi><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><msub><mi>D</mi><mrow><mtext>c</mtext><mn>,</mn><mi>i</mi></mrow></msub><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>ρ</mi><mi>F</mi></msub><msub><mrow><mtext>Sc</mtext></mrow><mtext>t</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mtext>I</mtext></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>Σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mn>,</mn><msup><mi>g</mi><mo>′</mo></msup></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> In addition, the estimate for the effective multiplication factor k eff , may be computed from the Eq. (13): <math display="block" id="M24"><mrow><msubsup><mi>k</mi><mrow><mtext>eff</mtext></mrow><mi>n</mi></msubsup><mo>=</mo><mfrac><mrow><msubsup><mi>k</mi><mrow><mtext>eff</mtext></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mstyle displaystyle="true"><mrow><mo>∭</mo></mrow></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mtext>G</mtext></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub><mo stretchy="false">)</mo></mrow><mi>g</mi></msub><msubsup><mi>ϕ</mi><mi>g</mi><mi>n</mi></msubsup><mtext>d</mtext><mi>V</mi></mrow><mrow><mstyle displaystyle="true"><mrow><mo>∭</mo></mrow></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mtext>G</mtext></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub><mo stretchy="false">)</mo></mrow><mi>g</mi></msub><msubsup><mi>ϕ</mi><mi>g</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mtext>d</mtext><mi>V</mi></mrow></mfrac><mn>,</mn></mrow></math> where, the superscript n and n -1 represent the number of iterations. In this work, the steady neutronics model consisted of two-group (G =2) neutron diffusion equations for fast and thermal neutron fluxes in the fuel molten salt and reflector, transport equations for six-group (I =6) DNPs in the fuel salt. The neutron diffusion equations in the fuel salt: <math display="block" id="M25"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><msub><mi>u</mi><mn>1</mn></msub></mrow></mfrac><mo>∇</mo><mo>⋅</mo><mo stretchy="false">[</mo><mi>U</mi><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mrow><mtext>n</mtext><mn>,1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>∇</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>r</mtext><mn>,1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mfrac><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mn>2</mn></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>Σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mn>,</mn><msup><mi>g</mi><mo>′</mo></msup><mn>,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mn>,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mtext>I</mtext></munderover></mstyle><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M26"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><msub><mi>u</mi><mn>2</mn></msub></mrow></mfrac><mo>∇</mo><mo>⋅</mo><mo stretchy="false">[</mo><mi>U</mi><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mrow><mtext>n</mtext><mn>,2,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>∇</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>r</mtext><mn>,1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>Σ</mi><mrow><mn>1</mn><mo>→</mo><mn>2,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The neutron diffusion equations in the reflector: <math display="block" id="M27"><mrow><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mrow><mtext>n</mtext><mn>,</mn><mn>1,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>∇</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>r</mtext><mn>,</mn><mn>1,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>=</mo><mn>0</mn><mn>,</mn></mrow></math> <math display="block" id="M28"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mrow><mtext>n</mtext><mn>,</mn><mn>2</mn><mn>,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo>▽</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>r</mtext><mn>,</mn><mn>2</mn><mn>,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>Σ</mi><mrow><mn>1</mn><mo>→</mo><mn>2,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><mn>1,</mn><mtext>ref</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>=</mo><mn>0</mn><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The transport equations for six-group DNPs in the fuel salt: <math display="block" id="M29"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><mo stretchy="false">[</mo><mi>U</mi><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><msub><mi>D</mi><mrow><mtext>c</mtext><mn>,</mn><mi>i</mi></mrow></msub><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><msub><mrow><mtext>Sc</mtext></mrow><mtext>t</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mn>2</mn></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>Σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mn>,</mn><msup><mi>g</mi><mo>′</mo></msup><mn>,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mn>,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mover><mi>r</mi><mi>Δ</mi></mover><mo stretchy="false">)</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The effective multiplication factor k eff : <math display="block" id="M30"><mrow><msubsup><mi>k</mi><mrow><mtext>eff</mtext></mrow><mi>n</mi></msubsup><mo>=</mo><mfrac><mrow><msubsup><mi>k</mi><mrow><mtext>eff</mtext></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mstyle displaystyle="true"><mrow><mo>∭</mo></mrow></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mn>2</mn></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub><mo stretchy="false">)</mo></mrow><mrow><mi>g</mi><mn>,</mn><mtext>fs</mtext></mrow></msub><msubsup><mi>ϕ</mi><mrow><mi>g</mi><mn>,</mn><mtext>fs</mtext></mrow><mi>n</mi></msubsup><mtext>d</mtext><mi>V</mi></mrow><mrow><mstyle displaystyle="true"><mrow><mo>∭</mo></mrow></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mn>2</mn></munderover></mstyle><msub><mrow><mo stretchy="false">(</mo><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub><mo stretchy="false">)</mo></mrow><mrow><mi>g</mi><mn>,</mn><mtext>fs</mtext></mrow></msub><msubsup><mi>ϕ</mi><mrow><mi>g</mi><mn>,</mn><mtext>fs</mtext></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msubsup><mtext>d</mtext><mi>V</mi></mrow></mfrac><mn>,</mn></mrow></math> where, the subscript fs and ref respectively represent the fuel salt and the reflector; ∑ r,1,fs , ∑ r,2,fs , ∑ r,1,ref and ∑ r,2,ref represent the removal cross section. C. Boundary conditions 2 1. Thermohydraulics 3 (1) Inlet boundary: at the inlet the velocity, the turbulent kinetic energy, the turbulent dissipation rate and the temperature are given: <math display="block" id="M31"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><msub><mi>u</mi><mi>z</mi></msub><mo>=</mo><msub><mi>u</mi><mrow><mtext>inlet</mtext></mrow></msub><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mi>k</mi><mo>=</mo><mfrac><mn>3</mn><mn>2</mn></mfrac><msup><mrow><mo stretchy="false">(</mo><mi>I</mi><msub><mi>u</mi><mi>z</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mi>ε</mi><mo>=</mo><msubsup><mi>C</mi><mi>μ</mi><mrow><mn>0.75</mn></mrow></msubsup><mfrac><mrow><msup><mi>k</mi><mrow><mn>1.5</mn></mrow></msup></mrow><mi>l</mi></mfrac><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> where, u inlet is given as 0.5 m/s; the turbulence intensity I is 10%; C μ is a constant assigned 0.09; l is length scale. T inlet is imposed at 873.15 K. (2) Outlet boundary: at the outlet a free outflow is assumed. A Neumann boundary condition is used for the velocity, the turbulent kinetic energy, the turbulent dissipation and the temperature. (3) Symmetry planes: the symmetry boundary is set for the velocity, the turbulent kinetic energy, the turbulent dissipation and the temperature. (4) Wall boundary: the boundary condition at the inner wall is treated by wall function method for the velocity, the turbulent kinetic energy, the turbulent dissipation and the temperature. In the outer wall the temperature is set as a constant (680 K). 2. Neutronics 3 (1) Inlet boundary: the fast and thermal neutron fluxes are imposed to the vacuum boundaries. Due to the DNPs’s decay characteristics and resident time in external loop, the DNPs may return to the inlet. The DNPs at the inlet is: <math display="block" id="M32"><mrow><msub><mi>C</mi><mrow><mi>i</mi><mn>,</mn><mtext>inlet</mtext></mrow></msub><mo>=</mo><msub><mi>C</mi><mrow><mi>i</mi><mn>,</mn><mtext>outlet</mtext></mrow></msub><mo>⋅</mo><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><mi>τ</mi></mrow></msup><mn>,</mn></mrow></math> where C i ,outlet is the DNP group i in the outlet and τ is the resident time out of the core for fuel salt. (2) Outlet boundary: the vacuum boundary is for the neutron fluxes at the outlet, and the DNPs adopt the zero gradient boundary condition. (3) Symmetry plane: the neutron fluxes and the DNPs are all set symmetry boundary condition. (4) Wall boundary: in the inner wall the neutron fluxes use the coupling boundary condition and the DNPs are set to vacuum condition. In the outer wall the vacuum condition is for the neutron fluxes. D. Thermophysical properties and group constants 2 Generally, the thermophysical properties in the thermohydraulics and the group constant in the neutronics are both dependent on the temperature. Ignatiev et al . [ 6 6 ] gave the LiF/NaF/BeF 2 solvent system properties including the density, the thermal conductivity, the viscosity and the heat capacity. <math display="block" id="M33"><mrow><msub><mi>ρ</mi><mrow><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mtext>kg</mtext><mo>⋅</mo><msup><mtext>m</mtext><mrow><mo>−</mo><mn>3</mn></mrow></msup><mo stretchy="false">)</mo><mo>=</mo><mn>2518</mn><mo>−</mo><mn>0.406</mn><mi>T</mi><mn>,</mn></mrow></math> <math display="block" id="M34"><mrow><msub><mi>λ</mi><mrow><mi>T</mi><mn>,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mtext>W</mtext><mo>⋅</mo><msup><mtext>m</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>⋅</mo><msup><mtext>K</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo stretchy="false">)</mo><mo>=</mo><mn>0.0429</mn><mo>+</mo><mn>0.0009</mn><mi>T</mi><mn>,</mn></mrow></math> <math display="block" id="M35"><mrow><mi>μ</mi><mo stretchy="false">(</mo><mtext>Pa</mtext><mo>⋅</mo><mtext>s</mtext><mo stretchy="false">)</mo><mo>=</mo><mn>0.001</mn><msup><mi>e</mi><mrow><mo>−</mo><mn>0.9942</mn><mo>+</mo><mn>1603.2</mn><mo>/</mo><mi>T</mi></mrow></msup><mn>,</mn></mrow></math> <math display="block" id="M36"><mrow><msub><mi>C</mi><mrow><mtext>p</mtext><mn>,</mn><mtext>fs</mtext></mrow></msub><mo stretchy="false">(</mo><mtext>J</mtext><mo>⋅</mo><msup><mrow><mtext>kg</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>⋅</mo><msup><mtext>K</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo stretchy="false">)</mo><mo>=</mo><mn>2090.</mn></mrow></math> The thermophysical properties of the graphite reflector are considered as constants ( Table 4 Table 4 ). 1 1 Properties 1 1 Graphite 1 1 ρ ref (kg·m -3 ) 1 1 1840 1 1 λ T ,ref (W·m -1 ·K -1 ) 1 1 31.2 1 1 C p,ref (J·kg -1 ·K -1 ) 1 1 1760 The MSR3DS code system uses DRAGON4 [ 14 14 ] code based on the XMAS-172 format libraries produced from the ENDF-VII.1 nuclear data to generate the two group macroscopic cross-sections, the diffusion coefficients and the neutron group velocities under different temperatures. The relations between group constants and the temperature are fitted by mean of the 5 th order polynomial curve: <math display="block" id="M37"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>Y</mi><mi>c</mi></msub><mo stretchy="false">(</mo><mi>T</mi><mo stretchy="false">)</mo><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><msub><mi>A</mi><mn>0</mn></msub><mo>+</mo><msub><mi>A</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>T</mi><mo>−</mo><msub><mi>T</mi><mrow><mtext>avg</mtext></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><msub><mi>A</mi><mn>2</mn></msub><msup><mrow><mo stretchy="false">(</mo><mi>T</mi><mo>−</mo><msub><mi>T</mi><mrow><mtext>avg</mtext></mrow></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>A</mi><mn>3</mn></msub><msup><mrow><mo stretchy="false">(</mo><mi>T</mi><mo>−</mo><msub><mi>T</mi><mrow><mtext>avg</mtext></mrow></msub><mo stretchy="false">)</mo></mrow><mn>3</mn></msup><mo>+</mo><msub><mi>A</mi><mn>4</mn></msub><msup><mrow><mo stretchy="false">(</mo><mi>T</mi><mo>−</mo><msub><mi>T</mi><mrow><mtext>avg</mtext></mrow></msub><mo stretchy="false">)</mo></mrow><mn>4</mn></msup></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>A</mi><mn>5</mn></msub><msup><mrow><mo stretchy="false">(</mo><mi>T</mi><mo>−</mo><msub><mi>T</mi><mrow><mtext>avg</mtext></mrow></msub><mo stretchy="false">)</mo></mrow><mn>5</mn></msup></mrow></mtd></mtr></mtable></mrow></math> where Y c denotes the group constants and T avg is reference temperature. A 0 – A 5 are fitting constants. E. Numerical method 2 The finite volume method (FVM) is widely employed for solution of computational fluid dynamics (CFD) problems in engineering. The solution domain is subdivided into a finite number of small control volumes and the conservation equations are applied to each control volume. And FVM can handle complex geometries. In this paper, the FVM is used for spatial discretization of all equations. We have discussed methods of discretizing the governing equations. This process results in a system of linear algebraic equations which needs to be solved. For the discretised equations of the steady neutron diffusion and the heat transfer Gauss-Seidel iterative method is employed. The Preconditioned Bi-Conjugate Gradient Method (PBiCG) is used to solve the discretized equations of the DNPs, the momentum, the turbulent kinetic energy and the turbulent dissipation. And the Preconditioned Conjugate Gradient Method (PCG) is applied in the discretized equations of the pressure correction [ 15 15 ]. The SIMPLE [ 13 13 ] algorithm gives a method of calculating pressure and velocities. The acronym SIMPLE stands for Semi-Implicit Method for Pressure-Linked Equations. It was originally put forward by Patankar and Spalding and is essentially a guess-and-correct procedure for calculation of pressure. It is an iterative method, and when other scalars are coupled to the momentum equations the calculation shall be done sequentially. The sequence of operations in the SIMPLE algorithm is given in Fig. 2 Fig. 2 . Fig. 3 Fig. 3 is the program flow diagram of the coupled solver in the MSR3DS code. F. The description of the MSR3DS code 2 The entire program diagram of the MSR3DS code is shown in Fig. 4 Fig. 4 . The MSR3DS code system consists of the pre-processing, coupled solver and post-processing module. The pre-processing module contains the CAD 3D modeling, meshing and generation of the group constants. The CAD 3D modeling and meshing tool use respectively the SolidWorks and ANSYS ICEM CFD software. The DRAGON4 generates the group constants under different temperature for the MSR3DS code. The code DRAGON4 is open source and contains a multi-group iterator conceived to control a number of different algorithms for the solution of the neutron transport equation. The SYBIL option solves the integral transport equation using the collision probability method for simple one-dimensional (1D) geometries (either plane, cylindrical or spherical) and the interface current method for 2D Cartesian or hexagonal assemblies. The EXCELL option solves the integral transport equation using the collision probability method for general 2D geometries and for 3D assemblies. The MCCG option solves the integro-differential transport equation using the long characteristics method for general 2D and 3D geometries. The ParaView [ 16 16 ] code is as the pos t-processing module in the MSR3DS code. The ParaView is an open-source, multi-platform data analysis and visualization application. It can be used quickly to build visualizations and to analyze data using qualitative and quantitative techniques. The data exploration can be done interactively in 3D or programmatically using ParaView’s batch processing capabilities. Furthermore, it is easily integrated into the MSR3DS code. III. RESULTS AND DISCUSSION 1 A. Code validation 2 The 3D TWIGL Seed/Blanket problem [ 17 17 ] is adopted to benchmark the neutronics calculation in the MSR3DS code. The effective multiplication factors calculated by CITATION [ 18 18 ] and MSR3DS are listed in Table 5 Table 5 , and the relative error is only 0.00147%. The fluxes normalized to 1 MW th are shown in Figs. 5 Figs. 5 and 6 6 . The results verify validity of the model for neutronics presented in this study. 1 1 Benchmark 2 1 k eff 1 1 Relative error/% 1 1 3D TWIGL 1 1 CITATION 1 1 MSR3DS 1 1 1 1 Seed/Blanket 1 1 0.885905 1 1 0.885892 1 1 0.00147 In order to evaluate the flow and heat transfer calculation in the MSR3DS code, a simple pipe case simplified from MOSART is used ( Fig. 7 Fig. 7 ). The thermophysical properties of the molten salt and the graphite at 873:15K are applied in the process of validation. Similarly, we assume that the velocity in the inlet is 0:5 m/s in the validation of the flow calculation of the MSR3DS. The results are compared in Figs. 8 Figs. 8 and 9 9 . From the figures, the results calculated by the MSR3DS agree well with those calculated by the Fluent. Therefore, the flow calculation in the MSR3DS is valid. For validation of the heat transfer calculation in the MSR3DS, the temperature in the inlet is assumed as 873:15 K, the temperatures in the outlet and other outer boundaries are 300 K. Figs. 10 Figs. 10 and 11 11 show that the temperatures obtained by the MSR3DS are in accord with those obtained by the Fluent, indicating that the heat transfer calculation in the MSR3DS is acceptable in an engineering context. B. Distributions of calculated physical fields in the core 2 Using the multiphysical coupling code, we calculted under different conditions the 3D distributions of velocity, the turbulent kinematic viscosity, fast and thermal neutron fluxes, DNPs and the temperature in the core. Figures 12 Figures 12 - 13 13 show the 3D distributions of the velocity and the turbulent kinematic viscosity. Fig. 13 Fig. 13 , the turbulent kinematic viscosity in the outlet is far greater than that of other locations. Figures 14 Figures 14 - 15 15 show the 3D distributions of the fast and thermal neutron flux without flow (a), with convective term (b), and with convective and turbulent diffusion term at Sc t = 0.7 (c) in the balance equations for six-group DNPs. The turbulent Schmidt number is set to 0.7 as the default value in the ANSYS Fluent. The temperature in the core is set to 900K under no flow condition. Figs. 14 Figs. 14 - 15(a) and (b) 15(a) and (b) show that the fuel salt flow has little effect on the distribution of the fast and thermal neutron fluxes. And the turbulent diffusion term hardly affect the distribution of the fast and thermal neutron fluxes comparing Figs. 14 Figs. 14 - 15(b) 15(b) with Figs. 14 Figs. 14 - 15(c) 15(c) . Figures 16 Figures 16 , 17 17 , 18 18 , 19 19 , 20 20 , 21 21 display the 3D distributions of DNPs without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc t = 0.7 (c) in the balance equations for six-group DNPs. The convective term affects the distribution of the DNPs significantly as shown in in Figs. 16 Figs. 16 , 17 17 , 18 18 , 19 19 , 20 20 , 21(b) 21(b) , and the smaller the delay constant, the greater the influence of the flow. After considering the turbulent diffusion term, due to the transport effect of the turbulent diffusion term, the concentrations of the DNPs decrease and the distribution of the DNPs in Figs. 16 Figs. 16 , 17 17 , 18 18 , 19 19 , 20 20 , 21(c) 21(c) are also changed observably, likewise, the smaller the delay constant, the greater the influence of the turbulent diffusion. In order to calculate the temperature distribution in the core, in this paper, the turbulent Prandtl number is imposed at 0.85, which is the same as the default value in the ANSYS Fluent. So, the distributions of calculated temperature with convective term (a), with convective and turbulent diffusion term at Sc t = 0.7 (b) are shown in Fig. 22 Fig. 22 . The results in the figure indicate that the distribution of the temperature in the core remains unchanged, because the DNPs affect the neutron flux slightly under the steady condition. In addition, because the transport effect of the turbulent diffusion term of the energy conservation equation at the outlet becomes more and more strong with the distribution of the turbulent kinematic viscosity, the fuel temperature at the outlet of the core decreases slightly in both cases. IV. CONCLUSION 1 In order to investigate the complex behavior of the core in the CMSFR, in this present research, a new multi-physical coupling model including the turbulent diffusion is presented. The model physically describes the mutually dependence in the neutron flux, the delayed neutron precursor (DNP), the heat transfer and the turbulent flow. The neutronics model consists of two group neutron diffusion equations for fast and thermal neutron fluxes considering the flow effect of fuel salt, and balance equations for six-group DNPs considering the flow effect of fuel salt and the transport effect of the turbulent diffusion. In the thermohydrualics, the RANS equations with Boussiesq’s closure hypothesis and the standard k - ε turbulence model were adopted. Based on the model, integrating open source DRAGON4 and ParaView code, the CAD modeling, structured and unstructured mesh technology, a 3D multi-physical coupling steady state code system is developed and validated. With the purpose of demonstrating the ability of the code, the 3D distributions of the velocity, the temperature, the neutron flux and the DNPs are obtained for simplified MOSART under steady-state condition. The main results are as follows: (1) The fuel salt flow has little effect on the distribution of the fast and thermal neutron fluxes. And the turbulent diffusion term hardly affect the distribution of the fast and thermal neutron fluxes. (2) The convective term affects the distribution of the DNPs significantly, and the smaller the delay constant, the greater the influence of the flow. The turbulent diffusion reduces the concentrations of the DNPs and observably changes the distribution of the DNPs, likewise, the smaller the delay constant, the greater the influence of the turbulent diffusion. (3) The turbulent diffusion item in these balance equations for six-group DNPs doesn’t change the distribution of the temperature in the core. But the transport effect of the turbulent diffusion term of the energy conservation equation has strong effect on the distribution of the temperature. Hence, the MSR3DS code system can be applied to simulate main physical fields and describe multi-physical coupling phenomena in the core of molten salt fast reactor and can well reflect the flow effect of the convective term and the transport effect of the turbulent diffusion term, which is peculiar to cantype molten salt reactor. References Bettis E S , Cottrell W B , Mann E R , et al . Nucl Sci Eng , 1957 , 2 : 804 - 825 . Nucl Sci Eng Robertson R C . MSRE design and operations report, part I, description of reactor design. Technical report: ORNL–TM–0728 , 1965 . MSRE design and operations report, part I, description of reactor design. Technical report: ORNL–TM–0728 Rosenthal MW , Haubenreich P N , Briggs R B . The development status of molten–salt breeder reactors. Technical report: ORNL–4812 , 1972 . The development status of molten–salt breeder reactors. Technical report: ORNL–4812 Furukawa K , Lecocq A , Kato Y , et al . J. Nucl Sci Technol , 1990 , 27 : 1157 - 1178 . J. Nucl Sci Technol Vergnes J and Lecarpentier D . Nucl Sci Eng , 2002 , 216 : 43 - 67 . Nucl Sci Eng Ignatiev V , Feynberg O , Gnidoi I , et al . 7548: Progress in development of Li, Be, Na/F molten salt actinide recycler & transmuter concept . Proceedings of ICAPP 2007 , Nice, France , May 13–18, 2007 . 7548: Progress in development of Li, Be, Na/F molten salt actinide recycler & transmuter concept Holcomb D E , Flanagan G F , Patton B W , et al . Fast spectrum molten salt reactor options. Technical report ORNL/TM–2011/105 , 2011 . Fast spectrum molten salt reactor options. Technical report ORNL/TM–2011/105 Merle–Lucotte E , Heuer D , Allibert M , et al . Introduction to the Physics of Molten Salt Reactors, Materials Issues for Generation IV Systems, NATO Science for Peace and Security Series – B , Editions Springer , 2008 , 501 - 521 . Introduction to the Physics of Molten Salt Reactors, Materials Issues for Generation IV Systems, NATO Science for Peace and Security Series – B Casino W A , Sorensen K F , Whitener C A . A small mobile molten salt reactor (sm–msr) for underdeveloped countries and remote locations . The 2007 ANS Student Design Contest, 2007 . A small mobile molten salt reactor (sm–msr) for underdeveloped countries and remote locations Wang S , Rineiski A , Maschek W . Nucl Sci Eng , 2006 , 236 : 1580 - 1588 . Nucl Sci Eng Nicolino C , Lapenta G , Dulla S , et al . Ann Nucl Energy , 2008 , 35 : 314 - 322 . Ann Nucl Energy Zhang D L , Qiu S Z , Su G H , et al . Ann Nucl Energy , 2009 , 36 : 590 - 603 . Ann Nucl Energy Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics: the finite volume method . England : Longman Group Ltd , 1995 . An introduction to computational fluid dynamics: the finite volume method Marleau G , Hebert A , Roy R . A user guide for dragon version4 . Qu´ebec : Ecole Polytechnique de Montreal , 2013 . A user guide for dragon version4 Saad Y . Iterative methods for sparse linear systems , 2 nd edition. Pennsylvania : Society for Industrial and Applied Mathmatics , 2003 . Iterative methods for sparse linear systems ParaView Guide . New York : Kitware, Inc. , 2013 . ParaView Guide Taylor J B . The development of a three–dimensional nuclear reactor kinetics methodology based on the method of characteristics, PhD. Thesis ( Pennsylvania State University , 2007 ). The development of a three–dimensional nuclear reactor kinetics methodology based on the method of characteristics, PhD. Thesis Fowler T B , Vondy D R , Cunningham GW . Nuclear reactor core analysis code: CITATION, ORNL–TM–2496 . Tennessee : Oak Ridge National Laboratory , 1971 . Nuclear reactor core analysis code: CITATION, ORNL–TM–2496 10.13538/j.1001-8042/nst.25.010601 Supported by the "Strategic Priority Research Program" of the Chinese Academy of Science (No.XD02001004) 2013-10-10 2013-11-25 2014-02-20 2014-02-20 2021-04-11 2014 Nuclear Science and Techniques 1 25