CONTINUOUS ENERGY COMET SOLUTION TO A SMALL MODULAR ADVANCED HIGH-TEMPERATURE REACTOR BENCHMARK PROBLEM (SmAHTR)

The continuous energy coarse mesh transport (COMET) method is a hybrid stochasticdeterministic solver that provides transport solutions to heterogeneous reactor cores. In this paper, COMET is tested against continuous energy Monte Carlo in solving the recently developed stylized Small Modular Advanced High-Temperature Reactor (SmAHTR) Benchmark Problems based on the Oak Ridge National Laboratory pre-conceptual design (core configurations). These problems are well-suited to test the performance of advanced neutronics tools because of their unique neutronics characteristics such as the multiple heterogeneities. The COMET solutions for the three benchmark problems were found to agree very well with the continuous energy Monte Carlo reference solutions. The discrepancy in the core eigenvalue (k-eff) varied from 40 pcm to 51 pcm. The average and maximum relative differences in the pin fission densities were in the range of 0.20% to 0.21% and 0.77% to 0.94%, respectively. It was also found that COMET was more than 2,000 times fast than MCNP. It can be concluded that COMET can model the SmAHTR core configuration with high fidelity and significantly high computational efficiency.


INTRODUCTION
The coarse mesh transport code COMET [1][2][3] is a hybrid stochastic-deterministic solver that provides transport solution to heterogeneous reactor core problems. The method consists of local stochastic transport a priori calculations to generate a response function library for global deterministic calculations to obtain the core solution. COMET is well-suited to model advanced reactors since material heterogeneities down to the fuel and burnable absorber particle levels are fully modeled in COMET response function generation. The benchmarking calculations [4,5] against Monte Caro for the Advanced High-Temperature Reactor (AHTR) and VHTR benchmark problems have shown that it is highly accurate with significantly faster computational speed than Monte Carlo. In this work, a set of stylized Small Modular Advanced High-Temperature Reactor (SmAHTR) problems [6] was used to evaluate the performance of COMET by comparing its results with the corresponding continuous-energy Monte Carlo solutions. These benchmark problems are ideal to evaluate the COMET performance because of their unique neutronics characteristics `such as strong leakage due to the small size of the modular design and high heterogeneity resulting from the use of TRISO fuel particles as well as burnable absorber particles.
The rest of the paper is organized as follows. The COMET method is briefly described in section 2 and the benchmark problems are summarized in section 3. The COMET calculations and the comparison with the Monte Carlo reference solutions are presented in section 4. The concluding remarks are found in section 5.

CONTINUOUS ENERGY COMET METHOD
For the convenience of the reader, the COMET method is briefly described this section. The detailed theory can be found in [1][2][3]. The method starts with the decomposition of the spatial domain of the problem of interest into a number of non-overlapping regions (referred as coarse meshes) where = ⋃ . It is assumed that the incoming/outgoing flux on the bounding surfaces of each coarse mesh can be expanded in terms of a set of expansion functions.
where ± ( ⃗, Ω , ) represent the incoming ("+") or outgoing ("-") flux on surface of coarse mesh , denoted as , at position ⃗, direction Ω and energy , ( ) stands for the asymptotic spectrum [3], Γ represent the orthogonal expansion functions, and ±, are the corresponding incoming or outgoing partial current expansion coefficients (moments). Generally, the tensor product of continuous orthogonal (Legendre, Chebyshev, etc.) polynomials are used to expand the spatial and angular distribution of the angular flux on the interfaces between coarse meshes. The energy dependence is decomposed into two components: a fast-varying component representing the asymptotic spectrum ( ) for a typical single assembly problem(s), and a slow varying (shape) component that is expanded in terms of the 0 th order Bsplines with coarse energy bins.
Because of the linearity of the Boltzmann equation, the angular flux distribution in coarse mesh can be formulated as the superposition of the contribution due to each incoming flux moment incident on its boundary as: where the surface-to-volume response function is the transport solution to the following local fixedsource problem: with the boundary condition In the above equations, H and F represent the standard transport and fission operators. is the core eigenvalue, and denotes the outward normal at surface .
Similarly, the outgoing partial current moments can also be computed as the superposition of the response to each incoming flux imposed on the coarse mesh boundary.
where the surface-to-surface response coefficients have the following relation with the surface-tovolume response functions ⃗, Ω , ; .
As it can be seen from Equations (3) and (4), the response functions are needed only for unique coarse meshes since they are a function of the property (geometry and material composition) of coarse meshes.
Once the expansion functions are chosen, the continuous energy stochastic response function generator is first used to solve a set of local fixed-source problems to obtain the response coefficients for all unique coarse meshes. The deterministic COMET core solver is then used to obtain the global solution in a twolevel iteration process. In the inner iterations, the local current balance equation (Eq. 5) is used to converge on the partial currents crossing coarse meshes. In the outer iteration, the particle balance equation over the entire core is used to repeatedly update the core eigenvalue. When the partial current moments and core eigenvalue are converged, the flux (or pin fission density) distribution can be constructed by Equation (2) via the pre-computed response function library.

2D SMAHTR BENCHMARK PROBLEM
The detailed description of the SmAHTR benchmark problem can be found in a companion paper [6]. A brief description of these problems is presented in this section for readers' convenience.
(a) radial core (b) fuel assembly Figure 1. SmAHTR core and fuel assembly configurations.
The 2D SmAHRT core shown in Figure 1a consists of 19 hexagonal fuel assemblies and a surrounding radial graphite block with 18 cylindrical coolant channels. The outer radius of the graphite reflector is 150 cm. Note, that the circular/cylindrical external boundary of the problem is where the graphite reflector ends and the reactor vessel begins (i.e., the vessel is replaced by vacuum). The fuel assembly shown in Figure  1b is comprised of 66 fuel pins (dark blue) filled with TRISO UCO fuel at a 50% volumetric packing density, 6 burnable absorber (BA) pins (grey), 12 graphite pins (light blue), a cluster of six control rods (yellow), and one reserved shut down rod/pin (green), surrounding coolant FLiBe (magenta) and graphite (light blue). For simplicity, both TRISO fuel and BA particles are modeled as a body-centered hexagonal lattice.
The three core configurations are considered in this paper:  All rods out (ARO): all control rods in 19 assemblies are withdrawn;  All rods in (ARI): all control rods in 19 assemblies are fully inserted;  Critical (CRT) core: It is a core configuration with an eigenvalue close to 1. In this configuration, all control rods in the six outer corner assemblies are withdrawn, while all control rods in the other assemblies are fully inserted.

NUMERICAL RESULTS
Both the continuous energy COMET and MCNP [7] were used to model the three 2D core configurations. The eigenvalues and pin fission density (PFD) distributions predicted by COMET were compared with the MCNP reference solutions. Given the large number of fuel pins, the following four statistical parameters are used to compare the agreement/disagreement between the PFDs computed by COMET and the reference (MCNP) solutions.  AVE: the average relative difference between the COMET and reference PFDs;  MRE: the mean (pin power weighted average) relative difference between the COMET and reference PFDs;  RMS: the root mean square of the relative differences between the COMET and reference PFDs;  MAX: the maximum relative differences between the COMET and reference PFDs; Similarly, those parameters are also used to summarize the uncertainties in the MCNP PFDs.

MCNP Calculations
MCNP was used to model the three (ARO, ARI, CRT) full core configurations to obtain their continuousenergy reference solutions. The fission source was first converged by tracking 40 million particles (1,000 cycles with 40,000 particles per cycle), and it took about 8 hours on a 56-core cluster. Then additional 60 million particle histories (300 inactive and 1,200 active cycles with 40,000 particles per cycles) were followed to compute the eigenvalue and tally the pin fission density distribution in all the fuel pins for each core configuration, and each calculation took about 16.5 hours on the same cluster.

COMET Calculations
A total of eight unique coarse meshes were used by COMET to model the 2D SmAHTR benchmark problems. They include three (uncontrolled without BA, uncontrolled with BA and controlled with BA) fuel assembly coarse meshes shown in Figure 1b, one full reflector block (Figure 2a), and four partial reflector blocks shown in Figs. 2b-2e. In those partial reflector blocks, the white region is filled with infinitely black absorber. This allows one to match the actual geometry (circular core boundary containing partial hexagonal reflector block) while retaining consistent hexagonal coarse meshes. The dimension and material property of each unique coarse mesh are shown in Table I. Figure 2. Configuration of five reflector coarse meshes. The asymptotic spectrum was first obtained by a continuous energy MCNP single assembly calculation for coarse mesh 2 (Table I)

Comparison
The comparison of the COMET and MCNP eigenvalues as well as their computation time is presented in Table II. 0.39 * MCNP calculations were performed on a 56-core cluster. ** COMET calculations were performed on a single core of the cluster.
It was found that the eigenvalues predicted by COMET agree very well with those computed by MCNP, with a difference of 40 to 51 pcm. Not counting the run time for source convergence, each MCNP calculation took about 16.5 hours on a 56-core cluster, while the COMET core calculation took about 0.39 to 0.47 minutes on a single processor.
The PFD distributions predicted by COMET are shown in Figs. 3a -3b for the ARO, ARI, and CRT cores, respectively. As expected, the CRT core has a relatively flat PDF distribution. Since the core has a 60º azimuthal rotation symmetry, the MCNP PFDs averaged over the corresponding six symmetric pins were used as the reference solutions. The comparison of the COMET PFDs and the reference solutions (symmetric average) is presented in Table III. For comparison, the deviations of the MCNP PFDs from the symmetric averages and the MCNP PFD uncertainties are also summarized in Table III.  for the three core configurations, respectively. It is also found that the COMET PFDs are in better agreement with the reference results than those predicted by MCNP. In order to achieve the same accuracy as the COMET PFDs, at least four times more active particle histories must be tracked for MCNP calculations. However, this also increases the MCNP computation time by about a factor four.
(a) Deviation of MCNP PFD from the symmetry (b) Relative difference between COMET and reference PFD   Figure 6. Comparison of the PDF distribution for the CRT core.

SUMMARY AND CONCLUSIONS
The recently developed 2D SmAHTR benchmark problems were used to test the performance of the COMET method. The three SmAHTR core configurations were modeled as a combination of eight unique coarse meshes. The continuous energy COMET response function generator was first used to solve the fixed-source problems for each unique coarse mesh to generate the response coefficient library, and then the COMET whole core solution was obtained by converging the core eigenvalue and partial currents crossing coarse meshes. It was found the eigenvalues predicted by COMET agreed very well with the MCNP results, with a difference 40 pcm to 51 pcm. The pin power distributions computed by COMET were in excellent agreement with the reference solutions, with the average and maximum relative differences of 0.20% to 0.21% and 0.77% to 0.94%, respectively. It can be concluded that COMET can model the SmAHTR core configuration with high fidelity and significantly high computational efficiency.

ACKNOWLEDGMENTS
This research is being performed using funding received from the DOE Office of Nuclear Energy's Nuclear Energy University Program.

DISCLOSURE
The second author owns equity in a company that has licensed the COMET technologies from Georgia Tech. This study which is a demonstration of COMET could affect his personal financial status. The terms of this arrangement have been reviewed and approved by Georgia Tech in accordance with its conflict of interest policies.