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 α-k 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 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 α-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.
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.
The MC α calculations are performed for 1000 active iterations on 10,000 sources per iteration. 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 (=) while the neutron is absorbed with the probability of (=). Considering that the time source generation is dominated by the second-group neutron by for a small value of , let us calculate the expected number of time source neutrons, , during the MC neutron transport simulations starting from a second-group source until it is absorbed. Because as given in Table 1, can be written as
When 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), is found to be closer to 1 with values of 0.11, 0.34, 0.57, 0.80, and 0.97 when and 0.28, 0.61, 0.80, 0.92, and 0.99 when as the subcriticality is deeper from 0.90 to 0.15. Because the α estimate from the α-k iteration method is updated iteratively and fluctuated due to a finite number of MC history simulations, a stability condition can be expressed from (18) as where is the fundamental mode alpha-eigenvalue estimated at the th fission source iteration. From Table 1, for the two-group infinite homogeneous problems is determined to be 231,967. Figure 1 shows for the cases of of 0.3 with different values of 0 and 2 on 107 histories per iteration starting from an initially guessed α, , of 185,568. From the figure, we can observe that becomes close to 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 of for 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.
3.2. 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 detectors used for Th-PE, Th-Gr, Th-Be, Th-HEU-PE, and NU-PE cores and one with two detectors for Th-HEU-5PE and Th-HEU-Gr-PE cores having the 14 MeV pulsed neutron generator. The effective multiplication factors, ’s, of the Th-PE, Th-Gr, Th-Be, Th-HEU-PE, NU-PE, Th-HEU-5PE, and Th-HEU-Gr-PE configurations are reported to be 0.00613, 0.00952, 0.00765, 0.58754, 0.50867, 0.85121, and 0.35473, respectively [16].
(a) Th-PE, Th-Gr, Th-Be, Th-HEU-PE, and NU-PE
(b) Th-HEU-5PE and Th-HEU-Gr-PE
The McCARD calculations are performed with continuous-energy cross section libraries produced from JENDL-4.0 for the seven cores with the 14 MeV pulsed neutron generator. The alpha-static calculations by the α iteration algorithm are conducted with 5000 active iterations on 10,000 sources per iteration. For the comparison, McCARD simulations for the PNS experiments are performed on 20 billion of 14 MeV neutron sources. Figure 3 shows flux tally results of the three detectors for the Th-Gr core and the two detectors for the Th-HEU-5PE core in 1 μs time bins from 0.0 ms to 4.0 ms. A time decay constant α can be obtained by fitting the tally results to the exponential function, , in the time range when the higher-mode fluxes sufficiently decay as [2]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.
(a) Th-Gr core
(b) Th-HEU-5PE core
(a) Th-Gr core
(b) Th-HEU-5PE core
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 detectors 2 and 3 for the Th-PE, Th-Gr, Th-Be, Th-HEU-PE, and NU-PE cores and 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.
4. 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was supported by the National Nuclear R&D Program through the National Research Foundation of Korea (NRF) funded by MSIP (Ministry of Science, ICT and Future Planning) (no. 2010-0018758).