2.1 Neutronics model

OpenMC is an open source Monte Carlo particle transport simulation code that can perform fixed source, k-eigenvalue, and subcritical multiplication calculations on three-dimensional solid geometry models [21]. It was originally developed by the Massachusetts Institute of Technology and is widely used for calculations of the criticality, depletion, multigroup cross section generation, etc. The Direct Accelerated Geometry technology has been integrated into OpenMC (DAG-OpenMC hereafter) to perform neutron transport simulation directly on a computer-aided design model [22]. In this work, this technology contributes to the solution of the problem of curved wall modeling for an MSFR owing to its advantages for precise modeling of an arbitrary complex geometry.

Based on the ENDF/B-VII.1 nuclear data, a microscopic cross section library between 500 K and 1600 K with 100K internals was produced by NJOY2016 to support the OpenMC continuous-energy transport calculation; a 20 K interval was used between 800 and 1400 K (the operating temperature range of MSFRs). The cross sections are given exactly for only a few points in the table; the temperature values for other materials can be obtained by statistical interpolation [23]. Suppose that cross sections are available at temperatures T₁, T₂,..., T_n, and a material is assigned a temperature T, where T_i < T < T_i+1. Statistical interpolation is applied as follows. A uniformly distributed random number (ξ) of unit intervals (0–1) is sampled. If ξ < (T − T_i)/(T_i+1 − T_i), the cross sections at temperature T_i+1 are used. Otherwise, the cross sections at T_i are used. To obtain the power distribution of the MSFR, the calculation model is divided into numerous regions, and each region is labeled with a unique material number.

2.2 Thermal-hydraulic model

The open source code OpenFOAM features automatic matrix construction and solution capabilities for scalar and vector equations [24]. Each fluid flow can be represented by a series of partial differential equations; the solvers for these fluid mechanics equations have been developed and implemented in OpenFOAM.

The fuel salt flow in the MSFR core without a moderator is considered to be turbulent [25] because the Reynolds number is estimated to be 2.5 × 10⁶ in the core region and 5.3 × 10⁵ in the out-of-core region [11]. The fuel and fertile salt flows can be assumed to be incompressible because the effect of fuel density variation in the gravity force term is taken into account by using the Boussinesq approximation. Moreover, the ReNormalization Group (RNG) k–ε turbulence model can provide a better description of the flow on the curved wall than the standard k–ε turbulence model [26]. Therefore, considering that the fuel salt flow is incompressible, the Reynolds-averaged Navier–Stokes equations with the Boussinesq closure hypothesis and the RNG k–ε turbulence model are adopted, and buoyantBoussinesqSimpleFoam [24] is selected as the OpenFOAM solver. The steady-state continuity, momentum, and energy conservation equations for the fuel and fertile salts can be written as

where $\overrightarrow{\text{u}}$ represents the Reynolds mean flow velocity vector, m∙s⁻¹; p represents the Reynolds average pressure, which is divided by the reference density, m²∙s⁻²; ν and ν_t represent the dynamic viscosity and turbulent dynamic viscosity, respectively, m²∙s⁻¹; $\overrightarrow{\text{g}}$ is the gravitational acceleration, m∙s⁻²; T represents the Reynolds-averaged temperature, K; β is the thermal expansion coefficient in K⁻¹ at a given reference temperature T_ref in K; S_T is the heat source term, K∙s⁻¹; and I, Pr, and Pr_t represent the identity matrix, Prandtl number, and turbulent Prandtl number, respectively. The turbulent kinetic energy k and turbulent dissipation rate $\epsilon$ are modeled by the RNG k–ε turbulence model, which is given by

with

$R=\frac{\eta \left(1+\eta /{\eta }_0\right)}{1+\alpha {\eta }^3}$; ${\nu }_{\text{t}}=C_{\mu }\frac{k^2}{\epsilon }$; $G=2{\nu }_{\text{t}}{S}^2$; $S=\frac{1}{2}\left(\nabla \overrightarrow{u}+{\left(\nabla \overrightarrow{u}\right)}^{\text{T}}\right)$; $\eta =\sqrt{2}Sk/\epsilon$; $C_{\mu }=0.085$; $C_{\epsilon 1}=1.42$; $C_{\epsilon 2}=1.68$; ${\sigma }_k=0.71942$; ${\sigma }_{\epsilon }=0.71942$; ${\eta }_0=4.38$; $\alpha =0.012$.

2.3 Coupling scheme

The coupling calculation flowchart of DAG-OpenMC and OpenFOAM is presented in Fig. 1, illustrating the exchange of core power distribution, temperature, and density distribution of the fuel and fertile salt using a Python script (pyMF). As each region of the neutronics calculation model contains many meshes for the TH calculation by OpenFOAM, the average fuel temperature and density of each region are assigned according to the material number (where each region has a unique material number) used in DAG-OpenMC. Similarly, the OpenFOAM meshes located in the same region are assigned the same power density from the neutronics calculation.

Fig. 1Full-screen View

Flowchart of the coupled code

The key points of the coupling code pyMF are summarized as follows:

1) The temperature and density of the fuel salt and structural material are first assumed to be uniform in the DAG-OpenMC input file to start the neutronics calculation.

2) The power distribution is obtained by transforming the kappa-fission tally of each region from OpenMC to the power density as follows:

where p_i is the power density of region i, J∙m⁻³∙s⁻¹; ν represents the neutron yield per fission, neutron/fission; P is the thermal power of the MSFR, W; Q represents the released energy per fission in the MSFR, eV/fission; V_i is the volume of region i, m³; k_eff represents the effective multiplication factor obtained from the DAG-OpenMC criticality calculation, neutron/source neutron; and T_i is the kappa-fission tally result for region i, eV/source neutron.

3) The heat source term in the energy conservation equation [Eq. (3)] is updated from the power distribution, and OpenFOAM is launched to start the TH calculation. The calculation is terminated when the temperature and velocity fields meet the residual convergence criteria set by users; then the temperature and density distributions can be obtained.

4) The above iteration is performed until either k_eff and the power density converge or the maximum number of iterations is reached [27–29].

In our previous work [30], to verify the reliability of the coupled code, a benchmark model of the MSFR was established, and the effects of the region number in the neutronics calculation and the initial conditions on the k_eff, fuel salt velocity, and temperature distributions were studied. The TH results simulated using pyMF were also compared to the reference result for the MSFR [31], and good agreement was obtained for the velocity and temperature fields (Table S1 and Figs. S1 and S2). It can be concluded that the pyMF code can provide reliable results for the steady-state neutronics and TH coupling calculations for an MSFR.

3.1 Description of the molten salt fast reactor

A schematic view of the MSFR system is shown in Fig. 2. The MSFR is a 3 GWth liquid-fueled reactor with a core fuel salt volume of 9 m³ [9, 19, 20]. In the primary loop, a molten fluoride salt system functioning as both fuel and coolant works at approximately 923 and 1023 K at the inlet and outlet of the core, respectively. A Hastelloy [19, 20] is selected as the axial reflector, and a radial fertile blanket is adopted to absorb the neutrons escaping from the core and thus improve the breeding capability of the reactor. The primary loop of the MSFR has 16 fuel salt loops, and each loop is equipped with a pump and a heat exchanger. The fuel salt circulates through the core and the 16 external loops, and the nuclear heat generated by fission is transferred to the secondary loop salt by the heat exchangers. The main parameters of the MSFR reactor core are listed in Table 1 [12, 20, 31].

Fig. 2Full-screen View

(Color online) Schematic view of the MSFR primary circuit. HX: heat exchanger

3.2 Core thermal-hydraulic analysis

As illustrated in Fig. 3, three core geometries are proposed in this work to analyze the TH behavior of the MSFR. Geometry I [Fig. 3 (a)] is described in the literature [17, 18]. Geometry II [Fig. 3 (b)] is obtained by further curving the upper wall of the core, and Geometry III [Fig. 3 (c)] is obtained by curving both the upper and lower walls of the core for comparison. All the arcs of the blanket and absorber in the three geometries are half-ellipses.

Fig. 3Full-screen View

(Color online) Schematic maps of optimized MSFR geometries (All dimensions are in centimeters.) a. Geometry I; b. Geometry II; c. Geometry III

In the OpenFOAM simulation, only the core, fertile blanket, and primary loop are modeled, and the solid regions, including the reflector and absorber, are not considered, because the coupling calculation is performed only in the above three regions. To obtain a good balance between the calculation efficiency and accuracy, the maximum bulk cell size is set to 3 cm [16, 17]. The RNG k–ε turbulence model and a wall function are adopted to model the flow behavior near the wall. In fluid dynamics, a dimensionless distance to the wall (y⁺) is used as an evaluation criterion of near-wall mesh quality. Therefore, a suitable y⁺ value is needed to obtain accurate results when the k–ε turbulence model and wall function are employed. To ensure that the y⁺ values of all the volume elements adjacent to the wall are distributed within a range of 11.63–300 [32], the cell size is reduced to approximately 3 mm near the core wall. Considering that the MSFR core configuration has good symmetry, only 1/16 of the core geometry (Fig. 4) is simulated in the TH calculation to conserve computational resources. The pump and heat exchanger are treated as a momentum source in the momentum conservation equation and as a heat sink in the energy conservation equation. The mesh numbers of these three geometries are 124,540 [Fig. 4 (a)], 132,784 [Fig. 4 (b)], and 138,562 [Fig. 4 (c)], respectively.

Fig. 4Full-screen View

(Color online) CFD meshes adopted in core TH calculation for Geometry I (a), Geometry II (b), and Geometry III (c)

A continuous-energy neutronics simulation is carried out using the DAG-OpenMC code; the calculation models (Fig. 5) include the core, fertile blanket, and solid regions. All of the components are simulated during the neutron transport simulation; a reflective boundary condition is set for the symmetry planes, and a vacuum boundary condition is set for the walls. No heat exchange is modeled between the core and the fertile blanket or between the fertile blanket and the solid region, and the temperatures of the fertile blanket and solid region are assigned as the average temperature of the fuel salt in the core TH analysis. To obtain the power distribution for the TH calculation, the core is divided into numerous regions, and each region is labeled with a unique material number during the neutronics calculation. Considering the efficiency and accuracy of the neutronics calculation, the size of these regions should be approximately 20 cm, as demonstrated in our previous work [30]. Only a small fraction of the power is released in the external loops, and the strong coupling between the neutronics and thermal hydraulics is significant in the core region. Therefore, no detailed regional division of the external loops is performed. For the neutron transport simulation, 300 generations (the first 50 iterations are skipped) with 50,000 neutrons per generation are adopted, and the obtained standard deviation of k_eff is approximately 0.00023.

Fig. 5Full-screen View

(Color online) Geometric models for neutron transport simulation for Geometry I (a), Geometry II (b), and Geometry III (c)

The fuel velocity fields of the three geometries are shown in Fig. 6, which shows that a low flow velocity occurs near the outer wall (close to the fertile blanket), whereas the large recirculation region near the fertile blanket in the reference MSFR model (see Fig.S1) is eliminated, and a relatively flat velocity distribution is obtained for all three geometries. Although a small vortex or recirculation pattern still occurs at the core inlet, no hot spot forms owing to the low power density in the area of interest and the low fuel salt temperature at the core inlet. Therefore, further optimization to eliminate this small vortex and recirculation pattern is not needed. The connecting pipe is much smaller than the core; thus, the flow area is widened or narrowed when the fuel salt is flowing into or out of the core. Furthermore, the connecting pipe is curved. Therefore, a large local pressure drop is obtained, which maximizes the flow velocity in the connecting pipe in the three geometries. The maximum fuel salt velocity obtained for Geometry I [Fig. 6 (a)] is approximately 5.05 m∙s⁻¹, which is slightly higher than those for Geometry II [Fig. 6 (b)] and Geometry III [Fig. 6 (c)]. The reason is that the connecting pipe in Geometry I is smaller than that in the other two geometries. As shown in Fig. 6 (a), a small stagnation zone appears near the upper and lower reflectors in the velocity field of Geometry I, whereas a stagnation zone appears in Geometry II only at the core bottom, as shown in Fig. 6 (b). For Geometry III, the stagnation zone is eliminated [Fig. 6 (c)], which indicates that Geometry III allows a more uniform velocity distribution.

Fig. 6Full-screen View

(Color online) Fuel salt velocity field results for Geometry I (a), Geometry II (b), and Geometry III (c)

The fuel salt temperature field in each geometry is presented in Fig. 7. The fuel salt temperature increases along the axial direction from the bottom to the top for all three geometries, and a small hot spot appears at the top of the central core in Geometry I [Fig. 7 (a)], which might be caused by the low velocity in this region. The temperature profile for Geometry II [Fig. 7 (b)] does not differ significantly from that for Geometry I, although the peak value of the fuel temperature is slightly lower owing to the difference in the geometry. Therefore, the hot spot at the upper core cannot be effectively eliminated by curving only the upper wall of the core. Unlike these two geometries, Geometry III [Fig. 7 (c)] exhibits a more uniform temperature distribution throughout the core and a relatively small radial gradient of the fuel temperature. Moreover, the peak temperature drops sharply, and the difference between the maximum core temperature and average outlet temperature is only approximately 9 K.

Fig. 7Full-screen View

(Color online) Fuel salt temperature distribution for Geometry I (a), Geometry II (b), and Geometry III (c)

The turbulent dynamic viscosity field of each geometry is presented in Fig. 8. These results are also compared to the result for the reference MSFR model (see Fig.S3). The location of maximum turbulent dynamic viscosity is shifted from the upper core cavity to the lower core cavity. In addition, a zone of high turbulent dynamic viscosity is formed near the entire arc wall of the fertile blanket. Therefore, there is a strong turbulent transport effect in this region, which prevents the formation of the hot spot on the arc wall near the fertile blanket even though the core flow velocity is low in this region (see Figs. 6 and 7). The distributions of the turbulent dynamic viscosity for Geometry I [Fig. 8 (a)] and Geometry II [Fig. 8 (b)] are both relatively flat; hence, the heat transfer in the core due to turbulent transport is almost uniform. By contrast, the turbulent dynamic viscosity in the core center in Geometry III [Fig. 8 (c)] is low, and it increases gradually along the radial direction; therefore, the turbulent transport effect also increases gradually. As shown in Fig. 6, the mean flow velocity in the core center is higher than that near the fertile blanket; therefore, the cooling capacity resulting from fuel salt flow in the core center is better than that near the fertile blanket. Owing to the competing effects of fuel salt flow and turbulent transport, Geometry III has the smallest radial temperature gradient in the core and the most uniform radial temperature distribution among the three optimized geometries.

Fig. 8Full-screen View

(Color online) Turbulent dynamic viscosity field for Geometry I (a), Geometry II (b), and Geometry III (c)

Table 2 lists the average core inlet and outlet temperatures, the temperature difference across the core, the core volumes, the total volume flow rate in the primary loop, and the maximum core temperature of each geometry. The heat exchanges are modeled as suitable heat sinks in the energy conservation equations to obtain the nominal value (923 K) of the core inlet temperature. A different momentum source is added to the momentum conservation equation for each geometry to maintain the temperature difference across the core at the nominal level (100 K), because the core volume and shape of each geometry is different. The temperature difference across the core and the average core inlet and outlet temperatures for Geometry I, II, and III are all approximately 100, 923, and 1023 K, respectively. However, the maximum temperatures differ significantly. The maximum fuel salt temperature in Geometry III is 1032.5 K, which is much lower than those of Geometry I (1070.3 K) and Geometry III (1062.1 K). Geometry III exhibits the best TH performance.

3.3 Preliminary thermal-hydraulic analysis of the fertile blanket

The fertile blanket has a much smaller heat source than the core; hence, it is neglected in the TH studies of the core. However, the effect of the core heat on the fertile blanket is unknown, and the feasibility of modeling the fertile blanket separately should be analyzed. In this work, a conjugate heat transfer model of the core and fertile blanket is established; it is shown schematically in Fig. 9. The calculation model presented in Fig. 5 (c) is adopted for the neutronics simulation, and the corresponding TH calculation mesh is shown in Fig. 10. A 3 GWth heat source is loaded in the core region, and three different power loading cases for the fertile blanket are analyzed. An adiabatic boundary condition is set for the blanket wall. In the OpenFOAM code, a heat flux balance is maintained on the interface wall between the core and the fertile blanket, which can be described as

Fig. 9Full-screen View

Schematic illustration of the conjugate heat transfer model

Fig. 10Full-screen View

(Color online) TH calculation mesh for conjugate heat transfer

where q is the heat flux; T_B is the temperature of the interface wall; and λ_c and λ_b are the thermal conductivity of the fuel salt and fertile salt, respectively. T_c and T_b are the core and fertile blanket cell center temperature of the first control volume on the interface wall, respectively; δ_c and δ_b are the distance of the interface wall from the core and the fertile blanket cell center, respectively.

The temperature distribution including the conjugate heat transfer of the fuel and fertile salts is illustrated in Fig. 11. A total heat source of 0 MWth [Fig. 11 (a)], 0.1 MWth [Fig. 11 (b)], and 1 MWth [Fig. 11 (c)] is loaded in the fertile blanket. Fig. 11 (a) shows that the temperature of the fertile blanket matches that of the core in the axial direction even when no heat source is loaded in the fertile blanket, which reveals that the nuclear heat of the core affects the temperature of the fertile blanket. Fig. 11 (b) and 11 (c) indicate that the nuclear heat generated in the fertile blanket cannot be removed by the fuel salt flow through the interface wall between the core and the fertile blanket. Therefore, an external cooling system for the fertile blanket should be considered. ²³³U will be generated in the fertile blanket because the thorium loaded in this region is exposed to the core neutron flux, and the heat source (including fission and neutron capture by thorium) generated in the fertile salt is estimated to be approximately 37 MW [19, 20]. However, the temperature of the fertile blanket greatly exceeds that of the core when a heat source of only 1 MWth is loaded in the fertile blanket [Fig. 11 (c)]. Therefore, from Fig. 11 (a) and 11 (c), one can deduce that the effect of the core heat on the fertile blanket is much smaller than that of the heat source in the fertile blanket itself, which demonstrates that modeling the fertile blanket separately is reasonable for the MSFR TH analysis.

Fig. 11Full-screen View

(Color online) Temperature distribution considering the conjugate heat transfer calculation for the fuel and fertile salts for (a) 0 MWth, (b) 0.1 MWth, and (c) 1 MWth heat sources in the blanket

Three cases are designed to study the required flow rate for obtaining an acceptable temperature distribution in the fertile blanket; the external circuit is ignored for convenience. A sketch map of the three cases is presented in Fig. 12. (case 1 [Fig. 12 (a)], case 2 [Fig. 12 (b)] and case 3 [Fig. 12 (c)]). To reduce the computational demands, a minimal symmetrical unit is adopted for the TH calculation; the CFD meshes for the three cases are shown in Fig. 13. According to the grid independence verification, the mesh numbers for case 1 [Fig. 13 (a)], case 2 [Fig. 13 (b)], and case 3 [Fig. 13 (c)] are 108,912, 56,576, and 29,772, respectively.

Fig. 12Full-screen View

(Color online) Sketch map of the geometry models of the fertile blanket for cases 1 (a), case 2 (b), and case 3 (c)

Fig. 13Full-screen View

(Color online) CFD meshes of fertile blanket adopted for TH calculations for cases 1 (a), case 2 (b), and case 3 (c)

The core region of Geometry III [Fig. 3 (c)] is selected to consider the effect of the core neutronics on the neutronics behavior of the fertile blanket. Consequently, the fertile blanket is also divided into several regions to obtain its power distribution, as illustrated in Fig. 14. All of the parts are simulated during the neutron transport simulation; a reflective boundary condition is set for the symmetry planes, and a vacuum boundary condition is set for the walls. The heat source of the fertile blanket is estimated to be approximately 37 MW, as mentioned above. Therefore, the proportion of ²³³U loaded in the fertile blanket is set to match the power deposition.

Fig. 14Full-screen View

(Color online) Neutronics calculation model used for fertile blanket

Considering the temperature tolerance of the structural material, the temperature difference (temperature gradient) between the core and the fertile blanket should be as small as possible. Therefore, the inlet temperature of the fertile blanket is set to 923 K, which is consistent with the average core inlet temperature, by setting the inlet boundary condition to 923 K. In addition, the volume flow rate is adjusted so that the average outlet temperature of the fertile blanket is consistent with the average core outlet temperature (1023 K).

The velocity vector profile of the fertile blanket for the three cases is shown in Fig. 15. The inlet jet for case 1 [Fig. 15 (a)] is larger than those for case 2 [Fig. 15 (b)] and case 3 [Fig. 15 (c)]. The main reason that the inlet velocity decreases as the number of external circuits increases is that when there are more external circuits, a lower inlet velocity for each pipe can maintain the same volume flow rate of the fertile blanket. A strong inlet jet causes a large reverse flow phenomenon in the entire fertile blanket, as illustrated in Fig. 15 (a). A small vortex appears in Fig. 15 (b) and 15 (c), and no reverse flow occurs near the arc wall owing to the relatively low inlet flow velocity.

Fig. 15Full-screen View

(Color online) Velocity vector profile of the fertile blanket for cases 1 (a), case 2 (b), and case 3 (c)

Fig. 16 shows the temperature distribution of the fertile blanket for the three cases. The temperature distribution of the fertile blanket becomes more uniform as the number of external circuits increases. To maintain a constant outlet temperature of approximately 1023 K, the total volume flow rates of the fertile blanket for the three cases are almost identical to each other owing to their similar geometries. As illustrated above, the required inlet flow velocity for case 1 is highest among the three cases to match the equivalent volume flow rate. As a result, the inlet jet in case 1 causes a highly nonuniform temperature distribution, and a relatively high temperature appears at the bottom arc, which is higher than the outlet temperature of the fertile blanket. Furthermore, a local high temperature, which is much higher than the outlet temperature, appears in the top cavity of the fertile blanket. By contrast, because no reverse flow occurs near the arc wall, the temperature distribution of the top cavity for cases 2 and 3 is relatively flat, and no hot spot forms.

Fig. 16Full-screen View

(Color online) Temperature distribution of the fertile blanket for cases 1 (a), cases 2 (b), and cases 3 (c)

Table 3 shows the maximum temperature, average outlet temperature, and required total volume flow rate of the fertile blanket. Case 3 has the lowest maximum temperature and the most uniform temperature distribution. Therefore, it has the best TH performance among the three cases. However, the calculated maximum temperature and temperature distribution for case 2 are also acceptable; this result indicates that eight cooling circuits with a total volume flow rate of approximately 0.0555 m³∙s⁻¹ can also produce an acceptable temperature distribution. Considering the engineering requirements, the geometry model for case 2 is recommended.