IMPLEMENTATION OF QUADRATIC AXIAL TRIAL FUNCTIONS IN THE HIGH-FIDELITY TRANSPORT CODE PROTEUS-MOC

PROTEUS-MOC is a pin-resolved high-fidelity transport code, in which the axial variation of angular flux is represented in terms of orthogonal polynomials. Currently, PROTEUS-MOC employs linear functions and requires relatively fine axial meshes to achieve high accuracy, which increases the number of axial meshes and hence the memory requirement. In this study, aiming to reduce the memory requirement and potentially the computational time by allowing larger axial meshes, we have extended the PROTEUS-MOC transport solution method to quadratic trial functions. Preliminary tests for the performance of quadratic trial functions have been performed using the 3-D C5G7 benchmark problem. Test results showed that for the same axial mesh configuration with relatively large sizes, the quadratic approximation yields about 2 to 5 times more accurate pin powers than the linear approximation, depending on the degree of axial variation of angular fluxes. The quadratic approximation also allows the use of about 3 times coarser axial meshes than the linear approximation for comparable pin power accuracy, which consequently reduces the memory requirement by about 2 times. The memory reduction is not proportional because of the increased number of coefficients in each element from 2 to 3. However, the quadratic approximation did not reduce the computational time as expected because of the deteriorated performance of the pCMFD acceleration scheme due to large axial mesh sizes.


INTRODUCTION
Three-dimensional (3-D) heterogeneous whole-core deterministic transport calculation is gaining increasing interest due to its capability to provide high-fidelity, pin-resolved flux solutions. The 2-D/1-D methods were first developed with successful applications to pressurized water reactor (PWR) problems [1,2]. In these methods, the 3-D problem domain is divided axially into multiple slices and the 3-D transport equation is converted into a system of 2-D transport equations coupled through the axial transverse leakages at the bottom and the top of each slice. The 2-D transport equation for each slice is solved using the 2-D method of characteristics (MOC) with assumed axial transverse leakages. The axial 1-D transport equations are typically solved with a lower fidelity method such as the diffusion theory nodal expansion method (NEM) or the simplified P 3 (SP3) NEM method [3]. The initial 2-D/1-D method showed numerical instabilities when applied to problems with small axial meshes, strong axial heterogeneities, or very low-density regions [5]. To eliminate the instability problem, an alternative 2-D/1-D approximation that preserves the exact transport physics in the radial directions and the approximate diffusion theory in the axial direction was proposed with a stable iteration scheme [4]. Full 3-D MOC codes with 3-D ray tracing techniques [5,6] were also developed but computationally expensive to be applied to large practical core problems.
An alternative 3-D transport method was developed by combining the discontinuous Galerkin weighted residual method with the 2-D MOC and implemented in the PROTEUS-MOC code [7]. PROTEUS-MOC solves the first-order Boltzmann transport equation for axially extruded geometries of 2-D unstructured meshes by combining MOC for the radial direction and the discontinuous Galerkin finite element method (DFEM) in the axial direction. It provides a fully consistent and accurate solution of the 3-D transport equation by overcoming the limitations of the 2-D/1-D method without any noticeable increase in computational time. A good computational efficiency is achieved through the single-and two-level consistent partial current-based coarse-mesh finite difference (pCMFD) acceleration schemes [8].
In PROTEUS-MOC, the axial variation of angular flux in each element is represented in terms of orthogonal polynomials. By the use of the Galerkin weighted residual technique, a system of 2-D partial differential equations for the expansion coefficients is obtained for each axial slice. This system of equations has a similar form to the 2-D transport equation, and thus it can be solved using a 2-D MOC solver with exponentials of matrices instead of numbers. However, the equations for two neighboring axial slices are coupled explicitly through the expansion coefficient vector at the interface. Therefore, they are solved sequentially from the bottom to the top of the problem for a positive polar direction and from the top to the bottom for a negative polar direction. Currently, the PROTEUS-MOC code uses linear trial functions and requires relatively fine axial meshes to achieve high accuracy, which increases the number of axial slices and hence the memory requirement and computational time. In this study, aiming to reduce the memory requirement and potentially the computational time by allowing larger axial meshes, we have extended the PROTEUS-MOC transport solution method to quadratic trial functions.
The purpose of this paper is to present the mathematical formulations of the quadratic expansion method and the numerical test results. The rest of this paper is organized as follows. In Section 2, the transport solution methods of PROTEUS-MOC with the linear and quadratic trial functions are discussed. Section 3 presents the numerical test results for the 3-D C5G7 benchmark problem [9]. Conclusions are given in Section 4.

METHODOLOGIES OF PROTEUS-MOC
In the 3-D transport solver of PROTEUS-MOC, the axial variation of the angular flux in each element is represented in terms of orthogonal polynomials. Omitting the element index for simplicity, the angular flux and source in the m-th direction in axial slice k can be represented as where J is the expansion order, ( ) where I is the J×J identity matrix, s is the segment length of the characteristic line, and m θ is the polar angle of the m-th direction. The solution in Eq. (4)

Basis Functions and Stiff Matrices of Linear Approximation
For the linear approximation, the basis functions where cos

Basis Functions and Stiff Matrices of Quadratic Approximation
For the quadratic approximation, the basis functions { 1, 2,3} ( ), and Eq. (6) and the following quadratic function: The explicit expression of the stiffness matrix , k m Σ for axial slice k and angular direction m is given by , ,

NUMERICAL EXCERCISES
The quadratic basis and trial functions have been implemented into PROTEUS-MOC. The performance and effectiveness of the quadratic approximation have been tested using the 3-D C5G7 benchmark problem [8]. As shown in Fig. 1, the 3-D C5G7 benchmark problem consists of two MOX and two UO2 assemblies surrounded by water reflectors. Each fuel assembly consists of 264 fuel pins, 24 guide tubes, and 1 fission chamber. The active core is 42.84 cm high, and an additional 21.42 cm reflector is added on the top. Reflective and vacuum boundary conditions are imposed on the bottom and top faces, respectively. Detailed descriptions can be found in the C5G7 benchmark report [9].

Figure 1. 3-D C5G7 Problem: (a) Radial Configuration, (b) Pin Arrangements in Four Assemblies, and (c) Axial Configuration.
Using the PROTEUS-MOC code, the eigenvalue and pin powers were calculated with the linear and quadratic approximations. Each fuel pin cell was discretized into 160 elements by dividing the fuel and moderator zones into 3 and 2 annular rings, respectively, and 32 azimuthal sectors. Each guide tube was discretized into 224 elements by dividing the absorber region into 5 annular rings to capture the steep flux gradient. The angular domain was discretized using 8 polar angles for π and 48 azimuthal angles for 2π. A two-norm residual error criterion of 1.0e-7 was used for convergence of the eigenvalue, fission source, and flux. In order to compare the performance of the quadratic approximation relative to the linear approximation, two different axial mesh configurations were prepared. In the first one, the core and reflector regions were divided into 18 and 12 meshes, respectively, so that the linear approximation yields a maximum pin power error less than 2%. The second one was prepared by coarsening the meshes by a factor of three so that the quadratic approximation would yield a similar accuracy to the linear approximation solution with the first axial mesh set.
The calculation results were compared with the reference MCNP solution [9]. Pin powers were calculated in each of three axial segments shown in Fig. 1(c) as in the reference solution. The solution accuracy was measured by the following four distinct errors: the maximum error (MAX), the average error (AVG), the root mean square (RMS) error, and the mean relative error (MRE), which are defined by where N is the number of fuel pins, n P is the n-th fuel pin power, n e is the percentage error of the n-th fuel pin, and avg P is the average pin power. Reflector reference MCNP eigenvalue was 1.143080. The linear approximation solution with 30 axial meshes yields an eigenvalue error of -0.0038% and a maximum pin-segment power error of 1.7%, which occurs in the segment 3 next to the reflector. The reduction of the number of axial meshes from 30 to 10 increases the eigenvalue error of the linear approximation solution from -0.0038% to 0.035% and the maximum pin segment power error from 1.7% to 4.0%. It can be seen that the accuracy of the quadratic approximation solution for 10 axial meshes is comparable to that of the linear approximation solution for 30 axial meshes. The eigenvalue error of the quadratic solution is only about 1 pcm compared to reference result. The maximum pin-segment power error is around 1% in the axial segments 1 and 2. The axial segment 3 next to the reflector shows slightly higher pin-segment power errors, yielding a maximum pinsegment power error of 1.7%. The AVG and RMS errors are much smaller than the MAX error, indicating that only a few pin-segment power errors are larger than 1% and the other pin-segment power errors are significantly smaller than 1%. The MRE error of segment 3 is 0.08% and comparable to that of segment 1 or 2, suggesting that all the relatively large pin-segment power errors in segment 3 occurred in the fuel pins of relatively small powers. These results indicate that the quadratic approximation allows using about 3 times coarser axial meshes than the linear approximation for comparable pin-segment power accuracy. With the about 3 times coarser axial meshes, the quadratic approximation reduces the memory requirement by about 2 times. By comparing the linear and quadratic approximation solutions for 10 axial meshes, it can be seen that relative to the linear approximation, the quadratic approximation reduces the eigenvalue error from 0.035% to 0.0014% and the maximum pin-segment power error from 4.0% to 1.7%. In the segment 1 and 2, where the neutron flux exhibits relatively smooth variations, the quadratic approximation yielded about 2 ~ 3 times more accurate results in terms of the AVG, RMS and MRE errors. In the segment 3 where the neutron flux shows steeper axial variations, the quadratic approximation improved the pin power accuracies by about 5 times compared to the linear approximation. Table II compares the computational performance of the linear approximation for 30 axial meshes and the quadratic approximation for 10 axial meshes. The same number of processors was used for both cases, but with different domain decomposition configurations because of the different axial mesh numbers. It can be seen that the quadratic approximation reduced the time for a single transport sweeping operation by 39%. The time reduction for a single transport sweep is smaller than the reduction of the axial mesh numbers (66%) due to the additional overhead of the floating-point operations in matrix inversion and exponential calculations. However, compared to the linear approximation, the quadratic approximation increased the number of outer iterations by 80% as shown in Figure 2, and the total number of transport sweeps by 68%. Consequently, the quadratic axial expansion method did not reduce the overall computational time as expected. The increased number of outer iterations of the quadratic approximation suggests that the larger axial meshes in the quadratic approximation deteriorate the performance of the pCMFD acceleration because the pCMFD acceleration scheme updates only the mesh-averaged flux values. In order to accelerate the iterative solution of the quadratic approximation effectively, the linear and quadratic expansion coefficients need to be updated properly, especially when large axial meshes are used.

CONCLUSIONS
The potential benefits of higher order axial trial functions in the discontinuous Galerkin finite element method of PROTEUS-MOC due to increased mesh size has been investigated. The trial function space of PROTEUS-MOC was extended from the current linear polynomial to the quadratic polynomial. Preliminary tests for the performance of the quadratic approximation were performed using the C5G7 benchmark problem. The solution accuracy was examined by comparing the calculated eigenvalue and pin powers in three axial segments with the reference Monte Carlo solutions. Numerical results showed that compared to the linear approximation, the quadratic approximation can reduce the memory requirement by a factor of two by allowing coarser axial meshes while maintaining the same accuracy. This makes the quadratic approximation an attractive option for practical large core problems. However, the quadratic approximation with coarser axial meshes did not reduce the computational time as expected. This is believed to be due to the deteriorated performance of the pCMFD acceleration method for large axial meshes. More effective acceleration schemes that can update the higher order expansion coefficients will be investigated in future.