A SEMI-ANALYTIC EIGENVALUE EXTENSION TO THE DOPPLER SLAB ANALYTIC BENCHMARK1

Advancement in multiphysics simulation has motivated interest in availability of analytic and semi-analytic benchmark solutions. These solutions are sought because they can be used to assess the accuracy of complicated numerical schemes necessary to simulate coupled physics systems. While there exist analytic solutions for fixed-source problems, benchmark-quality eigenvalue solutions are of interest because eigenvalue problems more closely align with analyses undertaken with coupled solvers. This paper extends a fixed-source benchmark, the Doppler Slab benchmark, to the eigenvalue case. A novel solution for this benchmark is derived. Numerical implementation of the benchmark is demonstrated through verification of numerical computation of the power reactivity coefficient.


INTRODUCTION
Interest in the characterization of the stability and the convergence characteristics of coupled-physics simulations for nuclear engineering applications has motivated the development of analytic, semi-analytic, and method of manufactured solution benchmarks, such as [1][2][3][4][5][6]. One such benchmark recently developed is the Doppler Slab benchmark [7], which couples neutron transport with thermal conduction physics via Doppler broadening. However, like most multiphysics analytic benchmarks, the radiation transport in the original Doppler Slab benchmark was driven by a fixed source. While the utility of these benchmarks is in their ability to highlight convergence behavior and provide verification of code implementations of established methods, eigenvalue benchmarks are desirable because they more accurately represent multiphysics problems of interest. In particular, certain quantities in eigenvalue calculations, such as reactivity coefficients, are sought that are fundamentally different than those computed in fixed-source problems.
There are few previous works that extend analytic benchmarks to the eigenvalue case. One example of such a benchmark is the eigenvalue Candlestick depletion benchmark by Kooreman and Griesheimer [3], which coupled neutron transport and depletion. Another example is Gonzales et al. [6], which derived analytic expressions for eigenvalue and for an infinite medium problem. While each of these benchmarks' applicability to their relevant physics considerations has been demonstrated, the utility of Kooreman and Griesheimer's work is twofold. Their methodology for extending a fixed-source problem to an eigenvalue problem can be extended to other benchmark applications. By posing an eigenvalue problem with the same assumptions presented in [3], an eigenvalue problem can be made to "look" like a fixed-source problem. If the solution to this fixed source problem is known, it can be used to construct a solution to the eigenvalue problem.
Using this approach, this paper extends the Doppler Slab benchmark to an eigenvalue problem. The solution can best be characterized as semi-analytic rather than fully analytic because expressions for flux and eigenvalue are dependent upon a parameter that must be found via root-finding algorithms. However, since these techniques can give results to any desired precision, the eigenvalue extension of the Doppler Slab problem is a high-precision benchmark applicable to neutron transport codes coupled with thermal conduction feedback. Further, the utility of extending this benchmark to the eigenvalue case is demonstrated through the calculation of reactivity coefficients. A coupled code system with the in-house Monte Carlo code MC21 [8] and an in-house thermal conduction solver is used to numerically evaluate the benchmark through evaluation of the power reactivity coefficient via a brute-force method.

BENCHMARK SPECIFICATION
To extend the Doppler Slab benchmark to the eigenvalue case, we modify the physical conditions of the benchmark given in [7], by making similar assumptions to those made in [3]. Consider the two energy group neutron transport in a one-dimensional homogeneous slab of purely-absorbing fuel occupying = [0, ] with cross section Σ = Σ in the thermal group, no interactions in the fast group, and thermal conductivity . To the left of the slab is a perfect moderator material with constant temperature , such that all exiting fast neutrons on the left side of the slab re-enter the fuel slab as thermal neutrons. The right of the slab is a vacuum, and the interface between the fuel and vacuum is a perfect thermal insulator with = 0. The geometry for the problem is given in Figure 1.  While Kooreman and Griesheimer were interested in a coupled neutron transport-depletion problem, we analyze the case of neutron transport with thermal feedback. In the fixed-source version of the Doppler Slab benchmark, two Doppler feedback mechanisms were considered. However, only inverse root feedback (in the language of the previous paper) is considered here. This feedback is given by Eq. (1): In Eq. (1), Σ is macroscopic cross section as a function of spatial coordinate , Σ is the unperturbed (reference) macroscopic cross section, and is some reference temperature (for simplicity, is assumed to be the temperature at the left-hand boundary of the slab). If neutrons are restricted to travel only on the -axis and scatter isotropically -what we call the "bisotropic" approximation -and we assume that the boundary between the fuel and the vacuum is perfectly insulated, then the governing equations for this problem are given by In Eq. (2), ! is the fast left-moving flux, & ' is the thermal right-moving flux, is thermal conductivity, ( is energy released at a fission event, and % is the slab eigenvalue. To obtain a solution, we proceed as in [3] by observing that the bottom two equations in Eq. (2) are similar to those in the original Doppler Slab benchmark specification with the exception of the boundary conditions. We require that the flux boundary condition at the left-hand side of the slab to be a constant and well-characterized; the right-hand temperature boundary condition will be addressed later. To characterize the left-hand boundary condition, first we note that the fast flux can be found by direct integration to be We also note that the power . of the slab is given by the equation With Eqs. (4) and (5), the left-hand boundary condition becomes well-specified: We are then to solve Equation (7) is a fixed-source problem rather than an eigenvalue problem. Further, if the solution to Eq. (7) is known, then the solution to the eigenvalue problem can be constructed. The derivation for the semi-analytic solution for this problem is similar to the solution to the similar fixed-source problem in [7] until boundary conditions are considered. We begin with the following substitutions: With these substitutions as well as the utilization of Eq. (1) for Doppler feedback, the dimensionless equation for neutronics is The dimensionless equation for the thermal physics is Using arguments similar to those given in [7], Eqs. (9) and (10) can be combined to yield a nonlinear integro-differential equation for 3, where > is an integration constant. Applying the variable substitution to Eq. (11) results in the nonlinear initial value problem which has the implicit solution The integration constant > can be found by applying the Neumann boundary condition, which gives With the solution to Eq. (17), > is known, and the benchmark solution can be constructed. It is convenient to say that the solution to Eq. (13) is ?(0), where ?(0) is defined as the inverse of the function By combining the expression for & ' and Eq. (3), the solution for ! is given as By utilizing the boundary condition & ' (0) = ! (0), the slab eigenvalue is given by With Eqs. (20)-(22), a solution to the benchmark is obtained.
One should note that this benchmark is semi-analytic rather than fully analytic. Indeed, the numerical implementation of the V-function represents a semi-analytic result. Furthermore, the V-function depends on the value of 4 through the solution of Eq. (17), while 4 is dependent upon the eigenvalue through However, by combining Eqs. (17), (22), and (23), an equation for H that explicitly only depends on problem inputs Σ , ., , , and is obtained: Once H is obtained, % can be determined via Eq. (22) and 4 can be determined via Eq. (23).
It is also possible to generate expressions for derivatives of % with respect to problem input parameters.
While for this problem, the input parameters are ., , , Σ , and , the focus in this paper is The expression in Eq. (25) is the power reactivity coefficient for the eigenvalue problem.

NUMERICAL IMPLEMENTATION
The solution presented in Section 2 is semi-analytic and is applicable for any physically meaningful input parameters. However, it is useful for the purposes of benchmarking numerical tools to dictate "canonical" values for the benchmark. These values were chosen such that desired thermal feedback of the benchmark is observed and that common behavior between numerical implementations can be easily verified. The canonical values for the benchmark are given in Table 1. For the values given in Table 1, we observe about a 20% change in cross section over the length of the slab. To demonstrate the numerical implementation of the benchmark, a computational model was constructed using MC21 and an in-house thermal conduction solver. MC21 was used to compute the flux distribution, and the thermal conduction solver determined the temperature distribution. The MC21 model consisted of contiguous fuel and moderator slabs with dimensions 3 × 1 × 1 cm and 5 × 1 × 1 cm, respectively. The fuel material consisted of a fictitious purely-fissioning nuclide with p For this exercise, the power reactivity coefficient is computed via a finite-difference approximation: where % ] ƒ is the eigenvalue evaluated at power . to be made, . and . & should be sufficiently close together. Therefore, one should expect the eigenvalues % ] and % ]& to be relatively close numerically as well. Because of these expectations, a large number of particle histories should be run for calculations to estimate power reactivity coefficient so that we may ensure that the estimates statistically distinct, that is, their 95% confidence intervals (CIs) do not overlap.
A running strategy of 50 discard and 10,000 active batches of 2 million particle histories per batch was used for all MC21 calculations. When running simulations, it was observed that apparent eigenvalue convergence is reached within 5 iterations, so only 5 feedback iterations were run in the calculations for this exercise. For all calculations, Robbins-Monro relaxation was used, as in [7]. To estimate power reactivity coefficient, . was 2 W, and . & was 2.01 W. The results for this exercise are given in Table 2. For all cases, the value for power at which the reactivity coefficient was calculated was . & . As seen in Table 2, the estimates for numerical eigenvalues for each power agree within less than 1 pcm with the benchmark eigenvalues and are within the 95% CI for each calculation. The numerical estimation for the power reactivity coefficient also agrees with the benchmark value to within its 95% CI. However, of note is that even though 20 billion active particle histories were run at each feedback iteration and the uncertainties of % ] ƒ and % ] ‚ are 2 pcm, the relative 95% CI of the numerical power reactivity coefficient is greater than 50% of its absolute value. This large relative uncertainty on the estimated power reactivity coefficient highlights the fact that accurate calculations of reactivity coefficients via brute force Monte Carlo methods are notoriously difficult due to the statistical uncertainty in the individual eigenvalue results. The reference power reactivity result for the Doppler slab eigenvalue benchmark can also be used to verify alternative reactivity estimation techniques, such as Monte Carlo perturbation methods. Unfortunately, a detailed examination and comparison of reactivity estimation methods is beyond the scope of this work.

CONCLUSIONS
The Doppler Slab benchmark initially described in [7] has been extended to an eigenvalue problem. The eigenvalue extension was solved semi-analytically with inverse-root Doppler broadening to yield an expression for eigenvalue (% " ) as a function of the boundary and initial conditions of the problem. Like the fixed-source version of the benchmark, the problem is described with simple geometry and is straightforward for numerical evaluation. Because of this, the eigenvalue extension to the Doppler Slab benchmark is an effective tool for benchmarking coupled thermal conduction/neutron transport solvers for quasistatic multiphysics calculations. In addition, an analytical expression for the power reactivity coefficient has been established for the Doppler slab eigenvalue benchmark problem. This analytical reactivity solution is valuable for verifying reactivity estimation techniques such as Monte Carlo perturbation methods. Numerical results produced with the MC21 Monte Carlo transport solver coupled with a simple thermal conductivity solver show agreement with the analytical benchmark solution for eigenvalue and power reactivity, as expected.