Two alternative approaches to the solution of cyclic chains in transmutation and decay problems

Journal Pre-proof Two alternative approaches to the solution of cyclic chains in transmutation and decay problems Carlos Antonio Cruz López, Juan Luis François

PII: DOI: Reference:

S0010-4655(20)30064-3 https://doi.org/10.1016/j.cpc.2020.107225 COMPHY 107225

To appear in:

Computer Physics Communications

Received date : 16 September 2019 Revised date : 22 January 2020 Accepted date : 18 February 2020 Please cite this article as: C.A. Cruz López and J.L. François, Two alternative approaches to the solution of cyclic chains in transmutation and decay problems, Computer Physics Communications (2020), doi: https://doi.org/10.1016/j.cpc.2020.107225. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier B.V.

Journal Pre-proof

lP repro of

Two alternative approaches to the solution of cyclic chains in transmutation and decay problems Carlos Antonio Cruz López and Juan Luis François Universidad Nacional Autónoma de México, Facultad de Ingeniería, Departamento de Sistemas Energéticos Paseo Cuauhnáhuac 8532, Col. Progreso, Jiutepec, Morelos, México [email protected], [email protected]

Abstract

rna

In several burnup and activation problems there is a recurrent issue related to the singularities in the Bateman’s solution due to its inability to solve linear transmutation schemes where there are repeated isotopes, or where there are two different isotopes with the same removal coefficients. Most of these types of transmutation schemes are called cyclic chains. In these cases, the Bateman’s solution fails due to the presence of subtractions between the mentioned coefficients in some denominators, which eventually become zero and undefined. In order to overcome this problem, two methodologies have been reported in the open literature. The first one consists in introducing small modifications in the repeated removal coefficients, preventing the presence of zeros in certain denominators. The second is to develop more general equations for the Bateman’s solution. Nevertheless, both methodologies are based on the approximation of the cyclic chains with the linear chain method, whose error has not been studied until now. The study of this error is fundamental to omit cyclic chains and to reduce the execution time of the algorithm. In the present work, a more general approach to the cyclic chains was studied, starting with the description and classification of the transmutation and decay networks that generate them. Afterward, two different approaches to solve some of these structures were proposed, which are not based on the linear chain method. One of them is based on a power series analysis, and the other one is related to a numerical analysis of the roots of a polynomial in the Laplace transform space. Additionally, computer algorithms were developed for each approach to facilitate their implementation in a burnup or activation code, and a numerical comparison with the linear chain approximation was carried out. Through the present work, it was possible to compute the actual error involved when the linear chain is used for approximating cyclic chains, and to conclude if a cyclic chain can be ignored in a burnup problem. 1. Introduction.

Jou

In several problems of nuclear engineering it is necessary to compute the concentration of a set of nuclides which undergo a process of successive transformations. The basic mathematical formulation that is used to model this problem consists of a system of coupled linear differential equations. In 1910, Harry Bateman Bateman 1910 proposed a solution that gained notoriety due to the methodology that the author used to solve this problem: the Laplace method. Originally, the Bateman’s solution was used to model successive transformations due to radioactive decay processes. Nevertheless, over the years, it was adapted to solve cases where some of the transformations were caused by the presence of the neutron flux. Additionally, through a process of break-down and linearization, the Bateman’s solution was extended to solve more complex problems with networks and ramifications. These two extensions of the application of the Bateman’s solution resulted in the methodology known as “linear chain 1  

Journal Pre-proof

lP repro of

method”, which was formally developed by T. S. England in 1962 England, 1962 . Among the recent research topics that are related to the linear chain method, it is the development of general solutions, which overcome one of the most remarkable limitations of the Bateman’s solution: the presence of repeated elements in the transmutation structures. One of the conditions in the original Bateman’s formulation was the assumption that all the isotopes in a linear chain were different. Such condition guarantees that several subtractions in the denominator do not become zero, and therefore it allows that the solution is well defined. Nevertheless, in several practical cases, transmutation and decay networks appear that lead to the generation of linear chains with repeated elements, and therefore the Bateman´s solution fails. The structures where this occurs are called cyclic chains. Two approaches were developed in order to solve the difficulty of the singularities. The first one is to introduce small modifications in the repeated removal coefficients, which avoid the divisions by zero. The second approach is to develop more general solutions. Nevertheless, all general solutions can model cyclic chains, but strictly speaking they do not solve the transmutation and decay networks that originate them. This is because the process of linearization only approximates the decay and transmutations networks with loops or cycles, and therefore an error is introduced when linear chains are used to model them. In the present work, the decay and transmutations schemes that contains loops or cycles are studied from a general approach without using linear chains to model them. A classification of these structures is proposed, identifying the simplest one and developing two methodologies to solve them. The first methodology is based on a solution expressed in terms of the equation of a linear chain, using a power series expansion. The second method consists of building a solution to the cyclic chains using the roots of a polynomial in the Laplace transform space.

rna

Through these approaches it is possible to determine the actual error of approximating the cyclic chains through the linear chain method and to conclude if a cyclic chain can be ignored. This is very important because it allows reducing the execution time of a burnup problem. As an example of the application of the developed solutions, a numerical comparison analysis using a cyclic chain composed of heavy metals that begins with U-235 and ends with Pu-239 is proposed.

Jou

The structure of this paper is as follows: Section 2 contains a brief discussion of the Bateman’s solution using the Laplace transform and a brief description of the linearization process. In Section 3 an analysis of the cyclic chains is carried out, starting from their structure and ending with a power series solution. Section 4 contains a study of how the power series solution can be extended. Section 5 presents the development of a solution of pure cyclic chains through the roots of a polynomial in the Laplace transform space, as well as the use of an integral method to extend its application. In Section 6 a numerical analysis is depicted and the conclusions are in Section 7. 2. The Bateman’s solution and the linearization process. 

2.1 Bateman’s solution and the Laplace transform. In several fields of nuclear engineering the following successive transformations appears: 𝑋

,

𝑋

,

,

… ⎯⎯ 𝑋 2

 

1

Journal Pre-proof

lP repro of

Where 𝑋 , 1 𝑖 𝑛, denotes the isotope in the 𝑖 position, which is originated due to the 𝑟 , reaction. This structure, where except for the first and the last isotopes all the other elements have only one father and only one daughter, is known as a linear chain. For the case of the radioactive decay, it is possible to model this process using the following set of coupled linear differential equations1: 𝑑𝑋 ⎧ 𝑑𝑡 ⎪ ⎪ 𝑑𝑋 𝑑𝑡 ⎨ ⋮ ⎪ ⎪𝑑𝑋 ⎩ 𝑑𝑡 Where 𝜆 , 1

𝑖

𝜆 𝑋

𝜆 𝑋

𝜆 𝑋

2

⋮

𝜆

𝑋

𝜆 𝑋

𝑛 are the decay constants. Considering the following initial conditions: 𝑋 𝑡

0

𝑋 0

0 with 2

𝑖

𝑛, 𝑋 0

0

3

The solution of system 2 is equal to: 𝑋 𝑡

𝑋 0

𝜆

𝑒

1

𝜆

𝜆

4

rna

The last expression is known as the Bateman´s solution, and originally it was set to model structures as the one described in (1), i.e. linear chains. Nevertheless, linear chains are not common in burnup and activation problems, and it is more probable to find more complex structures called transmutation and decay networks, like the one shown in Figure 1. Fortunately, the Bateman solution can be applied using the linearization process, a method that essentially consists of searching all the possible paths in a network, which will be the linear chains of it. In other words, the linearization process transforms a transmutation and decay network in a set of linear chains. 3. Cyclic chains

3.1 Definition of cyclic chains and the approximation using the linear chain.

Jou

There is a special type of networks where the linearization process fails. Such networks contain structures called “cyclic chains”, which appear commonly in several burnup and activation problems. Tasaka Tasaka, 1980 defined it as “a loop chain where a nuclide transforms certain times and then transforms to the original nuclide”. Based in Tasaka’s definition, a cyclic chain can be represented in a network diagram through a closed loop between two or more nuclides. For example, if the term 𝑋 is replaced by the term 𝑋 in the Figure 1, as result it will appear a cyclic chain given by the succession of isotopes 𝑋 → 𝑋 → 𝑋 → 𝑋 → 𝑋 → 𝑋 . Such a cyclic chain is shown in Figure 2, where these isotopes have been highlighted in yellow.

 

 Henceforth, the symbols 𝑋 will denote the concentration functions of the isotopes when they appear in balance equations. 

1

3  

lP repro of

Journal Pre-proof

Figure 1. Example of a transmutation and decay network.

Figure 2. Example of cyclic chain.

rna

It is not possible to use equation 4 to solve this loop or cycle due to a missing gain term given by 𝜆 𝑋 for the balance equation of 𝑋 . Except for this term, all the other balance equations of the loop can be perfectly described by the equations of the form given in 2 . Therefore, a first approach to find the concentration of the isotopes that are included in this loop is to assume 0, which is reduced to treat these isotopes as if they belong to a linear chain. It is that 𝜆 𝑋 possible to improve this approximation building an artificial succession of isotopes given by 𝑋 →𝑋 →𝑋 →𝑋

→𝑋

→𝑋 →⋯

5

Jou

In other words, the element 𝑋 can be added at the end of the linear chain in order to “simulate” the contribution 𝜆 𝑋 to 𝑋 .This methodology, where a repeated element is added to a linear chain, is equivalent to build artificial linear chains, and it is known as the approximation of a cyclic chain through the linear chain method. 3.2 Implications of simulating the loops through the linear chain method. When the isotope 𝑋 is added at the end in (5), there is a duplication of all the elements that are originated from the decay and transmutation of this isotope. Particularly, since this element is the beginning of the structure in the example showed in Figure 2, it implies that all the transmutation and decay networks appearing after it will be repeated as it is represented in Figure 3, where only three times of this process of duplication are showed. 4  

lP repro of

Journal Pre-proof

Figure 3. Duplication of the decay and transmutation network of Figure 2. This scheme is the result of building an artificial linear chain. In other words, the extension of the linear chain in (5) will be infinite and the transmutation and decay network will have a finite number of isotopes that appear an infinite number of times. Therefore, the use of the linearization process cannot be applied in such structures, because such procedure will not have an end.

rna

Wilson Wilson, 1999 proposed a method to overcome, to some extension, this problem. Essentially it is possible to approximate the solution of the linear chain, limiting the process of duplication to an appropriate number of times. Nevertheless, as it was showed in a past work Cruz-López and Francois, 2018 , the number of generated linear chains grows in an exponential way for the cases of cyclic chain, and therefore it is important to study alternative ways to approach this problem, in order to provide a more adequate methodology. As it will be showed later, when the equation and system related to a cyclic chain is solved without approximations, it is possible to reduce the number of generated linear chains, and therefore decreasing the execution’s time.

Jou

In addition to the problems with the number of linear chains, there are several issues related to the solution of the artificial linear chains described in 5 . Essentially, equation 4 does not admit repeated lambda coefficients and therefore it is not possible to use it. This limitation has its origin in the assumption that all the lambdas were different when the process of partial fraction decomposition was carried out. There are two ways to overcome this difficulty: the first one consists in introducing small modifications on the lambda coefficients Isotalo, 2013 , and the second one is related to develop more general solutions which do not use the mentioned assumption Cetnar, 2006 Dreher, 2013 . 0, which Nevertheless, these two approaches assumed conditions similar to 𝜆 𝑋 represented a kind of approximation, which until now has not been studied or analyzed in literature. Even more, such methodology produces a large number of linear chains through the linearization process, which becomes an advantage in terms of computational time. Considering the last discussion, the objectives of the present work are: 1 To find what is the error involved in such methodology, 2) regarding this last error, to determine the cases where this methodology can be applied, 5  

Journal Pre-proof

lP repro of

3) to reduce the number of generated linear chains, ignoring cyclic chains whose contribution can be neglected, 4) to develop an algorithm that facilitates the implementation of the cyclic chain’s solution. 3.3 A first approach to the cyclic chains’ classification.

As for its complexity, it is possible to divide the cyclic chains into two classes: the pure and the compound. The first one can be defined as a loop, where all the nuclides belong to the same cyclic chain. The compound type is the structure where at least one nuclide belongs to at least two different loops. Figure 4 shows some examples of this classification. Among the structures of the first type, there are two additional classes regarding the place where the loop is present: the “first position” and the “on-path position” class. In the first two diagrams it can be observed that all the elements in the cyclic chain have only one daughter and one father, except for the isotope denoted by the gray rectangle. The gray rectangles in the compound’s type correspond to isotopes that belongs to more than one cyclic chain, in other words, isotopes that have more than one daughter. The straight arrows correspond to linear chains, while the curved lines represent a succession of elements in a cyclic chain. As a first approach to the problem of cyclic chains, it will be found the solution for the pure’s type and, starting from here, a methodology that allows treating more general cases will be developed. 3.4 Analysis for the pure type of the first position class.

rna

A pure cyclic chain can be denoted by the succession of isotopes given by 𝑋 → 𝑋 → ⋯ → 𝑋 → 𝑋 , with 𝑛 ∈ ℕ. In a strictly sense, there is a limit value for the number 𝑛, which corresponds to the largest pure cyclic chain that can be found according to all the possible transformations of the isotopes. In other words, 𝑛 is limited by physically means. As it will be discussed later, for the case of practical nuclear engineering problems, it follows that 𝑛 7. For the cyclic chain’s structure, it can be set the following balance equation system: 𝑑𝑋 ⎧ 𝑑𝑡 ⎪ ⎪ 𝑑𝑋 𝑑𝑡 ⎨ ⋮ ⎪ ⎪𝑑𝑋 ⎩ 𝑑𝑡

𝜆 𝑋

𝜆 𝑋

𝜆 𝑋

𝜆 𝑋

ℒ

⋮

𝑋

Jou

𝜆

𝜆 𝑋

⎧𝑥 ⎪ ⎪ 𝑥 ⎨ ⋮ ⎪ ⎪𝑥 ⎩

⋮

𝑋 0 𝜆 𝑥 𝑠 𝜆 𝜆 𝑥 𝑠 𝜆 ⋮ 𝜆 𝑥 𝑠 𝜆

6

After a process of substitutions, beginning with 𝑥 and ending with 𝑥 , the last equation of the system in 6 can be written as: 𝑥

𝑠

𝜆

𝛼 𝑛 𝑠 𝜆

1 𝑋 0 … 𝑠 𝜆

𝛼 𝑛

6  

, with 𝛼 𝑛

𝜆

7

lP repro of

Journal Pre-proof

Figure 4. Classification of the cyclic chains.

From this point forward, the solution obtained through the linear chain method can be denoted by 𝑋 , and the solution of the cyclic chain will be represented as 𝑋 . It is possible to show that if 𝛼 𝑛 0 then 𝑋 → 𝑋 . Equation 7 can be arranged as: 𝑥

𝑠

𝛼 𝑛 𝜆 𝑠

1 𝑋 0 𝜆 … 𝑠

𝜆

𝑠

Using the assumption 𝛼 𝑛 | 𝑠 𝜆 expansion, it is possible to write 8 as: 𝑥

𝑋 0 𝜆

𝑠

𝜆

1

𝑠

𝑠

𝜆

1 𝜆

1 𝛼 𝑛 𝑠 𝜆 … 𝑠

𝜆

… 𝑠

8

𝜆

𝜆 |, and though a power series

… 𝑠

𝛼 𝑛

𝜆

𝑥

rna

The term inside the brackets in the last equation can be written as a sum, using partial-fraction decomposition: 𝑋 0 𝜆

𝑎

𝑠

𝜆

,

𝛼 𝑛

𝑎

1 𝜆

𝜆

Assuming the series is absolutely convergent it is possible to introduce the inverse Laplace operator into the summation and carry out the inverse Laplace transform: 𝑋 0 𝜆

Jou 𝑋

ℒ

𝑎 𝑠

𝜆

𝛼 𝑛

Using the following identity Cruz-López, 2019 : ℒ

𝑎 𝑠

𝜆

𝛼 𝑛

It is possible to write 9 as:

7  

𝑋 1 𝑋 0

9

Journal Pre-proof

𝑋 0 𝜆

𝜆 𝑋 𝑋 0

𝑋

𝑋 0 𝑋 0 𝜆 𝑋

10

lP repro of

𝑋

As it can be observed, the last equation solves the cyclic chain under the assumption that 𝛼 𝑛 | 𝑠 𝜆 𝑠 𝜆 … 𝑠 𝜆 |. This first solution, based on power series, shows that for small values of 𝜆 𝑋 , the linear chain solution is suitable for modeling cyclic chains. 3.5 The limit value of 𝒏.

In the previous section it was mentioned that the value 𝑛 in a pure cyclic chain is limited due to physical means. In other words, such value has as an upper bound given by the length of the longest pure cyclic chain, which is finite and is determined by the relationships between the isotopes and the decay and transmutation reactions. Tasaka proposed the largest pure cyclic chain, considering the following structure Tasaka, 1980 :

As it can be observed, there are seven different isotopes which are linked through the 𝑛, 𝛾 reaction and the 𝛽 and 𝛼 decay. A generalization of this structure can be represented as: A ZX

→

,

→

11

A 4 A Z 2 X → ZX

rna

Where denotes all the succession of five different isotopes where the mass and atomic number increased in four unities, which can be built through the permutations of the following reactions: 𝑛, 𝛾 , 𝑛, 𝛾 , 𝑛, 𝛾 , 𝑛, 𝛾 , 𝛽 , 𝛽

Jou

An example of this type of pure cyclic chain will be used in the numerical analysis of Section 6. It is worth mentioning that the upper bound of 7 is only valid to the pure cyclic chains, because for the compound’s type this value is bigger. For example, the following structure contains eight heavy isotopes that commonly appear in decay networks in nuclear engineering: ,

U

⎯ ⎯⎯⎯ ,

,

U

⎯ ⎯⎯⎯ ,

U

,

⎯ ⎯⎯⎯ ,

,

U

,

⎯ ⎯⎯⎯

U

,

⎯ ⎯⎯⎯ ,

,

U

⎯ ⎯⎯⎯ ,

,

U

⎯ ⎯⎯⎯ ,

,

U

⎯ ⎯⎯⎯

U

,

Nevertheless, such structure cannot be considered as a pure-cyclic chain, because it was built as a union of eight different pure cyclic chain, as it can be noted in Figure 5, and therefore it is a compound cyclic chain. Even more, only the reactions that involve isotopes of uranium have been represented in the last sequence, but the fission reaction, the alpha and beta decay, and all the network related to plutonium and neptunium isotopes has been ignored. In other words, the full network is even more complex. 8  

lP repro of

Journal Pre-proof

Figure 5. Example of a compound cyclic chain related to isotopes of uranium. Such structure was built as a union of eight pure cyclic chains. As it was mentioned before, the present work only considers the pure cyclic chain, but as it will be discussed in Section 5.3, it is possible to extend the results to other types of compound cyclic chains. Therefore, in the present study, the limit value for 𝑛 will be 7. 3.6 Forward and Backward method.

rna

A remarkable disadvantage of the system in 6 is the lack of symmetry between the equations of the terms corresponding to the isotope 𝑘, with 2 𝑘 𝑛 1 and the isotope 𝑛. Due to this, it cannot be extended the application of equation 10 to the other isotopes of another index 𝑘. Nevertheless, through the balance equations it can be found more straightforward methods for the solution of 𝑋 . From the balance equation it can be written: 𝑑𝑋 𝑑𝑡

𝜆

𝑋

𝜆 𝑋

12

Jou

From which it is possible to obtain the following recursive formula: 1

𝑋

𝜆

𝑑𝑋 𝑑𝑡

13

𝜆 𝑋

Since this equation used the solution of the isotope 𝑘 to obtain the solution of the isotope 𝑘 1, it will be named the backward method equation. On the other hand, if the equation 12 is multiplied by the integral factor given by 𝑒 , the differential equation can be solved as: 𝑋

𝑒

𝜆

𝑋

𝑡 𝑒

9  

𝑑𝑡

𝑋 0

14

Journal Pre-proof

lP repro of

In this last method it is necessary to know the solution for the isotope 𝑋 to obtain the solution of 𝑋 , therefore it will be called the forward method equation. Both methods will be used in the following sections to extend the solution to the other isotopes in a cyclic chain. 4. A backward formulation of the Bateman equation.

It is possible to transform the recursive expression given in 13 to a general equation, which will be known in this work as the backward formulation of the Bateman equation. This generalization will be very useful to extend the pure cyclic chain solution. It is only necessary to perform multiple substitutions in the recursive formula. For the case 𝑘 2, equation 13 is equal to: 1

𝑋

𝑑𝑋 𝑑𝑡

𝜆

𝜆

𝑋

in 13 is replaced in the last expression, 𝑋

Now, if the term 𝑋 𝑋

𝜔 𝑘

2

𝑑 𝑋 𝑑𝑡

𝜆

𝑑 𝑋 𝑑𝑡

𝜆

𝜆

can be written as:

𝜆 𝑋

Where the function 𝜔 𝑘 is defined as:

1 𝜆

𝜆 ∙

𝜔 𝑘

Following a similar procedure, the equation for 𝑘 𝑋

𝜔 𝑘

𝜆

𝜆

𝜆

𝑑 𝑋 𝑑𝑡

𝜆

rna

𝜆

3

15

3 will be:

𝑑 𝑋 𝑑𝑡

𝜆

𝜆

𝜆

𝜆

𝑑 𝑋 𝑑𝑡

𝜆 𝜆

𝜆

𝑋

Finally, through induction it is possible to show Cruz-López, 2019 that for the case 𝑘 𝑋

𝑚

Jou

With:

𝜔 𝑘

𝑑 𝑑𝑡

𝒮 𝑖, 𝑘

𝒮 𝑖, 𝑘

𝜆

𝑋 ,

1

𝑘

𝑚

𝑛 16

…

𝜆

𝑑 𝑑𝑡

 

4.1 Derivatives in the Bateman backward equation.  In order to use equation 16 it is necessary to compute the high order derivatives given by ,1

𝑖

𝑛, 1

𝑘

𝑛 . It is possible to divide such task in two parts. First, it is possible to

find the derivative for the case 𝑘 other indexes using 13 .

𝑛. Then, starting from it, compute the derivatives for the

10  

Journal Pre-proof

4.2 Case 𝒌

𝒏

lP repro of

It is convenient to write equation 10 as: 𝑋

𝑋 0 𝜆

17

𝜆 𝑋 𝑋 0 𝜆 𝑋

Then, the following functions are defined: 𝑈 𝑡 For the case of 𝑎

𝜆 𝑋 𝑡 , 𝑉 𝑡

𝑋 0 , it is clear that 𝑋

𝑈 𝑡 𝑎 𝑈 𝑡

𝑉 𝑡 . Therefore, the problem of computing the

high order derivatives of 𝑋 is reduced to find a formula for the derivatives of 𝑉 𝑡 and for the derivatives of 𝑈 𝑡 . 4.2.1 Derivatives of 𝑽 𝒕 .

In the present work, it was found that the derivatives of 𝑉 𝑡 are related to integer’s partitions, which are all the possible combinations in which an integer 𝑛 can be represented as a sum of positive integers. In order to show that relationship, in the following discussion the fifth order derivative of 𝑉 𝑡 will be built. First, such derivatives will have as many terms as partitions of the number to which the derivative corresponds. For 𝑉 , it is necessary to know the partitions of 5, which are: 5, 4 1, 3 2, 3 1 1, 2 2 1, 2 1 1 1, and 1 1 1 1 1. Then, there are seven partitions, and therefore the derivative of 𝑉 will have seven terms. Each term will have a derivative of 𝑈 or a product of derivatives of 𝑈, whose order is related, again, to the partitions. For example, for the partition 3 1 1, the corresponding product of derivatives will be: 1

∙𝑈 ∙𝑈

1→𝑈

rna

3

𝑈

∙ 𝑈

Additionally, each term will be a fraction, whose denominator is given as a subtraction raised to a power. Therefore, the preliminary shape of 𝑉 will be: order 5

order 5

𝑉

𝑐 𝑎𝑈 𝑎 𝑈

order 3

order 1

order 4 order 1 order 5

order 3 order 2 order 5

𝑐 𝑎𝑈 𝑈′ 𝑎 𝑈

𝑐 𝑎𝑈 𝑈′′ 𝑎 𝑈

order 1 order 5

order 2

Jou

𝑐 𝑎𝑈 𝑈 𝑎 𝑈

order 2

order 1

order 1 order 1 order 5

order 2

18

order 1 order 5

𝑐 𝑎 𝑈 𝑈′ 𝑎 𝑈

order 1 order 1

𝑐 𝑎𝑈 𝑎 𝑈

order 1

order 1 order 1 order 5

𝑐 𝑎 𝑈 𝑎 𝑈

The following step is to find the coefficients 𝑐 , … , 𝑐 and 𝑘 , … , 𝑘 in equation 18 . The first set of coefficients can be computed through a sequence of numbers also related to partitions, which 11  

Journal Pre-proof

𝑧 , where 𝑧

lP repro of

is denoted by the OEIS’s code2 as A049019. Such sequence will be denoted as 𝑝 is a not negative integer. The first 29 terms of this sequence are the following: 1, 1, 2, 1, 6, 6, 1, 8, 6, 36, 24, 1, 10, 20, 60, 90, 240, 120,

1, 12, 30, 20, 90, 360, 90, 480, 1080, 1800, 720, … In order to use this sequence, it is only necessary to fragment it using the number 1 as a division point or delimiter: 1 , 1, 2 , 1, 6, 6 , 1, 8, 6, 36, 24 , 1, 10, 20, 60, 90, 240, 120 , 1, 12, 30, 20, 90, 360, 90, 480, 1080, 1800, 720 , …

Then, each resulted block will correspond to the set of coefficients of a given derivative of 𝑉. In the present example, for the derivative 𝑉 , the set 𝑐 , 𝑐 , … , 𝑐 corresponds to the set 1, 10, 20, 60, 90, 240, 120 , which is the block in the fifth position. In most programming languages is possible to carry out the fragmentation of the sequence 𝑝 𝑧 . It is only necessary to use a “split” function of strings, using as a character separator or delimiter the number 1. A procedure of fragmentation related to the format used in Python and MATLAB is described in Algorithm 1 on the Appendix A. On the other hand, the coefficients 𝑘 , … , 𝑘 in 18 also have a relationship with the partitions. Their value can be computed considering the number of elements that are involved in a particular partition, which will be denoted by the function 𝑛 . For the case of 𝑉 such function is given by:

And, therefore: 𝑘

𝑛 𝑃

1, 𝑛 4

1

𝑛 2

1

1

𝑛 3

2

2, 𝑛 3

1

4,

𝑛 1

1

1

1

1

𝑛 2

2

1

5

1

1

3

rna

𝑛 5

1, 1

𝑙

Number total of terms in a partition of the number 𝑛

Where 𝑃 is one element of the partition of the integer 𝑛. The general procedure used to build the derivative 𝑉 is described in Algorithm 2, which is also in the Appendix A. It is worth mentioning that in this procedure only the equation is built, in a symbolical sense, but the evaluation is not carried out.

Jou

4.2.2 Derivatives of 𝑼 𝒕 .

For the case of 𝑉 𝑡 , it is necessary to find high order derivatives of 𝑈 𝑡 . For the case of 𝑈′ the following equation is valid: 𝑈 𝑡

2

𝜆

𝑑 𝑋 𝑡 𝑑𝑡

𝜆 𝜆

𝑋

𝜆 𝑋

𝜆 𝜆

𝑋

𝜆 𝑈 𝑡

 

On-Line Encyclopedia of Integer Sequences, https://oeis.org/

12  

𝜆 𝜆

𝑋

𝜆 𝑋

Journal Pre-proof

𝑋 𝑡

The equality

𝑋

𝜆

𝜆 𝑋 was obtained from the balance equation given in

𝑡

𝑈

𝑋

𝜆

lP repro of

12 . The general formula, through multiple replacements, is equal to Cruz-López, 2019 : 1 𝐶 𝐴

,𝑘 𝑋

1 𝜆 𝑈 𝑡

19

𝜆 ,𝜆 ,…,𝜆 and the function 𝐶 𝐴 , 𝑘 is equal to the sum of the Where the set 𝐴 products related to the combinations with repetitions of the set 𝐴 . For example: 𝐶 𝐴 ,4 𝐶 𝐴 ,5

𝜆 𝜆

𝜆

𝜆 𝜆

𝜆

5. Root‐based solution.

𝜆

𝜆 𝜆

𝜆

𝜆 𝜆

𝜆

𝜆 𝜆

𝜆

𝜆 𝜆

𝜆

𝜆

𝜆

𝜆 𝜆

𝜆

𝜆

In addition to the solution developed in 10 , it is possible to develop an alternative method based on the roots of the following polynomial: 𝑠

𝜆

𝑠

𝜆

… 𝑠

𝜆

𝛼 𝑛

20

If such roots are computed, equation 7 can be expressed as partial fractions, and therefore it is only necessary to find the inverse Laplace transform to solve the problem. Under the assumption that all the roots are real, it is possible to rewrite the last polynomial as: 𝑠

𝛽 ,

𝛽 is a root of 20

rna

In order to use symmetry in the analytic solution, the following variables 𝛽 ∗ 𝛽 , with 1 𝑗 𝑛 will be used. Using a similar procedure, like the one described in 8 - 9 , it is possible to solve the system in 7 through these roots. Such solution will be known as the root-based equation, being equal to 𝑋 : 𝑋

𝑡

𝑋 0

𝜆

𝑒

1

∗

𝛽∗

21

𝛽∗

Jou

As it can be observed, this last equation has a similar shape to equation 4 . Unfortunately, like equation 10 , equation 21 cannot be applied to other isotopes in the cyclic chain due to the lack of symmetry of the system. Nevertheless, through the backward equation 13 , the solution for the isotope 1 𝑘 𝑛 will be equal to: 𝑋

𝑋 0

𝜆

𝑒

∗

𝜆

𝛽∗

1 𝛽

∗

𝛽∗

22

5.1 Forward method and the extension to more complex cyclic chains. Until this point, two methodologies for solving cyclic chains of the pure’s type have been developed. Such methodologies depend on the backward method described in Section 3.6, without which it will be not possible to extend their application to all the isotopes in the cyclic 13  

Journal Pre-proof

𝑋

𝑒

lP repro of

chain. In principle, the forward method can also be used to carry out this task, nevertheless, for the power series solution this procedure is very difficult, because it is necessary to solve integrals as the following: 𝜆

𝑋 𝑡 𝑒

𝑑𝑡

𝑋 0

23

𝑒

𝜆 𝑋 0

𝑋 𝑡 𝑒 𝜆 𝑋 𝑡 𝑋 0

𝑑𝑡

𝑋 0

This type of integrals will be studied as a part of a future research, but for the moment the forward method cannot be applied in a direct way to the power series solution. Furthermore, for the root-based solution, the forward method is very useful because it allows computing the error involved when the cyclic chain is approximated by the linear chain method, and also because through this procedure it is possible to extend the solution of a pure cyclic chain to more complex structures. 5.2 Ignoring cyclic chains based on a relative error.

rna

Isotalo and Aarnio carried out a comparison between three versions of the linear chain methods (Isotalo and Aarnio, 2011): the basic methodology that ignores all cyclic chains (which will be called as “the linear chain method”), the variation type where small modifications are introduced in the repeated lambda coefficients, and the generalized method that uses a solution of the system (2) that does not have the restriction about the lambda coefficients. They analyzed two study cases: the pressurized light water reactor (PWR) and the liquid metal fast breeder reactor (LMFBR). They found, among other findings, that the cyclic chains are not particularly important in terms of the relative error of the results. Regarding the running times, they also found that the slowest method is the generalized one, followed by the variation’s type and being the basic linear chain the fastest. Another important conclusion of the Isotalo and Aarnio’s work, is that the relative errors are concentrated around certain levels, because they are not independent. In other words, there is a set of components that dominates the error of most of the isotopes in a burnup problem. Clearly, these components are related to the isotopes that are repeated in a cyclic chain.

Jou

It would be useful to generalize the results to all nuclear reactor problems, in order to ignore the cyclic chains and to reduce the execution’s time of the burnup codes. Nevertheless, it is necessary to have a practical methodology who computes the involved error when the basic linear chain method is used, and it is in this point where the developed solutions of the present work become important. Instead of performing the whole calculations, as Isotalo and Aarnio carried out, it is possible to find the exact solution of a cyclic chain and, starting from it, to compute the relative error related to a given isotope. 5.2.1 Deduction of the error. In order to compute the error involved in such procedure, it is possible to use the forward method and the root-based solution. For example, in a pure cyclic chain, denoted by the sequence 𝑋 → 𝑋 → ⋯ → 𝑋 → 𝑋 → ⋯, the solution for 𝑋 can be computed as: 14  

Journal Pre-proof

𝑋

𝜆

Through equation (22):

𝑋

𝑑𝑡′

𝑡 𝑒

𝑑𝑡

𝑡 𝑒

𝑋

𝑋 0

24

lP repro of

𝑒

⎡ 𝜆 ⎢⎢ ⎢ ⎣

𝑋 0

∗

𝑒

1

1

𝛽∗

𝜆

𝛽∗

⎤ ⎥ 𝛽 ∗ ⎥⎥ ⎦

25

Defining the variable 𝛽 ∗ 𝜆 , the term 𝜆 𝛽 ∗ can be included inside the product under the index 𝑗, and the expression inside the brackets can be written as: 1

∗

𝑒

1

𝛽∗

𝛽∗

𝛽∗

𝛽∗

For the second term of the last equation, the following relationship is valid (Slodička and Balážová, 2008 : 1

1

𝛽∗𝑖

𝛽∗𝑗

𝜆

Using this, the term 𝑒

𝑋 0

𝑋

⎛ 𝜆 ⎜

1

𝑑𝑡 in (24) will be equal to:

𝑡 𝑒

𝑒

𝛽𝑛∗

𝛽∗𝑗

1

∗

𝛽

∗

1

𝛽

𝑒 𝛽∗𝑛 1

𝛽𝑗∗

∗

𝑛

𝑗 1

1

𝛽

∗

𝛽

∗

𝑛 1

𝜆1 𝑡

𝑒

26

⎠

rna

⎝ By notation, it is possible to show that:

⎞ ⎟

𝑗 1 𝑗 𝑛 1

1

𝛽

∗

𝛽

𝑒

∗

𝛽∗𝑛 1 𝑡

Then, using the last expression, the equation (26) can be rewritten as:

Jou

𝑋 0

𝑒

𝜆

1

∗

𝛽∗𝑖

𝛽∗𝑗

Finally, substituting the last equation in (24), the solution of 𝑋 becomes: 𝑋

𝑋 0 𝑒

𝑋 0

𝜆

Term equal to

1

∗

𝛽∗𝑗

𝛽∗𝑖

Additional term due to cyclic structure

15  

𝑒

27

Journal Pre-proof

𝑋

lP repro of

As it can be observed, through the forward method it is possible to express the root-based solution as the superposition of the linear chain solution, denoted by 𝑋 , plus an additional term, which is related to the cyclic structure. This last term is very useful to compute the error involved with any approximation of the cyclic chain. Particularly for the general case, from (27) it can be deduced that: 𝑋

𝛿 𝑡 ,

𝛿 𝑡

𝑋

𝑋

Where the error function can be found through: 𝛿 𝑡

𝑋 0

1

∗

𝑒

𝜆

𝛽∗𝑖

𝛽∗𝑗

, 𝛽𝑛

𝜆𝑘 ,

𝑘

1

𝑘

𝑛,

28

And the absolute relative error will be given by: 𝐸

|𝛿 𝑡 | 𝑋

29

5.2.2 Reduction of the number of linear chains.

In Section 3.2 it was discussed that a direct implication of building artificial linear chains is the generation of an infinite patron of repeated networks. As it was mentioned, it is possible to apply the linearization process in a cyclic chain to a certain degree. Wilson compared this procedure to the process of converting a connected graph into a n-array tree and called it the “tree straightening” procedure Wilson, 1999 . Unfortunately, this process generates a large number of linear chains, which is dependent of their extension i.e. the number of isotopes that contain , that increases in an exponential way. As an example, it will be considered the following compound cyclic chain: ,

⎯ ⎯⎯⎯

rna

X

,

,

X

⎯ ⎯⎯⎯

X

30

,

X and X have only the reactions that appears in 30 , Considering that the isotopes X, and that only the isotope X has an initial concentration different from zero, it is possible to build the following linear chains for the initial time step: ,

⎧ X ⎯ ⎪

,

X ⎯⎯⎯

,

X ⎯

Jou

⎨ ⎪ ⎩

,

X ⎯

,

X ⎯

,

X ⎯

,

X ⎯

,

,

X ⎯⎯⎯ ,

X ⎯⎯⎯

,

,

X ⎯⎯⎯ X ⎯⎯⎯

X ⎯⎯⎯

X

31

X

X

These three linear chains are part of those produced by applying the linearization process, with 7, 5 and 3 isotopes. In all the cases the final isotope is X. If such linear chains are ignored, the concentration of X will be underestimated to a certain degree. The equations developed before can help to compute this error. Now, it will be considered that X undergoes fission in 30 , and that the fission products originate a set of linear chains given by:

16  

Journal Pre-proof

𝑋 𝑋 … 𝑋

… … … …

𝑋 𝑋 … 𝑋

(32)

lP repro of

𝐴

𝑋 𝑋 … 𝑋

Where the elements 𝑋 represent the descendent 𝑗 of the fission product 𝑖. In (32) it has been assumed that the linear chains have the same length, but this assumption is not necessary. Then, the set 𝐴 must be inserted in the linear chains that were generated in (31). Therefore, the number of linear chains that need to be solved to compute the concentration of the fission products is equal to three times the cardinality of 𝐴. If this cyclic chain is ignored, the concentration of the fission products will be underestimated proportionality to the underestimation of the concentration of X. In other words, the relative error due to this procedure is not independent of the isotopes. In this example, the number of linear chains was reduced three times. Nevertheless, in a compound cyclic chain, like the one showed in Figure 5, the number of linear chains decreased considerably. Therefore, if the cyclic chain structures are ignored the execution time of a burnup code will be reduced. 5.2.3 Considerations about the accuracy and running time in the linear chain method. A recurrent issue with the Bateman equations is related to precision and round-off difficulties when very small lambda coefficients are used as well as with the time’s execution. These topics were identified since the beginning of the development of the linear chain method in the 50’s and 60’s. The first author that addressed them in a very detailed work was Vondy Vondy, 1962 , who carried out a complete analysis about the way in which the modified Bateman solution can be implemented in digital machine calculations, particularly in the IBM-7090 computer on Fortran programming language.

rna

Such effort continues since then, being notorious the works of (Tobias, 1978), Tasaka (Tasaka 1980), Henderson (Henderson, 1986), Vukadin (Vukadin, 1991) among others, who focused mainly in the precision and round-off errors, proposing several ways in which the Bateman solution can be expressed in order to decrease the mentioned difficulties. Regarding the time execution, the first important work was published by Cetnar (2006) who developed a full theory related to the cut-off termination criteria, named Transmutation Trajectory Analysis (TTA), which allows reducing the number and the length of the linear chains that were used in burnup problems and therefore increasing the speed of the algorithm.

Jou

After these works, a new generation of linear chain codes appearing, who proposed new techniques or improved the existing ones, being worth mentioning three among them: Harr (Harr, 2007) developed an advanced method based in an integral formulation of the Bateman equation, that allows dealing with cases where the lambda coefficients differ in several orders of magnitude. Huang and his colleagues (Huang et. al ,2015) carried out two important improvements to the TTA, that are related to the termination criteria as well as to a recursive formula. More recently Chaitanya Tadepalli et. al. (Chaitanya Tadepalli et. al., 2017) developed an advanced activation code named ACTYS, where they proposed a new methodology for the chain termination, based in recurrence relationships of a passage function, as well as a technique to overcome the difficulties when there are linear chains with long-lived isotopes.

17  

Journal Pre-proof

lP repro of

This last method is based on separating the long-lived isotopes in the Bateman solution and carried out a Taylor series expansion of the exponential functions. It is important to note that the present work proposed an improvement in execution’s time terms from a new approach: analyzing the cyclic chains, being complementary to the listed works because in them there is not a depth discussion about this type of decay and transmutation networks. Additionally, as it was discussed in Section 3.1 and 3.2, the development of solution of cyclic chains can be used as reference solutions, because they do not involve assumptions related to building artificial linear chains as the standard methodology does. Therefore, it is important to make this distinction about the approach that is used in the present work in comparison with the existing ones. 5.3 Extension of the root‐based solution.

Besides to compute the error involved with the solution of cyclic chains with the linear chain method, the forward method can be used to extend the solution for more complex cyclic chains. Particularly it is possible to include structures with branches like the one showed in Figure 6, where are “branches” that emerge from a pure-cyclic chain. With the purpose to show how these structures can be solved, the following development will correspond to the solution of the branch denoted by the sequence of the 𝑌 isotopes that appears in the figure. The balance equations for the isotopes of the pure cyclic chain, denoted by: 𝑋 → 𝑋 → ⋯ → 𝑋 are the same that those described in system (6). The branch related with the isotopes 𝑌 can be described as the linear chain: 𝑌

𝑌

𝑌

…

Whose balance equation are given by:

𝜆

𝑌 𝜆𝑋

rna

𝑑 𝑌 𝑑𝑡

𝜆 𝑌 𝜆 𝑌

2 𝑝

𝑝 , 1

For the solution of 𝑌 the forward method can be used as follows: 𝑌

⎛ ⎜

𝑒

∗

𝜆

𝜆

𝑋

𝑌

𝑒 𝛽∗

𝑌

𝑋 0 𝜆

𝜆

0 33

𝛽∗

𝑛

18

𝑘, after which such term is

𝜆

34 1

𝛽

∗

⎞ 𝛽 ⎟ ∗

⎠

⎝

 

𝑑𝑡

𝑡 𝑒

𝑡 though equation (22) with 𝑗

Jou

It is possible to build 𝑋 replaced in (33):

𝑒

𝑌 0 𝑒

Journal Pre-proof

lP repro of



Figure 6. Scheme of cyclic chain where some of its elements have branches, from which emerge linear chains. Such scheme can be solved through the forward method. For the general case, the solution for the rest of isotopes 𝑌 in the branch was found in the same way see Appendix B :

𝑌

𝑋 0 𝜆

𝜆

⎛ 𝜆 ⎜

𝐹 𝛽 ,𝑌 ,…,𝑌

𝜆

1

𝛽∗

𝛽∗

⎝

Where

6. Numerical Analysis.

⎞ 𝛽∗ ⎟

35

⎠

𝐹 𝑎 ,𝑎 ,…,𝑎

𝑒

1

𝑎

𝑎



rna

A numerical analysis between the developed solutions in this work and the linear chain method has been carried out, with the purpose to show how the developed equations can be used. It is important to observe that such developed solutions are not equivalents, because the first one was developed using certain assumptions related to the lambda coefficients, while the rootbased solution does not use any, except the related to the numerical computation of the roots. In such sense, this last solution will be considered as the reference solution. The present analysis was based in the work of Blaauw Blaauw, 1993 .

Jou

The cyclic chain that will be used for the comparison is defined by the sequence of isotopes 235 U → 236U → 237U → 237Np → 238Np → 238Pu → 239Pu → 235U, as it is showed in left side of Figure 7. The model will be fed with data obtained from the simulation of a typical light water reactor unit cell in infinite medium, in a thermal neutron spectrum. Calculations were done with the well-known Monte Carlo code SERPENT Leppanenet al., 2015 . The description and main features of the unit cell are shown in the right side of Figure 5 and in Table I, Table II and Table III. For the decay constants, the ENDF/VII.1 Library was used, and the numerical operations were carried in the software Mathematica 11.3.0, using precision 40.

19  

lP repro of

Journal Pre-proof

Figure 7 Left side. Pure cyclic chain given by the sequence of isotopes that begins with 235U and ends with 239Pu. Right side. Geometry of the unit cell in an infinite medium used in the comparison scheme3. As it will be shown in the following sections, through one of the developed solutions it is possible to conclude if the linear chain method can approximate the solution of a cyclic chain, and using the other developed solution it is possible to compute the error involved in such approximation. 6.1 Generalization of the equations.

rna

Until now all the developments were made considering only the decay process. Nevertheless, in order to include nuclear reactions, it is necessary to define the following terms Isotalo, 2013 : 𝜆eff

𝜆

𝜙 𝜎

,

36

Where 𝜙 is the neutron flux for the group energy 𝑔. 𝜎 , is the microscopic cross section of the isotope 𝑖, for the reaction type 𝑚 and the energy-group 𝑔. Besides these elements, it is necessary to define the effective branch, as the following ratio: 𝑏,

Jou

𝑏 eff ,

𝜆

𝜆eff

𝜔

,

𝜔

𝛾,

,

𝜙 𝜎

37 ,

Where the constant 𝑏𝑛,𝑛 1 is the standard branching ratio, i.e., the fraction of decays of the 𝑔 isotope 𝑛 that produces the isotope 𝑛 1. 𝛾𝑖,𝑖 1,𝑚 is the yield of the isotope in the position 𝑖, related to the production of the isotope in the position 𝑖 1, for the reaction type 𝑚 and with the energy group 𝑔. The constants 𝜆eff will be called “the effective removal coefficients”.

3

 

→ In such scheme the symbols 𝐽 and 𝐽 denote the magnitude of the neutron current density vector 𝐽 .

20  

Journal Pre-proof

Table I. Main features of the unit cell used in the comparison scheme



Isotope

Inner radius cm Outer radius cm Square’s apothem cm Coolant Cladding neutrons Neutron Flux cm ∙s Volume m



         

U U 236 U 238 U 235

5.33267x10

Table II. Neutronic parameters of the isotopes belonging to the unit cell used in the comparison scheme 236 237 237 235 Parameter U Np U U 1.73591 𝜎 𝑏 5.51954x10 3.0372x10 5.05637𝑥10 𝜎 , 𝑏 8.65019 1.17859x10 6.29918𝑥10 3.75076𝑥10 𝜎 , 𝑏 3.52596x10 2.18561x10 1.12769𝑥10 9.28667𝑥10 1.503275x10 8.27197x10 4.72784𝑥10 1.37712𝑥10 𝜎 ∙𝜙 𝑠 𝜎 , ∙𝜙 𝑠 3.20995026x10 2.35592x10 1.71561𝑥10 1.02154𝑥10 𝜎 , ∙𝜙 𝑠 9.60313x10 5.95262x10 3.07132𝑥10 2.52927x10 3.12298x10 9.38495x10 1.18852x10 1.02516𝑥10 𝜆 𝑠 1.82437x10 2.43924x10 1.20615x10 1.03533𝑥10 𝜆eff eff 1.75949x10 9.65843𝑥10 9.85381x10 9.86673x10 𝑏,

rna

 

234

0.412 0.475 0.665 Water Zircalloy 2.72x10

Table III. Neutronic parameters of the isotopes belonging to the unit cell used in the comparison scheme. 238 238 239 Parameter Pu Pu Np 2.88213 𝜎 𝑏 1.94989x10 1.41846x10 𝜎 , 𝑏 1.97274x10 4.11423x10 7.96234x10 𝜎 , 𝑏 4.46679x10 1.08646x10 1.39825x10 5.31063x10 7.84963x10 3.86325x10 𝜎 ∙𝜙 𝑠 𝜎 , ∙𝜙 𝑠 5.37286x10 1.12053𝑥10 2.16858x10 𝜎 , ∙𝜙 𝑠 1.21655x10 2.95903x10 3.80821x10 3.78958x10 2.50622x10 9.11636x10 𝜆 𝑠 3.84806x10 1.22412x10 6.03196x10 𝜆eff 9.84802x10 9.15377x10 1.51134x10 𝑏 eff,

Jou



Using these terms, it is possible to generalize the equation 4 : 21  

Initial concentration atoms/ b ∙ 𝑐𝑚 6.15138x10 6.89185x10 3.16249x10 2.17092x10

lP repro of

Parameter





Journal Pre-proof

𝑏 eff,

𝑋 0

𝜆eff

𝑒

1

eff

𝜆eff

𝜆eff

38

lP repro of

𝑋 𝑡

Due to these modifications, from now on, equation 38 will be called the Modified Bateman Equation. For equation 10 the equivalent expression will be: 𝑋

𝑋 0

𝑋

39 𝑋 0 𝑏 eff, 𝜆eff 𝑋 Where the super script 𝑃 will denote “the power series solution”. For the equation 22 , the corresponding generalized equation is: 𝑋

𝑋 0

𝑏 eff,

And finally, for the function 𝛼:

𝜆eff

𝑒

∗

𝜆eff

1

𝛽∗

𝛽∗

𝛽∗

40

𝜆eff 𝑏 eff ,

𝛼 𝑛

6.2 Comparison with the Power series solution.

In order to compute the error when the linear chain solution is used in the structure given in the left side of Figure 7, it is necessary to analyze the term 𝑏 eff, 𝜆eff of equation 39 . From the definition provided in the past section it is known that: 𝑏 eff, 𝑏 239Pu, 235U 𝜆eff 239 Pu

𝜆 239Pu

𝜆 239Pu 𝜆eff 239 Pu

9.11636x10

rna

Then, 𝑏 eff, 𝜆eff

𝑏 239Pu, 235U

1.51134x10

.

5

On the other hand, Figure 8 contains

the graph of 𝑋239Pu . Through this graph, it is possible to conclude that: 10

atoms b ∙ cm

𝑋239Pu

10

atoms b ∙ cm

From the last equations:

atoms b ∙ cm

𝑏 239Pu, 235U 𝜆eff 239 𝑋239 Pu Pu

Jou

10

10

atoms b ∙ cm

Since the denominator of 39 contains a subtraction between 𝑋 0 and 𝑏 239Pu, 235U 𝜆eff 239 𝑋239 , Pu Pu it is possible to analyze through such operation if the solution of the cyclic chain can be approximated by the linear chain solution. In this particular case, it is clear that this subtraction is reduced to 𝑋 0 , and therefore the equation 38 is equal to the linear chain solution given in 37 . In other words, using this brief analysis it is possible to determine if the linear chain method is a very good approximation to the cyclic chain solution.

22  

Journal Pre-proof

Table IV. Roots and their comparison with the effective removal coefficients. 𝛽𝑠 -1.82437 x10 -2.439230 x10 -1.20616x10 -1.035350x10 -3.848060x10 -1.22411x10 -6.03196 x10

𝐸 %

𝛽∗ 𝑠 1.82437 x10 2.439230 x10 1.20616x10

𝜆eff 𝑠 1.82436x10 2.43923947 x10 1.2061559x10

2.15839 x10 3.88253 x10 3.34371x10

1.035350x10

1.035336x10

1.322642x10

3.848060x10 1.22411x10 6.03196 x10

3.848058x10

4.115806x10

1.224119 x10 6.03196184 x10

7.927498 x10 3.04933 x10

lP repro of

Isotope 235 U 236 U 237 U 237 Np 238 Np 238 Pu 239 Pu

As it can be observed, the main utility of equation 39 is to conclude if the linear chain method can be used to model a cyclic chain. Such conclusion can be obtained through the analysis about the subtraction that is the denominator of the equation. Nevertheless, from equation 39 it is not possible to compute the error involved in such approximation. From this last task it is more convenient to use equation 22 and 29 , which will be carried out in the next section. 6.3 Analysis using the root-based solution.

In order to build the root-based solution it is necessary to compute the roots of the following polynomial: 𝑠

𝜆eff 235 U

𝑠

𝜆eff 236 U

𝑠

𝜆eff … 𝑠 237 U

𝜆eff 239 Pu

𝛼 7

41

Such computations were carried out with the software Mathematica 11.3.0, and they are listed in Table IV, with their respective values of 𝛽 ∗ , as well as the absolute relative error given by their comparison with the effective removal coefficients, in other words: absolute relative error of 𝛽 ∗

rna

𝐸

100

|𝛽 ∗

𝜆eff |

eff

𝜆

%

42

Jou

As it can be observed from Table IV, the negative of the roots of equation (40) are closer to the effective removal coefficients. This is an important fact, because it implies that the cyclic structure modifies all the removal coefficients of a standard linear chain. Considering the data from Table II, Table III and Table IV, the 239Pu concentration was computed using equations (39) and (40) for a range of 5 days to 1000 days. Figure 8 shows both graphs, which overlap. A precision of 50 digits was used, because several of the isotopes of the cyclic chain have a very small effective removal coefficients, which can produce that the equations yield negative values of the concentration. Figure 9 shows the absolute relative error between the Modified Bateman equation and the root-based equation, which was computed using the first one as the reference value: 𝐸 %

100

𝑋

23  

𝑋 |

|𝑋

%

43

Journal Pre-proof

Pu-239

lP repro of

1.00E-06

Concentration [atoms/(barn-cm)]

1.00E-07 1.00E-08 1.00E-09 1.00E-10 1.00E-11 1.00E-12 1.00E-13 1.00E-14 1.00E-15 1.00E-16 1.00E-17 0

200

400

600

800

1000

Time (Days)

Root-based Equation

Modified Bateman Equation

Figure 8. Concentration for 239Pu compute with the Modified Bateman Equation, and the Root Based Equation. The two curves overlap. Absolute relative error (%)

rna

0.1

0.01

Jou



Absolute relative error (%)

1

0.001

0



100

200

300

400

500

600

700

800

900

1000

Time (Days)

Figure 9. Absolute relative error between the Modified Bateman Equation and the Root-based equation for the concentration of 239Pu. 24  

Journal Pre-proof

lP repro of

Except for the first point, such value is close to 1%, and for the rest of data show an increment that starts in an order of 10 and then it grows slowly to a value of the order of 10 . The source of the error in the first two steps of Figure 9 can be traced to the precision used. In fact, if such numerical precision is increased until 150 digits, it is possible to show that the error in the first two steps is of the order of 1𝑥10 %. As it can be observed from this analysis, the root-based equation can be useful to compute the actual error when the linear chain method is used to approximate the cyclic chains. A detailed graph of the dependence on the numerical precision of error is showed in Appendix C. 6.4 Results for U‐235.

It was found that the results for the concentration of the isotopes 238Pu, 238Np, 237Np, 237U and U are very similar to the results of 239Pu. In fact, for all these isotopes the graphs overlap, and the error has the same behavior. Therefore, it was considered convenient to discuss only the 235U of the cyclic chain, because the forward method can be used to solve it. Using equation (27), the forward Bateman solution for the 235U is given by 236

8

𝑋235U

𝑋235U

𝑋 235U 0

𝑏 eff,

eff

𝜆

𝑒

1

∗

𝛽∗𝑗

𝛽𝑖∗

44

From the last equation it is possible to identify the contribution that is associated with the structure of the cyclic chain, and particularly with the reaction 239Pu → 235U. Figure 10 shows the graph of equation (44) and for the modified Bateman solution for this isotope, as in the past case the curves overlap. The absolute relative error is shown in Figure 11; again, a precision of 50 digits was used to compute this error. The contribution due to the cyclic chain in this case is practically neglected, in fact, for a range of 1 Day 𝑡 1000 Days:

rna

atoms b ∙ cm

1.64547x10

8

𝑏 eff,

𝑋235U 0

4.0733x10

12

eff

𝜆

𝑒

1

∗

𝛽∗𝑗

𝛽∗𝑖

45

atoms b ∙ cm

Jou

As it was mentioned before, this contribution is due to the alpha decay of 239Pu to 235U, which involves a half-life of 24,110 years and, therefore, it has a very small value. It is important to note that (45) is an upper bound of the error for approximating a cyclic chain through the linear chain method for the 235U concentration. This outcome, along with the analysis carried out in Section 6.2, allows to compute the error of using the linear chain method. Since this error is very small, it is possible to omit the cyclic chain given in Figure 7 and reducing the number of linear chains that are generated with the linearization process. Therefore, the solutions that were developed in this work provide a complete methodology to analysis the pure cyclic chains. 7. Conclusions

In the present wok two alternative solutions to the cyclic chains were developed. Through them it is possible to answer: 1) In which cases the linear chain method is a good approximation to the pure cyclic chains? and 2) What is the error involved in such approximation? 25  

Journal Pre-proof

U-235

lP repro of

Concentration [atoms/(barn-cm)]

0.0008 0.0007 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0

200

400

600

800

1000

Time (Days)

Root-based equation

Modified Bateman Equation

Figure 10. Concentration for 235U compute with the Modified Bateman Equation, and the Root Based Equation. The two curves overlap. Absolute relative error (%)

1.00E-04



rna



1.00E-08 1.00E-10 1.00E-12 1.00E-14 1.00E-16 1.00E-18 1.00E-20

Jou



Absolute relative error (%)

1.00E-06



1.00E-22 1.00E-24

0

200

400

600

800

1000

Time (Days)

Figure 11. Absolute relative error between the Modified Bateman Equation and the Root-Based Equation, for the isotope U-235. 26  

Journal Pre-proof

lP repro of

The developed solutions have several advantages over the linear chain method. Firstly, using them it is possible to reduce the computational time because it is possible to conclude if a cyclic chain can be ignored, or if it is necessary to linearize it, reducing the number of equations in a considerably way. Secondly, the root-based solution developed in this work allows finding more accurate results, being adequate to find reference solutions, because its formulation does not involve approximations as in the case of the linear chain methodology. Using typical data of a light water reactor, it was possible to analyze a pure cyclic chain of heavy isotopes with both the developed solutions and the linear chain method. As it was expected in such example, the standard Bateman equation was a very good approximation to the cyclic chain, but this time it was possible to obtain as additional conclusion: the error involved when the linear chain is used. This last represents an important result over the standard analysis reported in literature. If the error involved in such approximation is very small, it is possible to use a cyclic chain solution instead a set of linear chains or even ignoring the cyclic chain structure, and therefore reducing the computational time. Even when the solutions proposed in the present work only solve pure cyclic chains, such study opens the possibility to extend the analysis to more complex cyclic chains through the backward and forward methods, as it was showed for the case of pure cyclic chains with branches emerging of them. Since the proposed solutions require a considerable number of algebraic operations as well as the calculation of high order derivatives, two computational algorithms were developed that reduce and simplify such operations in a considerable way. Finally, for the above discussion it is possible to conclude that the study of cyclic chains represents a potential field where some aspects related to the linear chain, from the algorithmic to the mathematical level, can be improved.

rna

Acknowledgments

Jou

Thanks to the Consejo Nacional de Ciencia y Tecnología for the scholarship provided to C.A. Cruz for the preparation of the present work, which is part of his doctoral research. To UNAM for its support through the project PAPIIT-IN115517. Also, the authors acknowledge the funding of the strategic Project No. 212602, AZTLAN Platform: ‘‘Desarrollo de una plataforma mexicana para el análisis y diseño de reactores nucleares’’, of the Sectorial Sustainability Energy Fund CONACYT — SENER.

27  

Journal Pre-proof

References.

lP repro of

Bateman, H., 1910. Solution of a System of Differential Equations Occurring in the Theory of Radioactive Transformations. Proceedings of the Cambridge Philosophical Society, 15, 423– 427. Blaauw, M., 1993. A Versatile Computer Algorithm for Linear First-order Equations Describing Compartmental Models with Backward Branching. Applied Radiation and Isotopes, 44 (9), 1225-1229. Cetnar, J., 2006. General Solution of Bateman equations for Nuclear Transmutations. Annals of Nuclear Energy 33, 640–645. Chaitanya Tadepalli, S., Kanth, P., Indauliya, G., Saikia, I., Deshpande, S. P., 2017. Development and Validation of ACTYS, an Activation Analysis Code. Annals of Nuclear Energy, Vol. 107, pp. 71-81. Cruz-López, C., 2019. Contributions to the Solution of the Bateman Equations for the Isotopes Transmutation in Fission Nuclear Reactors. Doctoral Dissertation in Engineering. National Autonomous University of Mexico. Mexico. Cruz-López, C., François, J., 2018. An alternative algorithm for the linearization process for transmutation and decay networks. Computer Physics Communications, 231, 122-139. Dreher, R., 2013. Modified Bateman Solution for Identical Eigenvalues. Annals of Nuclear Energy 53, 427-438. England, T. R., 1962. Time-Dependent Fission-product Thermal and Resonance Absorption Cross Section. In AEC Research and Development Report. Bettis Atomic Power Laboratory, Pittsburgh. U.S. Atomic Energy Commission. Westinghouse Electric Corporation.

rna

Harr, L. J., 2007. Precise Calculation of Complex Radioactive Decay Chains. Master Dissertation. Department of the Air Force. Air University. Wright-Patterson Air Force Base, Ohio. Henderson, D.L., 1986. DKR-1100: A Univac 1100 version of DKR Radioactivity code. University of Wisconsin, Fusion Technology Institute, UWFDM-671. Huang, K., Wu, H., Cao, L., Li, Y., Shen, W. 2015. Improvements to the Transmutation Trajectory Analysis of Depletion Evaluation. Annals of Nuclear Energy, 87, 637-647.

Jou

Isotalo, A., 2013. Computational Methods for Burnup Calculation with Monte Carlo Neutronics. Doctoral Dissertation. Aalto University Publication Series. Helsinki, Finland. Isotalo, A., Aarnio, P. A. 2011. Substep Methods for Burnup Calculations with Bateman Solutions. Annals of Nuclear Energy, 38, 2509-2514. Leppänen, J., Pusa, M., Viitanen, T., Valtavirta, V., Kaltiaisenaho, T., 2015. The Serpent Monte Carlo Code: Status, Development and Applications in 2013. Annals of Nuclear Energy, 82, 142150. Meyer, S., Schweildler, E., 1927. Radioaktivität. Springer Fachmedein Wiesbaden GmbH, 56-57.

28  

Journal Pre-proof

lP repro of

Slodička, M., Balážová, A., 2008. Singular Value Decomposition Method for Multi-species Firstorder Reactive Transport with Identical Decay Rates. 73 2 , 161-172. Tasaka, K., 1980. DCHAIN 2: A Computer Code for Calculation of Transmutation of Nuclides. Japan Atomic Energy Research Institute. Tokyo, Japan. Tobias A. 1978. FISP5-an Extended and Improved Version of the Fission Product Inventory code FISP, CEGB Report RD/B/N4303. Vondy, D. R., 1962. Development of a General Method of Explicit Solution to the Nuclide Chain Equation for Digital Machine Calculations. Oak Ridge National Laboratory. Vukadin, Z. 1991. Recurrence Formulas for Evaluating Expansion Series of Depletion Functions. Kerntechnik, 56, 395-397.

Jou

rna

Wilson, P. P. H., 1999. ALARA: Analytic and Laplacian Adaptive Radioactivity Analysis. Doctoral dissertation. Fusion Technology Institute. University of Wisconsin. Madison Wisconsin. USA.

29  

Journal Pre-proof











Appendix A

lP repro of

Algorithm 1 Fragmentation and search procedure Input: Order of the derivative of 𝑉 Output: Coefficients of the derivative of 𝑉

STEP 1: Select a subsequence of A049019 that contains its first element, and whose value is 𝑛. Store it as a string variable 𝑥: A049019 :𝑛 → 𝑥

STEP 2: Split the string variable 𝑥 using the number 1 as a separator or delimiter , and store the generated list in the variable 𝑦: 𝑠𝑝𝑙𝑖𝑡 𝑥, 1 → 𝑦

IF length 𝑦

order of derivative :

IF the split function of STEP 2 removed the separator: █ For

𝑖

1, 2, … , length of 𝑦: the element 1 at the beginning of the stored string in 𝑦 𝑖

█ Add

█ Find the element in 𝑦 whose position is equal to the order of the derivative of 𝑉, and store it in the variable 𝑧:

𝑦 order of the derivative of V → 𝑧

ELSE:

a larger value of 𝑛 and repeat the algorithm.

rna

█ Choose

Algorithm 2 Equation for the derivative of order 𝑖 of the function 𝑉 Input: 𝑖 ∈ ℕ Output: Equation for the derivative of order 𝑖 of the function 𝑉

Jou

STEP 1: Compute all the partitions of 𝑖, and list them in a descending order. STEP 2: Compute the set of coefficients of 𝑉 using Algorithm 1. STEP 3: For each partition: i

ii

iii iv

Determine how many numbers belong to the partition and call this quantity as 𝑛 . For each number 𝑠 that belongs to the partition: a Write the derivative of 𝑈, whose order is equal to 𝑠. b Store the term built in a . Multiply all the derivatives that were built in ii . Compute the value 𝑛 using 𝑛 . 30

 

Journal Pre-proof

Write the division of the product in iii by the expression 𝑎 Store the final term in the set 𝑇.

𝑈

.

lP repro of

v vi

STEP 4: Write an equation given by the sum of all the terms contained in the set 𝑇. STEP 5: Multiply each term of the equation built in STEP 3 by the corresponding coefficient computed in STEP 2.

Appendix B.

Proof of equation (34)

eff 𝑋 0 𝑏 eff , 𝜆

𝑌

⎛ ⎜

𝑏 eff,

𝜆eff

𝐹 𝛽 , 𝜆eff , … , 𝜆eff

⎝

𝜆eff

𝑏 eff,

𝜆eff

1

𝛽∗

𝛽∗

⎞ 𝛽∗ ⎟ ⎠

Proof

We will use induction over 𝑝. Considering that

𝐹 𝑎 ,𝑎 ,…,𝑎

𝑎

𝑎

1 was obtained in (33), being equal to

rna

The case 𝑝

1

𝑒

eff 𝑋 0 𝑏 eff , 𝜆

𝑌

⎛ ⎜

𝑒

∗

eff

𝜆

𝜆eff

∗

𝛽∗

⎝

Jou

𝜆eff

eff

𝑒

𝛽

𝑏 eff,

1 𝛽

∗

⎞ 𝛽 ⎟ ∗

⎠

It is clear that:

𝑒

∗

𝜆eff

𝑒

eff

𝐹 𝛽∗𝑖 , 𝑌1

𝛽∗

Therefore, case 𝑝 1 confirms equation (33). The next step is assuming that this equation is valid for 𝑝 𝑢, then it must be also valid for 𝑝 𝑢 1. Using the forward method, we have that:

31  

Journal Pre-proof

eff

𝑏 eff,

𝜆eff

𝑌 𝑡 𝑒

eff

𝑑𝑡′

lP repro of

𝑒

𝑌

Using the induction assumption, we have: Λ𝑒

𝑌

eff

With

𝛽∗, 𝑎

𝜆eff , 𝑎

𝐹 𝛽 , 𝜆eff , … , 𝜆eff 𝑒 1 𝑎

𝛽∗

eff 𝑋 0 𝑏 eff , 𝜆

Λ If we consider that 𝑎

1

𝛽∗

𝜆eff

𝑎

𝑒

𝑏 eff,

𝑑𝑡′

𝑏 eff,

𝜆eff , 𝑎

𝑎

𝑑𝑡′

𝑒

𝑒

𝑎

1

𝑎

𝑎

𝑒

𝑎

1

𝑎

𝑎

𝑎

𝑎

𝑎

1 𝑎

𝑎

𝑎

Considering that (Slodička and Balážová, 2008 :

rna

1

𝛽∗𝑗

1

𝛽∗𝑖

𝛽∗𝑗

𝛽∗𝑛

1

It is possible to write:

Jou

𝐹 𝛽 , 𝜆eff , … , 𝜆eff 𝑒

eff

𝑒

𝑑𝑡′

𝑎

𝑎

𝑒 𝑎

𝑎

Therefore

𝑌

Λ

𝜆eff

𝛽∗

1 𝛽

∗

32  

𝛽

∗

𝑒

eff

𝑑𝑡′

𝜆eff , then:

1

𝑑𝑡′

𝑒

eff

𝜆eff

𝐹 𝑎 ,𝑎 ,…,𝑎

1

𝑎

𝜆eff

𝜆eff , … , 𝑎

eff

𝐹 𝛽 , 𝜆eff , … , 𝜆eff 𝑒

𝛽∗

𝑒

Journal Pre-proof

𝜆eff

𝐹 𝑎 ,𝑎 ,…,𝑎

𝛽∗

1 𝛽∗

𝛽∗

lP repro of

Λ

Which is that we wanted to prove.

Appendix C

 

 

Absolute relative error

1000

Absolute relative error (%)

100 10 1 -10 0.1

50

0.01 0.001 0.0001

P30

0.00001

170

P50

rna

0.000001 0.0000001

110

230

P100

290

P150

Time (Days)

Absolute relative error error between the Modified Bateman Equation and the Root-Based Equation, for the isotope 239Pu, as a function of time and with different values of precision.

Jou

 

33  

Journal Pre-proof Declaration of Interest Statement

Declaration of interests

lP repro of

Manuscript: Two alternative approaches to the solution of cyclic chains in transmutation and decay problems. Authors: Carlos-Antonio Cruz-López, Juan-Luis François

Jou

rna

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.