CODE VERIFICATION OF MPACT USING GANAPOL CRITICAL ROD BENCHMARK

This work is dedicated to the code verification of MPACT, which is developed under the Consortium for Advanced Simulation of Light Water Reactors by the University of Michigan and Oak Ridge National Laboratory, where the numerical solution is compared to the reference solution of a benchmark problem with a known analytical solution. In this work, Benchmark Problem 3.4 in Barry Ganapol’s benchmark book was chosen as an MOC code verification test problem. Problem 3.4 is a bare cylinder of infinite height, which is an excellent benchmark problem for 2D MOC. To ensure that this benchmark problem exercised the same code as typically used by MPACT, the bare rod configuration was surrounded by a bounding box filled with a non-scattering material. To avoid implementing a critical rod search in the MPACT code, the MPACT analysis was performed using cross sections that yielded the given c, the average number of secondary neutrons per collision, and a rod radius that was the corresponding critical rod radius. MPACT agreed with all cases to within a few pcm. The convergence behavior was studied. The results show a 2nd order radial convergence, consistent with flat-source approximation. The convergence curves with respect to ray spacing and polar angle quadrature set order were also obtained. The other quantity of interest tabulated for Problem 3.4 was the radial distribution of the scalar flux. Two configurations were analyzed, and the resultant radial flux profiles agreed very well with the tabulated results. The verification of the production neutronics code MPACT has been augmented by the addition of the analytical solutions for an infinite cylinder from the Ganapol benchmark book. These test cases can be included in the regression suite for MPACT.


INTRODUCTION
The MPACT code [1] is developed under the Consortium for Advanced Simulation of Light Water Reactors (CASL) by the University of Michigan and Oak Ridge National Laboratory. MPACT uses 2D/1D methodology, where the 2D method of characteristics (MOC) is used to compute the radial solution for each plane and the 1D method is used to compute the axial leakage terms that couple the radial planes. A key component of the development and deployment effort is the verification and validation (V&V) of MPACT. Verification is an important step in the software development cycle, which includes code verification and solution verification. Code verification tests the ability of the code to converge to the correct answer, while solution verification tests its ability to converge for the intended application.
Previous work published by the authors with Method of Manufactured Solutions (MMS), which is a rigorous code verification method, and Method of Characteristics (MOC) has focused on monoenergetic cross sections and 1D plane geometry [2], [3]. MMS applied to two-group, multiple region problem was also demonstrated and documented in [4] and some other unpublished reports. Excellent results were obtained for both fixed source and eigenvalue problems, where observed error convergence rates are consistent with theory. The addition of this MMS verification problem is a significant step for code verification of MPACT. However, the assumed solutions are by definition not realistic transport solutions. Hence, an actual neutron transport benchmark problem was needed for code verification.
This work is dedicated to the code verification of the 2D MOC solver in MPACT using a benchmark problem with a known analytical solution. In this work, Benchmark Problem 3.4 in Barry Ganapol's benchmark book [5] was chosen as a MOC code verification test problem. In Section 2, a description of the benchmark problem is provided together with the MPACT benchmark model. In Section 3, the numerical results are presented. In section 4, the conclusions are drawn.

Ganapol benchmark description
Benchmark Problem 3.4 in Barry Ganapol's analytical benchmark book [5] is shown to be an excellent code verification test problem that can be modeled using MPACT without special coding. It can also be included in the MPACT regression suite. This benchmark problem is an analytic (or "semi-analytic" as Ganapol explains on page xxii) benchmark that is based on the exact solution of the singular integral equation that describes the cylindrical transport problem. This solution methodology is a complex sequence of steps that are described in detail in Ganapol's book.
The configuration is a homogeneous right circular cylinder that is infinite in height. The homogeneous cylinder has radius r R and height h o f with a total cross section t 6 and ܿ secondary neutrons per collision, where / s f t c Q 6 6 6 . The radius r is given in terms of mean free paths (mfp), so the physical radius is / t r 6 . Benchmark results are given for several cases: (a) Uniform isotropic source -the scalar flux ( ) r I is tabulated for selected values of ܿ < 1 (b) Critical rod -the critical radius is tabulated as a function of ܿ > 1 (c) Critical rod -the scalar flux ( ) r I is tabulated for critical rods as a function of ܿ > 1 The following tables are the benchmark results for the critical rod problem that were extracted from [5]. Figure 1 is  for several of the critical rod cases. The first column should be corrected such that ‫ܴ/ݎ‬ goes from 0.00 to 1.00. Table 3

MPACT benchmark problem
Since Problem 3.4 is a bare cylinder of infinite height, it is an excellent benchmark problem for 2D MOC. However, MPACT analyzes pin cells in an LWR lattice, special input processing would be needed to model the isolated cylinder. To ensure that this benchmark problem exercised the same code as typically used by MPACT, the bare rod configuration was surrounded by a bounding box filled with a non-scattering material. Vacuum boundary conditions were imposed. Since the incoming angular flux is zero at the surface of the box and stays zero until the ray intersects the rod, the cylinder surface is an effective vacuum boundary.
The critical rod cases give the critical rod size as a function of ܿ. To avoid adding special coding in MPACT to include an outer iteration to converge on the critical rod radius, MPACT solved for the eigenvalue ݇ for each of the critical rod radii given in Table 3 and ݇ = 1 is equivalent to finding the critical value of ܿ given the rod radius.

Bounding box assumption
While it seems physically plausible that the solution to the rod in a bounding box would be equivalent to solving the bare rod for the case of vacuum boundaries, there may be small effects due to the fact that the tracks incident on the rod are not uniformly and isotropically (U&I) distributed on the rod surface. They are U&I incident on the bounding box but will enter the rod with a non-U&I distribution (in azimuthal angle) depending on the size of the box relative to the rod. This should not be a practical issue for zero incoming fluxes, but this could change the numerical integration along the tracks, possibly causing disagreement with the benchmark results, perhaps in the 6 th or 7 th digit, for example. To answer this question, all the cases were rerun with bounding boxes that preserved the ratio of the bounding box side to the cylinder radius, which was ~ 17, as opposed to cases with fixed bounding box size 30 cm. The resulting ݇'s will be presented in Section 3, demonstrating that solving the critical rod problem inside a bounding box was equivalent to solving the bare rod.

NUMERICAL RESULTS
For MPACT, cases (b) and (c) were chosen for the MPACT benchmark cases. These problems were chosen because they exercise the 2D MOC solver and the eigenvalue solver. Results beyond case (b) are provided, including the convergence behavior with respect to certain discretization variables, e.g., the radial discretization, ray spacing, and polar angular discretization. Results for case (c), comparing the critical flux distribution, will follow.

Eigenvalue results
Different ܿ can be achieved by fixing the scattering cross section at 0.4 ܿ݉ ିଵ and absorption cross section at 0.6 ܿ݉ ିଵ , with a varying "fission yield cross section", which results in a varying critical rod radius.  Table 3.4.3(a). All cases were run with the following data: x Bounding box side length = 30 cm x Ray spacing = .0005 cm x # radial rings = 160 x # azimuthal slices = 32 x Quadrature set = CHEBYSHEV-GAUSS 32 24 x Convergence criterion = 1e-7 for both ݇ and ߶ It is assumed that the reference value of ݇ is known to within .01 pcm. The resulting ݇'s are all within a few pcm, demonstrating excellent agreement with the benchmark results for the infinite rod. This is a stringent code verification problem because it is an actual transport solution that is being computed by MPACT, not a manufactured solution. These cases, or a subset of them, could easily be included in the MPACT regression test suite. To validate the bounding box assumption, all the cases were rerun with bounding boxes that preserved the ratio of the bounding box side to the cylinder radius, which was ~ 17, as opposed to cases with fixed bounding box size 30 cm. The resulting ݇'s show that solving the critical rod problem inside a bounding box was equivalent to solving the bare rod.

Mesh convergence analysis
Next, we include the mesh convergence analysis for three discretization parameters, including the radial discretization, ray spacing, and polar angular discretization.

Radial convergence
We use ܿ=1.01 critical rod case as an example. The same set of discretization parameters as in the previous part are used for the radial convergence analysis.
The results are tabulated in Table 3.2. Errors are calculated by taking the difference between the eigenvalue ݇ for each individual case and unity. The last column is the assessed rate of convergence. The convergence curve is plotted in Figure 3.   Figure 3 shows a nice 2nd order radial convergence (up to 160 rings) for ܿ = 1.01. However, further study shows a larger ܿ comes with a narrower asymptotic region, e.g., when ܿ = 1.3, the 2nd order convergence can be observed only up to 30 rings. Past 30 rings, it is evident that other error components start to contaminate the 2nd order convergence, flattening out the convergence curve earlier. This is understandable since a larger ܿ comes with a smaller critical radius, where a smaller number of radial discretization is needed to achieve a certain accuracy level. If we assume a constant -0.7 pm error compensation for the other error components, the 2nd order convergence can be observed up to 100 rings.

Ray spacing convergence
We use ܿ = 1.3 critical rod case to demonstrate the solution convergence with respect to ray spacing. The following set of discretization parameters are used for the radial convergence analysis. The convergence curve is plotted in Figure 4.
We do not have a constant rate of convergence, but MOC never guarantees that. The convergence curve is generally monotonic. Local oscillation exists, which is consistent with the solution verification study in MPACT using the Progression Problem 1a. Generally, the convergence curve is nicely converging to the theoretical benchmark.

Polar angle convergence
The convergence curve is plotted in Figure 5. We do not observe a constant rate of convergence with respect to polar angular refinement. This should not be surprising since the Gaussian quadrature set does not have a constant order.

MPACT flux distribution for critical rods
The following set of discretization parameters are used for different critical rod configurations.
x Bounding box side length = 30 cm x # radial rings = 80 x # azimuthal slices = 32 x Quadrature set = CHEBYSHEV-GAUSS 32 24 x Convergence criterion = 1e-7 for both ݇ and ߶ The reference scalar fluxes against the relative radial location from [5] are plotted in Figure 5. The numerical results can be normalized such that the inner-most ring has cell-averaged scalar flux as 1, as in Figure 7. Since the benchmark results are pointwise results while the MPACT results are radially-averaged. In order to integrate the benchmark results over the radial mesh to allow a direct comparison, two solid curves in Figure 6, c=1.05 and 1.4 are digitized and fitted using a 7 th order polynomial. The polynomials introduce less than 1% error everywhere compared to the discrete data points, as shown in Figures 8. The cell-averaged critical flux and the numerical results are plotted as in Figure 9, where the Ganapol results come from the cell-averaging of the polynomial fit. The error is plotted separately in Figure 10. The cellaveraged flux agrees well with the Ganapol benchmark results, with a relative error within 1%. For c=1.4, the analysis is repeated. Again, the cell-averaged flux distribution agrees well with the Ganapol benchmark results, with a relative error within 1%.

CONCLUSIONS
Ganapol Benchmark Problem 3.4 in [5] has been used as a code verification test for MPACT. The bare rod configuration of Problem 3.4 was changed to mimic a square lattice by surrounding the rod with a bounding box with a non-scattering material. Vacuum boundary conditions were imposed on the surface of the bounding box, and several critical rod cases were analyzed with MPACT using the tabulated critical rod radius as a function of ܿ, the mean number of secondaries per collision. MPACT agreed with all cases to within a few pcm. The convergence behavior was studied. The radial rate of convergence is shown to be 2nd order, which is consistent with 2nd order flat-source approximation. The convergence curve with respect to ray spacing and polar angle quadrature set order were obtained, both of which converge to the analytic solution nicely. The other quantity of interest tabulated for Problem 3.4 was the radial distribution of the scalar flux. Two cylinders were analyzed, and the resultant radial flux profiles agreed very well with the tabulated results.
The verification of the production neutronics code MPACT has been augmented by the addition of the analytical solutions for an infinite cylinder from the Ganapol benchmark book. These test cases can be included in the regression suite for MPACT.