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

Authors: A. Stankovskiy; G. Van den Eynde
Journal: Science and Technology of Nuclear Installations
 (2012)

---

## 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.


---

## Body

## 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,(n,2n),(n,3n),(n,p), and (n,α)). 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 nuclidesyi at the end of irradiation step can be obtained by solving the system of first-order linear differential equations with constant coefficients:(1)dyi(t)dt=∑jλjitryj(t)-λitryi(t),
or, in compact matrix form,(2)dy⃗dt=Λy⃗,
with the matrix coefficients being λij=λjitr-λitrδij. Here λjitr is the transmutation (production) rate of nuclide i from nuclide j, and λitr is the disappearance rate of nuclide i. These rates are defined as follows:(3)λjitr=∑mσjimφm+λjid,(4)λitr=∑mσR,imφm+λid.
Here σjim is the energy-averaged production cross section of nuclide i from nuclide j in the reaction with incident particle m, φm is the average flux of type m particles, λjid stands for the decay rate of nuclide j into the nuclide i, σR,im represents the energy-averaged reaction cross section of nuclide i transforming it into another nuclide, and λid is the decay constant of nuclide i. In most cases m=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)dy(t)dt=f(t,y(t));y(t0)=y0.y is unknown N-size vector of values (concentrations in case of depletion equations), and f(t,y(t)) is the N-size vector of functions of known form. In case of depletion equations, it is simply f(t,y(t))=Λy(t), and Λis the effective reaction rate matrix.If the time step is subdivided intom substeps with a fixed step size(6)h=t1-t0m
with the mesh points tn,n=0,1,…,m,tm=t1, the implicit Runge-Kutta method is defined by(7)yn+1=yn+h∑i=1sbif(tn+cih,gi).
Here(8)gi=yn+h∑j=1saijf(tn+cjh,gj),i=1,…,s.
The stage number is denoted by s and aij,bi, and ci are parameters satisfying the conditions:(9)aij≠0,∑j=1saij=cifori=1,2,…,s,∑i=1sbi=1.
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 ci for i=1,2,…,s are selected as the zeros of Radau polynomial (ds-1/dts-1)(ts-1(t-1)s). Parameters aij are obtained by solving the linear system(10)∑j=1saijcjq-1=ciqq,i=1,…,s;q=1,…,s
and bi=as,i for i=1,2,…,s. Stage number s can be set to 1,3,5, or 7. The RADAU5 solver uses s=3 thus its accuracy order is 2s-1=5 and it is considered as absolutely stable.At each sub-step, the system (3) of rank sN has to be solved iteratively for the stage values {gi}i=1s. To reduce round-off errors, it is often reformulated as(11)zi=gi-yn,
where(12)zi=h∑j=1saijf(t0+cjh,yn+zj),i=1,…,s.
Since the matrix of coefficients A=(aij) is nonsingular (as it can be seen from (10)), (12) can be rewritten in vector form(13)(z1...zs)=A(hf(tn+c1h,yn+z1)...hf(tn+csh,yn+zs)),
so that (7) becomes(14)yn+1=yn+∑i=1sdizi,
where the coefficients di are defined as(15)(d1,…,ds)=(b1,…,bs)AT.
For s=3, since bi=as,i, the vector d is (0,0,1).The linear system with matrix(16)(I-ha11∂f∂y(tn+c1h,yn+z1)⋯-ha1s∂f∂y(tn+csh,yn+zs)⋮⋱⋮-has1∂f∂y(tn+c1h,yn+z1)⋯I-hass∂f∂y(tn+csh,yn+zs))must be solved to get the solution at the end of time step. The Jacobians in (12) are approximated as follows:(17)J=∂f∂y(tn+cih,yn+zi)≈∂f∂y(tn,yn).
It has to be noted that for the case of depletion equations the calculation of the Jacobians is straightforward since ∂f/∂yi=λii. The simplified Newton iterations become(18)(I-hA⊗J)ΔZk=-Zk+h(A⊗I)F(Zk)Zk+1=Zk+ΔZk.
Here the following notations are used:(i)
kth approximation to the solution Zk=(z1k,…,zsk)T;(ii)
incrementsΔZk=(Δz1k,…,Δzsk)T;(iii)
functionF(Zk)=(f(tn+c1h,yn+z1k),…, f(tn+csh,yn+zsk))T.Thus each iteration requiress=3 evaluations of f and the solution of a linear system with the dimension 3N. Each linear system of dimension sN can be split into one real and one complex linear system of dimension N 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 matrixI-hA⊗J 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]:(19)hnew=hold⋅θ(Tol‖Err‖)1/(s+1),
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 the210Po inventory accumulated in lead-bismuth eutectic coolant of MYRRHA after 90 days of irradiation has been investigated. Neglecting 210Po production chains from lead isotopes (whose contribution into total 210Po inventory is rather small), we focus on 210Po production from neutron capture on 209Bi through the chain:
(20)209Bi→(n,γ)Bi→5.013dPo→138.376dPb206210210.
The system (1) applied to this chain can be solved analytically for the given time step and initial concentration of  209Bi using the classic Bateman solutions [38]
(21)yn=y1(0)λntr∑i=1nλitrαie-λitrt,
where
(22)αi=∏j=1j≠inλjtr(λjtr-λitr).
The indexing here is according to (20) so that i=1 is assigned to 209Bi.
The numerical data used to obtain concentrations are summarized in Table1. The values of average capture cross section on 209Bi leading to the ground state of 210Bi, neutron flux, and initial concentration of 209Bi are typical for MYRRHA core. It has to be noted that the energy dependence of the branching ratio of (n,γ) reaction leading to ground and metastable states of 210Bi has been taken into account, since only the ground state of 210Bi decays to 210Po.
For the calculations with RADAU5,λ21=-λ21,λ32=-λ22,λ43=-λ33, all other λji=0.
For the calculations with ORIGEN-2.2, all the cross sections in ORIGEN-2.2 library were zeroed, except neutron capture cross section on209Bi. 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-4 which was selected from numerous tests as ensuring the best quality at the reasonable speed when handling full matrix with N=3851, and, for the sake of comparison, Tol=10-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 of210Po production.
iNuclidey(0), at/(b·cm)λii, s−11209Bi0.0165368-2.09·10-11*2210Bi0.0-1.60043·10-63210Po0.0-5.79554·10-84206Pb0.00.0*Obtained by multiplyingσ=0.0209 b with φ=1·1015 n/(cm2s).Table 2
Calculated nuclide inventories fort=90 days and relative differences ((yicalc-yianal)/yianal) from analytical solution.
NuclideAnalytical solutionORIGEN-2.2RADAU5,Tol=10-4RADAU5,Tol=10-10at/(b·cm)at/(b·cm)diff, %at/(b·cm)diff, %at/(b·cm)diff, %209Bi1.65341e–21.65341e–2-7.68e–51.65341e–20.01.65341e–20.0210Bi2.15920e–72.16954e–74.79e–12.15847e–7-3.38e–22.15920e–7-4.63e–11210Po2.02059e–61.99902e–61.072.02067e–63.75e–32.02059e–60.0206Pb4.50804e–74.74951e–75.364.50801e–7-5.75e–44.50804e–7-3.55e–10(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 N. 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 toN~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 Table3.
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 Table3 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 forN=3851 nuclides is used.
CPU wall time, sFraction, %MCNPX fluxes and spectra, rsd 0.1%3.58e+699.988ALEPH generation of average cross sections19.35.4e–4ALEPH solving ODE by RADAU5312.58.7e–3(or ALEPH solving ODE by ORIGEN-2.2)(14.7)(4.1e–4)Total100Figure 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 functionf(t,y(t)) in the right-hand side of (5) can be of arbitrary form, the coefficients λij 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 (t=0,t=1/2τ and t=τ, 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. (α,n) Neutron Source

Another major improvement of the code concerns the possibility to calculate the neutron source due to (α,n) 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 (α,n) 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.


## 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,(n,2n),(n,3n),(n,p), and (n,α)). 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 nuclidesyi at the end of irradiation step can be obtained by solving the system of first-order linear differential equations with constant coefficients:(1)dyi(t)dt=∑jλjitryj(t)-λitryi(t),
or, in compact matrix form,(2)dy⃗dt=Λy⃗,
with the matrix coefficients being λij=λjitr-λitrδij. Here λjitr is the transmutation (production) rate of nuclide i from nuclide j, and λitr is the disappearance rate of nuclide i. These rates are defined as follows:(3)λjitr=∑mσjimφm+λjid,(4)λitr=∑mσR,imφm+λid.
Here σjim is the energy-averaged production cross section of nuclide i from nuclide j in the reaction with incident particle m, φm is the average flux of type m particles, λjid stands for the decay rate of nuclide j into the nuclide i, σR,im represents the energy-averaged reaction cross section of nuclide i transforming it into another nuclide, and λid is the decay constant of nuclide i. In most cases m=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)dy(t)dt=f(t,y(t));y(t0)=y0.y is unknown N-size vector of values (concentrations in case of depletion equations), and f(t,y(t)) is the N-size vector of functions of known form. In case of depletion equations, it is simply f(t,y(t))=Λy(t), and Λis the effective reaction rate matrix.If the time step is subdivided intom substeps with a fixed step size(6)h=t1-t0m
with the mesh points tn,n=0,1,…,m,tm=t1, the implicit Runge-Kutta method is defined by(7)yn+1=yn+h∑i=1sbif(tn+cih,gi).
Here(8)gi=yn+h∑j=1saijf(tn+cjh,gj),i=1,…,s.
The stage number is denoted by s and aij,bi, and ci are parameters satisfying the conditions:(9)aij≠0,∑j=1saij=cifori=1,2,…,s,∑i=1sbi=1.
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 ci for i=1,2,…,s are selected as the zeros of Radau polynomial (ds-1/dts-1)(ts-1(t-1)s). Parameters aij are obtained by solving the linear system(10)∑j=1saijcjq-1=ciqq,i=1,…,s;q=1,…,s
and bi=as,i for i=1,2,…,s. Stage number s can be set to 1,3,5, or 7. The RADAU5 solver uses s=3 thus its accuracy order is 2s-1=5 and it is considered as absolutely stable.At each sub-step, the system (3) of rank sN has to be solved iteratively for the stage values {gi}i=1s. To reduce round-off errors, it is often reformulated as(11)zi=gi-yn,
where(12)zi=h∑j=1saijf(t0+cjh,yn+zj),i=1,…,s.
Since the matrix of coefficients A=(aij) is nonsingular (as it can be seen from (10)), (12) can be rewritten in vector form(13)(z1...zs)=A(hf(tn+c1h,yn+z1)...hf(tn+csh,yn+zs)),
so that (7) becomes(14)yn+1=yn+∑i=1sdizi,
where the coefficients di are defined as(15)(d1,…,ds)=(b1,…,bs)AT.
For s=3, since bi=as,i, the vector d is (0,0,1).The linear system with matrix(16)(I-ha11∂f∂y(tn+c1h,yn+z1)⋯-ha1s∂f∂y(tn+csh,yn+zs)⋮⋱⋮-has1∂f∂y(tn+c1h,yn+z1)⋯I-hass∂f∂y(tn+csh,yn+zs))must be solved to get the solution at the end of time step. The Jacobians in (12) are approximated as follows:(17)J=∂f∂y(tn+cih,yn+zi)≈∂f∂y(tn,yn).
It has to be noted that for the case of depletion equations the calculation of the Jacobians is straightforward since ∂f/∂yi=λii. The simplified Newton iterations become(18)(I-hA⊗J)ΔZk=-Zk+h(A⊗I)F(Zk)Zk+1=Zk+ΔZk.
Here the following notations are used:(i)
kth approximation to the solution Zk=(z1k,…,zsk)T;(ii)
incrementsΔZk=(Δz1k,…,Δzsk)T;(iii)
functionF(Zk)=(f(tn+c1h,yn+z1k),…, f(tn+csh,yn+zsk))T.Thus each iteration requiress=3 evaluations of f and the solution of a linear system with the dimension 3N. Each linear system of dimension sN can be split into one real and one complex linear system of dimension N 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 matrixI-hA⊗J 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]:(19)hnew=hold⋅θ(Tol‖Err‖)1/(s+1),
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 the210Po inventory accumulated in lead-bismuth eutectic coolant of MYRRHA after 90 days of irradiation has been investigated. Neglecting 210Po production chains from lead isotopes (whose contribution into total 210Po inventory is rather small), we focus on 210Po production from neutron capture on 209Bi through the chain:
(20)209Bi→(n,γ)Bi→5.013dPo→138.376dPb206210210.
The system (1) applied to this chain can be solved analytically for the given time step and initial concentration of  209Bi using the classic Bateman solutions [38]
(21)yn=y1(0)λntr∑i=1nλitrαie-λitrt,
where
(22)αi=∏j=1j≠inλjtr(λjtr-λitr).
The indexing here is according to (20) so that i=1 is assigned to 209Bi.
The numerical data used to obtain concentrations are summarized in Table1. The values of average capture cross section on 209Bi leading to the ground state of 210Bi, neutron flux, and initial concentration of 209Bi are typical for MYRRHA core. It has to be noted that the energy dependence of the branching ratio of (n,γ) reaction leading to ground and metastable states of 210Bi has been taken into account, since only the ground state of 210Bi decays to 210Po.
For the calculations with RADAU5,λ21=-λ21,λ32=-λ22,λ43=-λ33, all other λji=0.
For the calculations with ORIGEN-2.2, all the cross sections in ORIGEN-2.2 library were zeroed, except neutron capture cross section on209Bi. 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-4 which was selected from numerous tests as ensuring the best quality at the reasonable speed when handling full matrix with N=3851, and, for the sake of comparison, Tol=10-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 of210Po production.
iNuclidey(0), at/(b·cm)λii, s−11209Bi0.0165368-2.09·10-11*2210Bi0.0-1.60043·10-63210Po0.0-5.79554·10-84206Pb0.00.0*Obtained by multiplyingσ=0.0209 b with φ=1·1015 n/(cm2s).Table 2
Calculated nuclide inventories fort=90 days and relative differences ((yicalc-yianal)/yianal) from analytical solution.
NuclideAnalytical solutionORIGEN-2.2RADAU5,Tol=10-4RADAU5,Tol=10-10at/(b·cm)at/(b·cm)diff, %at/(b·cm)diff, %at/(b·cm)diff, %209Bi1.65341e–21.65341e–2-7.68e–51.65341e–20.01.65341e–20.0210Bi2.15920e–72.16954e–74.79e–12.15847e–7-3.38e–22.15920e–7-4.63e–11210Po2.02059e–61.99902e–61.072.02067e–63.75e–32.02059e–60.0206Pb4.50804e–74.74951e–75.364.50801e–7-5.75e–44.50804e–7-3.55e–10(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 N. 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 toN~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 Table3.
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 Table3 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 forN=3851 nuclides is used.
CPU wall time, sFraction, %MCNPX fluxes and spectra, rsd 0.1%3.58e+699.988ALEPH generation of average cross sections19.35.4e–4ALEPH solving ODE by RADAU5312.58.7e–3(or ALEPH solving ODE by ORIGEN-2.2)(14.7)(4.1e–4)Total100Figure 1
Elapsed CPU time as a function of the number of ODE equations.


## 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)dy(t)dt=f(t,y(t));y(t0)=y0.y is unknown N-size vector of values (concentrations in case of depletion equations), and f(t,y(t)) is the N-size vector of functions of known form. In case of depletion equations, it is simply f(t,y(t))=Λy(t), and Λis the effective reaction rate matrix.If the time step is subdivided intom substeps with a fixed step size(6)h=t1-t0m
with the mesh points tn,n=0,1,…,m,tm=t1, the implicit Runge-Kutta method is defined by(7)yn+1=yn+h∑i=1sbif(tn+cih,gi).
Here(8)gi=yn+h∑j=1saijf(tn+cjh,gj),i=1,…,s.
The stage number is denoted by s and aij,bi, and ci are parameters satisfying the conditions:(9)aij≠0,∑j=1saij=cifori=1,2,…,s,∑i=1sbi=1.
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 ci for i=1,2,…,s are selected as the zeros of Radau polynomial (ds-1/dts-1)(ts-1(t-1)s). Parameters aij are obtained by solving the linear system(10)∑j=1saijcjq-1=ciqq,i=1,…,s;q=1,…,s
and bi=as,i for i=1,2,…,s. Stage number s can be set to 1,3,5, or 7. The RADAU5 solver uses s=3 thus its accuracy order is 2s-1=5 and it is considered as absolutely stable.At each sub-step, the system (3) of rank sN has to be solved iteratively for the stage values {gi}i=1s. To reduce round-off errors, it is often reformulated as(11)zi=gi-yn,
where(12)zi=h∑j=1saijf(t0+cjh,yn+zj),i=1,…,s.
Since the matrix of coefficients A=(aij) is nonsingular (as it can be seen from (10)), (12) can be rewritten in vector form(13)(z1...zs)=A(hf(tn+c1h,yn+z1)...hf(tn+csh,yn+zs)),
so that (7) becomes(14)yn+1=yn+∑i=1sdizi,
where the coefficients di are defined as(15)(d1,…,ds)=(b1,…,bs)AT.
For s=3, since bi=as,i, the vector d is (0,0,1).The linear system with matrix(16)(I-ha11∂f∂y(tn+c1h,yn+z1)⋯-ha1s∂f∂y(tn+csh,yn+zs)⋮⋱⋮-has1∂f∂y(tn+c1h,yn+z1)⋯I-hass∂f∂y(tn+csh,yn+zs))must be solved to get the solution at the end of time step. The Jacobians in (12) are approximated as follows:(17)J=∂f∂y(tn+cih,yn+zi)≈∂f∂y(tn,yn).
It has to be noted that for the case of depletion equations the calculation of the Jacobians is straightforward since ∂f/∂yi=λii. The simplified Newton iterations become(18)(I-hA⊗J)ΔZk=-Zk+h(A⊗I)F(Zk)Zk+1=Zk+ΔZk.
Here the following notations are used:(i)
kth approximation to the solution Zk=(z1k,…,zsk)T;(ii)
incrementsΔZk=(Δz1k,…,Δzsk)T;(iii)
functionF(Zk)=(f(tn+c1h,yn+z1k),…, f(tn+csh,yn+zsk))T.Thus each iteration requiress=3 evaluations of f and the solution of a linear system with the dimension 3N. Each linear system of dimension sN can be split into one real and one complex linear system of dimension N 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 matrixI-hA⊗J 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]:(19)hnew=hold⋅θ(Tol‖Err‖)1/(s+1),
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 the210Po inventory accumulated in lead-bismuth eutectic coolant of MYRRHA after 90 days of irradiation has been investigated. Neglecting 210Po production chains from lead isotopes (whose contribution into total 210Po inventory is rather small), we focus on 210Po production from neutron capture on 209Bi through the chain:
(20)209Bi→(n,γ)Bi→5.013dPo→138.376dPb206210210.
The system (1) applied to this chain can be solved analytically for the given time step and initial concentration of  209Bi using the classic Bateman solutions [38]
(21)yn=y1(0)λntr∑i=1nλitrαie-λitrt,
where
(22)αi=∏j=1j≠inλjtr(λjtr-λitr).
The indexing here is according to (20) so that i=1 is assigned to 209Bi.
The numerical data used to obtain concentrations are summarized in Table1. The values of average capture cross section on 209Bi leading to the ground state of 210Bi, neutron flux, and initial concentration of 209Bi are typical for MYRRHA core. It has to be noted that the energy dependence of the branching ratio of (n,γ) reaction leading to ground and metastable states of 210Bi has been taken into account, since only the ground state of 210Bi decays to 210Po.
For the calculations with RADAU5,λ21=-λ21,λ32=-λ22,λ43=-λ33, all other λji=0.
For the calculations with ORIGEN-2.2, all the cross sections in ORIGEN-2.2 library were zeroed, except neutron capture cross section on209Bi. 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-4 which was selected from numerous tests as ensuring the best quality at the reasonable speed when handling full matrix with N=3851, and, for the sake of comparison, Tol=10-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 of210Po production.
iNuclidey(0), at/(b·cm)λii, s−11209Bi0.0165368-2.09·10-11*2210Bi0.0-1.60043·10-63210Po0.0-5.79554·10-84206Pb0.00.0*Obtained by multiplyingσ=0.0209 b with φ=1·1015 n/(cm2s).Table 2
Calculated nuclide inventories fort=90 days and relative differences ((yicalc-yianal)/yianal) from analytical solution.
NuclideAnalytical solutionORIGEN-2.2RADAU5,Tol=10-4RADAU5,Tol=10-10at/(b·cm)at/(b·cm)diff, %at/(b·cm)diff, %at/(b·cm)diff, %209Bi1.65341e–21.65341e–2-7.68e–51.65341e–20.01.65341e–20.0210Bi2.15920e–72.16954e–74.79e–12.15847e–7-3.38e–22.15920e–7-4.63e–11210Po2.02059e–61.99902e–61.072.02067e–63.75e–32.02059e–60.0206Pb4.50804e–74.74951e–75.364.50801e–7-5.75e–44.50804e–7-3.55e–10(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 N. 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 toN~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 Table3.
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 Table3 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 forN=3851 nuclides is used.
CPU wall time, sFraction, %MCNPX fluxes and spectra, rsd 0.1%3.58e+699.988ALEPH generation of average cross sections19.35.4e–4ALEPH solving ODE by RADAU5312.58.7e–3(or ALEPH solving ODE by ORIGEN-2.2)(14.7)(4.1e–4)Total100Figure 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 functionf(t,y(t)) in the right-hand side of (5) can be of arbitrary form, the coefficients λij 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 (t=0,t=1/2τ and t=τ, 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. (α,n) Neutron Source

Another major improvement of the code concerns the possibility to calculate the neutron source due to (α,n) 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 (α,n) 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 factorkeff~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 Figure3.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 Figures4, 5, and 6.Figure 4
Evolution of effective neutron multiplication factorKeff and external source neutron multiplication factor Ksrc.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 isotope194Hg) 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 Table4.Table 4
Hg yield in atoms per source proton.
Spallation target volumeCore volumeHg yield due to protons (at/p)0.0320.036Hg yield due to neutrons (at/p)0.1020.110Total Hg yield (at/p)0.1340.146Since ALEPH2 calculates neutron- and proton-induced production cross sections by averaging over the corresponding particle spectra, the yields shown in Table4 were derived from ALEPH2 output using the relation (23)Ym=V⋅γLBEσmφMCNPXm,
where Ym is the given nuclide yield due to reactions caused by particle type m, expressed in atoms per source proton, V is the scoring volume, γLBE represents the atomic density of LBE in atoms/(b·cm), σm stands for the energy-averaged production cross section (b) of the given nuclide by particle type m, and φMCNPXm(E) is the MCNPX flux tally of particles m, expressed in particles per cm2 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).Table5 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.
Spallation target volumeCore volumeTotalHg as impurity (g)0.5625.6226.18Mass of Hg due to spallation (g)5.285.7411.02Combined mass initial and spallation (g)5.8431.3637.20The evolution of activity with distance to the core center is shown in Figure7 for several radioactive mercury isotopes, ranging from short-lived to long-lived 194Hg 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 Figure2.Figure7 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 194Hg 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.Neutron (a) and proton (b) energy spectra averaged over spallation target and core volumes. Neutron- and proton-induced194Hg production cross section from lead and bismuth isotopes (c).
(a)(b)(c)


## 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 factorkeff~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 Figure3.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 Figures4, 5, and 6.Figure 4
Evolution of effective neutron multiplication factorKeff and external source neutron multiplication factor Ksrc.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 isotope194Hg) 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 Table4.Table 4
Hg yield in atoms per source proton.
Spallation target volumeCore volumeHg yield due to protons (at/p)0.0320.036Hg yield due to neutrons (at/p)0.1020.110Total Hg yield (at/p)0.1340.146Since ALEPH2 calculates neutron- and proton-induced production cross sections by averaging over the corresponding particle spectra, the yields shown in Table4 were derived from ALEPH2 output using the relation (23)Ym=V⋅γLBEσmφMCNPXm,
where Ym is the given nuclide yield due to reactions caused by particle type m, expressed in atoms per source proton, V is the scoring volume, γLBE represents the atomic density of LBE in atoms/(b·cm), σm stands for the energy-averaged production cross section (b) of the given nuclide by particle type m, and φMCNPXm(E) is the MCNPX flux tally of particles m, expressed in particles per cm2 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).Table5 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.
Spallation target volumeCore volumeTotalHg as impurity (g)0.5625.6226.18Mass of Hg due to spallation (g)5.285.7411.02Combined mass initial and spallation (g)5.8431.3637.20The evolution of activity with distance to the core center is shown in Figure7 for several radioactive mercury isotopes, ranging from short-lived to long-lived 194Hg 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 Figure2.Figure7 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 194Hg 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.Neutron (a) and proton (b) energy spectra averaged over spallation target and core volumes. Neutron- and proton-induced194Hg production cross section from lead and bismuth isotopes (c).
(a)(b)(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, (α,n) 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.



---

*Source: 545103-2012-04-09.xml*