Evaluation of RAPID for a UNF cask benchmark problem

This paper examines the accuracy and performance of the RAPID (Real-time Analysis for Particle transport and In-situ Detection) code system for the simulation of a used nuclear fuel (UNF) cask. RAPID is capable of determining eigenvalue, subcritical multiplication, and pin-wise, axially-dependent fission density throughout a UNF cask. We study the source convergence based on the analysis of the different parameters used in an eigenvalue calculation in the MCNP Monte Carlo code. For this study, we consider a single assembly surrounded by absorbing plates with reflective boundary conditions. Based on the best combination of eigenvalue parameters, a reference MCNP solution for the single assembly is obtained. RAPID results are in excellent agreement with the reference MCNP solutions, while requiring significantly less computation time (i.e., minutes vs. days). A similar set of eigenvalue parameters is used to obtain a reference MCNP solution for the whole UNF cask. Because of time limitation, the MCNP results near the cask boundaries have significant uncertainties. Except for these, the RAPID results are in excellent agreement with the MCNP predictions, and its computation time is significantly lower, 35 second on 1 core versus 9.5 days on 16 cores.


Introduction
To ensure criticality safety and maintain a high level of security of Used Nuclear Fuel (UNF), it is necessary to perform detailed Particle Transport simulations. The commonly used technique is the 3-D, continuous energy Monte Carlo method that is time-consuming and, for eigenvalue calculations, may have difficulties with the source convergence [1,2]. Wenner and Haghighat [2] and Dufek [3] demonstrated that in criticality problems with high dominance ratio and/or under-sampling possibilities, the fission matrix (FM) methodology [1] should be used in order to be able to obtain unbiased solution.
The RAPID (Real-time Analysis for Particle transport and In-situ Detection) code system [4] is based on the Multi-stage Response-function Transport (MRT) methodology [5], and uses the FM formulation to determine eigenvalue, subcritical multiplication and fission density distribution of nuclear systems. This paper uses a benchmark problem similar to the GBC-32 Cask [6] to evaluate the accuracy and performance of RAPID against the results of a reference MCNP [7] solution. We performed a detailed sensitivity analysis on a single-assembly model in order to determine the MCNP eigenvalue parameters necessary to obtain a reliable reference solution. These results provided a guide for the selection of an appropriate set of parameters for the eigenvalue calculation for the whole cask.
The paper is organized as follows. Section 2 discusses the RAPID code system and its MRT approach. Section 3 a e-mail: val@vt.edu b e-mail: haghighat@vt.edu c e-mail: roskofnj@vt.edu introduces the benchmark problem. Section 4 introduces the methodology used to analyze results. Section 5 discuss the results of different cases to identify the best set of eigenvalue parameters, and compares the selected case with the RAPID calculation. Section 6 gives the conclusions and future work perspectives.

RAPID Code System
To overcome the source convergence difficulty and computational effort of straightforward eigenvalue Monte Carlo calculations, the RAPID code system has been developed based on the MRT methodology. For a UNF pool or cask, RAPID is capable of quickly and accurately calculating eigenvalue, subcritical multiplication, pin-wise, axially dependent fission densities, and detector responses for monitoring a pool or cask. RAPID uses the Fission Matrix (FM) method for eigenvalue calculations, and the adjoint function methodology for the detector response calculations. RAPID partitions a problem into stages in which FM coefficients and importance functions are determined for different values of the key parameters such as enrichment, burnup, geometry, and cooling time. Then, it uses a system of equations to solve for eigenvalue and eigenfunction and the adjoint function methodology for determining a detector response. The methodology has been succesfully demonstrated for spent fuel pools [4].

The Fission Matrix approach
The FM method can take two forms, depending on the type of problem. For a sub-critical multiplication problem, in which the fission source is driven by an independent source in the spent fuel (i.e., spontaneous fission and (α, n) reactions), the induced fission source in cell i is given by Equation 1.
Where F j is the induced fission source strength in fuel pin j, S j is the intrinsic (or independent) source strength in fuel pin j, a i, j is the number of fission neutrons directly produced in fuel pin i due to a fission neutron born in fuel pin j, b i, j is the same as a i, j except for intrinsic source neutrons. These values are different because S and F have different spatial and spectral distributions within each cell. N is the total number of computational cells (which could be a whole assembly, a single fuel pin, or a fraction thereof). We also consider the eigenvalue problem as in Equation 2.
The method results in a set of N linear equations, which can be solved for F and k given the a i, j coefficients. The main task is to determine the coefficients a i, j and b i, j in the case of subcritical multiplication. Since not all the regions are coupled, i.e., have non-zero coefficients, we have developed a multi-layered approach to determine the necessary coefficients. Coefficients a i, j are calculated through fixed-source Monte Carlo calculations, in which neutrons are originated in cell j and fission neutrons born in cell i's are tallied. The source neutrons are sampled over energy according to the Watt fission spectrum, isotropically. For b i, j the same applies except that neutrons are sampled according to the spontaneous fission and (α, n) reactions spectra.

Benchmark problem
The GBC-32 cask is loaded with 32 Westinghouse 17x17 OFA/V5 fuel assemblies (FA) of different burnups [6]. For the purpose of this work, the cask has been loaded with fresh fuel with 4% wt. enrichment only. The model consists of the 32 FAs contained in a stainless steel canister. The canister is equipped with Boral plates encased in Aluminum clads between the FAs, and it is placed in a Stainless Steel cylinder. The system is flooded with water. The model is depicted in Figures 1a and 1b. Figure 2 shows one assembly of the GBC-32 Cask.
For the purpose of analysis, two representative pins are highlighted (circled black) in Figure 2: pin (1,9) near the absorber panel on the left of the assembly, and pin (10,10) top-right of the central guide tube. These two pins have been chosen to plot fission density axial distributions because of being subject to the lowest and the highest values of fission density, respectively.

Methodology
The aim of this work is to compare the RAPID calculated k e f f and pin-wise axial fission densities to those of a reference MCNP calculation. To obtain an accurate reference solution, an analysis is performed to obtain the best combination of the eigenvalue parameters used in the MCNP's KCODE calculations.

Determination of eigenvalue parameters in MCNP
In any Monte Carlo eigenvalue calculation, a major issue is achieving an accurate fission source distribution. This requires determination of three MCNP parameters, including number of skipped cycles (NSK), number of active cycles (NAC), and number of histories per cycle (NPS). Additionally, the variance of the fission neutron source may be significantly biased by the presence of cycle-tocycle correlation, which is difficult to detect and may cause MCNP reported variance to be underestimated [8]. This phenomenon arises from the fact that the power method [1] (used for eigenvalue Monte Carlo calculations) computes the fission source at iteration k + 1 based on the distribution at iteration k. The random numbers generated by the pseudo-random number generator (PRNG) for a certain iteration will be affected by the previous cycles, i.e., correlated. This effectively causes the fission density tallies to oscillate less than their actual statistical uncertainty, leading to an underprediction of the computed uncertainty. This effect usually arises in the fission density distribution tallies and not in the k e f f since the latter is an integral quantity.
In an attempt to detect the potential underprediction of the MCNP reported tally uncertainty, a replication study is performed for a fixed set of the eigenvalue parameters by changing the seed of the PRNG.

Analysis approach for the establishment of a reference MCNP single-assembly model
The Monte Carlo eigenvalue parameters, NPS, NSK, and NAC, have been investigated by analyzing the eigenvalue (k e f f ), the normalized fission density(i.e., source distribution) and the associated uncertainties for the singleassembly model. The fission density is determined in 24 axial intervals for every fuel pin in the model, resulting in 6336 tally regions. Convergence of the source is examined through the following techniques for different combinations of the eigenvalue parameters: • The Shannon entropy stabilization as a function of generations [9]. It is worth noting that, as demonstrated by Wenner and Haghighat in [2], this method may fail to correctly diagnose the source convergence in loosely coupled problems.
• The L ∞ , L 1 , and L 2 norms of the relative differences of the fission neutron density tallies with respect to the most accurate calculation are computed. These norms are compared to the norms of the uncertainties associated with each tally for the particular calculation, and are defined as follows: where X is a vector composed of N t elements. X will be either the vector of relative differences or uncertainties. In principle, the norm of the relative difference vector should be lower than the norm of the statistical uncertainty if the source has converged. This would indicate that, on average, the variation of the results from case to case is within the statistical uncertainty. Relative differences are computed for each of the 6336 tally with respect to the most detailed calculation as follows: where R k,l is the value of the relative difference of the k-th tally for the l-th case, S k,l is the value of the k-th tally for the l-th case, and S k,md is the value of the k-th tally for the most detailed (md) case.
• The Center of Mass (COM) [1] is calculated by: where m i is fission density at region i, r i is the distance between region i and geometric center, M is the total fission density of the assembly, andr is the COM distance to the geometric center. The COM position evaluated in Equation 5 should oscillate around the geometric center of the model for a symmetric problem.
In addition to the aforementioned analyses, based on a select parameters combination we performed repetition of the same run with different PRNG seeds in order to identify the presence of possible cycle-to-cycle correlation through the approach showed in [8]. From the repetition of the run, the average MCNP calculated variance is calculated as follows: where N s is the number of different seeds used and σ i is the standard deviation obtained with MCNP for each of the runs.
Having N s different values for the same quantity, the "real" standard deviation of the sample can be defined as: where x i is the value of a certain quantity calculated for the i-th seed, andx is the average of all these values. By calculating the value of these two standard deviations, an indication of the underestimation of the MCNP variance due to cycle-to-cycle correlation is provided: The f s parameter is evaluated using Equation 8 for the system eigenvalue, k e f f , and for each of the axially dependent pin-wise fission density tallies. Based on the above described evaluation process, a reference MCNP calculation is established for the singleassembly model to be compared with RAPID results.

Establishment of the full-cask MCNP reference model
Based on the results of Section 4.2, the set of parameters selected for the single-assembly are scaled to represent a full-cask model. Since the assemblies are separated by absorbing racks, thus assemblies are loosely coupled, the whole-cask model is established by keeping NAC and NSK constant with respect to the single-assembly case.
The NPS parameter, in principle, should be scaled proportionally to the increase in the number of tally regions, resulting in a scale-up factor of 32. However, since a NPS value of 10 6 (with 144 particles per tally region) has been used for the single assembly model, the directscaling approach would lead to prohibitive computation times. Therefore, we have used a reasonable NPS of 10 5 per assembly (demonstrated in Section 5.1), i.e., 3.2 · 10 6 for the full cask.

Comparison of MCNP reference models to RAPID calculations
The RAPID system calculation is compared to the MCNP reference calculation for both the single-assembly and the full-cask models. Note that the full-cask model is only partially converged because of time limitations. The quantities analyzed are the criticality eigenvalue parameter of the system k e f f and the pin-wise axial fission density distributions. The relative differences between the two methods are calculated and plotted, and maximum and average values are tabulated.

Results
This section is divided into four main parts: in Section 5.1 the parametric study of the single-assembly problem to assess adequate eigenvalue calculation parameters to establish a referene MCNP solution is presented; Section 5.2 addresses the cycle-to-cycle correlation issue; while in Sections 5.3 and 5.4 the results of the RAPID code are compared to the MCNP reference calculation for the single assembly and full cask models, respectively.

Eigenvalue Parameters
The parametric analyses are made by independently varying the three MCNP criticality parameters (NAC, NPS, and NSK) in order to isolate the effects of the variation of each of them. The NAC study has been conducted fixing NSK to 1000 and NPS to 10 6 , varying NAC from 500 to 3500 in steps of 500. The Shannon entropy behavior is presented in Figure 3 as a function of the cycle.  Figure 3 shows that the value of the Shannon entropy oscillates around a stable value for NAC greater than 300.
Such a behavior is usually interpreted as a sign of source convergence [9].  Above figure shows that k e f f has converged, and the various runs differ by less than one standard deviation. Figure 5 instead shows the behavior of the statistical uncertainty σ k as a function of NAC. The statistical uncertainty shown in Figure 5 correctly drops as 1/ √ N, as expected based on the Central Limit Theorem [1]. Figure 6 shows the percentage distance of the COM from the expected geometric center, relative to half the active fuel height, for varying NAC: where r i is the distance of the i th region from the center of the assembly, S i,n is the fission density value in the i t h region for the n th run, and H is the active fuel height.  Figure 6 shows that the COM oscillates of less than 0.5% (with respect to half the active fuel height along the z-direction). Such an oscillation is reasonable and indicates that the source has converged.
For the fission density distribution, the norms of the R k,l vector (|R| l , R l , and R ∞,l ) have been computed as described in Section 4.2 for every run l. Additionally, the norms of the vector U k,l , containing the relative uncertainties of each tally k as output from MCNP (|U| l , U l , and U ∞,l ) have been calculated. Figures 7, 8, and 9 show the behavior of these norms as a function of NAC.    Figures 7, 8, and 9 indicate that the relative differences are on average higher than the statistical uncertainty associated to the calculation. Figures 10 and 11 show the fission density axial distribution in the representative pins (10,10) and (1,9) shown in Figure 2. . Fission density axial distribution in pin near absorber - (1,9) Above figures show that the changes in the axiallydependent fission density distributions are on average higher than the MCNP predicted uncertainty. However, the shape is basically unaffected by the increasing number of active cycles, except for the NAC=500 (i.e., the lowest number considered) curve being slightly offset. Similar analyses have been done for the NSK and NPS parameters independently, keeping fixed the other two. In all the cases, simlar results are obtained, leading to the following conclusions: • The criticality eigenvalue, k e f f , converges and varies within the uncertainty for all the various runs; • The shape of the fission density distribution is converged while on average the variation of the tallies is above the statistical uncertainty, suggesting that these could be underestimated; • The Shannon entropy and COM behaviors are stabilized and would suggest that the fission source is not varying much and is symmetric as expected.
• The norms of the relative differences are higher than the statistical uncertainties ones.
The first three points suggest that the fission source has converged. However, the norms show an unexpected behavior: this can be caused either by a bias in the sample averages calculated by MCNP (i.e., a non-converged source, in contrast with the first three analyses) or an underestimation in the statistical uncertainties (e.g., due to cycle-to-cycle correlation). Section 5.2 analyzes the latter phenomenon.

Cycle-to-cycle correlation investigation
In order to assess the effect of the cycle-to-cycle correlation on the computed variance, we performed a replication study keeping the same set of MCNP criticality parameters and changing the PRNG seed. In this way, a "real" standard deviation is calculated and compared to the average MCNP uncertainty, as described in Section 4.2 and in [8]. This study was performed with smaller criticality parameters (NPS=1E6, NSK=300, and NAC=500) with respect to the study done in Section 5.1. This is necessary because significant computation is needed for running the 50 replication cases.
The approach is applied to both the criticality eigenvalue, k e f f , and each of the 6336 fission density tallies. Table 1 give k e f f and the associated standard deviations ("real" and MCNP calculated). The analysis reveals that the "real" standard deviation, calculated through this replication study, is actually lower than the MCNP calculated one. This is mostly related to the fact that MCNP standard deviations for k e f f are provided with a resolution of 1 pcm. This result is in line with what is shown in Figure 4: the MCNP k e f f is converged and the estimation of the statistical uncertainty is accurate.
For the fission density values, the following quantities are analyzed: • The maximum of the real to MCNP standard deviation ratios for all the k-th tallies: • The average of the real to MCNP standard deviation ratios for all the k-th tallies: • The weighted average of the real to MCNP standard deviation ratios using the fission density source for all the k-th tallies: Table 2 summarizes the real and the MCNP standard deviations and these ratios. The statistical uncertainties for the fission density distribution are underestimated by the MCNP calculations by a factor greater than 2, identifying a strong cycle-to-cycle correlation. A 3D plot of the f σ ratios for these tallies is reported in Figure 12.  Figure 12 demonstrates that there is, on average, a significant underestimation of the uncertainties and that this underestimation changes significantly in different regions of the model.
As discussed in References [8,9], cycle-to-cycle correlation is related to the method itself and it is not affected by the use of a "better" set of MCNP parameters. Therefore, we can compute the adjusted relative uncertainties for the NAC, NSK, and NPS parametric analyses (Section 5.1), multiplying each of the tally uncertainties by the weighted average ratios, f σ,wgt . This would take into account the presence of cycle-to-cycle correlation in the comparison of the relative differences and relative uncertainties for the various cases. The relative differences L 1norms, |R| l , are plotted against the adjusted relative uncertainties L 1 -norms, |R| l,ad j , for each of the parametric studies in Figures 13, 14, and 15.   The f σ should hold similar results even for higher NAC, NSK, and NPS, since cycle-to-cycle correlation is independent of these parameters, as long as the source has converged, i.e., NSK is high enough [9]. Figures 13,  14, and 15 suggest that the tally variations are on average within the statistical uncertainty associated to the method for all the analyzed cases. Therefore, source convergence is achieved for reasonably lower values of the MCNP criticality parameters although the uncertainty calculated by the code is underestimated due to the cycle-to-cycle correlation. From this latter analysis we conclude that the NSK=500, NAC=1000, NPS=10 6 constitutes a suitable reference model for comparison with the RAPID code system.

Comparison of RAPID and MCNP results for a Single Assembly
The RAPID code was used to model the single assembly reflected in the x−y directions. Table 3 presents the results of the eigenvalue calculation. Table 3. k e f f comparison of RAPID compared to MCNP reference calculation.
Case Above table shows that the RAPID k e f f calculation (NSK=500, NAC=1000, NPS=10 6 ) is in very good agreement with the MCNP reference calculation, with a relative difference of 52 pcm.   Table 4 indicates that the RAPID calculation on 1 computer core is 6500 times faster than the MCNP reference calculation on 16 computer cores. Figure 16 shows the axial fission distribution of the entire assembly, i.e., integrated over the x−y dimensions) for RAPID and the MCNP reference calculation. Note that all uncertainties on the MCNP reference calculation are 1−σ. Uncertainties are within 0.05% for the x − y integrated axially-dependent fission density and within 0.59% for the pin-wise axially-dependent one.   The relative differences for all the regions are below 5%, with regions above 2% being those in the top and bottom layers. This behavior at the fuel rod tips is due to the assumed axial symmetry in the way FM coefficients are calculated. Maximum, average, fission-density weighted average and adjusted-weighted average relative it uncertainties for the MCNP single assembly calculation are given in Table 5. Table 6 gives maximum, average, and weighted-average relative differences of RAPID vs. MCNP.

Comparison of RAPID and MCNP results for the full-cask
As already discussed in Section 4.3, the set of parameters used for the full-cask reference MCNP calculation are NPS=3.2 · 10 6 , NSK=500, NAC=1000 . k e f f for RAPID and MCNP reference calculations are given in Table 7 for the full cask. These results show excellent agreement. Table 7. k e f f comparison of RAPID compared to MCNP reference calculation -Full cask.
Case k e f f Rel. Diff.
[pcm] The maximum, average, and weighted average uncertainty values for the full-cask MCNP calculated fission density distribution tallies are reported in Table 8. Similar to the single-assembly case, the RAPID to MCNP relative differences for the pin-wise axiallydependent fission density distribution have been calculated. Figure 18 shows the 3D distribution of these relative differences. Figure 18. Relative differences of RAPID vs MCNP calculated pin-wise axially-dependent fission density for the full cask. Figure 18 shows how the relative differences are lower than 2% the inner regions, where the MCNP statistical uncertainties are smaller. The highest relative differences are on the boundaries: it is important to note that, due to under-sampling [10], the peripheral fission densities calculated by MCNP are far from convergence and might be significantly biased. Accurate estimation of the fission density in these regions is required for detector responses estimation and dose calculations on the outside surface of the cask. Hence, the fission matrix approach should be used instead of straight-forward Monte Carlo calculations to achieve accurate unbiased detector responses. Cask dose calculation has been recently implemented in RAPID and demonstrated in [11]. The average relative difference, weighted on the fission density distribution, amounts to 1.56%, comparable to the statistical uncertainty reported in Table 8. Table 9 shows the computational time requirements for both the methods, demonstrating how a large nuclear system as the full cask can be simulated in real time with RAPID, with a significant speedup with respect to MCNP.  Figure 19 shows RAPID calculated 3-D fission density distribution.

Conclusions
This paper examines the performance of the RAPID code system for the simulation of a UNF cask benchmark, GBC-32. Because of the presence of absorbing plates in between the assemblies in a cask, any eigenvalue Monte Carlo simulation faces with the difficulty of particle undersampling. Hence, it is necessary to examine carefully the source convergence. To do so, a detailed statistical analysis is performed using a single assembly model with reflective boundary conditions. This analysis showed the importance of carefully selecting the criticality run parameters when using MCNP, and underlined that substantial bias in the standard deviation of the tallies might be caused by cycle-to-cycle correlation, related to the nature of criticality calculations and impossible to foresee. The RAPID code system, since it is based on a set of physically based transport coefficients that do not require any eigenvalue calculation, doesn't suffer from any of these difficulties. Based on the parametric analysis, an appropriate set of eigenvalue parameters (NAC, NSK, NPS) have been selected to establish both a single-assembly and full-cask reference model for comparison with RAPID. It is demonstrated that the calculated eigenvalue and fission density distribution by the RAPID and MCNP code systems are in excellent agreement with the reference MCNP results for both the models, leading to average relative differences of the order of a few percents mostly related to boundary effects that will be addressed in future work.
In addition, under-sampling may cause biased boundary MCNP tallies even when large criticality parameters are used. This issue is avoided using the FM approach, since no Monte Carlo criticality calculations are necessary. The calculation of external doses and detector responses is highly affected by the accuracy of the outer regions, nearby the detector location, hence making the use of codes based on the FM methodology, such as RAPID, necessary to calculate these quantities.
Moreover, it is shown that the RAPID computation time is significantly (several orders of magnitude) smaller than the MCNP one, i.e., 35 seconds vs. 9.5 days on 16 cores. This speedup increases as the model becomes larger.