3D Real-Time Simulation of Reactor Core based on Discontinuity Factor Ratio

During the development of full scope simulator (FSS) in nuclear power plant, it is difficult to satisfy with high precision and efficiency of the 3D distribution calculation of reactor power. The accuracy of conventional coarse mesh finite difference method (CMFD) cannot meet the requirement if the nodal size is large. According to the generalized equivalence theory (GET), the finite difference solution of diffusion equation based on discontinuity factor correction was deduced, the calculation method of the discontinuity factor and boundary albedo are presented through the process of assembly homogenization. Through the calculation of typical LMW benchmark problem, the simulation results showed that using CMFD with discontinuity factor ratio (CMFD-DFR) to solve the time-dependent neutron kinetic equation, can not only improve the simulation precision, and can keep the calculation speed, and can meet the requirement of reactor core real-time model of FSS.


Introduction
During the past three decades, full scope simulators (FSS) for nuclear power plant (NPP) became very important tool used for training of personnel. In general, the transient neutron flux and power distribution in three-dimensional space can be obtained by solving the time-dependent neutron kinetic equation. The calculation accuracy and efficiency of real-time model are equally important. In spatial domain, the reactor core is generally divided into a certain number of coarse meshes, which radially each fuel assembly or quarter assembly as a node, and axially divided into a certain number of layers (including reflective layer).
The method of space-time separation of neutron flux is often used in the FSS. That is, neutron flux is decomposed into the product of amplitude function which changes rapidly with time and shape function which changes slowly with time [1]. The speed of calculation can meet the real-time requirements. But the point kinetics model is usually used in the amplitude function; the spatial variation of power at different positions in the reactor cannot be accurately represented.
With the functional requirements of the FSS of NPP shifting from a simple operator-training simulator to an engineering analytic simulator [2], advanced nodal methods for core design and safety analysis, such as the Nodal Expansion Method (NEM) and Analytic Nodal Method (ANM), have been gradually applied to solve the three-dimensional transient diffusion equations in real-time [3]. But because of the large resource consumption and long computation time, the equations must be solved by parallel computing strategy [4].
Using the coarse mesh finite difference method (CMFD) to solve the few group three dimensional time-dependent neutron kinetics equations is one of the effective methods to obtain the dynamic and static characteristics of neutron. The considerable error will be caused when the mesh spacing is large. By contrast, the computation resources will be lost large when it is solved by fine mesh. Then the realtime performance of the simulator will not meet the requirements.
The theory of discontinuity factor correction has been widely used in static diffusion calculation, especially in the area of strong absorber and the control rod tip effect [5]. The three-dimensional power distribution real-time calculation method by CMFD with discontinuity factor ratio correction (CMFD-DFR) is presented in this paper. Firstly, the heterogeneous parameter distributions of the typical condition are calculated with fine mesh. Then, the homogenous assembly group parameters (AXSs), the assembly discontinuity factors (ADFs) between the adjacent nodes and the outer boundary albedos (BAs) are calculated by the homogenization process. At last, they can be adjusted in real-time according to the thermal hydraulic parameters, such as rod position, boron concentration, fuel temperature, moderator density and so on. Transient simulation result confirm that CMFD-DFR can meet the two performances of accuracy and real-time with large mesh spacing.

Coarse Mesh Finite Difference (CMFD)
In general, the propagation and motion of neutrons are described by multidimensional neutron diffusion equation. The generation and decay of delayed neutron precursor group should be considered in the transient process. The nodal neutron balance equation for node ijk in Cartesian coordinate space can be written as follows Where, i j k is the index of a node in Cartesian coordinate space, g is the index of G energy groups; d is the index of D delayed neutron precursor groups; g  is the neutron velocity of g group; g  is the average neutron flux of g group; gx J  , gy J  , gz J  is the right and left neutron current of g group in the x, y, z directions; i x h , j y h , k z h is the width of the node in the x, y, z directions; rg  is the removing cross section of the g group, g Q is the neutron production term of the g group, d C is the density of the delayed neutron precursor of the d group, d  is the decay constant of the d group precursor, d  is the fraction of the d group precursor produced by fission, and fg v is the fission cross section of the g group. Among them, the neutron generation term is: Where, p g  and d g  is the fraction of prompt neutrons and delayed neutrons of g group, s g g '  is the cross section of g' group scatter to g group, eff k is the effective multiplication coefficient, the total fraction of delayed neutrons produced by fission, and the time coordinate t is omitted at the following description.
In the conventional CMFD method, the additional coupling equations are obtained by writing Fick's law into the approximate form of finite difference. The continuous equation of neutron net current at the interface between adjacent node i-1 and i in the x direction can be written as follows:  is the g group neutron flux at the right interface of node i-1; i gx   is the neutron flux at the left interface of node i. Because interface neutron flux is continuous, coupling relation for the average flux with respect to the net current between node i-1 and node i can be obtained as By substituting the coupling relation (5) into the nodal balance equation (1) and extending the coordinates from the x direction to the y and z directions, six supplementary equations of the node ijk are obtained. These three-dimensional diffusion equations can be solved by CMFD method. This method is suitable for the fine mesh, and the error of calculation is large for the coarse mesh.

CMFD with Discontinuity Factors Ratio correction (CMFD-DFR)
The reason for the large error in coarse mesh is that the coupling relation (5) These correction coefficients i gx f  and 1 i gx f   are called discontinuity factors (DFs). Put equation (6) into (4), we can get net current equation which is more accurate in x direction interface. The net current at the right side of node i-1 and the left side of node i can be calculated by average neutron flux of these two nodes: By the continuity of the net current, the heterogeneous interface neutron fluxes 1 From equation (8), the discontinuity factors in the coupled equation can be expressed as a ratio, so the ratio rather than the coefficient can be used to reduce the storage in the algorithm and to improve the efficiency. The discontinuity factor ratios (DFRs) to the left and right interface of node i is defined as follow, and is inversely related to adjacent nodes.
The coupling equation (8) and the definition equation (9) are incorporated into nodal equation (1), and the nodal neutron balance equation corrected by discontinuity factors ratio is obtained.
From equation (10), it can be seen that if the DFRs are given, the set of coupled equations will be solved by the implicit Euler method.  Table 1, it can be seen that the calculation time increases with the decreasing of node size and the power distribution accuracy increases with the decreasing of node size. In addition, from the T 1 and T 2 , the conventional CMFD method exhibits a very fast convergence rate when the neutron flux has a relative stable value, which is suitable for the transient iterative computation.

Homogenization and calculation of discontinuity factor
The whole core is divided into 11x11 coarse meshes in radial and 10 layers in axial. Homogenization is the process of calculating the AXSs, ADFs, and BAs of each assembly. Several fine mesh cases in Table 1 are homogenized into 11x11x10 coarse mesh respectively. The static simulation results of coarse mesh using CMFD_DFR are shown in Table 2.  As can be seen from Table 2, using the AXSs, ADFs and BAs after homogenization, the CMFD_DFR solving speed is very fast and is basically controlled at the level of milliseconds. Although the static accuracy using CMFD_DFR is slightly lower than the corresponding original fine mesh, the accuracy is greatly improved compared with conventional CMFD method. The mean square error of the power distribution is within 2%.

Coarse mesh transient calculation
In the LMW problem, the transient process is the movement of two group control rods at constant speed. According to reference [6], set total simulation time as 60s, the G1 group control rod is increased from 100cm to 180cm in 0s ~ 26.67s with 3cm/s, and the G2 group control rod is decreased from 180cm to 60cm in 7.5S ~ 47.5s, the moving speed is also 3cm/s. The whole reactor core is divided into 11x11x10 coarse nodes, and the power spatial distribution is simulated real-time by two different methods. The time step dt is 10ms, the transient curve of the average power density is recorded as shown in Figure 1. The calculation time is both about 35s. Compared CMFD and CMFD_DFR, the average power and 4 special space nodes (P1, P2, P3, P4) power and the corresponding transient reference described in literature [6] compared one by one moment, mean square error comparison results as shown in Table 3.  Table 3, it is shown that the mean square deviation of the average power density of the whole reactor core is less than 2.5%, and those of the special space nodes are less than 3% during the whole transient process. This shows the accuracy of transient calculation is greatly improved by the CMFD_DFR method, and the calculation time is basically the same.

Conclusion
This paper represents how to improve the calculation accuracy and efficiency at the same time during the development of the FSS of NPP. Based on the general equivalent theory, the solution formula of diffusion equation with correction of discontinuity factor is deduced, and the calculation method of discontinuity factor and boundary albedo is given through the homogenization process of fine mesh. The results of typical benchmark LMW problem show that the accuracy of conventional CMFD method is deficient at the coarse mesh partition. The CMFD_DFR method is applied to the real-time calculation of reactor core, which can not only achieve the corresponding accuracy, but also maintain the original calculation speed. This method has both real-time characteristics and high precision enough to meet the reactor core model requirements of nuclear power plant simulator.