Monte Carlo K algorithms

The source iteration and fission matrix algorithms for Monte Carlo criticality calculations are discussed below using the collision estimator. For simplicity of notation, we assume a homogeneous system made up of a single isotope material that has cross sections $\bar\nu\Sigma_f$, $\Sigma_t$ and $\Sigma_s$.

Source Iteration Algorithm:
The source iteration method is formulated by solving the source equation (equation (2.15)),

$$S^n(\vec r) = \int S^{(n-1)}(\vec r') A(\vec r,\vec r')dr'$$ (21)

and calculating K as the ratio of source neutrons of two successive neutron generations [Uss58, Men68] as follows:

$$K = \frac{\int S^n(\vec r)dr}{\int S^{(n-1)}(\vec r)dr}.$$ (22)

$S^n(\vec r)$ is the source distribution of the n^th generation of neutrons and $A(\vec r,\vec r')$ describes the probability that a neutron originating at $\vec r'$ produces a fission neutron at $\vec r$. The Monte Carlo source iteration algorithm is given below:

1.

N source neutrons are started with an initial guessed spatial distribution. Each neutron has a starting weight W = 1.0.

2.

At every collision point where fission is possible we perform the following operations.

(a)

Accumulate the tally for the collision estimator,

$$TALLY_j = TALLY_j + W_i {{\bar\nu\Sigma_f}\over {\Sigma_t}}$$

where,
i = each collision by the neutron, and
j = all neutrons, j = 1,2,....,N.

(b)

We record the location and the number (n'_i) of source neutrons to be started from that location (where collision occurred) for the next generation:

$$n'_i = INTEGER({{W_i\bar\nu\Sigma_f}\over {\Sigma_t}}\frac{1}{K} + \xi) ,$$

where $\xi$ is a random number between 0 and 1, and K is the last batch's eigenvalue or a guessed one for the first batch.

$$n'_j = \sum\limits_{i}n'_i.$$

(c)

For survival biasing, we set the new weight after the i^th collision, W_i+1, to,

$$W_{i+1} = W_i\frac{\Sigma_s}{\Sigma_t}.$$

(d)

We also apply the Russian roulette technique [Lew93] to eliminate low weight neutrons.

3.

After all the N neutron histories have been tracked, the eigenvalue for that generation or batch is calculated as;

$$K = {{\sum\limits_{j=1}^{N}{TALLY_j}}\over {N}}.$$

4.

The next batch is started with N' source neutrons, where

$$N' = \sum\limits_{j}n'_j.$$

For this next batch each source particle has a starting weight,

$$W = \frac{N}{N'}.$$

This completes the source iteration algorithm. It should be noted that the number of neutrons, started in each batch, times the weight of each neutron, for that batch, is constant for the all batches, i.e., W*N' is same for all the batches. N' is close to N due to the K normalization done in step 2(b). If the initial guess of K is too low or too high, the number of fission neutrons produced for the next batch will be too high or too low relative to the desired nominal number N. The number of fission neutrons, n_i' stored at each collision is rounded up or down (step 2(b)) to an integer (including zero) with a probability proportional to its deviation from that integer.

Fission Matrix Algorithm:
The fission matrix method [Kap58, Mih67] solves the eigenvalue equation

$$KS(\vec r) = \int S(\vec r') A(\vec r,\vec r') dr',$$ (23)

which holds for every generation. Equation (2.23) can be discretized and is then reduced to the matrix equation,

KS = AS ,

where S is the source eigenvector and A is the fission matrix. The element A(l,m) of the fission matrix describes the probability that a neutron starting in cell m will generate a fission neutron in cell n [Ave58]. Assuming that the physical space is discretized into M cells, the Monte Carlo fission matrix algorithm can be described as:

1.

N source neutrons are started with an initial guessed spatial distribution. Each neutron has a starting weight W = 1.0. An array C(m) keeps track of the number of neutrons that start from cell m where m = 1,2,...,M.

2.

At every collision point where fission is possible, the following operations are performed:

(a)

If the collision occurs in cell l, accumulate into matrix B(l,m) the contribution to K due to a neutron that started from cell m, i.e.,

$$B(l,m) = B(l,m) + W_i {{\bar\nu\Sigma_f}\over {\Sigma_t}}$$

where
i = every collision that a neutron encounters and $\bar\nu\Sigma_f$ and $\Sigma_t$ are macroscopic cross sections for cell l.

(b)

Record the location and the number (n'_i) of source neutrons to be started from that location (where collision occurred) for the next generation.

$$n'_i = INTEGER({{W_i\bar\nu\Sigma_f}\over {\Sigma_t}}\frac{1}{K} + \xi) ,$$

where $\xi$ is a random number between 0 and 1, and K is the previous batch eigenvalue (or a guessed one for the first batch).

$$n'_j = \sum\limits_{i}n'_i.$$

(c)

For survival biasing, adjust the new weight after the i^th collision,

$$W_{i+1} = W_i\frac{\Sigma_s}{\Sigma_t}.$$

(d)

Apply the Russian roulette technique [Lew93] to eliminate low weight neutrons.

3.

After all the N neutron histories have been tracked, compute the fission matrix A(l,m) element,

$$A(l,m) = \frac{B(l,m)}{C(m)}.$$

The element A(l,m) gives the expected number of fission neutrons produced in cell l due to a unit source neutron in cell m. The dominant eigenvalue, or K, of the fission matrix A can be calculated by using a matrix iterative algorithm.

4.

Start the next batch with N' source neutrons, where

$$N' = \sum\limits_{j}n'_j,$$

where each source particle has a starting weight,

$$W = \frac{N}{N'}.$$

This completes the fission matrix algorithm. As with source iteration, the quantity W*N' is preserved for all batches. There are two ways to formulate the fission matrix algorithm. These two formulations are explained in detail in section 2.5. The numerical examples of section 2.4 are solved using the cycle fission matrix approach, where the elements of the matrix B(l,m) and the array C(m) are zeroed out at the begining of every batch.

The following information is provided to a Monte Carlo eigenvalue simulation algorithm:

• the nominal number of source histories per batch;
• the initial guess of K;
• the number of initial batches, I_i, to be discarded before tally accumulation;
• the total number of batches, I_t, for the problem; and
• an initial spatial distribution for the source.

The number of active batches, I_a that are used to evaluate the mean and the standard deviation for K is then given by,

I_a = I_t - I_i.