Abstract
A fast numerical method for the calculation in a zero-dimensional approach of the equilibrium isotopic composition of an iteratively used transmutation system in an advanced fuel cycle, based on the Banach fixed point theorem, is described in this paper. The method divides the fuel cycle in successive stages: fuel fabrication, storage, irradiation inside the transmutation system, cooling, reprocessing, and incorporation of the external material into the new fresh fuel. The change of the fuel isotopic composition, represented by an isotope vector, is described in a matrix formulation. The resulting matrix equations are solved using direct methods with arbitrary precision arithmetic. The method has been successfully applied to a double-strata fuel cycle with light water reactors and accelerator-driven subcritical systems. After comparison to the results of the EVOLCODE 2.0 burn-up code, the observed differences are about a few percents in the mass estimations of the main actinides.
1. Introduction
The implementation of the advanced technologies of partitioning and transmutation (P&T) depends on the particular topics that energy policy makers of a country or region choose to optimize because of its singular constraints or motivations. Each of the different possible P&T configurations, that is, the advanced fuel cycle scenarios, will have different objectives that can impact technology choices and performance expected from such systems.
However, the full potential of any P&T policy can be exploited only if the advanced fuel cycle strategy is utilized for a minimum time period of about a hundred years [1, 2]. During this period of time, the irradiated fuels of the advanced reactors are reprocessed to obtain the useful elements for their later recycling as part of the new fresh fuel, in an iterative procedure along successive cycles. The fuel isotopic composition will evolve along a number of cycles, approaching progressively an equilibrium composition, defined as the long-term state reached after infinite cycle iterations.
The calculation of the equilibrium isotopic composition of a transmutation system in an advanced fuel cycle is a matter of special interest since this composition is often used as representative of the nuclear reactor through the whole fuel cycle [3, 4].
In this paper, we describe a fast numerical method for the calculation of the equilibrium isotopic composition in a zero-dimensional approach, taking advantage of the Banach fixed point theorem [5] for the demonstration of uniqueness and existence of a solution. The method is based on the division of the fuel cycle in successive stages: fuel fabrication, storage, irradiation inside the transmutation system, cooling, reprocessing, and incorporation of the external material into the new fresh fuel. The fuel isotopic composition is represented by an isotope vector and will be modified by each fuel cycle stage that is described in a matrix formulation. The fuel cycle is hence represented as a contractor mapping over the fuel composition vector.
With the aim of validating the methodology, a simulation of the fuel cycle with the advanced burn-up code EVOLCODE 2.0 [6, 7] has been performed, checking that the isotopic composition has successfully reached the equilibrium and discussing the range of the possible deviations.
2. Computational Tool
EVOLCODE 2.0 is an in-house development to solve the burn-up problem, that is, the coupled problem of neutron transport and isotopic evolution. Its cycle data flow is shown in Figure 1.
The neutron transport stage is solved by the MCNPX code [8], allowing an important degree of the heterogeneity description in the reactor core model. The isotope evolution stage can be solved by the ORIGEN code [9] or by ACAB [10]. The user selects some geometrical regions or cells inside which neutronic properties and material composition are considered constant. For each of these cells, EVOLCODE 2.0 creates one-group effective cross-section libraries using the neutron flux energy spectrum provided by MCNPX in such a way that the whole calculation is faster than allowing MCNPX to do the job. Then, one ORIGEN execution is made independently for each cell providing the isotopic evolution of the materials. The ORIGEN executions are made using the neutron flux and materials at the beginning of the time step. In case that the averaged reactor power (estimated by the ORIGEN simulation) is not equal to the desired value or if its variation along the time step is large (causing a large uncertainty in the final results [7]), then the predictor/corrector method is activated to fix it. This whole procedure is repeated for several successive irradiation steps (or cycles) to reach the total burn-up specified by the user.
3. Double-Strata Fuel Cycle Scenario
One of the proposed advanced fuel cycle scenarios for waste reduction is the so-called double-strata scenario, whose scheme is shown in Figure 2. This fuel cycle has been studied in detail in the frame of different projects [11, 12] and consists in two separated strata or stages, one dedicated to energy production and the other dedicated to waste consumption.
In this particular specification of the double-strata scenario, the first stratum of the fuel cycle consists in the irradiation of UO2 in light water reactors (LWRs) and the later advanced Purex reprocessing of the irradiated fuel for the Pu reutilization, only once, as MOX, again in LWR. Recovered MA coming from the first partitioning process and Pu and MA coming from the reprocessing of the MOX irradiated fuel are reutilized as fuel for the accelerator-driven subcritical system (ADS) in the second stratum. This second stratum is based on a fast spectrum ADS, which operates with continuous recycling of the main actinides (U, Pu, Np, Am, and Cm). This recycling is supposed to be a pyro-metallurgical process.
4. ADS Fuel Isotopic Composition
The isotopic composition of the ADS fuel depends on the LWR park since the LWR recovered material (from both UO2 and MOX reprocessed fuels) becomes a part of the ADS fuel content in each new ADS irradiation (a 15% approximately). The other part of the ADS fuel comes from the pyro-reprocessing of the preceding ADS irradiated fuel. In different cycles, the ADS fuel isotopic composition would change (and being consumed), but maintaining similar ADS operating conditions and a constant additional supply per cycle of actinides from the first LWR-stratum, the fuel approaches an equilibrium (as it is demonstrated below) in which the fuel isotopic composition of each cycle is equal to the following one.
The iterative second stratum of this fuel cycle scenario can be divided in a succession of stages: fuel fabrication, storage, irradiation inside the ADS, cooling, reprocessing, and incorporation of the external material (from the first stratum) into the new fresh fuel. Each process can be described in a matrix formulation as an operator over the fuel composition vector , obtained just before the th ADS irradiation. Vectors and matrices length is determined by the total number of different isotopes considered in the problem. This number of isotopes is equal to in this study (the detailed list is shown in Table 1), although a larger value of has been used in other problems. If is the fuel composition vector before the th ADS irradiation, it is related to following the expression where is the matrix operator representing the irradiation inside the ADS, represents the cooling decay time after the irradiation (two years in this scenario), is the reprocessing operator, represents the storage time before irradiation (one year), is the material incorporation from the first stratum (taken at the moment just before the beginning of the ADS irradiation), and . Since the fuel composition vector can be expressed as , the process fulfils the requirements of the Banach fixed point theorem, and therefore there is a solution vector under the requirement that is a contractor mapping, that is, if there is a constant such as where the symbol represents the norm of the operator. Applying , this inequality can be written as which is satisfied if . In order to prove that this condition is satisfied, we must analyze the norm of each of the matrix operators appearing in the definition of . The norm of the decay operators and is equal to the unity because they only change one isotope into another leaving the total amount of material unchanged (being precise, it might be smaller than the unity as some decay products are not included in the dimensions of the composition vector). The reprocessing operator, , is a diagonal matrix having a constant diagonal value smaller than the unity (reprocessing recovery factor of 0.999) for the reprocessing isotopes (only the reprocessing elements U, Pu, Np, Am, and Cm have been considered), hence . Finally, the norm of the irradiation matrix is strictly smaller than the unity because it represents consumption of heavy material (approximately a 15% for this ADS burn-up). This can be easily shown when the norm is defined as the maximum and demonstrated for other definitions. For all these reasons, , therefore this process fulfils the Banach fixed point theorem and there exists a unique solution of the equation of the type , for an arbitrary .
Coming back to the previous notation, we can find the solution directly from the equilibrium equation; hence the expression of the fuel composition vector can be written as
The solution of this equation is
In order to reach this solution, it is necessary to calculate , , , , and in a matrix formulation from the data of the problem. The operator is represented by a diagonal matrix with a constant value of 0.999, as mentioned above. is a vector containing the contribution of the LWR stratum so it can be provided after solving the first noniterative stratum of the fuel cycle. The operator can be represented by the generic evolution of a nuclear system under irradiation. The generic evolution of the isotopic composition vector can be described in the zero-dimensional approximation as where is the cross-section matrix, is the radioactive decay matrix, and is the averaged neutron flux. Operators and can be obtained in case that . The formal solution of this coupled linear system of equations of constant coefficients is
In the base of eigenvectors of the matrix , this solution can be written as where is the matrix of the eigenvectors associated with the eigenvalues . The calculation of the eigenvalues and eigenvectors needed to compute the exponential and the inversion of the matrix has been a difficult problem because of the large numerical instabilities involved in the manipulation of these large matrices. However, the easy access to arbitrary precision arithmetic, made recently available in popular computer codes as Maple [13] and Mathematica [14], allows using brute force and direct methods of solving these problems. For our case, it was necessary to set the arithmetic precision to 70 decimal digits to obtain the required stability. In spite of this condition, the eigenvectors computed had to be regularized by discarding very small imaginary components (). The process was implemented with the Maple program, where efficient procedures for the calculation of eigenvalues, eigenvectors, and matrix inversions were readily available. The typical calculation time is few minutes in a modern PC.
In order to estimate the flux level and the one-group zero-dimensional equivalent cross-sections, a fully detailed transport simulation was performed. In this simulation, the exact geometry, the required total power, and an initial guess of the equilibrium composition of the reactor materials were used. Careful weighting had to be used in the procedure of averaging the cross-sections over the whole reactor volume, which was divided in different cells with isotopic evolution. For every one of these cells, one set of one-group effective cross-sections were calculated using the EVOLCODE 2.0 simulation system. The final reactor-averaged cross-sections were calculated weighting the different one-group cross-sections of each fuel cell by its fuel mass. The initial guess of the fuel composition was set to the same vector as the vector of actinides from the external refuelling, .
Since the initial ADS fuel composition, chosen for the calculation of matrices and , is different from the final equilibrium composition, this method of calculation of the equilibrium isotopic composition of the ADS fuel introduces an error in the results. An iterative process can be carried out to increase the precision of the method. Once the equilibrium isotopic composition is calculated using the previous noniterative method, the resulting isotopic composition can be used as the initial ADS fuel content and a new set of cross-sections/neutron flux could be obtained with EVOLCODE 2.0. Then, the method is repeated in order to obtain a more precise value of the equilibrium isotopic composition. The new procedure can be successively performed until the isotopic composition of a certain iteration is deviated less than a limit value from the previous iteration. Nevertheless, one iteration is typically enough to achieve precisions acceptable for fuel cycle studies.
5. Results
The isotopic composition of the different stages of the iterative stratum of the fuel cycle has been calculated using the matrix equation method. Table 1 shows the initial (BOL, beginning of life, before the irradiation) ADS fuel isotopic composition (in terms of actinides), calculated as the solution of the equation for described in the previous section. Moreover, the isotopic composition of the ADS fuel at BOL has also been calculated for the following cycle by means of the advanced burn-up code EVOLCODE 2.0 for checking purposes.
The deviation between the equilibrium isotopic compositions obtained from the matrix equation method and from the irradiation (and the successive decay and reprocessing as indicated in the fuel cycle scenario) using EVOLCODE 2.0 can be seen in the fourth column in Table 1. These deviations, due to the different accuracy of both methods, are smaller than 3% for all actinides (with significant contributions to the total mass) except for Cm-243, with 5.7% difference.
The objective of the fuel cycle scenario will define the maximum allowed deviation between one cycle and the following one and also the final decision concerning the necessity of a second iteration on the matrix equation method if the equilibrium isotopic composition has not been reached. For instance, for a fuel cycle scenario where the objective is the study of a representative set of the different (primary and secondary) waste streams, the obtained equilibrium isotopic composition is a good approximation so no iterative calculations are needed.
6. Conclusions
A fast numerical method, based on an arbitrary precision arithmetic solution to the matrix equations, for the calculation in a zero-dimensional approach of the equilibrium isotopic composition of a transmutation system in an advanced fuel cycle has been developed and successfully applied to a double-strata fuel cycle with LWR and ADS.
The observed differences between the matrix equation method and the EVOLCODE 2.0 results for the successive cycle are smaller than a few percents in the mass estimations of the main actinides. The precision of the calculation can be improved, if needed, by performing a second iteration on the matrix equation method, using the one-group cross-sections of the ADS with the equilibrium composition obtained in the first iteration.