Monte Carlo Alpha Iteration Algorithm for a Subcritical System Analysis

The α-k iteration method which searches the fundamental mode alpha-eigenvalue via iterative updates of the fission source distribution has been successfully used for theMonte Carlo (MC) alpha-static calculations of supercritical systems. However, the αk iterationmethod 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 newMCalpha-static calculationmethod named as the α iteration algorithm.The newmethod 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.


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 alphastatic 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 as where Φ is the neutron angular flux and the subscript indicates ignorance of the delayed fission neutron. is named as the time source distribution. V( ) is a neutron speed corresponding to its energy . Other notations follow convention.
For the alpha-static MC calculations, the -k iteration methods [5][6][7][8][9] searching the fundamental mode alphaeigenvalue 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 − /V for prompt supercritical systems ( < 0). The -k iteration methods for prompt subcritical systems ( > 0) treat a negative value of the time absorption as the time production ( /V) 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 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 − /VΦ ( > 0) 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 -k iteration methods have been implemented in the Seoul National University MC code, McCARD [15]. The effectiveness of the new method is examined by comparing 's calculated by the iteration and the -k 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 's calculated by the iteration algorithm with those from McCARD dynamic simulations as well as measurements of the PNS experiments.

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: (r, , Ω) = ∫ r (r → r | , Ω) (r , , Ω) 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: (r → r | , Ω) = Σ (r, ) 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 Science and Technology of Nuclear Installations Multiplying V −1 Σ −1 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 R 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.

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 (= 1, 2, . . .), which is governed by (12) and (13), the time source neutrons for the next iteration can be sampled at the collision site as many as where 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 (0, 1). −1 is an alphaeigenvalue estimated at iteration − 1. When = 1, 0 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, (r , , Ω ), are set to the sampled sources' parameters. It is meaningful to compare (15) with a typical fission source site sampling denotes a -eigenvalue estimated at cycle −1. 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 −1 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 −1 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 + 1, the initial weight of the time source neutrons, +1 , is determined from the number of the sampled sources at iteration by where is a number of time source neutrons per iteration inputted by a code user.
To estimate ( = 1, 2, . . .) by (14), a collision estimator for an inverse of , that is, ∫ r ∫ ∫ ΩR −1 , sums up /{V( )Σ (r , )} at each collision. Then after finishing each iteration for the time source neutrons, ( = 1, 2, . . .) can be calculated from a mean of history results as In 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.

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  Table 2 shows comparisons of 's calculated by the new algorithm and the -k iteration method applying the pseudo absorption adjustment [6] with analytic solutions. In the table, the pseudo absorption adjustment parameter, , of zero means the conventional -k iteration method [1] without the pseudo absorption adjustment. From the table, one can see that the MC results from the iteration method agree well with the analytic references within 95% confidence intervals while the -k iteration method fails when are 0.3 and 0.15 due to abnormal terminations.
These abnormal terminations of the -k iteration method can be explained by counting the number of time source neutrons generated from a source neutron. In the -k iteration algorithm with the parameter, a time source neutron is generated at each collision site with the probability of (={(1 + ) /V }/(Σ + /V )) while the neutron is absorbed with the probability of (={Σ + /V −Σ }/(Σ + /V )). Considering that the time source generation is dominated by the second-group neutron by 2 ≫ 1 for a small value of , let us calculate the expected number of time source neutrons, time,2 , during the MC neutron transport simulations starting from a second-group source until it is absorbed. Because Σ 12 = 0 as given in Table 1, time,2 can be written as When time,2 is greater than or equal to 1, the total number of time source neutrons generated in a source history becomes infinity which leads to the abnormal termination. By inserting the analytic solution of the fundamental mode alpha-eigenvalue in Table 2 and the constants given in Table 1 into (18) number of MC history simulations, a stability condition can be expressed from (18) as whereis the fundamental mode alpha-eigenvalue estimated at the th fission source iteration. From Table 1, limit for the two-group infinite homogeneous problems is determined to be 231,967. Figure 1 showsfor the cases of of 0.3 with different values of 0 and 2 on 10 7 histories per iteration starting from an initially guessed , 0 -, of 185,568. From the figure, we can observe thatbecomes close to limit of 231,967 at cycle 8 in the run with of 0 and at cycle 6 with of 2, which yield the abnormal terminations at the very next iteration calculations.
The -k iteration method with the pseudo absorption term may reduce occurrence of the infinite number of time source generations by suppressing fluctuations offor MC calculations with continuous cross section data, as reported in [12]. However, it is noteworthy that the abnormal terminations in the -k iteration calculations can happen by inevitable statistical fluctuations of estimates.

Application to the Th-ADS Experimental Benchmarks.
The proposed iteration method is applied for the experimental benchmarks on thorium-loaded accelerator-driven system (Th-ADS) [16]. The Th-ADS experiments were performed using the solid-moderated and solid-reflected type-A core of KUCA for seven core configurations combined with a 14 MeV pulsed neutron generator or a synchrotron type proton accelerator. The highly enriched uranium (HEU), thorium (Th), and natural uranium (NU) fuel was loaded together with the reflectors, including polyethylene (PE), graphite (Gr), and beryllium (Be). Figure 2 shows a core configuration with three 3 He detectors used for Th-PE, Th-Gr, Th-Be, Th-HEU-PE, and NU-PE cores and one with two Ref.
where and are fitting constants and and are the time after the neutron burst and the starting time of fitting, respectively. Figure 4 shows the convergence of to a stable one with increasing for the fixed fitting intervals of 800 s. From these convergence diagnoses, for the Th-PE, Th-Gr, Th-Be, Th-HEU-PE, NU-PE, Th-HEU-5PE, and Th-HEU-Gr-PE cores is determined to be 1.3, 1.6, 1.6, 1.7, 1.7, 1.6, and 1.3 ms, respectively. Table 3 shows comparisons of 's estimated by the iteration method with the measurements given in [16] and the values from the McCARD PNS simulations which are calculated by averaging fitted values of 3 He detectors 2 and 3 for the Th-PE, Th-Gr, Th-Be, Th-HEU-PE, and NU-PE cores and 3 He detectors 1 and 2 for the Th-HEU-5PE and Th-HEU-Gr-PE cores. From the comparisons between the values calculated by the iteration method and the McCARD PNS simulations, we can see that the absolute values of the relative difference for the Th-PE, Th-Be, Th-HEU-5PE, and Th-HEU-Gr-PE cores are less than 1.0% while being 1.3%, 2.1%, and 3.4% for the Th-Gr, Th-HEU-PE, and NU-PE cores, respectively. From the comparisons with measurements, we can see that the iteration results are quite comparable with   the experimental results for core configurations where two estimates from different detectors show good agreements, that is, the Ph-HEU-PE, Th-HEU-5PE, and Ph-HEU-Gr-PE cores.

Conclusion
A new MC alpha-static calculation method is developed to stably estimate the fundamental-mode alpha-eigenvalue of prompt subcritical systems. While the fission source distribution is updated in the existing -k iteration method, the new method named as the iteration algorithm directly updates the time source distribution iteration-by-iteration. The calculation stability of the iteration algorithm is achieved by controlling the number of time source neutrons generated at each iteration by the normalization scheme for the time source distribution.
It is demonstrated that the iteration algorithm does not suffer from the instability problem in the two-group homogeneous problems with deep subcriticalities where the -k iteration method fails to obtain the value due to the huge number of neutron productions. From the comparisons with measurements in the Th-ADS experimental benchmarks [16], it is observed that the results calculated by the iteration method are quite comparable with each other for the cases where the experiments provide reliable estimates.