Probabilistic approach of solving burnup problems

This paper presents a new approach to solve nuclear burnup problems using a probabilistic method. Unlike the traditional methods that rely on complex matrix exponential calculations, the proposed method tracks the transformation of nuclei over a period of time, enabling the estimation of nuclide concentrations. The method is implemented in a C++ program CNUCTRAN and verified against the Chebyshev Rational Approximation Method (CRAM). Three sample calculations are presented, showing a detailed comparison of results obtained using CNUCTRAN and CRAM48. The relative error analysis emphasizes the accuracy of the probabilistic method, with the most significant error observed below 0.001%. The computational efficiency is discussed in terms of CPU time, indicating that CNUCTRAN requires more CPU time than CRAM but with improved accuracy as the time step decreases. The numerical results for various burnup problems demonstrate the accuracy and efficiency of the probabilistic method, making it a promising alternative for simulating realistic nuclear transmutations.


Introduction
Analyzing nuclear reactors requires careful consideration of how radioactive materials change over time and in response to nuclear reactions, for instance, decay, neutron capture and fission.The burnup equation describes the behaviour of the material composition over time in a nuclear reactor, which depends on the reactor's core parameters and material compositions [1].The equation is usually a set of first-order linear differential equations.These equations are stiff, especially when the nuclides included in the calculation combine short-lived and long-lived ones.The Chebyshev Rational Approximation Method (CRAM) [2] and Transmutation Trajectory Analysis (TTA) [3] are two methods that earlier academics used to address challenging situations [4], and both are computationally efficient but are mathematically abstract.On the other hand, [5] recently found an unconventional way of solving burnup problems by deriving a probability distribution that predicts whether a nucleus will transform into another species.This paper presents a new non-deterministic approach to solve burnup problems, and the method is implemented in a C++ program named CNUCTRAN to verify the method.

Burnup equation
The burnup equation can be expressed in Eq. (1) below.
where y i is the density of nuclide i, λ i is the decay constant of nuclide i, b ij is the branching ratio from nuclide j into nuclide i, f ik is the production yield from nuclide k into nuclide i and A ik is the production rate of nuclide i due to the removal of nuclide k from the system.Eq. ( 1) is a system of I linear differential equations dependent on one another.The solution of Eq. ( 1) is the concentration of I nuclides.The term on the left-hand side of Eq. ( 1) is known as the rate of change of the nuclide-i concentration.The left-hand side of Eq. ( 1) is balanced with the terms on the right-hand side.Here, each term describes the various transmutations of nuclides that lead to the production and removal of nuclide-i within the system.
In particular, the first term on the right-hand side of Eq. ( 1) represents the natural decay of various nuclides that transform into nuclide-i.In contrast, the second term of the right-hand side of Eq. ( 1) describes the formation of nuclides-i via different transmutation reactions of all nuclides, excluding the natural decay reaction.Finally, the last term on the right-hand side of Eq. ( 1) is a combination of natural decay and mixed transmutation reactions that transform nuclide-i into other nuclides, which removes it from the system.Alternatively, Eq. ( 1) can be expressed in its compact matrix form given by, where A is now known as the transmutation matrix.Also, A is a sparse I × I Square matrix, where most elements are zero, and y is a column vector containing the instantaneous concentration of all I nuclides defined in the calculation.Technically, the norm of A is in order of 10 21 , and includes both off-diagonal elements and diagonal elements.However, using the separation of variables Eq. ( 2) can also be expressed as given below: Alas, evaluating the matrix exponential in Eq. (3) proves to be a complex task.One might attempt to simplify the process by employing a Taylor series exponential or Padѐ approximation of the matrix exponential, which is unfortunately impractical and does not always lead to accurate and reliable results, especially for stiff burnup problems.
Pusa [6] observed that the eigenvalues of matrix A are situated in close proximity to the negative real axis.This allows for the estimation of the matrix exponential using Chebyshev approximants.This approach is commonly referred to as the Chebyshev Rational Approximation Method (CRAM), and it provides an approximation for the matrix exponential as follows: where a i and θ i are the complex numbers published by Pusa [6] for the approximation order k ∈ {16 ,32 , 48 }.

The probabilistic method
In the previous study, the establishment of π ij (∆t ) (probability distribution) allows us to estimate the probability of various transmutation events experienced by a nucleus.The primary objective of this work is to estimate the time evolution of nuclide concentrations without directly solving Eq. ( 1) mathematically.This endeavour is underpinned by the utilization of the probability distribution given by [7], where γ i is the normalization constant and δ jl is the Kronecker delta.Also, Λ ij is the reaction rate constant of reaction j of nuclide i, and it corresponds to λ and A in Eq. ( 1).The value of γ i is given by, Generally, the concentration of nuclides i at a given time t +∆t, is contributed by two factors.Firstly, it is contributed by the un-reacted nuclides i within the interval (t , t+∆ t ).Secondly, it is contributed by the transmutation of other nuclides that occur within the interval (t , t+∆ t ), which leads to the production of nuclide i.Therefore, the concentration of nuclide i after the next sub-step, t +∆t , can be expressed as, where π i →i is the probability of nuclide i not undergoing a transmutation event; π l →i is the probability of any nuclides l undergoing a transmutation event that transforms themselves into nuclide i .By definition, where R is a set of transmutation events that transform nuclide l into daughter nuclide i.If removal event j is a fission reaction, then π j in Eq. ( 8) must be scaled to the fission yield of the daughter nuclide i. Eq. ( 7) can be systematically written in matrix form via the definition given in Eq. ( 8).Here, by considering all I nuclides in a burnup problem, ( Or in a simple form, Notice that all elements in P are calculated based on the sub-step interval, ∆ t.Eq. ( 10) is a recursion relation where the future nuclide concentrations, y (t +∆ t ), can be obtained by transforming the nuclide concentrations of the previous time interval via the transfer matrix, P. Here, P serves as the transformation operator.The column indices of P correspond to the various parent nuclides, and the row indices of P correspond to the different daughter nuclides.Notice that the diagonal elements of square matrix P are probabilities of no transmutation events occurring.The P element can be calculated using Eq. ( 7), provided all reaction rates associated with the nuclides are known.
Suppose that the evolution of nuclides concentration is monitored from the initial time, t=0, to the final time, t=T .Then, a small-time interval, ∆ t, is chosen such that the total number of sub-steps is given by l=T / Δt.Therefore, y (t ) can be obtained by repeating Eq. ( 10), for l times, giving, where y (0)≡ y 0 is the initial nuclides concentration at t=0.Here, the smaller the value of ∆ t, the larger the number of sub-steps, thus the better the accuracy [7].However, such a computation with small ∆ t requires arbitrary precision arithmetic to avoid the floating-point error.The evaluation of sparse matrix exponentiation, P l , can be accomplished at an incredible speed via exponentiation by squaring.Evaluating P l requires l repeated multiplications of P. Such a naïve approach will incur an unreasonably long CPU time, especially when l is large.Gratefully, the exponentiation by squaring helps to evaluate P l with only log 2 l Multiplications [8].The trick is that P l can be efficiently calculated by repeated squaring of P. Since evaluating P l also involving repeated matrix multiplications, the conventional matrix multiplication algorithm will impose unnecessary operations on zero elements of the sparse matrix, P. Therefore, this work also implements the sparse matrix multiplication (SpMM) algorithm proposed by [9].
Remember that the probabilistic method requires l self-multiplications of the transfer matrix P. in the previous section, such a computation can be done at an incredible speed via exponentiation by squaring.To save computational time, the number of substeps is chosen as where k is the total number of self-multiplications of P. Our objective is to obtain the suitable value of k, such that the substeps interval Δ t is in order of 10 −n .Here, n is known as the approximation order, and the time step, T, of the calculation is given below: The floor operation in Eq. ( 13) is needed because k must be a positive integer number.Also, ∆ t must be corrected to compensate for the error due to floor operation:

Sample calculations
The proposed method is implemented in a C++ program named CNUCTRAN and has been verified against the state-of-the-art method, CRAM.The nuclides data, such as the half-lives and other transmutation chains of over 4000 nuclides, are fetched from an XML file containing the transmutation information of various nuclides derived from ENDFBVII.1.The results obtained were compared with CRAM of order 48 (CRAM48).Table 1 summarizes all the information of sample 1,2 and 3 calculations.

Sample Calculation I: decay of U238 and U235
The first sample calculation is a decay reaction with a total of 1489 nuclides that involve U-235 and U-238 with initial concentrations of 0.1 mol and 0.9 mol, respectively.The system of nuclides is under a constant burnup such that the (n, a), (n,γ) and fission reaction rates are kept constant at 10 − 4 sec −1 , 10 − 4 sec −1 and 10 −5 sec −1 .The sample monitored the concentrations of all nuclides involved for the time step of 8.64E4 sec.The results obtained from the proposed method exhibit negligible differences with the reference solution.Interestingly, 1249 nuclides exhibited zero concentrations in both the proposed methods, while 224 nuclides maintained non-zero concentrations.Additionally, 16 displayed negative CRAM concentrations, coinciding with zero concentrations in the proposed method.Notably, the relative error between solutions computed using CNUCTRAN and CRAM48 is illustrated in Fig. 1, offering a visual representation of the comparative accuracy of these computational methods.

Sample Calculation II: U238 and U235.
The second sample calculation is similar to the first sample which tracks 1603 total number of nuclides over the time step of 8.64 × 10 6 sec (100 days).The results obtained show 224 nuclides with non-zero concentration.25 nuclides with negative CRAM concentration, and these nuclides have zero concentration on the proposed method.1354 nuclides with zero concentration on both two methods.A comprehensive assessment of the relative error between solutions computed using CNUCTRAN and CRAM48 is presented in Fig. 2.This graphical representation illustrates the comparative accuracy of these computational methods and provides a visual insight into their performance across the spectrum of tracked nuclides.The results obtained from the three distinct samples exhibit no significant differences from the reference solution, and the largest error observed is less than 0.001%.It is worth highlighting that the overall accuracy of CNUCTRAN improves as the time step size decreases.Therefore, this phenomenon is a known limitation of the probabilistic method outlined in section 3. It is important to note that only nuclides with concentrations 10 −40 in three samples, respectively, were considered to minimize rounding errors.Additionally, the relative error increases with decreasing nuclide concentrations.Furthermore, 100% of the important nuclides have a relative error of less than 10 −5 , indicating the outstanding accuracy of CNUCTRAN.For the benefit of readers, the CPU times for CNUCTRAN and the CRAM method are in Table 2.It is noteworthy that CNUCTRAN requires more CPU time than the CRAM method; as discussed in section 3, the CPU time of CNUCTRAN increases as the time step length increases.

Conclusion
In summary, the probabilistic method used an unconventional way to simulate the time evolution of nuclide concentrations.It stems from the direct simulation of nuclides transformations, and the implementation requires elementary probabilistic theories that can be easily understood.The presented method enables the inclusion of short-lived and long-lived nuclides in a single calculation.Even though it can stimulate realistic nuclide transmutations, it has a few limitations that necessitate future development.For instance, dealing with slowly evolving concentrations of long-lived nuclides requires arbitrary precision arithmetic to account for such a subtle change.Following that, the method requires a huge number of sparse matrix multiplications.Here, the speed of these multiplications depends on the computing.Further research and implementation of this method should focus on improving the speed of arbitrary precision sparse matrix multiplication.

Figure 2 . 7 Figure 3 :
Figure 2. Relative error versus sample II after 8.64 × 10 6 sec.4.3.Sample Calculation III: pure decay of Np-237This sample calculation involves a pure decay reaction consisting of the linear decay chain of Np-237 with an initial concentration of 1 mol, and the nuclide concentrations were monitored over the time steps of 3.154 × 10 10 sec.The sample involved 1632 total number of nuclides starting from atomic mass number of A =137 to A = 237.The obtained results shed light on the behaviour of these nuclides, revealing 15 with non-zero concentrations.Interestingly, most 1617 nuclides exhibit zero concentrations in the CRAM and CNUCTRAN, underlining a consistent behaviour across these computational approaches.To provide a quantitative perspective on the precision of the solutions, the relative error between computation using CNUCTRAN and CRAM48 is shown in Fig.3.This graphical representation serves as a visual aid, offering a clear comparison of the accuracy of these computational methods in modelling the dynamic evolution of nuclides concentrations.

Table 2 .
CPU Time (in seconds) comparison between CNUCTRAN and CRAM48.Time step, t (s)