Research Article  Open Access
Monte Carlo Alpha Iteration Algorithm for a Subcritical System Analysis
Abstract
The αk iteration method which searches the fundamental mode alphaeigenvalue via iterative updates of the fission source distribution has been successfully used for the Monte Carlo (MC) alphastatic 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 alphastatic 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 twogroup 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 thoriumloaded acceleratordriven 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 timedependent neutron transport equation and the alphastatic 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 timedependent 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 timestepwise 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 higherorder 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 alphastatic MC calculations, the αk iteration methods [5–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 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 Tripoli4 [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 alphastatic 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 twogroup infinite homogeneous problems. The applicability of the proposed method is tested for the thoriumloaded acceleratordriven 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 alphaeigenvalue from (11), we directly apply the power method [14] to (11):where is the iteration index. Equation (13) implies that the fundamentalmode alphaeigenvalue 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 alphaeigenvalue 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 fundamentalmode alphaeigenvalue 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. TwoGroup 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 twogroup infinite homogeneous medium problems. Table 1 shows twogroup 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 secondgroup 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 secondgroup 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 alphaeigenvalue 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 alphaeigenvalue estimated at the th fission source iteration. From Table 1, for the twogroup 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 10^{7} 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 ThADS Experimental Benchmarks
The proposed α iteration method is applied for the experimental benchmarks on thoriumloaded acceleratordriven system (ThADS) [16]. The ThADS experiments were performed using the solidmoderated and solidreflected typeA 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 ThPE, ThGr, ThBe, ThHEUPE, and NUPE cores and one with two detectors for ThHEU5PE and ThHEUGrPE cores having the 14 MeV pulsed neutron generator. The effective multiplication factors, ’s, of the ThPE, ThGr, ThBe, ThHEUPE, NUPE, ThHEU5PE, and ThHEUGrPE configurations are reported to be 0.00613, 0.00952, 0.00765, 0.58754, 0.50867, 0.85121, and 0.35473, respectively [16].
(a) ThPE, ThGr, ThBe, ThHEUPE, and NUPE
(b) ThHEU5PE and ThHEUGrPE
The McCARD calculations are performed with continuousenergy cross section libraries produced from JENDL4.0 for the seven cores with the 14 MeV pulsed neutron generator. The alphastatic 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 ThGr core and the two detectors for the ThHEU5PE 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 highermode 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 ThPE, ThGr, ThBe, ThHEUPE, NUPE, ThHEU5PE, and ThHEUGrPE cores is determined to be 1.3, 1.6, 1.6, 1.7, 1.7, 1.6, and 1.3 ms, respectively.
(a) ThGr core
(b) ThHEU5PE core
(a) ThGr core
(b) ThHEU5PE 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 ThPE, ThGr, ThBe, ThHEUPE, and NUPE cores and detectors 1 and 2 for the ThHEU5PE and ThHEUGrPE 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 ThPE, ThBe, ThHEU5PE, and ThHEUGrPE cores are less than 1.0% while being 1.3%, 2.1%, and 3.4% for the ThGr, ThHEUPE, and NUPE 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 PhHEUPE, ThHEU5PE, and PhHEUGrPE cores.

4. Conclusion
A new MC alphastatic calculation method is developed to stably estimate the fundamentalmode alphaeigenvalue 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 iterationbyiteration. 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 twogroup 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 ThADS 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. 20100018758).
References
 D. E. Cullen, C. J. Clouse, R. Procassini, and R. C. Little, “Static and dynamic criticality: are they different?” Report UCRLTR201506, Lawrence Livermore National Laboratory, Livermore, Calif, USA, 2003. View at: Google Scholar
 T. Suzaki, “Subcritically determination of lowenriched UO_{2} lattices in water by exponential experiment,” Journal of Nuclear Science and Technology, vol. 28, no. 12, pp. 1067–1077, 1991. View at: Publisher Site  Google Scholar
 T. Yamamoto, K. Sakurai, T. Arakawa, and Y. Naito, “Accurate estimation of subcriticality using indirect bias estimation method, (I): theory,” Journal of Nuclear Science and Technology, vol. 34, no. 5, pp. 454–460, 1997. View at: Publisher Site  Google Scholar
 D. E. Cullen, “TART 2002: a coupled neutronphoton, 3D, combinatorial geometry, time dependent Monte Carlo transport code,” Tech. Rep. UCRLID126455, Rev. 4, Lawrence Livermore National Laboratory, 2002. View at: Google Scholar
 D. Brockway, P. Soran, and P. Whalen, “Monte Carlo α calculation,” Tech. Rep. LAUR851224, Los Alamos National Laboratory, 1985. View at: Google Scholar
 T. Ye, C. Chen, W. Sun, B. Zhang, and D. Tian, “Prompt time constants of a reflected reactor,” in Proceedings of the Symposium on Nuclear Data, T. Fukahori, Ed., Ibarakiken, Japan, January 2007. View at: Google Scholar
 S. D. Nolen, T. R. Adams, and J. E. Sweezy, “Integral criticality estimators in MCATK(U),” Tech. Rep. LAUR1225458, Los Alamos National Laboratory, Los Alamos, NM, USA, 2012. View at: Google Scholar
 A. Zoia, E. Brun, and F. Malvagi, “Alpha eigenvalue calculations with TRIPOLI4,” Annals of Nuclear Energy, vol. 63, pp. 276–284, 2014. View at: Publisher Site  Google Scholar
 A. Zoia, E. Brun, and F. Malvagi, “A Monte Carlo method for prompt and delayed alpha eigenvalue calculations,” in Proceedings of the International Conference on the Role of Reactor Physics Toward a Sustainable Future (PHYSOR '14), CDROM, Kyoto, Japan, SeptemberOctober 2014. View at: Google Scholar
 J. Lieberoth, “A Monte Carlo technique to solve the static eigenvalue problem of the boltzmann transport equation,” Nukleonik, vol. 11, article 213, 1968. View at: Google Scholar
 T. Yamamoto and Y. Miyoshi, “An algorithm of α and γmode eigenvalue calculations by Monte Carlo method,” in Proceedings of the 7th International Conference on Nuclear Criticality Safety (ICNC '03), JAERIConf 2003019, Japan Atomic Energy Research Institute, Tokai, Japan, October 2003. View at: Google Scholar
 T. Yamamoto, “Higher order α mode eigenvalue calculation by Monte Carlo power iteration,” Progress in Nuclear Science and Technology, vol. 2, pp. 826–835, 2011. View at: Google Scholar
 O. Petit, F.X. Hugot, Y.K. Lee, C. Jouanne, and A. Mazzolo, TRIPOLI4 Version 4 User Guide, CEAR6169, CEA Saclay, 2008.
 S. Nakamura, Computational Methods in Engineering and Science: With Applications to Fluid Dynamics and Nuclear Systems, WileyInterscience, 1977.
 H. J. Shim, B. S. Han, S. J. Jong, H. J. Park, and C. H. Kim, “McCard: Monte Carlo code for advanced reactor design and analysis,” Nuclear Engineering and Technology, vol. 44, no. 2, pp. 161–176, 2012. View at: Publisher Site  Google Scholar
 C. H. Pyeon, Experimental Benchmarks on ThoriumLoaded AcceleratorDriven System at Kyoto University Critical Assembly, KURRTR(CD)48, Research Reactor Institute, Kyoto University, Kyoto, Japan, 2015.
 G. I. Bell and S. Glasstone, Nuclear Reactor Theory, Van Nostrand Reinhold Company, New York, NY, USA, 1970.
 J. J. Duderstadt and L. J. Hamilton, Nuclear Reactor Analysis, John Wiley & Sons, New York, NY, USA, 1976.
 I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations, CRC Press, Boca Raton, Fla, USA, 1991.
 J. F. Briesmeister, “MCNP. A general Monte Carlo Nparticle transport code, version 4B,” Tech. Rep. LA13181, Los Alamos National Laboratory, Los Alamos, NM, USA, 1997. View at: Google Scholar
Copyright
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.