ANALYSIS OF TIME-EIGENVALUE AND EIGENFUNCTIONS IN THE CROCUS BENCHMARK

Time-dependent neutron transport in non-critical state can be expressed by the natural mode equation. In order to estimate the dominant eigenvalue and eigenfunction of the natural mode, CEA had extended the α-k method and developed the generalized iterated fission probability method (G-IFP) in the TRIPOLI-4® code. CRIEPI has chosen to compute those quantities by a time-dependent neutron transport calculation, and has thus developed a time-dependent neutron transport technique based on k-power iteration (TDPI) in MCNP-5. In this work, we compare the two approaches by computing the dominant eigenvalue and the direct and adjoint eigenfunctions for the CROCUS benchmark. The model has previously been qualified for keffs and kinetic parameters by TRIPOLI-4 and MCNP-5. The eigenvalues of the natural mode equations by α-k and TDPI are in good agreement with each other, and closely follow those predicted by the inhour equation. Neutron spectra and spatial distributions (flux and fission neutron emission) obtained by the two methods are also in good agreement. Similar results are also obtained for the adjoint fundamental eigenfunctions. These findings substantiate the coherence of both calculation strategies for natural mode.


INTRODUCTION
Over the last ten years, there has been a growing interest for bias and uncertainty quantification in reactor physics calculation based on Monte Carlo methods. In this respect, the discovery of the Iterated Fission Probability (IFP) method has been a major breakthrough for continuous-energy Monte Carlo transport [1][2], enabling the deployment of rigorous first-order reactivity perturbation theory and k-eigenvalue sensitivity analysis [3][4]. By means of these computational tools, it is now possible to quantify the uncertainty of a designed core due to nuclear data and nuclide densities in continuousenergy Monte Carlo calculations. However, it is well known that the reactivities are not directly and absolutely measurable, unless the reactor core is exactly critical. In non-critical state, the neutron asymptotic reactor period [5]. The quantity α is measurable parameter that can be obtained by integral experiments: as such, the determination of α and the associated sensitivity by Monte Carlo methods is of utmost importance for the quantification of the bias due to errors in cross sections data and/or nuclide densities [6]. For this purpose, CEA has first improved the α-k method to compute α and the corresponding direct neutron flux, explicitly taking delayed neutron emission into account [7][8]; later on, they have developed the Generalized IFP (G-IFP) method to compute the adjoint fundamental αeigenfunction [9]. Based on the G-IFP, novel techniques to estimate the perturbations/sensitivity of the α eigenvalue with respect to the nuclear data and the nuclide density have been introduced [10]. Similar efforts have been carried out at CRIEPI, by resorting to a different approach: the direct and adjoint αeigenfunction can be obtained by observing the neutron population for sufficiently long time after a source is placed at the initial position [11]: the distribution of the neutron population at the final time represents the fundamental direct mode; furthermore, the asymptotic population can be used to assess the importance of the initial point-source, and in this respect it acts as an estimator of the adjoint flux for the natural mode equation. In order to compare their newly developed methods, CEA and CRIEPI have jointly performed benchmark calculations of α and the adjoint flux distribution in a 1-D rod-type problem with 2 energy group cross sections. The neutron population obtained by a time-dependent analog transport calculation yields values of α and adjoint flux shapes that are entirely consistent with those given by G-IFP [12]. In this work, we further extend this preliminary investigation by considering full-scale three-dimensional and continuous-energy problems. For this purpose, CEA has implemented the α-k and the G-IFP methods in the development version of the TRIPOLI-4® code [13], whereas CRIEPI has developed a time-dependent neutron transport technique based on the k-power iteration (TDPI) [14] and implemented it into the MCNP-5 code [15]. For our code-to-code comparisons of α eigenvalues and eigenfunctions we have selected the CROCUS benchmark [16]. To the best of our knowledge, this is the first time that similar comparisons are attempted.

Generalized Iterated Fission Probability Method
G-IFP method has been extensively described in [9,12]; here it will be recalled briefly. It consists of two steps. At the first step, the natural mode equation [5] | | is solved by the α-k iteration (see references [7][8] for a thorough description of the α-k method). Here ℒ is the (net) disappearance operator, V is the neutron velocity, and ℱ and ℱ are the prompt and delayed fission production operators, respectively, the index j denoting the precursor family. The quantities λ i denote the precursor decay constants. At convergence, i.e., after a sufficiently large number of cycles has been discarded, the α-k iteration yields both the dominant eigenvalue α and the associated eigenfunction ϕ(r, V). In the second step, neutrons are initially placed at a position r and velocity V.
The G-IFP then proceeds iteratively. From the source, a neutron transport calculation is performed in accordance with the left hand side of eq. (1), and the source for the next step is obtained from the production terms (the right-hand-side of the equation). The delayed neutron emission ℱ ϕ is weighted by λ j /( λ j +α). During the G-IFP iterations, α is kept fixed to the value determined in the α-k iteration.
By iterating the neutron transport and source estimation, the distribution of the neutron population in a cycle corresponding to the initial source converges to the adjoint dominant mode ϕ † (r, V) satisfying after a sufficiently large number of cycles has been discarded.

Time-Dependent Neutron Transport
CRIEPI has recently investigated the neutron population balance for sufficiently long time after the neutrons are placed at the phase space position (r, V), and obtained eq. (2) by extending the derivation by Bell and Glasstone [17]. In the asymptotic regime (i.e., after a "sufficiently long time"), the neutron population varies exponentially with time and is distributed according to the dominant α-mode, with exp(-(α+λ i )t)→0 [11,12]. Following this strategy, the adjoint flux ϕ † (r, V) can be obtained by putting a neutron source at the position (r, V) and simulating the time-dependent neutron transport up to some long time. However, analog time-dependent transport is a very challenging task for Monte Carlo codes, due to population multiplication or attenuation, especially if the emission of delayed neutrons from precursors must be taken into account. In order to overcome this issue, CRIEPI has developed TDPI, whereupon time-dependent calculations are performed by adapting the k-power iteration (PI) [14]. In the conventional PI, the ratio of fission neutron emission in cycle i to that in cycle i-1 is estimated in order to obtain k i . Thus, the neutron emission per fission is artificially scaled by 1/ k i-1 for the fission source in cycle i: the neutron population is constantly adjusted by this procedure. The ratio of the real number of neutrons to that of the PI calculation is expressed as Π m=0 i-1 k m . This scheme can be used to "normalize" analog time-dependent calculations as follows. Consider a neutron flight in cycle i at time t, originated from the source at (r, V) at the initial time t 0 =0. The current time is quantified by summing over the flight times of all the ancestors (including the precursor decay times). Accordingly, the neutron flux at time t is tallied at every PI cycle and written on a file. Then, the tallied flux is multiplied by the scaling factor Π m=0 i-1 k m . By summing the scaled flux contributions, we obtain the time dependent neutron flux resulting from a point source [14]. This calculation scheme requires only minimal modifications in MCNP-5 code.

Description of the CROCUS Benchmark
CROCUS is an open-tank type zero-power reactor operated at the EPFL (Lausanne, Switzerland). As schematically shown in Fig. 1, two kinds of fuel rods are aligned in square lattices in inner and peripheral regions. In the inner region, UO 2 fuel at 1.806 wt% 235 U /U is loaded, with 1.837 cm pitch. In the peripheral region, metallic uranium at 0.947 wt% 235 U /U is loaded, with 2.917 cm pitch [16]. In the core, the criticality is attained by filling with light water from the bottom cadmium plate, up to the appropriate height. According to the benchmark specifications, the criticality is attained at a water level of 96.51 cm. CEA and CRIEPI have already performed code-to-code comparisons against a benchmark model of CROCUS, for the V&V of the reactor kinetics parameters by the IFP method implemented into TRIPOLI-4 and MCNP-5 [18]. In [18] k eff s were computed by TRIPOLI-4 and MCNP-5 for water levels of 96.51, 98.51, 99.00, and 99.51 cm, and the obtained values were within 5 pcm (for three different nuclear data libraries, JEFF3.1.1, ENDF/B-VII.0 and JENDL-4.0, specially processed for the two codes). These results show that the independent CROCUS models developed by the two institutions are coherent. Similar tests were performed on the effective delayed neutron fraction and the neutron generation time, and the agreement was also extremely good [18].
In the present work, we focus on code-to-code comparisons for α-eigenvalue and eigenfunctions, using the JEFF-3.1.1 library. The α eigenvalues are very sensitive to small reactivity changes, as discussed in reference [19]. Moreover, in the benchmark model, the k eff is known to be biased by ~ 250 pcm [18]: in other words, at the supposedly critical water height the codes yield a positive reactivity. The water levels for our calculations are listed in table 1 (the level that yields a unit k eff has been also determined). The α-k and G-IFP methods are not particularly sensitive to the reactivity level; on the contrary, for supercritical configurations TDPI calculations are limited to small positive reactivity values, up to 200 pcm, and larger statistical fluctuations are expected for delayed super-critical systems [14].   [18]; (**): deduced from k eff and the kinetic parameters by MCNP-5 [19].

Direct calculations: dominant α-eigenvalue and eigenfunctions
Time-dependent neuron transport was simulated by TDPI for a source of 2 MeV located in the fuel pin at the center of the core (the black dot surrounded a rectangle in Fig. 1), corresponding to the water levels listed in Table 1. Nominally, 5.0E+6 neutrons per cycle were simulated, along 20000 cycles. For water levels from 80.0 cm to 91.66 cm, the end of the last cycle corresponded to more than 780 s after the source. The estimated neutron flux integrated over the whole fuel pins is shown in Fig. 2. For sub-critical configurations, the flux decreases with time and after a while decays exponentially (a straight line in the semi-log graph). For lower water levels, statistical fluctuations become significant after 300 s, especially for the H=80.0 cm case. For the water levels from 80.0 cm to 91.66 cm, the dominant eigenvalue α is deduced by fitting the time decay slope from 300 to 500 s 1 . As for positive reactivity, TDPI is hindered by statistical uncertainty. For the water level H=96.51 cm, 8 calculations varying the random number sequences were performed. The standard deviation of the 8 runs exceeds the average value at 43 s. Thus, α was estimated by fitting from 10 to 42 s. For water levels from 98.61 to 99.51 the statistical uncertainty prevented TDPI from converging.
The dominant eigenvalues α were also estimated by the α-k method for all water levels listed in Table 1, with 5000 inactive and 5000 active cycles and 20000 neutrons per cycle. Furthermore, the asymptotic reactor period was also deduced by solving the inhour equation using the kinetics parameters at H=96.51 cm obtained in our previous work [18]. The simulation findings are compared in Fig. 3.
TDPI and the α-k method yield α values that are both consistent with each other and with the inhour equation. For H=96.51cm, the small discrepancy between TDPI and α-k is possibly due to the cut-off at 42 s for the TDPI, which might lead to an under-sampling of the delayed neutrons of the j=1 precursor family.  For the water level of H=90.80 cm, the spatial shape of the neutron flux and fission neutron emission was tallied by TDPI from 400 to 500 s after a source is placed. The tallied fuel pins are shown in Fig.  1. Also, the average neutron spectra over all fuel pins were tallied. In the TDPI calculation, the source of 2 MeVwas placed as the same central position as above. A total of 1E+06 neutrons per cycle were simulated, along 20000 cycles. The computed neutron flux ϕ and fission neutron emission ℱ ϕ are compared to those obtained by the α-k method in Figs. 4 and 5, respectively. In the α-k method, we have used the same number of cycles and neutrons per cycle as above. The results obtained by TDPI with MCNP-5 and the α-k iteration with TRIPOLI-4 are entirely consistent, both spatially and spectrally, for fluxes as well as reaction rates.

CONCLUSIONS
In this work, we have focused on the natural mode equation and its adjoint formulation, in view of performing code-to-code comparisons on benchmark configurations. For this purpose, CEA has chosen to solve the natural mode equation and its adjoint equation as an eigenvalue problem, by using the α-k iteration and the G-IFP method, respectively, implemented in the development version of TRIPOLI-4. On the other hand, CRIEPI has chosen to estimate the direct flux and the neutron importance function by using time-dependent transport methods and following the neutron population from a source to its asymptotic behavior. For this purpose, CRIEPI has developed the TDPI technique, which was implemented in MCNP-5. In order to verify these approaches in continuous-energy Monte Carlo transport problems, we have computed the dominant eigenvalue α, the neutron flux ϕ , the fission neutron emission ℱ ϕ, and the adjoint flux ϕ † for some significant configurations of the CROCUS benchmark. The benchmark models have been previously qualified for k eff s and kinetics parameters by TRIPOLI-4 and MCNP-5 [17]. In the range of reactivities from -850 to 250 pcm, for the α eigenvalues TDPI gives results entirely consistent with those of the α-k iteration (and both agree with the predictions of the inhour equation). For larger reactivities, TDPI is hindered by statistical fluctuations, but the α-k iteration is consistent with the inhour equation in the range explored in this work (up to 360 pcm). For a slightly sub-critical state, we have further compared the direct flux ϕ and reaction rates ℱ ϕ, as well as the adjoint flux ϕ † : again, the TDPI is in excellent agreement with the α-k and G-IFP methods.
Our findings provide a firm basis for the use of G-IFP, and at the same time confirm the TDPI method as an option to estimate both direct and adjoint tallies related to the natural mode equation. Generally speaking, G-IFP involves a smaller computational burden, but TDPI is extremely easy to implement in Monte Carlo codes already having a standard k-power iteration tool. The complementary use of these simulation techniques might contribute to shed new light on the sensitivity of α-eigenvalues to nuclear data, which would involve the knowledge of direct and adjoint eigenfunctions [6,10].

ACKNOWLEDGMENTS
TRIPOLI-4® is a registered trademark of CEA. A. Jinaphanh and A. Zoia acknowledge partial financial support from EdF.