Comparisons among Three Coupling Forms of Cylindrical Nodal Expansion Methods

Nodal expansion method (NEM) is a typical nodal method in solving neutron diffusion equation in coarse mesh spatial discretization. NEM has been extended to cylindrical geometry in previous studies. Cylindrical NEM can be realized by different forms, where the variables in the nodal coupling equations are classified as net current, flux or partial current at the nodal surface, respectively. The three types of coupling forms of NEM have been implemented in PANGU code for the high temperature gas-cooled reactor (HTGR) physics analysis. This paper derives the three types of coupling forms of NEM, and analyzes their performance in solving the HTGR model containing void region.


INTRODUCTION
Nodal expansion method (NEM) is one of typical nodal methods [1] in solving neutron diffusion equation in coarse mesh spatial discretization. In previous studies [2] [3] [4], NEM has been extended to cylindrical geometry. There are three types of NEM forms, where the difference lies in the variables of the nodal coupling equations. The first type is the NEM with surface net current coupling (NEM-SNC), where the nodal surface net current is chosen as the coupling variable. The second type is the NEM with surface flux coupling (NEM-S ), which uses the surface flux as the coupling variable. The third type is the NEM with surface partial current (NEM-SPC) coupling, which uses the node surface partial current as the coupling variable.
PANGU code [5] is a modern computer code, developed by INET, Tsinghua university, for the pebble-bed high temperature gas-cooled reactor (HTGR) physics analysis. The original diffusion solver for wholecalculation in PANGU code is based on two dimensional (2D) r-z geometry fine-mesh finite difference (FMFD) method. In order to enhance the calculation efficiency and pave the way for the 3D whole-core simulations in the near future, a new 2D cylindrical NEM solver has recently been developed in PANGU code. Originally, NEM-SNC is chosen as the coupling form of the NEM solver. However, the presence of void cavity in HTGR limits the applicability of NEM-SNC. Therefore, the NEM-SF and NEM-SPC are included into PANGU code, to improve the computation performance in the calculation of the models, containing void regions.
The remainder of the paper is organized as follows. Section 2 gives the derivation of NEM-SNC and discusses its weakness in the treatment of void region. Section 3 gives the derivation of NEM-SF and NEM-SPC, respectively, which have advantage in treating void region. In Section 4, the computation performances of the three coupling forms are compared by the numerical tests of HTGR models. The conclusion remarks are given in Section 5.

Derivation of NEM-SNC in r-z cylindrical geometry
The 2D diffusion equation in a cylindrical-geometry node can be written as: where ( , ) m g r z I and ( , ) m g Q r z is the flux and source density in the node m. Most of nodal methods adopt the transverse integration technique [6], which reduces the multi-dimensional diffusion equation to a set of one dimensional equations: where N and K is the expansion order of the partial flux and leakage term, respectively. The basis function set proposed by Bandini [7] has been widely used in Cartesian-geometry NEM, however, it is not applicable in cylindrical geometry. PANGU code employs an improved basis function set proposed by Wang D et al., [2] to treat the intra-node flux distribution in cylindrical geometry, and N and K are respectively set as 4 and 2.
The basic equation set of NEM is mainly constituted by neutron balance equation, nodal surface net current continuity equation and nodal surface flux continuity equation. The neutron balance equation is as follows: Apart from the aforementioned three basic conditions, in NEM, the residual weighted procedure [9] is generally introduced as extra conditions to eliminate the higher order expansion coefficients of the partial flux  Then substituting Eq. (5) into Eq. (8), the neutron balance equation can be written: The surface flux continuity equation The zero order transverse leakage expansion coefficient (the averaged transverse leakage): As such, Eqs. (11), (12), (13) and (9) constitute the NEM-SNC's coupling equation set.

The weakness of NEM-SNC in the treatment of void region
The equation set of the NEM-SNC can be solved by various iterative procedures, such as Gauss-Seidel (GS) method and Generalized minimal residual (GMRES) method. GMRES method always requires precondition to accelerate the convergence of the iteration. The incomplete LU (ILU) precondition is most widely implemented in GMRES method.
In pebble bed HTGRs, there is a large cylindrical cavity between the pebble bed and the top reflector. PANGU code treats the cavity by using an equivalent diffusion coefficient derived from transport theory, while the cross sections of the area are zero [10]. However, the presence of void regions causes the node averaged flux term in neutron balance equation Eq. (11) to vanish. This will lead to zero diagonal elements of the coefficient matrix in NEM-SNC, which limits the usability of iterative methods. Firstly, the GS method cannot be applied, because the iterative procedure will break off in the row, where the matrix diagonal element is zero. Secondly, incomplete LU decomposition cannot be implemented for the matrix, containing zero diagonal elements, which means GMRES method cannot be combined with ILU precondition. The other two forms of NEM, i.e. NEM-SF and NEM-SPC, are able to overcome the difficulty, caused by the void region.

Derivation and analysis of NEM-SF and NEM-SPC
Considering the weakness of NEM-SNC in the treatment of void region mentioned above, NEM-SF and NEM-SPC have been derived. The neutron balance equation of the two coupling forms can avoid the vanishing of averaged flux term in neutron balance equation.

The NEM with surface flux coupling (NEM-SF)
NEM-SF chooses the node averaged flux g I and surface flux u m g I r as the variables to be solved [9]. Then substituting Eq. (5) into Eq. (8), the neutron balance equation is given by the following:

The NEM with surface partial current coupling (NEM-SPC)
Based on Fick's law, the relation between the surface partial current Attention should be paid to the ordering of variables and the construction of the coefficient matrix in NEM-SPC. The red-black ordering is utilized to improve the convergence of the iterative methods, taking into account the fact that, the surface incoming currents of each node are equal to the outgoing currents of the adjacent nodes. The nodes are divided into two types, red and black. The ordering of the variables to be solved can be described by Figure 1.

Description of test cases
In order to test the performance of the three forms of NEM solvers in PANGU code, a simplified 2D pebblebed HTGR model is established for full-core diffusion calculation, as is shown in Fig. 1

Numerical results
Numerical results of the two cases are presented in Table I and Table II, respectively. It can be seen from the results that the three NEM forms are all able to solve the diffusion equation in cylindrical geometry for HTGR. The computation efficiency of the NEM solvers are generally superior to that of the FMFD solver.
The keff results of the three forms of NEM solvers are very close, which is due to the fact that the three coupling forms are mathematically equivalent.
In Case 1, when there is no cavity, all forms of NEM solvers show excellent accuracy and efficiency compared with the FMFD solver. In Case 2, it is found that the computational time of both FMFD and NEM solvers are increased. It indicates that the performance of diffusion solvers deteriorates to some extent when there exist void regions. In specific, NEM-SNC with GMRES iteration method gives very low efficiency in Case 2, because the ILU precondition is not applicable in this case.
By our tests, it is found that the NEM-SF and NEM-SPC both provide remarkable efficiency in Case 1 and Case 2. NEM-SPC with GS iteration method gives the best performance among all of the NEM solvers.

Conclusion
This paper derives the three types of coupling forms in nodal expansion method in cylindrical geometry. The three coupling forms are mathematically equivalent, but show different performance in numerical iterative solution. Among the three coupling forms, the NEM-SNC has ill conditioned coefficient matrix when the calculated model contains void region. This feature limits the application of some iterative method or precondition technique in solving the coupling equation, thus making the solving process less efficient. NEM-SF and NEM-SPC can both overcome the drawback of NEM-SNC, and thus show better performance in treating void region. Preliminary numerical tests suggest that NEM-SPC combined with the GS method show the best performance in this work, though NEM-SPC consumes more memory than the two other forms. That may be ascribed to red-black ordering's improvement on the convergence of GS method.