Fission matrix based Monte Carlo criticality calculations

Abstract

We have described a fission matrix based method that allows to cancel the inactive cycles in Monte Carlo criticality calculations. The fission matrix must be sampled in the course of the Monte Carlo calculation using a space mesh with sufficiently small zones as it causes the fission matrix be insensitive to errors in the initial fission source. The $k_{\text{eff}}$ and other quantities can be derived by means of the final fission matrix. The confidence interval for the $k_{\text{eff}}$ estimate can be conservatively determined via the variance in the fission matrix.

Introduction

The criticality calculations solve the eigenvalue equation

$$\mathit{ks}(\overrightarrow{r}\text{,}E)=\mathit{Hs}(\overrightarrow{r}\text{,}E)\text{,}$$

where k is the eigenvalue, $s(\overrightarrow{r}\text{,}E)$ is the concentration of fission neutrons with energy E at $\overrightarrow{r}$, and

$$\mathit{Hs}(\overrightarrow{r}\text{,}E)\equiv {\int }_0^{\infty }\text{d}{E}^{\prime }{\int }_{\text{S}}{\text{d}}^3{r}^{\prime }f({\overrightarrow{r}}^{\prime }\text{,}{E}^{\prime }\to \overrightarrow{r}\text{,}E)s({\overrightarrow{r}}^{\prime }\text{,}{E}^{\prime })\text{,}$$

where $f({\overrightarrow{r}}^{\prime }\text{,}{E}^{\prime }\to \overrightarrow{r}\text{,}E)\text{d}E{\text{d}}^{3}r$ is an expected number of first generation fission neutrons produced in the volume element ${\text{d}}^{3}r$ at $\overrightarrow{r}$, in the energy element $\text{d}E$ at E, resulting from a fission neutron born at ${\overrightarrow{r}}^{\prime }$ with an energy $E^{\prime }$, and S denotes the space of the system. A Monte Carlo fission $\mathcal{S}$ of neutrons with specific positions, energies, and statistical weights. Standard Monte Carlo calculations solve Eq. (1) by the power iteration

$${\mathcal{S}}^{(n)}\leftarrow \frac{H{\mathcal{S}}^{(n-1)}}{k_{{\text{eff}}_c}^{(n-1)}}\text{,}n=1\text{,}2\text{,}\dots \text{,}$$

where $H{\mathcal{S}}^{(n-1)}$ is randomly sampled, and $k_{{\text{eff}}_c}^{(n)}$ is an estimate of the fundamental mode eigenvalue $k_{\text{eff}}$ obtained in cycle (iteration) n as

$$k_{{\text{eff}}_c}^{(n)}=\frac{〈H{\mathcal{S}}^{(n-1)}〉}{〈{\mathcal{S}}^{(n-1)}〉}\text{,}$$

where $〈\mathcal{S}〉$ gives the total statistical weight of the neutrons in $\mathcal{S}$. At each cycle n, histories of the fission neutrons from ${\mathcal{S}}^{(n-1)}$ are simulated, and a new fission bank ${\mathcal{S}}^{(n)}$ is sampled. The initial fission source ${\mathcal{S}}^{(0)}$ and $k_{{\text{eff}}_c}^{(0)}$ must be specified by the user. A number of inactive cycles must be performed first to converge the fission source to the fundamental mode. Quantities of interest, like the $k_{\text{eff}}$, detector response, neutron flux, etc., are then scored and combined over the active cycles where the fission source is assumed to be converged. The number of inactive cycles is not known in advance, and must be either guessed (by the user) or diagnosed during the calculation. Nevertheless, in loosely-coupled systems or systems with dominance ratios close to unity the fission source may falsely appear converged due to its very slow convergence, and the inactive cycles may be stopped prematurely; consequently, the final results may get corrupted by the incorrect fission source in the active cycles (Whitesides, 1971).

The Monte Carlo criticality calculations would be more credible if the inactive cycles were not needed at all. We show that an appropriate application of the fission matrix (Carter and McCormick, 1969) can make this possible. The fission matrix H represents a space discretization of the operator H over a space mesh. The $(i\text{,}j)\text{th}$ element of H stands for the probability that a fission neutron born in zone j causes the subsequent birth of a fission neutron in zone i,

$$\mathbf{H}[i\text{,}j]=\frac{\int {\int }_0^{\infty }\text{d}E\text{d}{E}^{\prime }{\int }_{{\text{Z}}_i}{\text{d}}^{3}r{\int }_{{\text{Z}}_j}{\text{d}}^3{r}^{\prime }f({\overrightarrow{r}}^{\prime }\text{,}{E}^{\prime }\to \overrightarrow{r}\text{,}E)s_0({\overrightarrow{r}}^{\prime }\text{,}{E}^{\prime })}{{\int }_0^{\infty }\text{d}{E}^{\prime }{\int }_{{\text{Z}}_j}{\text{d}}^3{r}^{\prime }{s}_0({\overrightarrow{r}}^{\prime }\text{,}{E}^{\prime })}\text{,}$$

where $s_0(\overrightarrow{r}\text{,}E)$ is the fundamental mode fission source. The fundamental mode eigenvalue of H equals $k_{\text{eff}}$, and the corresponding eigenvector equals the discretized fundamental mode fission source. Some Monte Carlo codes, e.g. TRIPOLI-4 (OECD/NEA, 2008) and KENO V.a (RSICC, 2006), can optionally calculate the fission matrix during the standard Monte Carlo calculations. Moreover, these codes can alternatively estimate $k_{\text{eff}}$ by the fundamental mode eigenvalue of the fission matrix. The computed fission matrix can, eventually, be also used for deriving the higher mode k-eigenvalues and eigensources, and for computing the dominance ratio of the modelled system.

The method presented in this paper allows to cancel the inactive cycles; it should not be mistaken with the fission matrix acceleration methods (Carter and McCormick, 1969, Kitada and Takeda, 2001) that aim at accelerating the convergence of the Monte Carlo fission source during the inactive cycles. To prevent possible confusions, we suggest the calculations made by the presented method be referred to as the “fission matrix based Monte Carlo criticality calculations”.

Section 2 describes the theory of the fission matrix based Monte Carlo criticality calculations; Section 2.1 provides the basic principles, Section 2.2 describes computation of the fission matrix and its variance, and Section 2.3 suggests sampling the semi-fixed fission source in the criticality safety calculations. Section 3 numerically demonstrates the method on the “k-effective of the world” model problem (Whitesides, 1971). Our conclusions are drawn in Section 4.

In the following text, the operator $〈·〉$ is also applied on continuous functions; e.g.

$$〈s(\overrightarrow{r}\text{,}E)〉\equiv {\int }_0^{\infty }\text{d}E{\int }_{\text{S}}{\text{d}}^3\mathit{rs}(\overrightarrow{r}\text{,}E)\text{.}$$

Next, let the operator $〈·{〉}_j$ be defined only on zone j; e.g.

$$〈s(\overrightarrow{r}\text{,}E){〉}_j\equiv {\int }_0^{\infty }\text{d}E{\int }_{{\text{Z}}_j}{\text{d}}^3\mathit{rs}(\overrightarrow{r}\text{,}E)\text{.}$$

Similarly, let $〈\mathcal{S}{〉}_j$ give the total statistical weight of fission neutrons from $\mathcal{S}$ that are located in zone j. We refer $〈\mathcal{S}{〉}_j\text{and}〈s(\overrightarrow{r}\text{,}E){〉}_j$ to as the intensity of $\mathcal{S}$ and $s(\overrightarrow{r}\text{,}E)$ in zone j, respectively.