Advanced Method for Calculations of Core Burn-Up, Activation of Structural Materials, and Spallation Products Accumulation in Accelerator-Driven Systems

Abstract

The ALEPH2 Monte Carlo depletion code has two principal features that make it a flexible and powerful tool for reactor analysis. First of all, it uses a nuclear data library covering neutron- and proton-induced reactions, neutron and proton fission product yields, spontaneous fission product yields, radioactive decay data, and total recoverable energies per fission. Secondly, it uses a state-of-the-art numerical solver for the first-order ordinary differential equations describing the isotope balances, namely, a Radau IIA implicit Runge-Kutta method. The versatility of the code allows using it for time behavior simulation of various systems ranging from single pin model to full-scale reactor model, including such specific facilities as accelerator-driven systems. The core burn-up, activation of the structural materials, irradiation of samples, and, in addition, accumulation of spallation products in accelerator-driven systems can be calculated in a single ALEPH2 run. The code is extensively used for the neutronics design of the MYRRHA research facility which will operate in both critical and subcritical modes.

1. Introduction

SCK·CEN, the Belgian Nuclear Research Centre in Mol, is designing a Multipurpose Hybrid Research Reactor for High-tech Applications (MYRRHA) [1, 2]. The Accelerator-Driven System (ADS) concept has been chosen as a basis for this reactor, assuming that it can operate in both subcritical and critical modes (clearly without the accelerator running). The 600 MeV, 4 mA proton accelerator is coupled with the MOX-fueled core. The lead-bismuth eutectic (LBE) serves as coolant and spallation target to generate the source neutrons. To perform the neutronics design of such a system, dedicated codes must be used. The majority of neutronics codes used in the analysis of critical systems cannot be applied to calculations of ADS since multiparticle physics extending far beyond the energy range typical for nuclear reactors is involved.

The growth of computing power caused sustainable increase of the share of Monte Carlo codes in nuclear reactor and nuclear criticality research and development. These Monte Carlo codes can provide the most accurate locally dependent neutronics characteristics in realistic 3D geometries of any complexity. Among them, only the general-purpose radiation transport code MCNPX [3, 4] is fully capable to treat ADS-related problems, since it tracks almost all particle types of nearly all energies. In the upper energy region (above 20 MeV) it relies on the model calculations using different intranuclear cascade, preequilibrium and equilibrium model combinations. However, recent progress in nuclear data library extensions to higher energies (up to 200 MeV) allows implementing evaluated nuclear data in neutronics analysis of ADS systems. These evaluated data are generally more precise than model calculations. Thus MCNPX coupled with nuclear data libraries is a powerful tool to perform steady-state particle transport calculations.

However, it has limited applicability to follow the evolution of the system neutronics characteristics in time. Its depletion module based on the CINDER90 code [5] can only be applied to fuel burn-up calculations in criticality (KCODE) mode and not in the fixed source mode necessary to model ADS. In addition, the material activation and spallation products accumulation cannot be treated at all thus external depletion codes are required. Among such depletion codes oriented to calculate time behavior of nuclide inventories, the SNT code [6] and EASY-2007 code system [7] can be mentioned. The ADS core burn-up is reported to be studied with EVOLCODE [8, 9] and MCB [10] codes. However, to our knowledge, none of these codes is capable to simulate the behavior of ADS neutronics in time.

There is a noticeable progress in the development of reactor-oriented Monte Carlo depletion codes such as those mentioned above: MCNPX with CINDER’90 depletion capability, MCB, as well as Serpent [11], MURE [12], ACAB [13], Monteburns [14], Casmo/Simulate [15], and other codes. The common principle of these codes is that a steady-state Monte Carlo calculation of neutron fluxes and spectra is followed by solving the system of first-order ordinary differential equations (ODEs) for a given time step using the reaction rates calculated from these fluxes and spectra. Then the nuclide concentrations are updated and passed to the Monte Carlo code for recalculation of fluxes and spectra, and so forth. This procedure repeats up to the end of the irradiation and/or decay history. As a rule, these codes wrap around general-purpose Monte Carlo radiation transport codes (MCNP [16], MCNPX) and depletion codes (ORIGEN [17, 18], CINDER’90). However, some of them use built-in depletion modules (MCB, Serpent). The method of coupling is often realized as a set of scripts which complicates its utilization by nonexperienced users.

The general-purpose burn-up code ALEPH merging MCNP(X) (any version of MCNP or MCNPX) Monte Carlo radiation transport and ORIGEN-2.2 [17] depletion codes is being developed at SCK·CEN since 2004 [19, 20]. Belonging to the same category of shells coupling Monte Carlo transport and “deterministic” depletion codes, ALEPH possess some unique features that distinguish it from other codes. The most important feature is full consistency of cross-section data. The same unionized cross section tables (i.e., cross sections are linearized on the same energy grid) for a given nuclide are used for Monte Carlo transport and subsequent depletion calculations. The reaction rates are calculated at the beginning of each time step using the fluxes and spectra provided by MCNP(X). ALEPH allows the user to change materials, temperature, and geometry after each time step to reflect the irradiation conditions. Finally, the code is easy to use since it requires only several extra cards in the MCNP(X) input deck. The validity of ALEPH has been confirmed for different types of problems [21–23].

However, since ALEPH invokes ORIGEN-2.2 to perform time evolution calculations, it retains all limitations inherent to ORIGEN-2.2. This concerns the accuracy of matrix exponential method to solve ODEs and the limited number of nuclear data involved. The application of the code is limited to (classical) reactor systems.

The new version of the code, ALEPH2 [24], is free from these limits and possess some other outstanding features which altogether make it a powerful and flexible tool to perform time evolution analysis of various nuclear systems. It is fully applicable to ADS problems: core burn-up calculations, activation of structural materials, spallation product accumulation, and associated radiation source terms (decay heat release, heat generation by delayed photons, absorbed doses, neutron sources due to decay) can be calculated in a single code run which requires a few extra cards to MCNPX input deck. The features of the code especially important for ADS-related problems are described in Section 2 with the main focus on the method to perform depletion calculations. Section 3 describes the application of the code to the neutronics design of MYRRHA system and the conclusions are drawn in Section 4.

2. Basic Working Principles of ALEPH2

Two principal modifications have been made to ALEPH code: extending the nuclear data set and replacing ORIGEN-2.2 by built-in depletion module. Besides that, the predictor-corrector mechanism, calculation of heating during irradiation and decay, calculation of neutron sources due to decay, fuel management as well as some minor corrections to improve the code performance and facility of use were made.

2.1. Nuclear Data

One of two major improvements concerns nuclear data treatment. ORIGEN-2.2 is capable to treat only a limited number of reactions (radiative capture, fission, (𝑛,2𝑛),(𝑛,3𝑛),(𝑛,𝑝), and (𝑛,𝛼)). However, the ADS problems require more extensive nuclear data treatment. The number of open reaction channels can be significantly higher and increases drastically for the problems involving energies beyond those characteristic to nuclear reactors. The auxiliary utility, ALEPH-DLG [20] wrapping around the NJOY code [25], was modified to be able to build the nuclear data library containing nuclide production cross sections up to 1 GeV. The library can be constructed for a set of temperatures from basic general-purpose evaluated data files like JEFF-3.1.1 [26, 27], ENDF/B-VII [28], or JENDL-4 [29]. Then it is extended to a higher number of nuclides and reactions per target nuclide by adding the data from the TENDL-2010 library [30]. TENDL-2010 contains ~2000 neutron files for the nuclides with half-lives greater than one second. The upper energy limit in these files is either 200 MeV for the majority of stable and important unstable nuclides or 60 MeV for short-lived isotopes. Actinides, however, have upper limit 20 MeV. The final step is an extension to higher energies using the HEAD-2009 activation data library [31]. The HEAD-2009 library contains neutron- and proton-induced data for 684 stable and unstable nuclides from C to Po in the primary particle energy range 150 MeV to 1 GeV.

Proton-induced reactions can play an important role in investigating the time behavior of spallation and activation products in an ADS. The capability to treat proton data was added to ALEPH. To summarize, the new library constructed by ALEPH-DLG consists of ~2000 neutron files and ~1200 proton files for individual nuclides with half-lives greater than a second. Files contain up to ~2000 reaction cross sections with upper energy limit 1 GeV.

Since the number of reactions for an individual nuclide can be rather high, significant efforts were directed to accelerate the cross section averaging over the spectrum produced by MCNP(X). The format of activation data storing (PENDF tapes generated by NJOY code) was found to be ineffective. It was changed to ACE (A Compact ENDF) dosimetry type [3, 4, 16] where only reaction cross sections are recorded. Tests have shown that the total time to get one-group averaged cross sections for most complex problem (involving all reactions for all nuclides in the library, i.e., about 460,000 neutron and proton induced reactions) is comparable to the time to process ~3000 neutron reactions for most complex problem handled by the previous version of ALEPH.

The previous version of ALEPH was updating only cross sections in the ORIGEN library at the beginning of each time step; moreover, the cross sections of reactions leading to metastable states remained unchanged. All other information, such as neutron fission product yields and radioactive decay data, was not updated. The ALEPH2 code is capable to use almost all special-purpose data supplied in basic libraries, namely, radioactive decay data (e.g., JEFF-3.1.1 radioactive decay file contains information for 3851 nuclides), total recoverable energy per fission of fissile nuclides, spontaneous neutron fission product yields, and direct neutron and proton fission product yields.

2.2. Depletion Algorithm

Assuming that the particle fluxes and spectra remain constant during the time step and locally in selected volume, the concentrations (atom densities) of nuclides 𝑦_𝑖 at the end of irradiation step can be obtained by solving the system of first-order linear differential equations with constant coefficients:𝑑𝑦_𝑖(𝑡)𝑑𝑡=_𝑗𝜆^𝑡𝑟_𝑗𝑖𝑦_𝑗(𝑡)−𝜆_𝑖^𝑡𝑟𝑦_𝑖(𝑡), (1) or, in compact matrix form,𝑑𝑦⃗𝑑𝑡=𝚲𝑦,⃗ (2) with the matrix coefficients being 𝜆_𝑖𝑗=𝜆^𝑡𝑟_𝑗𝑖−𝜆_𝑖^𝑡𝑟𝛿_𝑖𝑗. Here 𝜆^𝑡𝑟_𝑗𝑖 is the transmutation (production) rate of nuclide 𝑖 from nuclide 𝑗, and 𝜆_𝑖^𝑡𝑟 is the disappearance rate of nuclide 𝑖. These rates are defined as follows:𝜆^𝑡𝑟_𝑗𝑖=_𝑚𝜎^𝑚_𝑗𝑖𝜑^𝑚+𝜆^𝑑_𝑗𝑖,𝜆_𝑖^𝑡𝑟=_𝑚𝜎^𝑚_𝑅,𝑖𝜑^𝑚+𝜆^𝑑_𝑖. (3)(4) Here 𝜎^𝑚_𝑗𝑖 is the energy-averaged production cross section of nuclide 𝑖 from nuclide 𝑗 in the reaction with incident particle 𝑚, 𝜑^𝑚 is the average flux of type 𝑚 particles, 𝜆^𝑑_𝑗𝑖 stands for the decay rate of nuclide 𝑗 into the nuclide 𝑖, 𝜎^𝑚_𝑅,𝑖 represents the energy-averaged reaction cross section of nuclide 𝑖 transforming it into another nuclide, and 𝜆^𝑑_𝑖 is the decay constant of nuclide 𝑖. In most cases 𝑚=neutrons since only neutron reactions are considered, but when dealing with spallation products accumulation and activation of structural materials close enough to spallation target of accelerator driven system, the contribution of proton-induced reactions has to be taken into account.

As it was mentioned before, the previous version of ALEPH used the ORIGEN-2.2 depletion code [17] to obtain the nuclide concentrations at the end of irradiation step. The essence of the algorithm used by ORIGEN-2.2 is the matrix exponential method where the solution of (2) is obtained by truncating the power series for the exponential. The short-lived nuclides are assumed to decay instantly and are removed from the system. They are treated separately by a Gauss-Seidel iterative technique under the assumption of secular equilibrium of decay and production chains. The same method is used by the successor of ORIGEN-2.2, ORIGEN-S which is a part of the SCALE code system [18]. It is clear that the applied series expansion has a limited accuracy due to truncation and round-off errors. Moreover, due to some inherent workings of ORIGEN-2.2 it cannot be applied to the problems involving particle types other than neutron.

Modern computers, however, allow applying more precise methods of solving the systems of stiff first-order differential equations with constant coefficients. The common practice is to use the backward differentiation formulas (BDF methods) or implicit Runge-Kutta (IRK) methods. The IRK methods offer some important advantages over BDF. Being single step, the IRK methods do not suffer from systems discontinuities and changes in the order or time step. The IRK methods are also self-starting which means that the only starting values required are the initial conditions. However, the drawback is leading after the time discretization to larger and more complex sets of linear equations. Numerous tests have shown that the IRK method Radau IIA (3 stages, accuracy order 5) as implemented in RADAU5 solver [32] is more stable and faster than the BDF method applied in DLSODA solver [33, 34]. As a consequence the Radau IIA algorithm has been incorporated in ALEPH2. The details of implementation are given below.

2.2.1. Implementation

Implicit Runge-Kutta methods have been known to have a high order of solution accuracy and excellent stability properties in solving stiff ordinary differential equations (ODEs) [32, 35]:𝑑𝐲(𝑡)((𝑡𝑡𝑑𝑡=𝐟𝑡,𝐲));𝐲₀=𝐲_𝟎. (5)𝐲 is unknown 𝑁-size vector of values (concentrations in case of depletion equations), and 𝐟(𝑡,𝐲(𝑡)) is the 𝑁-size vector of functions of known form. In case of depletion equations, it is simply 𝐟(𝑡,𝐲(𝑡))=Λ𝐲(𝑡), and Λis the effective reaction rate matrix.

If the time step is subdivided into 𝑚 substeps with a fixed step size𝑡ℎ=¹−𝑡₀𝑚 (6) with the mesh points 𝑡_𝑛,𝑛=0,1,…,𝑚,𝑡_𝑚=𝑡₁, the implicit Runge-Kutta method is defined by𝑦_𝑛+1=𝑦_𝑛+ℎ^𝑠_𝑖=1𝑏_𝑖𝑓𝑡_𝑛+𝑐_𝑖ℎ,𝑔_𝑖. (7) Here𝑔_𝑖=𝑦_𝑛+ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡_𝑛+𝑐_𝑗ℎ,𝑔_𝑗,𝑖=1,…,𝑠. (8) The stage number is denoted by 𝑠 and 𝑎_𝑖𝑗,𝑏_𝑖, and 𝑐_𝑖 are parameters satisfying the conditions:𝑎_𝑖𝑗≠0,^𝑠_𝑗=1𝑎_𝑖𝑗=𝑐_𝑖for𝑖=1,2,…,𝑠,^𝑠_𝑖=1𝑏_𝑖=1. (9) There are different IRK methods based on different sets of above parameters, the Radau IIA method being one of them. In this case, the parameters 𝑐_𝑖 for 𝑖=1,2,…,𝑠 are selected as the zeros of Radau polynomial (𝑑^𝑠−1/𝑑𝑡^𝑠−1)(𝑡^𝑠−1(𝑡−1)^𝑠). Parameters 𝑎_𝑖𝑗 are obtained by solving the linear system^𝑠_𝑗=1𝑎_𝑖𝑗𝑐_𝑗^𝑞−1=𝑐^𝑞_𝑖𝑞,𝑖=1,…,𝑠;𝑞=1,…,𝑠 (10) and 𝑏_𝑖=𝑎_𝑠,𝑖 for 𝑖=1,2,…,𝑠. Stage number 𝑠 can be set to 1,3,5, or 7. The RADAU5 solver uses 𝑠=3 thus its accuracy order is 2𝑠−1=5 and it is considered as absolutely stable.

At each sub-step, the system (3) of rank 𝑠𝑁 has to be solved iteratively for the stage values {𝑔_𝑖}^𝑠_𝑖=1. To reduce round-off errors, it is often reformulated as𝑧_𝑖=𝑔_𝑖−𝑦_𝑛, (11) where𝑧_𝑖=ℎ^𝑠_𝑗=1𝑎_𝑖𝑗𝑓𝑡₀+𝑐_𝑗ℎ,𝑦_𝑛+𝑧_𝑗,𝑖=1,…,𝑠. (12) Since the matrix of coefficients 𝐴=(𝑎_𝑖𝑗) is nonsingular (as it can be seen from (10)), (12) can be rewritten in vector form𝑧₁...𝑧_𝑠𝑡=𝐴ℎ𝑓_𝑛+𝑐₁ℎ,𝑦_𝑛+𝑧₁...ℎ𝑓𝑡_𝑛+𝑐_𝑠ℎ,𝑦_𝑛+𝑧_𝑠, (13) so that (7) becomes𝑦_𝑛+1=𝑦_𝑛+^𝑠_𝑖=1𝑑_𝑖𝑧_𝑖, (14) where the coefficients 𝑑_𝑖 are defined as𝑑₁,…,𝑑_𝑠=𝑏₁,…,𝑏_𝑠𝐴^𝑇. (15) For 𝑠=3, since 𝑏_𝑖=𝑎_𝑠,𝑖, the vector 𝑑 is (0,0,1).

The linear system with matrix𝐼−ℎ𝑎₁₁𝜕𝑓𝑡𝜕𝑦^𝑛+𝑐₁ℎ,𝑦_𝑛+𝑧₁⋯−ℎ𝑎_1𝑠𝜕𝑓𝑡𝜕𝑦^𝑛+𝑐_𝑠ℎ,𝑦_𝑛+𝑧_𝑠⋮⋱⋮−ℎ𝑎_𝑠1𝜕𝑓𝑡𝜕𝑦^𝑛+𝑐₁ℎ,𝑦_𝑛+𝑧₁⋯𝐼−ℎ𝑎_𝑠𝑠𝜕𝑓𝑡𝜕𝑦^𝑛+𝑐_𝑠ℎ,𝑦_𝑛+𝑧_𝑠 (16)must be solved to get the solution at the end of time step. The Jacobians in (12) are approximated as follows:𝐽=𝜕𝑓𝑡𝜕𝑦^𝑛+𝑐_𝑖ℎ,𝑦_𝑛+𝑧_𝑖≈𝜕𝑓𝑡𝜕𝑦^𝑛,𝑦_𝑛. (17) It has to be noted that for the case of depletion equations the calculation of the Jacobians is straightforward since 𝜕𝑓/𝜕𝑦_𝑖=𝜆_𝑖𝑖. The simplified Newton iterations become(𝐼−ℎ𝐴⊗𝐽)Δ𝑍^𝑘=−𝑍^𝑘+ℎ(𝐴⊗𝐼)𝐹𝑍^𝑘𝑍^𝑘+1=𝑍^𝑘+Δ𝑍^𝑘. (18) Here the following notations are used:(i)𝑘th approximation to the solution 𝑍^𝑘=(𝑧^𝑘₁,…,𝑧^𝑘_𝑠)^𝑇;(ii)increments Δ𝑍^𝑘=(Δ𝑧^𝑘₁,…,Δ𝑧^𝑘_𝑠)^𝑇;(iii)function F(𝑍^𝑘)=(𝑓(𝑡_𝑛+𝑐₁ℎ,𝑦_𝑛+𝑧^𝑘₁),…, 𝑓(𝑡_𝑛+𝑐_𝑠ℎ,𝑦_𝑛+𝑧^𝑘_𝑠))^𝑇.

Thus each iteration requires 𝑠=3 evaluations of 𝑓 and the solution of a linear system with the dimension 3𝑁. Each linear system of dimension 𝑠𝑁 can be split into one real and one complex linear system of dimension 𝑁 through a dimensionality reduction procedure [36]. As a result, the computing time of solving each related linear system is sharply reduced. This dimensionality reduction makes the simplified Newton method efficient in the implementation of IRK methods. However, this method has poor convergence properties since it is only an approximation of the classic Newton-chord method [37]. Reducing the step size can accelerate the convergence rate of the simplified Newton method. RADAU5 automatically selects a sufficiently small step size of IRK to yield a fast rate of convergence.

The matrix 𝐼−ℎ𝐴⊗𝐽 is computed only once for all iterations for linear problems with constant coefficients (depletion problems belong to this category).

The step size is not a constant and changes according to the following relation [32]:ℎ_new=ℎ_old⋅𝜃Tol‖Err‖^1/(𝑠+1), (19) where ‖Err‖ is error estimate, Tol is prescribed tolerance, and 𝜃 is model parameter ranging from 0.6 to 0.9 and usually being set to 0.9 if not too many step rejections occur.

2.2.2. Tests

(1) Accuracy
To test the accuracy of the RADAU5 solver applied to depletion calculations, a simple problem of calculation the ²¹⁰Po inventory accumulated in lead-bismuth eutectic coolant of MYRRHA after 90 days of irradiation has been investigated. Neglecting ²¹⁰Po production chains from lead isotopes (whose contribution into total ²¹⁰Po inventory is rather small), we focus on ²¹⁰Po production from neutron capture on ²⁰⁹Bi through the chain: ²⁰⁹Bi−−→^{(𝑛,𝛾)210}Bi−−−→^5.013𝑑210Po−−−−−→^{138.376𝑑206}Pb. (20) The system (1) applied to this chain can be solved analytically for the given time step and initial concentration of  ²⁰⁹Bi using the classic Bateman solutions [38] 𝑦_𝑛=𝑦₁(0)𝜆_𝑛^𝑛𝑡𝑟_𝑖=1𝜆_𝑖^𝑡𝑟𝛼_𝑖𝑒^{−𝜆𝑖𝑡𝑟𝑡}, (21) where 𝛼_𝑖=^𝑛_{𝑗=1𝑗≠𝑖}𝜆_𝑗^𝑡𝑟𝜆_𝑗^𝑡𝑟−𝜆_𝑖^𝑡𝑟. (22) The indexing here is according to (20) so that 𝑖=1 is assigned to ²⁰⁹Bi.
The numerical data used to obtain concentrations are summarized in Table 1. The values of average capture cross section on ²⁰⁹Bi leading to the ground state of ²¹⁰Bi, neutron flux, and initial concentration of ²⁰⁹Bi are typical for MYRRHA core. It has to be noted that the energy dependence of the branching ratio of (𝑛,𝛾) reaction leading to ground and metastable states of ²¹⁰Bi has been taken into account, since only the ground state of ²¹⁰Bi decays to ²¹⁰Po.
For the calculations with RADAU5, 𝜆₂₁=−𝜆₂₁,𝜆₃₂=−𝜆₂₂,𝜆₄₃=−𝜆₃₃, all other 𝜆_𝑗𝑖=0.
For the calculations with ORIGEN-2.2, all the cross sections in ORIGEN-2.2 library were zeroed, except neutron capture cross section on ²⁰⁹Bi. The considered decay chain was also separated from other chains in the decay library by zeroing decay constants to ensure full consistency with analytical solution at the level of input data. The results of comparison are shown in Table 2.
The calculations with RADAU5 were performed with two different tolerance values, Tol=10⁻⁴ which was selected from numerous tests as ensuring the best quality at the reasonable speed when handling full matrix with 𝑁=3851, and, for the sake of comparison, Tol=10⁻¹⁰ which requires longer calculation time. As it can be seen from Table 2, ORIGEN-2.2 code cannot even provide accurate enough results for the small system of 4 equations. In contrast, RADAU5 has shown the outstanding quality. Thus, the obvious conclusion is that matrix exponential method of ORIGEN-2.2 cannot solve concentrations as accurately as RADAU5 does.

Table 1 

Input parameters for calculation of ²¹⁰Po production.

Table 2 

Calculated nuclide inventories for 𝑡=90 days and relative differences ((𝑦^calc_𝑖−𝑦^anal_𝑖)/𝑦^anal_𝑖) from analytical solution.

(2) Speed
RADAU5 has shown an excellent accuracy compared to matrix exponential method of ORIGEN-2.2, but at the cost of computational time which strongly depends on the size of matrix Λ with effective reaction rates. The comparison of RADAU5 (IRK method) and DLSODA (BDF method) is shown in Figure 1 as a function of matrix dimension 𝑁. The results were normalized to the CPU time required for ORIGEN-2.2 to solve the same number of equations. It has to be noted that ORIGEN is capable to treat only ~1700 nuclides.
The matrix exponential method of ORIGEN-2.2 has shown weak dependence on the size of the ODE system. Two precise numerical methods are competitive with ORIGEN only up to 𝑁∼500. In terms of nuclide inventory calculations, the calculations of structural materials activation can be done even faster than by ORIGEN. RADAU5 runs almost 2 times faster than BDF. This is an additional justification of our choice of proper method to solve system of ODE.
As it is mentioned in [35], poor convergence properties of the simplified Newton method may force the code to generate too small step sizes, which reduces the performance. Use of larger step sizes to solve an ODE system is one of main advantages of IRK. To retain this advantage, more effective iterative algorithms were proposed in [35] with better convergence properties than the simplified Newton method. The gain, as reported in [35], could be a reduction of the CPU time needed by a factor of two.
But this improvement, once implemented, will not significantly influence the total calculation time for typical depletion problems. Usually, Monte Carlo steady-state particle transport calculation at the beginning of time step requires much more CPU performance than subsequent generation of matrix with effective reaction rates and solving the ODEs. The typical calculation time distribution among ALEPH modules required to complete one burn-up step of MYRRHA sub-critical core is shown in Table 3.
To reduce the impact of the RADAU5 solver on the total calculation time, especially for the problems with a huge number of irradiated materials, a simple parallelization was recently implemented when each material is treated separately in its own parallel task. The CPU wall time share shown in Table 3 is reduced to ~30 sec.
To summarize, the speed-up of RADAU5 is subject for future work, in view of implementation of nuclear data uncertainty propagation algorithms which are also in our high-priority to-do list for future ALEPH2 development. However, even now for typical ADS problems involving numerous materials (fuel, clad, structural materials, LBE spallation target, and coolant, etc.) ALEPH2 thanks to the RADAU5 solver provides very accurate nuclide inventories at an almost negligible CPU time compared to the total CPU time needed for the Monte Carlo simulation. Moreover, the uncertainties on the nuclide concentrations at the end of time step are associated with nuclear data uncertainties and Monte Carlo statistical uncertainties only and not with the depletion calculations.

Table 3 

CPU time share between ALEPH modules to calculate one irradiation step (90 days) of MYRRHA subcritical core. The MCNP(X) calculations have been performed on 64 CPU cluster using Open MPI. The number of irradiated materials is 11 (so that the reaction rates generator and the RADAU5 solver were called 11 times). Full matrix for 𝑁=3851 nuclides is used.

Figure 1 

Elapsed CPU time as a function of the number of ODE equations.

2.3. Predictor-Corrector Algorithm

The previous version of ALEPH did not use any predictor-corrector mechanism. A predictor-corrector mechanism aims to increase the macro time steps by using intermediate information. It should be noted that for the same macro time steps taken, the predictor-corrector method will obviously slow down the code since additional steady-state Monte Carlo calculations of fluxes and spectra are required. So if the macro time steps are kept the same, the option with predictor/corrector will increase the accuracy of the calculation at the cost of more CPU time. In case the required accuracy is fixed, the predictor/corrector will allow increasing the macro time steps, reducing the total number of macro time steps for the full irradiation history simulation. Of course, the best option from the quality/performance point of view remains the fine binning of irradiation history; however it is not always possible to properly define it. This situation is somehow similar to the variance reduction technique inside a Monte Carlo code itself: the best variance reduction is simply higher statistics, but for some problems such increase is not feasible.

In order not to severely overburden memory consumption, the number of Monte Carlo calculations per step using predictor-corrector is limited to 2. This means that the fluxes and spectra are calculated at the beginning of time step, and then the concentrations are obtained at the half of the time step. Using these concentrations, the fluxes and spectra are recalculated and finally the concentrations at the end of time step are computed on the basis of new fluxes and spectra. The ALEPH2 code decides whether to apply predictor-corrector or not by comparing the nuclide concentrations and reaction rates at the beginning and end of time step.

A new and unique feature of ALEPH2 is the possibility of using the time-dependent matrix coefficients when solving the system of ODEs by the RADAU5 method. Since the function 𝐟(𝑡,𝐲(𝑡)) in the right-hand side of (5) can be of arbitrary form, the coefficients 𝜆_𝑖𝑗 can in general be time-dependent. This allows to better reflect the realistic irradiation conditions. Starting from the second macro step, when the matrix coefficients are computed for 3 time points (𝑡=0,𝑡=1/2𝜏 and 𝑡=𝜏, here 𝜏 denotes time step length), ALEPH2 uses the linear extrapolation of the matrix coefficients and RADAU5 solver computes the nuclide concentrations using these linearly time-dependent reaction rates. The implementation of this feature is currently in the test phase.

2.4. Calculation of Photon Heating

The calculation of heating caused by photons is divided in two steps. When ALEPH2 calculates regular tallies of particle flux and spectra, it appends the photons to the list of transported particles and computes the flux and spectrum of prompt gammas. By folding the heating numbers (kerma factors) from nuclear data files with photon spectra, it computes the local heating caused by prompt photons. After that, ALEPH2 creates an additional MCNP(X) photon-only problem with the source of photons built from the information on nuclide concentrations at the beginning of time step and associated decay gammas. Nuclides are assumed to be uniformly distributed inside the material of interest. If the user has provided the cylindrical mesh covering the problem geometry with irradiated materials, the photon source is distributed similarly to the MCNP(X) power law “−21”. The additional calculation of photon-only problem usually does not increase significantly the total calculation time required for the given irradiation step. Using again the photon heating numbers, ALEPH2 computes the delayed photon heating and then sums up the two contributions—prompt and delayed—to get total photon heating.

2.5. (𝛼,𝑛) Neutron Source

Another major improvement of the code concerns the possibility to calculate the neutron source due to (𝛼,𝑛) reaction. It has to be noted, however, that this option can work with the MCNPX code only since it requires the simulation of alpha particle transport which is not compatible with MCNP code family.

The code generates the dependent source distribution of alpha particles at the beginning of each decay step (if requested by the user) in a similar way as it does for delayed photons described above. The energy distribution and intensity of this source depends on the material composition in the cells containing materials to be followed. Then the calculation of alpha particle-only problem is launched at the beginning of each decay step. The alpha particle flux and spectrum is tallied. It has to be noted again that ALEPH2 optimizes the source and physics parameters in such way that the MCNPX calculation completes quickly so that the total ALEPH2 calculation time remains almost unaffected. After the completion of spectrum calculation, the neutron source is computed using (𝛼,𝑛) reaction cross section library generated from TENDL-2010 library data [30].

This feature allows calculating the decay neutron source in a much more precise way than it is realized in ORIGEN-2.2 and ORIGEN-S.

2.6. Termination by Burn-Up Tracers Concentration

In the majority of problems the irradiation history cannot be explicitly known. Experimentally, one uses so-called burn-up tracers like Neodymium to estimate the burn-up in energy produced per heavy metal mass. For researchers working with this experimental data, it is important to be able to end the irradiation history, or better simulation history, when one or more of these tracers have reached a certain concentration corresponding to the experimentally measured values. This option has been added to ALEPH2. The user has either the ability to define a weighted sum of selected isotopes which has to converge with a certain precision (minimization of the deviation of the weighted sum) or to create a “virtual” tracer consisting of a sum of different isotopes.

2.7. Other Improvements

Besides major improvements some minor corrections were done, such as the possibility to simulate fuel management by allowing material, temperature, density replacement, and fuel assembly permutation at each desired time step.

The code possess the very useful capability to perform criticality calculations as complementary to the main fixed source calculation in the problems such as core burn-up of an ADS. This allows the user to better optimize the core shuffling, for instance. However, for these applications the calculation time approximately doubles since an additional Monte Carlo calculation is required.

Next, the improved version of the code contains significantly lower number of commands which have to be added to the standard MCNP(X) input deck than its predecessor. This reduces the complexity of running ALEPH2 and the probability for users to make errors.

3. Myrrha Neutronics Calculations

ALEPH2 is currently one of the main tools to perform neutronics design of the MYRRHA fast spectrum research facility. As it was mentioned throughout the paper, the sub-critical (and critical) core burn-up, test materials irradiation performance, activation of structural materials, spallation products accumulation, radiation source terms—all these quantities—can be calculated (and are being calculated) with ALEPH2. In this section, the examples of sub-critical core burn-up calculations and spallation products analysis are shown.

3.1. Core Burn-Up with Mixed Shuffling

A possible irradiation scheme envisages batches of three 90-day irradiations with two 30-day beam-off periods for maintenance between irradiations and long 90-day maintenance period after three irradiations. Current facility design implies only a limited number of interventions to increase the proton beam current in order to compensate reactivity losses. In the ALEPH2 model, the beam current was increased in the middle of the cycle (45 days) to compensate the reactivity loss. However, the limitations on the beam current and on the maximum clad temperature do not allow to achieve high burn-up values without fuel reshuffling. That’s why several reshuffling strategies have been investigated with ALEPH2 code. Since ALEPH2 allows permuting or changing materials, the modeling of the reshuffling procedure becomes very simple in ALEPH2. The beginning-of-life (BoL) core consists of 58 fresh fuel assemblies (FAs) and is characterized by effective neutron multiplication factor 𝑘_eff∼0.965. The fuel assemblies are grouped in six groups in order to minimize the flux variations between assemblies in the group. During first five maintenance periods, 2 to 4 assemblies with fresh fuel are added so that the core shown in Figure 2 is built with 72 FA (6 groups of 12 FA each). The central assembly hosts the spallation target structures. Besides that, there are shown 6 in-pile sections (IPSs) for test materials irradiations and one assembly hosting the irradiation device for production of radioisotopes. Control rods are in upper positions and are not assumed to be inserted in the core during normal operation.

Figure 2 

MYRRHA sub-critical core with 72 FA collected in 6 groups 12 FA each according to average neutron flux values.

The reshuffling strategy aimed at minimizing the reactivity swing during the cycle and reducing the radial power factor because the sub-critical cores with the source in the center are characterized by a more steep flux gradiant in the radial direction. A sensitivity analysis allowed to select the reshuffling scheme shown in Figure 3.

Figure 3 

Mixed reshuffling scheme.

In this scheme, fresh fuel at the beginning of each cycle (BoC) is always inserted in Zone 2. Then, at each reshuffling, it subsequently goes to Zone 4, then Zone 3, Zone 1, Zone 6, and during the last irradiation cycle it resides in Zone 5. Results of main neutronics parameters calculations with ALEPH2 are shown in Figures 4, 5, and 6.

Figure 4 

Evolution of effective neutron multiplication factor 𝐾eff and external source neutron multiplication factor 𝐾src.

Figure 5 

Evolution of total thermal power release.

Figure 6 

Evolution of proton beam current.

The system reaches equilibrium in ~3 years. The reshuffling scheme chosen allows keeping the power almost the same at the beginning of each cycle. The reactivity swing during the equilibrium cycle, ~1100 pcm, appears to be much less than that for the core loaded with fresh fuel (~1700 pcm). This is clearly a result from the optimization in positioning the fuel assemblies with different burn-up level.

3.2. Spallation Products Accumulation—A Mercury Example

An important question to be answered is whether the confinement of the spallation source in a separate loop is worthwhile from the viewpoint of containing spallation and high-energy fission products. In order to analyze this, the specific case of mercury (with the very long-lived isotope ¹⁹⁴Hg) was studied.

The presence of Hg (all isotopes) has been evaluated for both the spallation volume and the volume of LBE circulating in the core. In the LBE composed of lead and bismuth, an impurity level of 1.734 mg/kg of mercury (natural occurring isotopes) has been assumed.

The different mercury yields, that is, atoms of Hg produced per proton coming from the accelerator, have been calculated. A distinction has been made in yields coming from proton-induced reactions and neutron induced reactions. The numerical values are given in Table 4.

Table 4 

Hg yield in atoms per source proton.

Since ALEPH2 calculates neutron- and proton-induced production cross sections by averaging over the corresponding particle spectra, the yields shown in Table 4 were derived from ALEPH2 output using the relation 𝑌^𝑚=𝑉⋅𝛾_LBE𝜎^𝑚𝜑^𝑚_MCNPX, (23) where 𝑌^𝑚 is the given nuclide yield due to reactions caused by particle type 𝑚, expressed in atoms per source proton, 𝑉 is the scoring volume, 𝛾_LBE represents the atomic density of LBE in atoms/(b⋅cm), 𝜎^𝑚 stands for the energy-averaged production cross section (b) of the given nuclide by particle type 𝑚, and 𝜑^𝑚_MCNPX(𝐸) is the MCNPX flux tally of particles 𝑚, expressed in particles per cm² per source proton. As it can be seen from this table, the relative contribution of proton-induced reactions is only about 25% both in the spallation target volume and in the core region. The rest is produced by neutrons (which in turn are produced in spallation target).

Table 5 summarizes the different origins of mercury present in the spallation target and core volume after 90 days of irradiation at 2.5 mA, which corresponds to 70 MW of thermal power release at the equilibrium cycle. It is clear that the majority of the mercury is produced outside the spallation target volume.

Table 5 

Origin of mercury in the LBE at the end of 90-day irradiation cycle.

The evolution of activity with distance to the core center is shown in Figure 7 for several radioactive mercury isotopes, ranging from short-lived to long-lived ¹⁹⁴Hg whose activity after 90 days of irradiation is rather low but will remain the main contributor to the activity of the irradiated LBE for thousands of years.

Figure 7 

The radial dependence of Hg isotopes activity. Rings 1 to 6 denote the fuel assembly groups corresponding to Figure 2.

Figure 7 clearly shows that the presence of high-energy tail in neutron spectrum will produce mercury isotopes even in the periphery of the core. The thresholds of most of the neutron- and proton-induced reactions leading to production of mercury isotopes from lead and bismuth isotopes are typically 60–80 MeV. Neutron- and proton-spectra in the spallation target and core as well as ¹⁹⁴Hg production cross sections from neutron and proton-induced reactions on LBE isotopes are plotted in Figure 8 for illustration purposes. The cross sections are those used by ALEPH2; they have been generated from TENDL-2010 and HEAD-2009 libraries. Figure 8 confirms that spallation reactions will occur elsewhere (although with less intensity than in spallation target) where particles with energies above production threshold are present.

(a)

(b)

(c)

Figure 8 

Neutron (a) and proton (b) energy spectra averaged over spallation target and core volumes. Neutron- and proton-induced 194Hg production cross section from lead and bismuth isotopes (c).

4. Conclusions

Two principal features of ALEPH2 Monte Carlo depletion code, an extensive nuclear data library and a state-of-the-art numerical solver for the Bateman equations, make this tool flexible and powerful. The nuclear data library covers neutron- and proton-induced reactions, neutron and proton fission product yields, spontaneous fission product yields, radioactive decay data and, total recoverable energies per fission. The numerical solver applied is an implicit Runge-Kutta method of the RADAU IIA family. The versatility of the code allows using it for time behavior simulation of various systems ranging from single pin models to full-scale reactor models, including specific facilities as accelerator-driven systems. The core burn-up, activation of the structural materials, irradiation of samples, and, in addition, accumulation of spallation products in accelerator-driven systems can be calculated in a single ALEPH2 run. The associated source term values, including decay heat release caused by particular radiation type, (𝛼,𝑛) and spontaneous fission neutron sources, absorbed doses, are calculated in much more precise way compared to other codes. All these features make ALEPH2 superior over other Monte Carlo depletion codes. The code is extensively used for the neutronics design of the MYRRHA research facility which will operate in both critical and sub-critical modes.