I. INTRODUCTION
The Generation IV International Forum (GIF) on nuclear energy aims at developing technologies to achieve safety performance, waste reduction and increased proliferation resistance, and to provide a nuclear energy option that is economically competitive [1]. Of the six future nuclear power plants under the framework of GIF, only the Molten Salt Reactor (MSR) is of the liquid fuel type. The concept of MSR was initially proposed by Oak Ridge National Laboratory (ORNL) [2, 3] where the feasibility was proved by both theories and experiments. In 2011, Chinese Academy of Sciences launched the project of Thorium Molten Salt Reactor Nuclear Energy System (TMSR). As part of the pre-study project, the molten salt reactor transient analysis considers key parameters (the power, flow rate, temperature and the controlling orders) and the physical properties of the core (flow rate of the initial zero-power and low-power mode, temperature and the primary loop pressure). The transient analysis software also judges status of the reactor via different parameters. Therefore, the simulation methods are of much concern.
Current researches on MSR mainly focus on finding different descriptive approaches of the reaction procedures. Lapenta et al. [4] analyzed neutronics in a fluid fuel system through a point-kinetics model. Dulla et al. [5, 6] used the quasi-static method, which considered a prescribed velocity field in one direction without taking the temperature cross-section feedback into account. Yamamoto et al. [7, 8] and Suzuki et al. [9] analyzed the steady-state and transient state for a small molten salt reactor by coupling the neutron diffusion equations with the heat transfer in fuel salt and graphite. Krepel et al. used the DYN1D–MSR [10] and DYN3D–MSR [11] codes for solving the two-group neutron diffusion equation with a point kinetics model, and models of delayed neutrons drift and heat release distribution were included. The prompt fission neutron source was calculated by means of the same nodal flux expansion method. Part of the neutrons was released with a time delay when the precursor core drifted with the liquid fuel along axial direction. An extended delayed neutrons model was developed. Transient simulating of 3D neutronics and parallel channel thermal hydraulics were used.
Lecarpentier et al. [12] developed the Cinsf1D program, which replaced MSRE at Oak Ridge to become the European benchmark of the transient simulation, and to study the AMSTER system [13]. Both steady-state and transient state could be simulated by the Cinsf1D program. After modeling the salt fuel and the graphite in the reactor, the axial and radial power distribution of the core were calculated by the DRAGON program [14]. Based on the diffusion theory with two energy groups, the neutron flux and the distribution of delayed neutron precursors were described in neutronics according to the cross sections. The thermal negative feedback depending on the cross sections was related to the temperatures of various mediums. Thermo hydraulics could be described by the energy balance between the salt fuel and graphite. A simplified geometric model represented the typical core channel in the molten salt reactor. Also, a point kinetic model was included, where graphite temperature was a single average value, and the relationship between the temperature of the molten salt and the height of the core was linear.
In this paper, we report our work on simulating neutron diffusion and transient analysis in an MSR. Optimization of the model based on Cinsf1D will be discussed in Sec. II, with a focus on the transient analysis, which is separated from the point kinetic model. The implementation of object-oriented programming is discussed in Sec. III. The software architecture is improved to meet the standards of software-engineering, high-performance computing and management via graphical user interface, and to acquire the simulation data and diagrams in visualized forms. The results are shown in Sec. IV. In Sec. V, we conclude that the neutron diffusion and transient analysis of the MSR are simulated in a precise and efficient way.
II. THEORETICAL MODEL OF THE MSR PROCESS
A vertical-section diagram of the primary loop is shown in Fig. 1. The molten reactor core consists of graphite reflector. The molten salt fuel contains beryllium fluoride, lithium fluoride, zirconium fluoride, and uranium or thorium fluoride. With the fissile materials and convertible molten salt at 600 ℃, the molten salt fuel exits from the core. The fuel salt flows sequentially from the left of the core to Pipe 1, next to the heat exchanger, then to Pipe 2, and finally back to the core. The molten salt reactor core is a cylindrical cavity with 199 parallel hexagonal graphite elements, each of which has a circular coolant channel. Each channel has a fuel passage in the middle. The calculation of the reactor core is divided into three zones, with h1, h2, and h3 describing the length of each zone, rg being the channel diameter and rf being the inner diameter. The liquid salt serves as both the coolant and fuel. The graphite moderator in the reactor is heated by the neutrons and gamma rays produced in the fission reaction process. The heat between the graphite moderator and coolant is exchanged through the heat exchanger. Controlling of the reaction is modeled as follows.
A. Neutron Diffusion Modeling
In the neutron diffusion modeling, the neutrons are of two energy groups (the fast and thermal neutron energy groups) and one spatial dimension (z). The neutron energy groups are cut at 1.8 eV. The flux of the molten salt reactor is calculated by the thermal spectrum. It has been shown that the axial direction (z) and radial direction (x) are relatively decoupled. The main object of the model is to compare the influence of liquidity for the delayed neutron precursors, which are produced only in the core, drifted along axial (z) with the MSR core movement. The life spans of these neutrons are much shorter than the molten salt fuel cycle, so the neutrons are not affected by the velocity field [15, 16].
The neutron and precursors dynamics equations in the fuel material of molten salt reactor are derived from the particle conservation equations as follows:
where, i is the delayed neutron group number. Vj=1,2 indicates the neutron speed of group j, with j=1 being the fast neutron energy group and j=2 the thermal neutron energy group. ϕj(z,t) is the neutron flux of group j, t is the time, Dj(z,t) is the diffusion coefficient of group j. βj is the proportion of emitted heat from the delayed neutron of group i, and β is the sum of all βj, ∑fj(z) and ∑αj(z) are the macroscopic neutron fission cross-section and the macroscopic neutron absorption cross-section of the group j, respectively; ∑r(z) is the neutron diffusion cross-section from the fast neutron group towards the thermal group; B2 is the radial leakage coefficient; λj is the delayed neutron decay constant of group i; Cj(z,t) is the linear precursor concentration of group i. ψ(z,t) is the number of neutrons emitted by fission per unit height of the reactor; and V(z,t) is the molten salt speed vector. Due to the long decay time of the delayed neutrons, they are supposed to drift with the molten salt, which means that the pioneer nuclear has the same speed with the molten salt without considering the interaction of them.
B. Thermal Negative Feedback
The diffusion cross-section can be written as a function of z and temperature. Thus, according to the diffusion theory, the cross-section as a function of z at a certain temperature can be calculated through the interpolation using standard cross-section data. Thus, Eq. (1) can be solved to obtain the flux and power released. The reactivity is affected by temperature and balance configuration of salt and graphite, and the control rod depth, but preset data can be prepared to evaluate the diffusion process simulation.
C. Thermal Hydraulics
The calculation is divided into two parts, the core and the heat exchanger, using lumped parameter model.
The core thermal hydraulics calculation is calculated in two steps: (1) regarding the core graphite and fuel salt as heat sources and calculating the amount of heat generated; and (2) calculating the heat transferred by the heat source. The power generated is calculated by
where, ξi = exp(-λiΔt); Cth,i(z,t) is the density per unit length for delayed precursors of group i; and λ, β and V are of the similar definitions to those in Eq. (1).
The heat exchange between the salts in the primary and secondary loops is described by thermal convection equations, taking into account the heat emission due to the decrease of delayed thermal precursors. The flow rate D is related with volume V, temperature T, z and time t as:
where footnote 1 stands for the primary loop and 2 for the secondary loop. ρ is salt density, and cp is the heating capacity at constant pressure, S is the exchange surface, and h is a coefficient of exchange per unit surface area. The total salt flow in the primary loop is calculated by D(t) = ∑jDj,1(t), (j=1→6).
Between the core outlet and the inlet of heat exchanger, the fuel salt passes a high-temperature zone without heat transfer, namely Pipe 1 in Fig. 1; similarly, between the bottom outlet of the heat exchanger and the core inlet the fuel salt passes a low-temperature zone, namely Pipe 2. Therefore, the heat exchanger is calculated in three parts: above the heat exchanger, the heat exchanger and below the heat exchanger. Each of the parts is calculated on the assumption that temperatures of the boundary and precursor concentrations in all parts are continuous.
III. THE REALIZATION OF THE SIMULATION SOFTWARE
According to the physical model and the geometric model of the MSR, simulation code Neutron Diffusion Program for 1D (NDP1D) for one-dimensional neutron diffusion and transient analysis was finished with the following featuress:
• To analyze and design the structure of the coupled program of neutron diffusion and transient analysis. Based on physical function of the neutron diffusion process, design the data structure and file operations in a hierarchical and modular way;
• To achieve the function of each module and to complete optimization of the program algorithm and computational efficiency;
• To transplant the original code into an architecture suitable for parallelization;
• To reserve open programming interface for further improvements the transient analysis program on multi-group and complicated geometry.
NDP1D is written in C++, the object-oriented language, combined with the linear algebra numerical library ublas (part of the library Boost, launched by the C++ Standards Committee Library Working Group members), shared memory parallel library OpenMP, as well as the cross-platform graphical user interface framework tools Qt.
A. Calculation Process of NDP1D
The NDP1D is divided into critical and transient calculations. Their flow charts are shown in Figs. 2 and 3, respectively.
The salt fuel temperature along the axis z is calculated. The macroscopic cross-section is linearly interpolated to calculate the thermal characteristics of molten salt and the graphite by adjusting the radial leakage and the uranium concentration and by using the input table of reaction cross-section data. At a fixed temperature, the neutron flux equations containing a concentration of delayed neutrons are resolved, and the flux is normalized by the core power. In transient calculations, the depth and speed of control rod is considered. Flow rate, coupled with heat exchange rate, is determined by the status of reactor. And the input table of reaction cross-section data is linearly interpolated to get the macroscopic cross-section at time t.
Then the neutron flux and delayed neutron precursors concentration at the next time step are calculated by the cross-section and flux at time t. By using the distribution of temperature and the power released in the time interval to calculate temperature of the core salt, temperatures of the graphite and heat exchanger are calculated at the next moment. Considering situations where some more vectors, matrix operations, file reading and writing operations are involved in the calculation and frequently repeated calls, a certain amount optimization and parallelization have been done to adapt the multi-dimensional data simulation.
B. Optimization of NDP1D
Based on the Cinsf1D code, the program is optimized in three aspects: functionality, syntax and data structures.
• The calculation of diffusion model has been separated according to functionality into several relatively independent modules, which can be tested separatly. Coupling of the program is reduced to improve the readability, being conducive to expansion and the second revision;
• With the unique pointer characteristics in C++, the data transfer procedure of the program is optimized to avoid the operations in Cinsf1D written in Fortran, which only supports transferring the real parameters instead of formal parameters;
• The stringstream object of C++ is used in NDP1D to simplify the type and upgrade file reading and writing to the memory-level and object-level, hence the reading speed and security is greatly improved.
C. Parallelization of NDP1D
For one-dimensional case, the code can be run on just a single computer or one node of a computing cluster, so it is parallelized using shared memory model, i.e. OpenMP. The parallel execution of the model is the Fork-Join. With time analysis, physical models and code analysis as well as the OpenMP parallelization, the performance of this program has been improved and CPU load can be equalized.
1. Parallel analysis of NDP1D
Combined with the above physical model, parameters of the model and calculation involved have been analyzed:
a. The mathematical equations resolution: MATHFONC :: SOLVETSYS, MATHFONC :: SOLVES_SYSTEM are respectively the LU decomposition and the main–element elimination method, solving the linear equations AX = Y. In both methods, the former can be paralleled in the process of back substitution; the latter can be in the elimination process of the primary line to other lines. Due to the usage of lower proportion of LU decomposition algorithm, smaller size of the one-dimensional program, and the larger cost of parallelization, we have considered the main-element elimination method to get effective parallelization.
b. File reading and writing: the program inputs data from the neutron, thermo, transient, geometry and other files; for outputing codes, results are written to some data files. Since the reading and writing of files does not interfere with each other, the CPU can exchange data with the cache in parallel, and then read and write the hard disk to reduce the waiting time between peripherals. On the other hand, file reading and writing is in parallel with other calculations in order to effectively use the waiting time of reading and writing.
c. The main cycle in transient simulation: pipeline parallelism is a better choice concerned for continuous cycle of code. Because the loop of the transient simulation in main program contains a large number of interdependent variables which are in serial execution and each cycle is highly dependent on the calculation results of the last cycle, it is very difficult to parallelize and the model and algorithm modified should be sought to accommodate parallelization.
2. Results of parallelization
The CPU time of the optimized code is 37.84 s, whereas it was 51.32 s before. In comparison, the speed-up ration is S ≈ 1.36. Although it does not reach linear acceleration, p/log2p ࣘ S ࣘ p, yet it has been greatly improved. For current scale of the code, there are two trends for the parallelization results: when only a small amount of parallelization is implemented, the calculation speed gets higher, but CPU load is unbalanced, as shown in Table 1; when the parallelized part is greater, the calculation speed is not significantly increased, or even gets slower, but CPU load is more balanced. The reason is that the large serial part is unable to be parallelized. Also, the small code scale is of high cost of parallelization. A possible way to improve the application performance of the code is to increase the number of threads based on Amdahl’s law [19], or to increase workload according to Gustafson’s Law [20]. Another way of parallelization using GPU [21], specifically CUDA [22], will also be considered in future work.
| Number of simultaneously active coresa | ||
|---|---|---|
| 1 | 2 | |
| Before | 96.1% | 3.9% |
| After | 79.8% | 20.2% |
IV. THE SIMULATING RESULTS OF NDP1D
The graphical interface of NDP1D is written in Qt, a C++ cross-platform graphical user interface development framework, as shown in Fig. 4. The diagrams displayed are generated by Root codes, which are integrated into NDP1D.
An important simulating process is reactor shutdown. When the reactor is shut down, the power of the reactor goes from the level of operation down to the neutron source level. For validation, a benchmark test, the shutdown of constant power, is performed, and the results are shown in Fig. 5, with both the molten salt flow rate and temperature mainly tested. Also, the control rod depth is shown in Fig. 4. The simulation results of NDP1D program consists with the calculation results of analytical method.
V. CONCLUSION
The simulation of neutron diffusion and transient analysis of the molten salt reactor has been discussed in this paper, which established an effective physical and geometric model and built a theoretical foundation for the development of the one-dimensional neutron diffusion and transient analysis calculation software NDP1D. The optimized program has reduced the degree of coupling among modules and the redundancy, in line with the software engineering. On the other hand, the program’s performance has been improved and the hardware load has been more balanced after parallelization. With the help of simulation software NDP1D, multiple results have been acquired, such as the simulation of constant power shutdown, including changes of molten salt flow rate, temperature and control rod depth. Although the software NDP1D we developed is just one-dimensional simulation, and some conditions are set close to the ideal states of transient analysis, the studying methods and framework of the software as well as the reactor physics and thermal coupling model can be the basis for follow-up projects to further explore transient analysis on molten salt reactors with higher dimension and more complicated geometry.
Ph.D. Thesis
,