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

NUCLEAR ENERGY SCIENCE AND ENGINEERING

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

CHENG Mao-Song，

DAI Zhi-Min

Nuclear Science and Techniques

Vol.25, No.1

Article number 010601

Published in print 20 Feb 2014

Available online 20 Feb 2014

DOI：10.13538/j.1001-8042/nst.25.010601

36201

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 reactorTurbulent modelDelayed neutron precursorNeutronicsThermo-hydraulicsTurbulent diffusion

I. INTRODUCTION

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] in 1954 operated successfully. In 1965, the Molten Salt Reactor Experiment (MSRE) [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] 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], Actinides Molten Salt TransmutER (AMSTER) [5], MOSART [6], Fast Spectrum Molten Salt Reactor (FS-MSR) [7], European Molten Salt Fast Reactor (MSFR) [8], and Small Mobile Molten Salt Reactor (SM-MSR) [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] 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] developed a new approach to describe the strong coupling 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] 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

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₂/(TRU)F₃ as a fuel salt, and an European MSFR concept using LiF/NaF/(TRU)F₃ or LiF/NaF/(U+Th)F₄. In this paper, in order to establish the theoretical model and develop the simulation code, a core configuration shown in Fig. 1 is adopted based on the core configuration of the MOSART. The main parameters of the core are list in Table 1. The mass proportions at equilibrium in finite critical core with 20 cm graphite reflector are given in Table 2. Table 3 lists six group delayed neutron fractions and the precursor decay constants.

Main parameters of the core [6]

Parameters	Values
Thermal power(MW)	2400
Reflector thickness(m)	0.2
Inlet temperature(K)	873.15
Inlet velocity(m⋅s^-1)	0.5
Resident time out of the core(s)	3.94

Mass proportion at equilibrium in critical core [6]

Compositions	Atomic number density/10^-24cm^-3
Li-6	3.856E-06
Li-7	4.301E-03
Na-23	1.774E-02
Be-9	8.030E-03
F-19	3.856E-02
Np-237	5.713E-06
Np-239	4.407E-11
U-234	3.993E-07
U-235	3.378E-07
U-236	2.522E-07
U-237	1.617E-09
U-238	2.172E-09
Pu-236	4.843E-12
Pu-238	2.088E-05
Pu-239	4.111E-05
Pu-240	7.585E-05
Pu-241	3.322E-05
Pu-242	4.954E-05
Pu-244	1.049E-11
Am-241	3.728E-06
Am-242m	1.266E-07
Am-243	1.789E-05
Cm-242	7.802E-07
Cm-243	1.359E-07
Cm-244	2.849E-05
Cm-245	1.286E-05
Cm-246	9.889E-06
Cm-247	3.473E-06
Cm-248	1.202E-06
Bk-249	2.336E-07
Cf-249	2.878E-07
Cf-250	3.215E-07
Cf-251	2.169E-07
Zr-93	1.121E-06
Nd-143	1.252E-06
Nd-145	8.960E-07
Nd-147	3.608E-08
Pm-147	3.011E-07
Sm-149	2.055E-07
Sm-150	2.402E-07
Sm-151	1.554E-07
Sm-152	1.736E-07
Eu-153	1.966E-07
Eu-154	9.447E-08
Eu-155	8.558E-08
Gd-157	2.378E-08
B-10	5.282E-07

Delayed neutron fractions and the precursor decay constants [12]

Groups	λ_i (s^-1)	βi(10^-5)
C1	0.0124	22.3
C2	0.0305	145.7
C3	0.111	130.7
C4	0.301	262.8
C5	1.14	76.6
C6	3.01	28.0

Fig. 1.

(Color online) Geometry of the reactor core (Unit/mm).

A. Thermohydraulics model

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⁵. 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] were adopted.

The continuity equation is:

\nabla \cdot U = 0 .

(1)

The momentum conservation equation is:

\begin{array}{l} ρ_{fs} (U \cdot \nabla) U = & - \nabla \cdot (p I + \frac{2}{3} ρ_{fs} k I) \\ + \nabla \cdot [(μ + μ_{t}) (\nabla U + {(\nabla U)}^{T})] . \end{array}

(2)

The transport equations for k and ε in the standard k-ε turbulence model are:

ρ_{fs} U \cdot \nabla k = \nabla \cdot (μ + \frac{μ_{t}}{σ_{k}}) \nabla k + G - ρ_{fs} ε,

(3)

and

\begin{array}{l} ρ_{fs} U \cdot \nabla ε = & \nabla \cdot (μ + \frac{μ_{t}}{σ_{ε}}) \nabla ε + C_{ε 1} \frac{ε}{k} G \\ - C_{ε 2} ρ_{fs} \frac{ε^{2}}{k}, \end{array}

(4)

$μ_{t} : C_{μ} \frac{k^{2}}{ε}$ , $G : ν_{t} S^{2}$ $S : \nabla U + {(\nabla U)}^{T}$

$C_{ε 1} = 1.44$ , $C_{ε 2} = 1.92$ , $C_{μ} = 0.09$

$σ_{k} = 1.0$ , $σ_{ε} = 1.3$ .

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:

\begin{array}{l} ρ_{fs} c_{p, fs} \nabla \cdot (U T_{fs}) = \nabla [(λ_{T, fs} + c_{p, fs} \frac{μ_{t}}{P r_{t}}) \nabla T_{fs}] + s_{fs} . \end{array}

(5)

The energy conservation equation of the reflector expressed by temperature:

\nabla \cdot λ_{T, ref} \nabla T_{ref} + s_{ref} = 0 .

(6)

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.

s_{fs} = γ E_{f} \sum_{g = 1}^{G} (ϕ_{g} \cdot Σ_{f . g}),

(7)

s_{ref} = (1 - γ) E_{f} \sum_{g = 1}^{G} (ϕ_{g} \cdot Σ_{f, g}),

(8)

where T_fs, T_ref, λ_T,fs, λ_{T, ref}, c_p,fs and Prt 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

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]:

\begin{array}{l} \frac{1}{u_{g}} \frac{\partial ϕ_{g}}{\partial t} & + \frac{1}{u_{g}} \nabla \cdot [U ϕ_{g} (\overset{Δ}{r}, t)] \\ = \nabla \cdot D_{n, g} (\overset{Δ}{r}) \nabla ϕ_{g} (\overset{Δ}{r}, t) - Σ_{r, g} (\overset{Λ}{r}) ϕ_{g} (\overset{Δ}{r}, t) \\ + χ_{p, g} (1 - β) \sum_{g^{'} = 1}^{G} {(ν Σ)}_{f, g^{'}} (\overset{Δ}{r}) ϕ_{g^{'}} (\overset{Δ}{r}, t) \\ + \sum_{i = 1}^{I} χ_{d, i, g} λ_{i} C_{i} (\overset{Δ}{r}, t) + \sum_{g^{'} = g}^{G} Σ_{g^{'} \to g} (\overset{Δ}{r}) ϕ_{g^{'}} (\overset{Δ}{r}, t) . \end{array}

(9)

The delayed neutron precursors are classified into six groups by half-life periods. The balance equation of the DNP concentration for group i:

\begin{array}{l} \frac{\partial C_{i} (\overset{Δ}{r}, t)}{\partial t} & + \nabla \cdot [U C_{i} (\overset{Δ}{r}, t)] \\ = \nabla \cdot [(D_{c, i} + \frac{μ_{t}}{ρ_{fs} {Sc}_{t}}) \nabla C_{i} (\overset{Δ}{r}, t)] \\ + β_{i} \sum_{g^{'} = 1}^{I} {(ν Σ)}_{f, g^{'}} (\overset{Δ}{r}) ϕ_{g^{'}} (\overset{Δ}{r}, t) - λ_{i} C_{i} (\overset{Δ}{r}, t) . \end{array}

(10)

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 $\nabla \cdot D_{c, i} \nabla C_{i} (\overset{Δ}{r}, t)$ part represents the molecular diffusion of the DNP group i and the $\nabla \cdot (\frac{μ_{t}}{ρ_{fs} {Sc}_{t}}) \nabla C_{i} (\overset{Δ}{r}, t)$ 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, ug, Ci, D_n,g, D_c,i, βi, β, ν, $\sum_{g ’ \to g}$ , χ_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.

\begin{array}{l} \frac{1}{u_{g}} \nabla \cdot [U ϕ_{g} (\overset{Δ}{r})] = & \nabla \cdot D_{n, g} (\overset{Δ}{r}) \nabla ϕ_{g} (\overset{Δ}{r}) - Σ_{r, g} (\overset{Δ}{r}) ϕ_{g} (\overset{Δ}{r}) \\ + \frac{χ_{p, g} (1 - β)}{k_{eff}} \sum_{g^{'} = 1}^{G} {(ν Σ)}_{f, g^{'}} (\overset{Δ}{r}) ϕ_{g^{'}} (\overset{Δ}{r}) \\ + \sum_{i = 1}^{I} χ_{d, i, g} λ_{i} C_{i} (\overset{Δ}{r}) + \sum_{g^{'} = g}^{G} Σ_{g^{'} \to g} (\overset{Δ}{r}) ϕ_{g^{'}} (\overset{Δ}{r}), \end{array}

(11)

\begin{array}{l} \nabla \cdot [U C_{i} (\overset{Δ}{r})] = & \nabla \cdot [(D_{c, i} + \frac{μ_{t}}{ρ_{F} {Sc}_{t}}) \nabla C_{i} (\overset{Δ}{r})] \\ + \frac{β_{i}}{k_{eff}} \sum_{g^{'} = 1}^{I} {(ν Σ)}_{f, g^{'}} (\overset{Δ}{r}) ϕ_{g^{'}} (\overset{Δ}{r}) - λ_{i} C_{i} (\overset{Δ}{r}) . \end{array}

(12)

In addition, the estimate for the effective multiplication factor k_eff, may be computed from the Eq. (13):

k_{eff}^{n} = \frac{k_{eff}^{n - 1} ∭ \sum_{g = 1}^{G} {(ν Σ_{f})}_{g} ϕ_{g}^{n} d V}{∭ \sum_{g = 1}^{G} {(ν Σ_{f})}_{g} ϕ_{g}^{n - 1} d V},

(13)

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:

\begin{array}{l} \frac{1}{u_{1}} \nabla \cdot [U ϕ_{1, fs} (\overset{Δ}{r})] = & \nabla \cdot D_{n,1, fs} (\overset{Δ}{r}) \nabla ϕ_{1, fs} (\overset{Δ}{r}) \\ - Σ_{r,1, fs} (\overset{Δ}{r}) ϕ_{1, fs} (\overset{Δ}{r}) \\ + \frac{(1 - β)}{k_{eff}} \sum_{g^{'} = 1}^{2} {(ν Σ)}_{f, g^{'}, fs} (\overset{Δ}{r}) ϕ_{g^{'}, fs} (\overset{Δ}{r}) \\ + \sum_{i = 1}^{I} λ_{i} C_{i} (\overset{Δ}{r}), \end{array}

(14)

\begin{array}{l} \frac{1}{u_{2}} \nabla \cdot [U ϕ_{1, fs} (\overset{Δ}{r})] = & \nabla \cdot D_{n,2, fs} (\overset{Δ}{r}) \nabla ϕ_{1, fs} (\overset{Δ}{r}) \\ - Σ_{r,1, fs} (\overset{Δ}{r}) ϕ_{1, fs} (\overset{Δ}{r}) \\ + Σ_{1 \to 2, fs} (\overset{Δ}{r}) ϕ_{1, fs} (\overset{Δ}{r}) . \end{array}

(15)

The neutron diffusion equations in the reflector:

\nabla \cdot D_{n, 1, ref} (\overset{Δ}{r}) \nabla ϕ_{1, ref} (\overset{Δ}{r}) - Σ_{r, 1, ref} (\overset{Δ}{r}) ϕ_{1, ref} (\overset{Δ}{r}) = 0,

(16)

\begin{array}{l} \nabla \cdot D_{n, 2, ref} (\overset{Δ}{r}) ▽ ϕ_{1, ref} (\overset{Δ}{r}) & - Σ_{r, 2, ref} (\overset{Δ}{r}) ϕ_{1, ref} (\overset{Δ}{r}) \\ + Σ_{1 \to 2, ref} (\overset{Δ}{r}) ϕ_{1, ref} (\overset{Δ}{r}) \\ = 0 . \end{array}

(17)

The transport equations for six-group DNPs in the fuel salt:

\begin{array}{l} \nabla \cdot [U C_{i} (\overset{Δ}{r})] = & \nabla \cdot [(D_{c, i} + \frac{μ_{t}}{ρ_{fs} {Sc}_{t}}) \nabla C_{i} (\overset{Δ}{r})] \\ + \frac{β_{i}}{k_{eff}} \sum_{g^{'} = 1}^{2} {(ν Σ)}_{f, g^{'}, fs} (\overset{Δ}{r}) ϕ_{g^{'}, fs} (\overset{Δ}{r}) \\ - λ_{i} C_{i} (\overset{Δ}{r}) . \end{array}

(18)

The effective multiplication factor k_eff:

k_{eff}^{n} = \frac{k_{eff}^{n - 1} ∭ \sum_{g = 1}^{2} {(ν Σ_{f})}_{g, fs} ϕ_{g, fs}^{n} d V}{∭ \sum_{g = 1}^{2} {(ν Σ_{f})}_{g, fs} ϕ_{g, fs}^{n - 1} d V},

(19)

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

1. Thermohydraulics

(1) Inlet boundary: at the inlet the velocity, the turbulent kinetic energy, the turbulent dissipation rate and the temperature are given:

\begin{array}{l} u_{z} = u_{inlet}, \\ k = \frac{3}{2} {(I u_{z})}^{2}, \\ ε = C_{μ}^{0.75} \frac{k^{1.5}}{l}, \end{array}

(20)

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

(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:

C_{i, inlet} = C_{i, outlet} \cdot e^{- λ_{i} τ},

(21)

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

Generally, the thermophysical properties in the thermohydraulics and the group constant in the neutronics are both dependent on the temperature. Ignatiev et al. [6] gave the LiF/NaF/BeF₂ solvent system properties including the density, the thermal conductivity, the viscosity and the heat capacity.

ρ_{fs} (kg \cdot m^{- 3}) = 2518 - 0.406 T,

(22)

λ_{T, fs} (W \cdot m^{- 1} \cdot K^{- 1}) = 0.0429 + 0.0009 T,

(23)

μ (Pa \cdot s) = 0.001 e^{- 0.9942 + 1603.2 / T},

(24)

C_{p, fs} (J \cdot {kg}^{- 1} \cdot K^{- 1}) = 2090.

(25)

The thermophysical properties of the graphite reflector are considered as constants (Table 4).

The thermophysical properties of the graphite reflector at 873.15 K [11]

Properties	Graphite
ρ_ref(kg·m^-3)	1840
λ_T,ref(W·m^-1·K^-1)	31.2
C_p,ref(J·kg^-1·K^-1)	1760

The MSR3DS code system uses DRAGON4 [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:

\begin{array}{l} Y_{c} (T) = & A_{0} + A_{1} (T - T_{avg}) + A_{2} {(T - T_{avg})}^{2} \\ + A_{3} {(T - T_{avg})}^{3} + A_{4} {(T - T_{avg})}^{4} \\ + A_{5} {(T - T_{avg})}^{5} \end{array}

(26)

where Yc denotes the group constants and T_avg is reference temperature. A₀–A₅ are fitting constants.

E. Numerical method

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].

The SIMPLE [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. 3 is the program flow diagram of the coupled solver in the MSR3DS code.

Fig. 2.

The SIMPLE algorithm.

Fig. 3.

The program flow diagram in the coupled solver.

F. The description of the MSR3DS code

The entire program diagram of the MSR3DS code is shown in 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.

Fig. 4.

The program diagram of the MSR3DS code.

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] code is as the post-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

A. Code validation

The 3D TWIGL Seed/Blanket problem [17] is adopted to benchmark the neutronics calculation in the MSR3DS code. The effective multiplication factors calculated by CITATION [18] and MSR3DS are listed in Table 5, and the relative error is only 0.00147%. The fluxes normalized to 1 MW_th are shown in Figs. 5 and 6. The results verify validity of the model for neutronics presented in this study.

The effective multiplication factors of benchmark

Benchmark	k_eff		Relative error/%
3D TWIGL	CITATION	MSR3DS
Seed/Blanket	0.885905	0.885892	0.00147

Fig. 5.

(Color online) Distributions of the neutron fluxes calculated by CITATION and MSR3DS at (Z=79:2 cm,Y=40 cm).

Fig. 6.

(Color online) Distributions of the neutron fluxes calculated by CITATION and MSR3DS at (X=40 cm,Y=40 cm).

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). 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 and 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.

Fig. 7.

(Color online) A simple pipe for validation of the flow and heat transfer (Unit/m).

Fig. 8.

(Color online) Velocities along the symmetry axis.

Fig. 9.

(Color online) Velocities along radial direction at the midplane.

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 and 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.

Fig. 10.

(Color online) Temperatures at the symmetry axis.

Fig. 11.

(Color online) Temperatures at the mid-plane.

B. Distributions of calculated physical fields in the core

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-13 show the 3D distributions of the velocity and the turbulent kinematic viscosity. Fig. 13, the turbulent kinematic viscosity in the outlet is far greater than that of other locations.

Fig. 12.

(Color online) Velocity field in the core (U/m·s^-1).

Fig. 13.

(Color online) Turbulent kinematic viscosity in the core (vt·s).

Figures 14-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-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-15(b) with Figs. 14-15(c).

Fig. 14.

(Color online) Fast neutron fluxs without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/neutron·m^-2·s^-1).

Fig. 15.

(Color online) Thermal neutron fluxs without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/neutron·m^-2·s^-1).

Figures 16, 17, 18, 19, 20, 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, 17, 18, 19, 20, 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, 17, 18, 19, 20, 21(c) are also changed observably, likewise, the smaller the delay constant, the greater the influence of the turbulent diffusion.

Fig. 16.

(Color online) C1 group precursor concentration without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/m^-3).

Fig. 17.

(Color online) C2 group precursor concentration without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/m^-3).

Fig. 18.

(Color online) C3 group precursor concentration without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/m^-3).

Fig. 19.

(Color online) C4 group precursor concentration without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/m^-3).

Fig. 20.

(Color online) C5 group precursor concentration without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/m^-3).

Fig. 21.

(Color online) C6 group precursor concentration without flow (a), with convective term (b), with convective and turbulent diffusion term at Sc_t = 0.7 (c) (Unit/m^-3).

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. 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.

Fig. 22.

(Color online) Temperatures with convective term (a), with convective and turbulent diffusion term at Sc_t = 0.7 (b).

IV. CONCLUSION

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

[1]

Bettis E S, Cottrell W B, Mann E R, et al. Nucl Sci Eng, 1957, 2: 804-825.

[2]

Robertson R C.

MSRE design and operations report, part I, description of reactor design. Technical report: ORNL–TM–0728

, 1965.

Baidu Scholar

Google Scholar

[3]

Rosenthal MW, Haubenreich P N, Briggs R B.

The development status of molten–salt breeder reactors. Technical report: ORNL–4812

, 1972.

Baidu Scholar

Google Scholar

[4]

Furukawa K, Lecocq A, Kato Y, et al. J. Nucl Sci Technol, 1990, 27: 1157-1178.

[5]

Vergnes J and Lecarpentier D. Nucl Sci Eng, 2002, 216: 43-67.

[6]

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.

Baidu Scholar

Google Scholar

[7]

Holcomb D E, Flanagan G F, Patton B W, et al.

Fast spectrum molten salt reactor options. Technical report ORNL/TM–2011/105

, 2011.

Baidu Scholar

Google Scholar

[8]

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.

[9]

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.

Baidu Scholar

Google Scholar

[10]

Wang S, Rineiski A, Maschek W. Nucl Sci Eng, 2006, 236: 1580-1588.

[11]

Nicolino C, Lapenta G, Dulla S, et al. Ann Nucl Energy, 2008, 35: 314-322.

[12]

Zhang D L, Qiu S Z, Su G H, et al. Ann Nucl Energy, 2009, 36: 590-603.

[13]

Versteeg H K, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method. England: Longman Group Ltd, 1995.

[14]

Marleau G, Hebert A, Roy R. A user guide for dragon version4. Qu´ebec: Ecole Polytechnique de Montreal, 2013.

[15]

Saad Y. Iterative methods for sparse linear systems, 2nd edition. Pennsylvania: Society for Industrial and Applied Mathmatics, 2003.

[16]

ParaView Guide. New York: Kitware, Inc., 2013.

[17]

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).

Baidu Scholar

Google Scholar

[18]

Fowler T B, Vondy D R, Cunningham GW. Nuclear reactor core analysis code: CITATION, ORNL–TM–2496. Tennessee: Oak Ridge National Laboratory, 1971.