Table of Contents Author Guidelines Submit a Manuscript
Science and Technology of Nuclear Installations
Volume 2012, Article ID 545103, 12 pages
Research Article

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

Institute of Advanced Nuclear Systems, SCK·CEN, Boeretang 200, 2400 Mol, Belgium

Received 14 November 2011; Accepted 26 January 2012

Academic Editor: Alberto Talamo

Copyright © 2012 A. Stankovskiy and G. Van den Eynde. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 [2123].

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, , 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: or, in compact matrix form, 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: 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 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]: 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 with the mesh points , the implicit Runge-Kutta method is defined by Here The stage number is denoted by and , and are parameters satisfying the conditions: 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 are selected as the zeros of Radau polynomial . Parameters are obtained by solving the linear system and for . Stage number can be set to , or 7. The RADAU5 solver uses thus its accuracy order is 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 . To reduce round-off errors, it is often reformulated as where Since the matrix of coefficients is nonsingular (as it can be seen from (10)), (12) can be rewritten in vector form so that (7) becomes where the coefficients are defined as For , since , the vector is ().

The linear system with matrixmust be solved to get the solution at the end of time step. The Jacobians in (12) are approximated as follows: 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 Here the following notations are used:(i)th approximation to the solution ;(ii)increments ;(iii)function ,, .

Thus each iteration requires evaluations of and the solution of a linear system with the dimension . 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]: where 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 210Po 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: 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] where The indexing here is according to (20) so that is assigned to 209Bi.
The numerical data used to obtain concentrations are summarized in Table 1. 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 () 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, , all other .
For the calculations with ORIGEN-2.2, all the cross sections in ORIGEN-2.2 library were zeroed, except neutron capture cross section on 209Bi. 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, which was selected from numerous tests as ensuring the best quality at the reasonable speed when handling full matrix with , and, for the sake of comparison, 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 210Po production.
Table 2: Calculated nuclide inventories for days and relative differences 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 . 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 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 ( 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 . 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 and external source neutron multiplication factor .
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 194Hg) 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 where is the given nuclide yield due to reactions caused by particle type , expressed in atoms per source proton, is the scoring volume, represents the atomic density of LBE in atoms/(bcm), stands for the energy-averaged production cross section (b) of the given nuclide by particle type , and is the MCNPX flux tally of particles , 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).

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

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.


  1. H. Aït Abderrahim, “MYRRHA, a new future for nuclear research. DRAFT-2 pre-design file,” Report R 4234, SCK·CEN, 2006. View at Google Scholar
  2. “Independent evaluation of the MYRRHA project,” MIRT Report, OECD-NEA-6881, 2009.
  3. D. B. Pelowitz, Ed., “MCNPX user's manual, version 2.6.0,” Tech. Rep. LA-CP-07-1473, 2008. View at Google Scholar
  4. D. B. Pelowitz et al., “MCNPX 2.7.C extensions,” Tech. Rep. LA-UR-10-00481, Los Alamos National Laboratory, 2010. View at Google Scholar
  5. W. B. Wilson, T. R. England, and P. Möller, “A Manual for CINDER’90, Version 06.1 Codes and Data,” Advanced Fuel Cycle Program, 2006.
  6. Y. A. Korovin, A. Y. Konobeyev, P. E. Pereslavtsev et al., “The code for the calculation of nuclide composition and activity of irradiated materials,” Journal Yadernye Konstanty, vol. 3-4, p. 117, 1992. View at Google Scholar
  7. R. A. Forrest, “The European Activation System: EASY-2007 Overview,” UKAEA FUS 533, 2007.
  8. E. Gonzalez, M. Embid, R. Fernandez et al., “EVOLCODE: ADS combined neutronics and isotopic evolution simulation system,” in Proceedings of the International Conference (MC '99), p. 963, Madrid, Spain, 1999.
  9. F. Álvarez-Velarde, P. T. León, and E. M. González-Romero, “EVOLCODE2, a combined neutronics and burn-up evolution simulation code,” in Proceedings of the 9th Information Exchange Meeting on Actinide and Fission Product P&T, NEA/OECD, Nîmes, France, 2007.
  10. J. Cetnar, “General solution of Bateman equations for nuclear transmutations,” Annals of Nuclear Energy, vol. 33, no. 7, pp. 640–645, 2006. View at Publisher · View at Google Scholar · View at Scopus
  11. J. Leppänen, “Serpent Progress Report 2009,” VTT-R-01296-10, VTT Technical Research Centre of Finland, 2010.
  12. O. Meplan et al., MURE, MCNP Utility for Reactor Evolution—User Guide—Version 1.0, IPNO-09-01, IPN Orsay, 2009.
  13. N. García-Herranz, O. Cabellos, J. Sanz, J. Juan, and J. C. Kuijper, “Propagation of statistical and nuclear data uncertainties in Monte Carlo burn-up calculations,” Annals of Nuclear Energy, vol. 35, no. 4, pp. 714–730, 2008. View at Publisher · View at Google Scholar · View at Scopus
  14. D. L. Poston and H. R. Trellue, “User's Manual, version 2.0 for MONTEBURNS version 1.0,” Tech. Rep. LA-UR-99-4999, Los Alamos National Laboratory, 1999. View at Google Scholar
  15. M. Edenius, K. Ekberg, B. H. Forssen, and D. Knott, “CASMO-4: A Fuel Assembly Burn-up Program User’s Manual,” Studsvik/SOA-95/1, Studsvik, 1995.
  16. T. E. Booth and X-5 Monte Carlo Team, “MCNP—a general monte carlo N-particle transport code, version 5, volume II: user's guide,” Tech. Rep. LA-CP-03-0245, Los Alamos National Laboratory, 2003. View at Google Scholar
  17. A. G. Croff, “A user's manual for the ORIGEN2 computer code,” Tech. Rep. ORNL/TM-7175, Oak Ridge National Laboratory, 1980. View at Google Scholar
  18. I. C. Gauld et al., “Origen-S: SCALE system module to calculate fuel depletion, actinide transmutation, fission product buildup and decay, and associated radiation source terms,” SCALE-6 User Manual, ORNL/TM-2005/39, vol.2, sect. F7, 2009.
  19. W. Haeck and B. Verboomen, “ALEPH 1.1.2—a monte carlo burn-up code,” BLG-1003, SCK·CEN, Brussels, Belgium, 2006. View at Google Scholar
  20. W. Haeck and B. Verboomen, “ALEPH-DLG (Data Library Generator)—creating a general purpose validated application library for MCNP(X) and ALEPH,” BLG 1002, SCK·CEN, Brussels, Belgium, 2005. View at Google Scholar
  21. W. Haeck and B. Verboomen, “An optimum approach to Monte Carlo burn-up,” Nuclear Science and Engineering, vol. 156, no. 2, pp. 180–196, 2007. View at Google Scholar · View at Scopus
  22. A. Schubert, P. Van Uffelen, J. van de Laar, C. T. Walker, and W. Haeck, “Extension of the TRANSURANUS burn-up model,” Journal of Nuclear Materials, vol. 376, no. 1, pp. 1–10, 2008. View at Publisher · View at Google Scholar · View at Scopus
  23. D. Chandler, R. T. Primm, and G. I. Maldonado, “Validation of a Monte Carlo based depletion methodology via High Flux Isotope Reactor HEU post-irradiation examination measurements,” Nuclear Engineering and Design, vol. 240, no. 5, pp. 1033–1042, 2010. View at Publisher · View at Google Scholar · View at Scopus
  24. A. Stankovskiy, G. Van den Eynde, and T. Vidmar, “Development and validation of ALEPH Monte Carlo burn-up code,” in Proceedings of the International Workshop (NEMEA-6 '10), pp. 161–169, Krakow, Poland, October 2010, NEA/NSC/DOC(2011)4.
  25. R. E. MacFarlane and D. W. Muir, “The NJOY nuclear data processing system, Version 91,” Tech. Rep. LA-12740-M, Los Alamos National Laboratory, 1994. View at Google Scholar
  26. A. Santamarina et al., “The JEFF-3.1.1 nuclear data library,” JEFF Report 22, NEA No. 6807, OECD, 2009. View at Google Scholar
  27. A. J. Koning et al., “The JEFF-3.1 nuclear data library,” JEFF Report 21, NEA No. 6190, OECD, 2006. View at Google Scholar
  28. M. B. Chadwick, P. Obložinský, M. Herman et al., “ENDF/B-VII.0: next generation evaluated nuclear data library for nuclear science and technology,” Nuclear Data Sheets, vol. 107, no. 12, pp. 2931–3060, 2006. View at Publisher · View at Google Scholar · View at Scopus
  29. K. Shibata, O. Iwamoto, T. Nakagawa et al., “JENDL-4.0: a new library for nuclear science and engineering,” Journal of Nuclear Science and Technology, vol. 48, no. 1, pp. 1–30, 2011. View at Publisher · View at Google Scholar · View at Scopus
  30. A. J. Koning and D. Rochman, “TENDL-2010: reaching completeness and accuracy,” JEFF-DOC 1349, 2010.
  31. Y. A. Korovin, A. A. Natalenko, A. Y. Stankovskiy, S. G. Mashnik, and A. Y. Konobeyev, “High energy activation data library (HEAD-2009),” Nuclear Instruments and Methods in Physics Research A, vol. 624, no. 1, pp. 20–26, 2010. View at Publisher · View at Google Scholar · View at Scopus
  32. E. Hairer and G. Wanner, Solving Ordinary Differential Equations. Stiff and Differential-Algebraic Problems, vol. 14 of Springer Series in Computational Mathematics, 2nd edition, 1996.
  33. A. C. Hindmarsh, “ODEPACK, a systematized collection of ODE solvers,” in Scientific Computing, R. S. Stepleman et al., Ed., pp. 55–64, North-Holland, Amsterdam, The Netherlands, 1983. View at Google Scholar
  34. L. R. Petzold, “Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations,” Journal on Scientific and Statistical Computing, vol. 4, pp. 136–148, 1983. View at Google Scholar
  35. D. Xie, “An improved approximate Newton method for implicit Runge-Kutta formulas,” Journal of Computational and Applied Mathematics, vol. 235, no. 17, pp. 5249–5258, 2011. View at Publisher · View at Google Scholar
  36. W. Liniger and R. Willoughby, “Efficient integration for stiff systems of ordinary differential equations,” SIAM Journal on Numerical Analysis, vol. 7, no. 1, pp. 47–66, 1970. View at Google Scholar
  37. S. González-Pinto, J. I. Montijano, and S. Pérez-Rodríguez, “Implementation of high-order implicit Runge-Kutta methods,” Computers and Mathematics with Applications, vol. 41, no. 7-8, pp. 1009–1024, 2001. View at Publisher · View at Google Scholar · View at Scopus
  38. H. Bateman, “Solution of a system of differential equations occurring in the theory of radioactive transformations,” Proceedings of the Cambridge Philosophical Society, Mathematical and Physical Sciences, vol. 15, pp. 423–427, 1910. View at Google Scholar