A SIMPLIFIED TWO-NODE COARSE-MESH FINITE DIFFERENCE METHOD FOR PIN-WISE CALCULATION WITH SP3

For accurate and efficient pin-by-pin core calculation of SP3 equations, a simplified twonode Coarse Mesh Finite Difference (CMFD) method with the nonlinear iterative strategy is proposed. In this study, the two-node method is only used for discretization of Laplace operator of the 0th moment in the first equation, while the fine mesh finite difference (FMFD) is used for the 2nd moment flux and the second equation. In the two-node problem, transverse flux is expanded to second-order Legendre polynomials. In addition, the associated transverse leakage is approximated with flat distribution. Then the current coupling coefficients are updated in nonlinear iterations. The generalized eigenvalue problem from CMFD is solved using Jacobi-Davidson method. A protype code CORCA-PIN is developed. FMFD scheme is implemented in CORCA-PIN as well. The 2D KAIST 3A benchmark problem and extended 3D problem, which are cell homogenized problems with strong absorber, are tested. Numerical results show that the solution of the simplified twonode method with 1×1 mesh per cell has comparable accuracy of FMFD with 4×4 meshes per cell, but cost less time. The method is suitable for whole core pin-wise calculation.


INTRODUCTION
With the design and calculation requirements of new and complex reactor types, such as MOX fuel, BWR, SCWR, and so on, traditional two-step framework based on single assembly homogenization is not accurate to capture the rapid changes in space and energy spectrum between adjacent assemblies and inner assembly. Compared to assembly homogenization, cell homogenization is expected to capture these rapid variations of flux. With the development of computer capabilities, pin-by-pin calculation has gained more research attention.
For a typical pressurized water reactor, the number of unknowns is more than 10 million, which is increased by three orders of magnitude compared to the traditional two-step method. Therefore, the efficiency of solving pin-by-pin problem is the key to the practical application. Several methods [1][2][3][4][5][6], EPJ Web of Conferences 247, 02023 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124702023 such as the response matrix method, finite element method, global/local methods, and so on, have been proposed to improve high efficiency.
In this paper, a simplified two-node method based on nonlinear iterative strategy is used to solve the SP3 transport equation. The method is described in section 2. The numerical results are shown in section 3. In section 4, we will give the conclusions.

SP3 Equations
The SP3 equations used in this paper is shown as follows, , ߮ is 0 th moment of angular flux ߮ ଶ is 2 nd moment, the subscript g means energy group, the rest are customary symbols.
The old form of SP3 equations is used in this paper, in which ߮ and ߮ ଶ considering associated discontinuity factors are continuous at the interface between adjacent meshes and ‫ܦ−‬ ∇߮ and ‫ܦ−‬ ଶ ∇߮ ଶ are continuous as well. The rigorous form of SP3 equations [7] with correct interface continuity conditions will be considered in future.

Two-node Method
In two-node method, the current at the interface in two-node method is written as follows, where h is mesh size, f is discontinuity factor, ‫ܦ‬ ା is current coupling coefficient which is updated in nonlinear iteration. If current coupling coefficient are equal to zero, the method degraded into fine mesh finite difference method (FMFD).
Then the coarse mesh finite difference (CMFD) discretization systems are established by substituting Eq.
(3) into Eq. (1) and Eq. (2), which is written as follows, The structures of matrixes A and B are the same as ones of FMFD.
Since the 2 nd moment flux is one order of magnitude less than 0 th moment flux, only the 1 st moment currents at interface in Eq. (1) are treated with nonzero current coupling coefficients. The rest of currents in Eq. (1) and Eq. (2) are treated by FMFD. On the other hand, more strictly criterion on flux convergence would be required to ensure the convergence of 3 rd moment currents coupling coefficients. Therefore a simplified two-node method is used in this paper. To update the current coupling coefficients, the transverse integration is employed to Eq. (1), and the transverse flux along x-direction satisfies where L is transverse leakage term.
The transverse leakage term L is approximated using flat distribution as follows The transverse 0 th moment is expanded to second-order Legendre polynomials satisfying where Pi is the i th Legendre polynomials.
According to orthogonal properties of Legendre polynomials, b0 equals to mesh average flux which is solved from Eq. (4) in the previous non-linear iteration. Therefore, there are four unknows, b1 and b2 of adjacent meshes, i and i+1, for each group. Multiplying the two sides of Eq. (5) by P0 and then integrating it over each mesh, we obtain two equations for each group. The continuity conditions of flux and current are used to get the rest two equations. Then four unknows are determined by these equations. Then currents ‫ܬ‬ ା ቀ ଶ ቁ at the interface are calculated using Eq. (7). Finally, the coupling coefficients are updated as follows, The iteration strategy of the coupling coefficients is the standard non-linear iteration strategy shown below, And we set the convergence criterion is equal to 0.001 in this paper.

Iteration Method of CMFD
The mainly computing burden is to solve Eq. (4). As well known, this k-eigenvalue problem is a generalized non-Hermitian eigenvalue problem in which the largest eigenvalue in real need to be solved. Jacobi-Davidson method, which a kind of subspace method, is used to solve this eigenvalue problem. The subspace is expanded by the approximate solution of the Jacobi orthogonal correction equation as follows, where ‫ݎ‬ = ‫ݑܣ‬ − ‫ݑܤߠ‬ is residual, t is correction, K is preconditioner of ‫ܣ‬ − ‫,ܤߠ‬ u is approximation eigenvector and ߠ is associated approximation eigenvalue. Since Eq. (9) is required to be solved approximately, Jacobi-Davidson method is more efficient than the Arnoldi method to solve generalized eigenvalue problem. In this paper, GMRES iterative solver with Jacobi preconditioner is used and the relative tolerance is set to 1E-3.
The iterative methods are implemented based on SLEPc v3.9 [8] and PETSc v3.9 [9] which are famous open source libraries providing a large suite of parallel eigenvalue and linear system solvers. Then CORCA-PIN code is developed.

2D 4 Groups ARI Benchmark problem
The original problem [10] is a three-dimensional 1/4 symmetry ARI core loaded with 157 assemblies with 4 groups cell homogenized cross sections. In order to obtain the mesh-independent solution with less calculation, we use two-dimensional geometry, and the reflective boundary conditions are used in the axial directions. The homogenized cells are divided into 1×1, 2×2, 4×4 meshes respectively. The reference solution is calculated by FMFD using 32×32 meshes per cell, and the mesh size is about 0.04 cm, which ensures that the mesh discretization error of the reference solution is negligible.
The Keff and power distribution results are given in Table 1 and Fig. 1. The error of Keff by two-node method is less than 202 pcm. Keff calculated by two-node method using 1×1 mesh is accurate than that of FMFD using 4×4 meshes. The maximum relative error of the power distribution of two-node method using 1×1 mesh per cell is 3.3%, 2×2 is 2.1%, and 4×4 is 0.6%. The power errors calculated by two-node method using 2×2 meshes per cell are less than that of FMFD using 4×4 significantly. And the accuracy of two-node method using 1×1 mesh per cell is similar with that of FMFD using 4×4 except on a few cells.
FMFD is of first-order accuracy for nonuniform material case, which explains the performance on this problem. It occurred frequently that the reactivity errors of two-node method using 1×1 mesh per cell is less than ones using 2×2 mesh per cell. But calculation with 2×2 mesh per cell have higher accuracy on flux and power distribution than calulations with 1×1 mesh per cell. "FD" means finite difference, "TN" means two-node "n by n" means mesh division of n n meshes per pin is used

2D KAIST 3A Benchmark problem
The 2D KAIST 3A benchmark problem [11] provides 7-group cell homogenized cross section, including UOX fuel rods, MOX fuel rods, control rods, and poison rods. The core layout is shown in Fig. 2. The discretization errors of two-node method are validated using different meshes and compared with FMFD. The homogenized cells are divided into 1×1, 2×2, 4×4 meshes respectively. The reference solution is calculated by FMFD using 32×32 meshes per cell.

Figure 2. Layout of KAIST 3A Benchmark Problem
The Keff comparison is given in Table 2. The error of Keff by two-node method is less than 100 pcm. Keff calculated by two-node method using 1×1 mesh is more accurate than that of FMFD using 2×2 meshes, even 4×4 meshes. In addition, the maximum relative error of the power distribution of two-node method using 1×1 mesh per cell is 3.9%, 2×2 is 2.6%, and 4×4 is 0.7%. The relative error of pin-wise power distribution is shown in Fig. 3. The comparison results of relative errors are similar to the previous section.   There are six types assembly in KAIST 3A benchmark problem. We use them to construct a threedimensional 1/4 reactor shown in Fig. 4. The active zone of the reactor is 365.76 cm high, and both the upper and lower reflection layers are 25.4 cm high. There are four kinds of grid, including 1.26cm×1.26cm×2.54cm, 0.63cm×0.63cm×2.54cm, 0.315cm×0.315cm×2.54 cm, and 0.315cm×0.315cm×1.27cm. The reference solution is calculated by two-node method using the finest grid. It is noted that the grid with 1.26cm×1.26cm mesh size is the same as the grid of 1×1 mesh per cell shown in section 3.1. For the coarsest mesh of 1.26cm×1.26cm×2.54cm, the number of grids is about 3 million, and the number of unknows is about 42 million.
The Keff comparison is shown in Table 3. The errors of Keff of two-node method are less than 20 pcm. Keff calculated by two-node method using 1.26cm×1.26cm mesh size is more accurate than that by FMFD using 0.63cm×0.63cm mesh size, even using 0.315cm×0.315cm mesh size.
Since the difference of Keff between axial mesh size of 1.27cm and 2.54cm can be ignored, axial mesh size of 2.54cm is sufficiently accurate in practical calculation. Therefore, from the perspective of engineering applications, two-node method using 1.26cm×1.26cm×2.54cm mesh size meets the accuracy requirements of practical engineering simulations.

CONCLUSIONS
In this paper, a simplified two-node method is proposed to discretize the SP3 equations, and the CORCA-PIN code is developed. The numerical results show that: 1) The calculation accuracy of two-node method using n×n meshes per cell is better than that of FMFD using 2n×2n meshes per cell significantly.
2) The efficiency of two-node method using 1×1 mesh per cell increases by one order of magnitude compared to that of FMFD using 4×4 meshes per cell. 3) Relative errors of pin power of two-node method using 1×1 mesh per cell are within 5%, errors of Keff are within 200 pcm. 4) It takes ~9 core hours to perform a pinby-pin 7 groups calculation of a practical three-dimensional 1/4 reactor using 1×1 mesh per cell.