Science and Technology of Nuclear Installations

Volume 2015, Article ID 859242, 7 pages

http://dx.doi.org/10.1155/2015/859242

## Monte Carlo Alpha Iteration Algorithm for a Subcritical System Analysis

Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 151-744, Republic of Korea

Received 14 April 2015; Revised 21 June 2015; Accepted 29 June 2015

Academic Editor: Valerio Giusti

Copyright © 2015 Hyung Jin Shim et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The *α*-*k* iteration method which searches the fundamental mode alpha-eigenvalue via iterative updates of the fission source distribution has been successfully used for the Monte Carlo (MC) alpha-static calculations of supercritical systems. However, the *α*-*k* iteration method for the deep subcritical system analysis suffers from a gigantic number of neutron generations or a huge neutron weight, which leads to an abnormal termination of the MC calculations. In order to stably estimate the prompt neutron decay constant (*α*) of prompt subcritical systems regardless of subcriticality, we propose a new MC alpha-static calculation method named as the *α* iteration algorithm. The new method is derived by directly applying the power method for the *α*-mode eigenvalue equation and its calculation stability is achieved by controlling the number of time source neutrons which are generated in proportion to *α* divided by neutron speed in MC neutron transport simulations. The effectiveness of the *α* iteration algorithm is demonstrated for two-group homogeneous problems with varying the subcriticality by comparisons with analytic solutions. The applicability of the proposed method is evaluated for an experimental benchmark of the thorium-loaded accelerator-driven system.

#### 1. Introduction

The prompt neutron decay constant (referred to as *α*) in subcritical systems is a basic neutronics parameter which can be directly measured from reactor experiments. In the Monte Carlo (MC) neutron transport analysis, *α* has been estimated by two approaches [1]: dynamic MC simulations to directly solve the time-dependent neutron transport equation and the alpha-static MC calculations to solve the *α*-mode eigenvalue equation. The dynamic MC calculations cover the *α* estimation from numerical simulations of the pulsed neutron source (PNS) experiment by fitting time-dependent tally results of a neutron detector to the exponential function. However, it is well known that a starting time of fitting should be carefully determined by ensuring the convergence of estimated values [2, 3]. The TART code [4] is equipped with a unique method to measure *α* in time-stepwise MC simulations with controlling the neutron population by the combing algorithm [1].

Differently from the MC dynamic simulations, the alpha-static MC methods calculate the fundamental mode or higher-order solutions of the *α*-mode eigenvalue equation for prompt neutron which can be expressed in operator notation aswhere is the neutron angular flux and the subscript indicates ignorance of the delayed fission neutron. is named as the time source distribution. is a neutron speed corresponding to its energy . Other notations follow convention.

For the alpha-static MC calculations, the * α-k* iteration methods [5–9] searching the fundamental mode alpha-eigenvalue in the middle of iterative fission source updates by the MC power iteration method (

*k*iteration method hereafter) [10] have been mainly used since Brockway et al. [5] developed an MC algorithm with considering the time absorption [1] of for prompt supercritical systems (). The

*iteration methods for prompt subcritical systems () treat a negative value of the time absorption as the time production ( by which the time source neutrons are generated in the MC neutron simulations to avoid negative absorption reactions. Yamamoto [11] suggested a neutron simulation algorithm with correcting the neutron weight in each track by the time production reaction. However, the time production strategy of the conventional*

*α*-k*-*

*α**k*iteration method [1] and Yamamoto’s weight correction method for highly subcritical systems are reported [6, 8, 12] to cause a gigantic number of source neutron generations or a huge value of the neutron weight which leads to an abnormal termination of the MC calculations. To reduce this instability problem, Ye et al. [6] introduced a pseudo neutron absorption term of into both sides of (1) and its slightly modified formulation [8] is used in Tripoli-4 [13]. However, an appropriate adjustment parameter, , is required to ensure a stable running without a halt in the pseudo neutron absorption approaches.

To overcome this instability problem in the * α-k* iteration methods, we propose a new alpha-static MC method named as the

*α*iteration algorithm for the prompt subcritical system analysis. In Section 2, the

*α*iteration algorithm is derived by applying the power method [14] for the

*α*-mode eigenvalue equation of (1) with a normalization scheme to stably control the number of time source neutrons at each iteration. The new

*α*iteration and the existing

*iteration methods have been implemented in the Seoul National University MC code, McCARD [15]. The effectiveness of the new method is examined by comparing*

*α*-k*α*’s calculated by the

*α*iteration and the

*iteration methods with analytic solutions for two-group infinite homogeneous problems. The applicability of the proposed method is tested for the thorium-loaded accelerator-driven system [16] at Kyoto University Critical Assembly (KUCA) by comparing*

*α*-k*α*’s calculated by the

*α*iteration algorithm with those from McCARD dynamic simulations as well as measurements of the PNS experiments.

#### 2.
*α* Iteration Algorithm

*α*

##### 2.1. Derivation

In order to obtain a MC neutron tracking algorithm from an integrodifferential form of the *α*-mode eigenvalue equation written by (1), it is required to transform (1) into its integral form. By integrating (1) along with a characteristic curve [17] in the same way to derive the integral form of the neutron transport equation [17, 18] and expressing the resulting integral equation for the collision density, , defined by [19], one can obtain the collision density equation for the *α*-mode eigenvalue equation:where is the transition kernel ignoring the delayed fission neutrons defined as a product of the collision kernel without the delayed fission neutron, , and the transport kernel, , as follows:where is used to index the neutron reaction except fission. is the average number of neutrons produced from a reaction type and is the probability that a collision of type by a neutron of direction and energy will produce a neutron in direction interval about with energy in about .

Then the Neumann series solution [19] to (5), which describes the MC neutron simulations, can be written by

Multiplying on both sides of (9), one can obtain the *α*-mode eigenvalue equation for the time source distribution as

To calculate the fundamental mode alpha-eigenvalue from (11), we directly apply the power method [14] to (11):where is the iteration index. Equation (13) implies that the fundamental-mode alpha-eigenvalue and are calculated by iterative updates of the times source distribution until they converge. It is notable that the transition kernels in operator defined by (12) are the same as the ones used in the fixed source mode MC calculations except that the delayed fission neutrons are not considered, which means that all the prompt fission neutrons should be simulated within an iteration. Note that an expected number of prompt fission neutrons generated within an iteration is less than 1 for prompt subcritical systems of which the prompt criticality is less than 1.

##### 2.2. Application to Monte Carlo Calculations

The derived *α* iteration method can be implemented by slightly modifying the iteration algorithm [10, 20] used in existing MC codes. The main difference between the *α* iteration and the iteration methods is that the time source distribution is updated in the *α* iteration algorithm while the fission source distribution in the iteration one. When a collision occurs in the course of the *α* iteration at iteration , which is governed by (12) and (13), the time source neutrons for the next iteration can be sampled at the collision site as many aswhere and are indices of time source neutrons and collisions, respectively. denotes the largest integer not exceeding . is the neutron weight for the th collision from the th time source neutron at iteration . is a uniform random number on the interval of (). is an alpha-eigenvalue estimated at iteration . When , indicates an initial *α* value guessed by a code user. When is greater than or equal to one, the location, energy, and direction of the collision neutron, , are set to the sampled sources’ parameters. It is meaningful to compare (15) with a typical fission source site sampling formulation of [20] where denotes a -eigenvalue estimated at cycle . Note that the direction of the collision neutron, , should be stored in the *α* iteration algorithm while the direction of a fission source neutron may not be banked in the iteration algorithm because of an assumption of its isotropic distribution. The multiplication of in the right hand side of (15) plays an important role in controlling the total number of time source neutrons per iteration like the division by in the iteration algorithm. As noted in the previous section, the fission reactions should be sampled in a routine of the collision type selection from the collision kernel defined by (7) in the *α* iteration algorithm, whereas the fission reaction is not allowed to occur in a common implementation of the iteration algorithm [20].

At the beginning of the next iteration , the initial weight of the time source neutrons, , is determined from the number of the sampled sources at iteration bywhere is a number of time source neutrons per iteration inputted by a code user.

To estimate by (14), a collision estimator for an inverse of , that is, , sums up at each collision. Then after finishing each iteration for the time source neutrons, can be calculated from a mean of history results asIn the same reason to apply the inactive cycle runs in the iteration method, the fundamental-mode alpha-eigenvalue should be estimated by averaging ’s at active iterations after proper number of inactive iterations to converge the time source distribution.

#### 3. Numerical Results

##### 3.1. Two-Group Infinite Homogeneous Problems

In order to investigate the effectiveness of the proposed method, the *α* and the * α-k *iteration algorithms implemented in McCARD [15] are tested for two-group infinite homogeneous medium problems. Table 1 shows two-group cross sections varying the prompt criticality . The differential scattering cross section of the first group, , is set to 0.265714, 0.197143, 0.128571, 0.060000, or 0.008571, which correspond to of 0.9, 0.7, 0.5, 0.3, or 0.15.