A simple formula for local burnup based on constant relative reaction rate per nuclei

A simple and analytical formula is suggested to solve the problems of the local burnup and the isotope distributions. The present method considers two extreme conditions of neutrons penetrating the fuel rod. Based on these considerations, the formula is obtained to calculate the reaction rates of $^{235}$U, $^{238}$U, and $^{239}$Pu and straightforward the local burnup and the isotope distributions. Starting from an initial burnup level, the parameters of the formula are fitted to the reaction rates given by a Monte Carlo (MC) calculation. Then the present formula independently gives very similar results as the MC calculation from the starting to high burnup level, but takes just a few minutes. The relative reaction rates are found to be almost independent on the radius (except $(n,\gamma)$ of $^{238}$U) and the burnup, providing a solid background for the present formula. A more realistic examination is also performed when the fuel rods locate in an assembly. A combination of the present formula and the MC calculation is expected to have a nice balance on the accuracy and the cost on time.


Introduction
To increase the efficiency of the fuel, one possible way is to increase the burnup of the fuel before discharge. The local burnup on the edge of the UO 2 fuel rod is much higher than the average burnup. Thus it is of great importance to investigate the properties of the rim of the fuel rod when considering the increment of the average burnup. Many investigations show that the mechanical structure close to the surface is rather different from that at the center of the fuel rod. In the high burnup range, a microstructure change was found on the rim of the fuel rod through transmission electron microscopy [1]. One explanation of the formation of the high burnup structure supposed that the bubbles in high burnup region are nucleated and stabilized by fission fragments, which depends on the fission rate [2]. The small pores at the high burnup region are calculated to be highly overpressurized [3]. Recently, the UO 2 fuel in both light and heavy water reactor are investigated to understand the high burnup structure of the fuel [4,5].
In a thermal reactor, the neutrons generated from fission need to be slowed down in the moderator to be able to induce next fission. When the low energy neutrons go from the moderator to the fuel rod, they are firstly absorbed by the fuel close to the surface. In general, the reaction rate is higher in the rim when induced by slowed neutrons, such as the (n, f ) reaction of 235 U and 239 Pu induced mainly by the thermal neutron. In contrast, the (n,γ) reaction of 238 U is mainly induced by the resonance neutrons. The local burnup phenomena is mainly caused by the (n, γ) reaction of 238 U, which has large cross section for resonance neutrons, especially some high peaks of cross section at certain energies. Such reaction produces much more 239 Pu near the surface of the fuel rod, thus giving origin to the increased local burnup. As a consequence of the different reaction rates along the radial direction, the U and Pu isotopes also have different radial distributions. In a fast reactor, the Pu isotopes are redistributed along the radial direction [6]. A large local burnup is not expected because the neutrons do not need to be slowed down in the moderator.
In nuclear science, theoretical models are very important and powerful to solve various of problems. Our previous works have shown that the mass, levels, electromagnetic moments and transitions, and other properties of the nuclei can be described with high precision in the frame of the nuclear shell model even for very unstable nuclei, such as 22 C and 21 Al [7]. In nuclear reactor, the simulation is even more important because the experimental data are relatively difficult to be obtained, sometimes with high risk. For the present problem, various models can be used to calculate the local burn up and the isotope distributions along the radial direction. The TRANSURANUS model is very successful in the description of the local burnup and related properties in various reactors including light and heavy water reactors [8][9][10]. In TRANSURANUS model, a fixed neutron flux is assumed except for the absorption reaction of 238 U and 240 Pu [10], while the RAPID model considers detailed properties of the neutron flux at each burnup [11]. Recently, an empirical formulation is suggested to describe the local burnup and high burnup properties based on DIONISIO code [12][13][14].
Here we suggest an analytical and simple formula to calculate local burnup and isotope distributions based on the constant relative reaction rates at the different radius (except for 238 U) and burnup. The present method considers the radial distribution of the neutron flux through two extreme conditions. The parameters are fitted to the reaction rates of the starting burnup level. From this level, the present method are shown to have independent and similar results as the Monte Carlo (MC) simulation, but just a few minutes is needed because its simplicity.

The description of the model
The rate of certain reaction type (x) between neutrons and a given nuclide (A) atan average burnup (bu) can be written as a function of the radius (r) and the energy of neutron (E).
where N(r, A) is the concentration of the nuclide A at the radial position r, φ(r, E) is the neutron flux at the energy E and the radius r, σ(A, x, E) is the microscopic cross section between neutron and A for reaction type x at the energy E. The (n, f ) and (n, γ) reactions are assumed to be the two more important reaction types in the fuel rod. In the present work, the main purpose is to examine the validity of the present model. Thus only the 235 U, 238 U, and 239 Pu isotopes are considered to constrain the model testing to a simple case.
In the present study, the Eq. (1) is assumed to be two terms for (n, f ) and (n, γ) of 235 U and 239 Pu, and the (n, f ) of 238 U: and three terms for (n, γ) of 238 U: where the coefficients c are specified for each nuclide and the reaction type. The coefficient b 1 is expected to be the same for all three nuclides. The coefficient b 2 in the (n, γ) of 238 U is much larger than b 1 .
The following discussion concentrates on the explanations of the above assumptions of the reaction rates. In a thermal reactor, two kinds of the neutron reactions are most important, the scattering reaction in the moderator and the absorption reaction in the fuel rod. When the slowed down neutrons penetrate into the fuel rod, we mainly concentrate on the neutron flux at the thermal and resonance region which are important to the absorption reaction in the fuel rod. The average energy of fission neutron is around 2 MeV, which is much larger than that of the thermal and resonance neutron. Thus no multiplication of the neutron flux is considered at the energy of thermal and resonance region in the following discussion. Two extreme conditions can be considered for certain energy of the neutron. The first one is that the neutrons are strongly absorbed by the fuel, which indicates that the mean free path (MFP) of the absorption, λ(a, E) = 1 Σa(E) = 1 Σ a,5 (E)+Σ a,8 (E)+Σ a,9 (E) , is much smaller than the size of the fuel rod at the energy E. The second one is that few neutrons are absorbed in the fuel at the energy E.
For the first extreme condition, one can further simplify that the velocity of the neutrons are all perpendicular to the surface of the fuel rod. Such assumption is acceptable because the velocity is almost symmetric around the radial directions in the real case. Then the transportation of the neutrons at certain energy E is the solution of the Boltzmann equation in the cylindrical coordinate without the terms of the scattering and the source. The radial part of the equation is: where the negative sign comes from the opposite direction between the velocity and the radial vector, and the solution is: where r 0 is the radius of the fuel rod. Because of the assumption of the strong absorption, the neutron flux is meaningful near the surface of the fuel rod. The divergence at r = 0 should not be considered. Because λ(a, E) is much smaller than the size of the fuel rod, r 0 , the neutron flux decreases very quickly to zero at small (r − r 0 ). The neutron flux φ(r, E) is approximately φ(r 0 , E)e (r−r 0 )/λ(E) near the surface and zero at the other radial region. The reaction rate per nuclei R(r, A, x, E)/N(r, A) is proportional to e (r−r 0 )/λ(a,E) , resulting the third term in Eq.(3). One example of such situation corresponds to the (n, γ) reaction of 238 U at certain energy. The magnitude of the atomic concentration of 238 U is around 10 22 /cm 3 in the fuel rod. For (n, γ) reaction of 238 U, some peaks of cross section at the resonance region can achieve 10 4 barn [15], resulting λ(a, E) = 1 σ a,8 (E)N 8 at the magnitude of 10 −2 cm, which is much smaller than the size of the fuel rod.
For the second extreme condition, the λ(a, E) is much larger than the size of the fuel rod r 0 . The φ(r, E) is almost unchanged in the fuel rod and R(r, A, x, E)/N(r, A) is also approximately constant, resulting the first term in Eq. (2) and (3). The magnitude of the concentrations of 235 U and 235 Pu are around or less than 10 20 /cm 3 in the fuel rod. If the σ a,5;a,9 (E) and the σ a,8 (E) are much smaller than 10 4 barn and 10 2 barn, respectively, the λ(a, E) is much larger than r 0 .
For the situations between these two limits, the second term in Eq.(2) and (3) are assumed with b 1 the same magnitude of 1 r 0 . Please note that the above discussion is restricted in the fuel rod. If the neutrons go out of the fuel rod, they may be slowed and re-enter the fuel rod with a different energy.
In principle, c, b in Eq. (2) and (3) can be calculated through the cross section data and the neutron flux at r 0 . But in some energy regions, the cross sections changes dramatically, and hence it is difficult to obtain the coefficients. Spatially Dependent Dancoff Method can be used to calculated the cross-sections in the resonance region [16]. In the present work, these coefficients are fitted using a MC calculation. An initial burnup level bu 1 is assumed and with it a MC calculation is done for one fuel rod to calculate the reaction rates at different radial positions. The coefficients are then fitted to these reaction rates and used to calculate the burnup levels after bu 1 . After bu 1 , the MC calculation and the present formula are performed independently for all burnup levels, bu 2 , bu 3 and so on. The concentrations can be calculated through: by assuming that the reaction rates do not change during the time T . The 5 time duration T corresponds to the burnup change between two levels: .
where Q is the average energy released by fission and M 0U (πr 2 0 ) is the initial mass of U isotopes in volume πr 2 0 . One can transform: ∆bu.
∆bu. (7) In MC calculation, the reaction rates in above equations are simulated with the concentration of all three isotopes at each burnup level. The present formula calculate all reaction rates through the coefficient fitted to the reaction rate at bu 1 . The R 9,γ (r) reaction rate is neglected in above equations in both MC and analytical calculation because the present model does not include 240 Pu. It is acceptable as present work concentrates on the examination of the formula not on a real burnup problem. The next section shows that the present formula can obtain results for bu n quite close to those obtained with the MC calculation.
The above equation looks similar to the formula in TRANSURANUS model [8,10], in which the relationship between N(r) and bu can be generally written as: Pu, respectively. The σ a and σ c are the one-group effective cross sections for the total neutron absorption and neutron captured, respectively. The cross sections are obtain differently for UO 2 and MOX fuel because of the very different initial concentrations in each fuel and the corresponding different neutron spectrum. A is a conversion constant. The f (r) is the radial form factor, which is 1 + p 1 exp(−p 2 (r 0 − r) p3 ) for 238 U and 240 Pu and unit  for all other nuclide. The f (r) comes from the resonance absorption and the parameters are determined by comparison with measurements [8,10]. The local burnup and isotope distributions can be calculated through Eq. (8). More details can be found in Ref. [8,10]. The present work considers not unit f (r) (actually the different radial neutron flux) for all nuclei.

Calculations and discussions
In the present work, the continuous energy Monte Carlo code TRIPOLI-4 [17] is used for the MC calculations as the starting point of the analytical formula and the reference for after calculations. The geometry of the fuel rod in the present investigation is set to be r 0 = 0.4127 cm, with cladding between 0.4127 and 0.4744 cm. The moderator is in an outside box with the length 1.2647 cm with all surface reflection all neutrons. The fuel rod is divided to seven parts in the MC simulation, with the dividing point located at r/r 0 = 1/2, 3/4, 7/8, 15/16, 31/32, 63/64. The corresponding center r ′ /r 0 of each part is listed in Table 1. For a normal UO 2 fuel, there is no Pu isotopes at the beginning. It is reasonable to start present calculation at a certain burnup bu 1 with Pu included. The MC calculation is done with 3.3% enrichment 235 U fuel in the fuel cell to obtain (n, f ) and (n, γ) reaction rates of 235 U and 238 U. With the reaction rates, the number of 235 U, 238 U, and 239 Pu in a certain burnup can be calculated by assuming the reaction rates do not change in this period. The concentration of 235 U, 238 U, and 239 Pu at bu 1 = 3.4 MWd/kgU is given in Table 1 as the starting point of the following calculations.
The (n, f ) and (n, γ) reaction rates of 235 U, 238 U, and 239 Pu can be simulated by using concentrations listed in Table 1. The corresponding relative reaction rates defined below are also listed in Table 1: The above equations firstly scale the reaction rate by the corresponding concentration. Because the concentration is different at different radial position.
The scaling cancels such difference in the consideration of reaction rate. And then calculate the reaction rates relative to that of (n, f ) reaction of 235 U. We find that the radial distribution of these relative reaction rates are almost constant except (n, γ) reaction of 238 U. Thus all the reaction rates except that of (n, γ) of 238 U can be treated in Eq (2) with the same coefficients. The radial distributions of these reaction rates are fitted as, (9) where c is the reaction rate per nuclei of (n, f ) of 235 U at r = r 0 /4 and canceled in Eq. (7). Please note that the same b 1 is used for all nuclei. From bu 2 , MC and analytical calculations are independent. It is straightforward to see the advantage of transforming the burnup problem to the reaction rate problem in the present work. The reaction rate actually is the most direct quantity to obtain the burnup and the concentration. In many models, the reaction rate must be obtained after solving the neutron flux problem. The present work treats the neutron flux problem in a simple way in section 2 and considers the relative reaction rate. The validity of these considerations will be discussed later. In the following discussion ∆bu is set to be 2.85 MWd/kgU and the total burnup is calculated to bu 15 = 43.30 MWd/kgU. The local burnup from two calculations are presented in FIG. 1. The analytical calculations can give almost the same results as MC calculations. The detailed data of three burnup leves in Table. 2 shows that the average errors are less than 1%.
The radial distributions of U and Pu isotopes at three burunp levels are presented in Table 3. The bu 2 is the first burnup level that MC calculation and present formula are performed independently. The error between two  Table 1) to bu 2 (in Table 3). The average error 0.05% indicates that the present formula well describes the isotope distribution of 235 U at the burnup level. The average error of the isotope distribution of 239 Pu is around 1%. It is partially from the large change of the concentration from bu 1 to bu 2 and partially from the coefficient of the (n, f ) reaction rate for 239 Pu in Eq. (9). The coefficient in Table 1 is decreasing from 2.50 to 2.47. The use of 2.49 in the Eq. 9 contributes to the 1% error between MC calculation and the present formula at bu 2 . At the burnup level bu 8 and bu 14 , the errors become larger, indicating that the reaction rates in Eq. (9) change slightly during burnup. The detail will be discussed later. It can be seen that the radial distribution of N 9 does not change much after a few periods because the production and reaction rates find a balance at such distribution. Both the comparisons of the local burnup and the isotope distributions imply that the present simple and analytical formula can give very nice description if a starting point is given.
As the present model is constrained in a single rod, it is worth to examine its validity in a more realistic case, such as in an assembly. Figure 2 shows the model of an assembly which has the same size for each fuel rod as the previous single one. There are 24 positions for the control rods and one postion in the center for the detectors. In the present calculation, the control rods are not inserted and the rod for detectors are filled with water. The reaction rates of the fuel rods in the assembly are simulated by setting all fuel rods the same   The relative reaction rates, defined by Eq. (9), in the single rod and four rods in the assembly are compared at the burnup level bu 2 in Table 4. The relative reaction rates of all five rods are generally similar, especially for the two important reactions, (n, f ) reaction of 235 U and 239 Pu. The relative reaction rate of (n, γ) reaction of 238 U in the single rod is a few percent less than those in the four rods in the assembly. A slight modification can be applied for the Eq. (9) to fit to the fuel rods in the assembly. In the case of assembly, it is clear to see that the relative reaction rates are almost constant for all four fuel rods, while the reaction rates are a little different depending on the position of the fuel rods, seen from Table 4. Table 5 presents the comparison of the relative reaction rates between the single rod and the rod1 in the assembly at the burnup level bu 9 and bu 15 . As relative reaction rates in the four rods in the assembly are similar to each other, only the data of rod1 are shown in Table 5. At both bu 9 and bu 15 , the (n, f ) reaction of 235 U and 239 Pu are similar in both rods, while the relative reaction rate of (n, γ) of reaction 238 U in the single rod is a few percents larger than that in the rod1. The situations are similar to that at bu 2 . The comparisons in Table 4 and 5 imply that the single rod model is a nice estimation of the fuel rods in the assembly.
It is worth to consider why such simple formula Eq. (9) works. Take R ′ 9,f (r) for example, which is the ratio of the one-group effective cross section. An approximately constant R ′ 9,f (r) indicate a constant neutron dynamics for the ratio between the fission reaction of 239 Pu and 235 U along the radial direction. Although the one-group effective cross section changes as the function of r, the ratio of the one-group effective cross section of each reaction included in the present study keeps almost the same except the (n, γ) reaction of 238 U. One of the advantage of the present method is the consideration of the reaction rate compared with the previous methods [10][11][12]. Because the local burnup and Table 4: The comparison of the relative reaction rates for a single fuel rod and four selected fuel rods in an assembly at bu 2 .  Table 5: The comparison of the relative reaction rates for a single fuel rod and the fuel rod in an assembly at bu 9 and bu 15 . the isotope distributions do not directly depend on the neutron flux but on the reaction rate, it is an alternative method which consider the properties of the latter one. The TRANSURANUS model considers the constant neutron flux and the one-group effective cross sections during ∆bu. The assumption in the present model is actually the constant ratio of the one-group effective cross sections, which is relatively easy to be achieved. From the Table 1, 4, and 5, it is seen that the relative reaction rates changes not much along the radial direction and among each burnup level. It should be mentioned that the relative (n, γ) reaction rates of 239 Pu are also close to constant. It is neglected in the present calculation because it links to other Pu isotopes which are not included in the present study. If the error of a few percents is acceptable, the Eq. (9) can be performed independently until bu 15 because the relative reaction rates are almost independent both on radius (except the (n, γ) reaction of 238 U) and burnup. More exact investigations need consider that the relative reaction rates are not exactly constant as the burnup increasing. Such as from bu 2 to bu 9 , R ′ 9,f decreases from 2.49 to 2.25, seen from Table 4 and 5. This is why the difference of isotope distributions between the MC and present calculations becomes larger as the increment of the burnup. One can use the present formula to obtain the concentration at a certain burnup level, such as bu 4 , re-simulate the reactions with the concentrations through MC, and refit the coefficients in the Eq. (9).
The present calculation is limited to 235 U, 238 U, and 239 Pu to examine its validity in a simple case. It is expected that it can be expanded to more complicated case, such as the consideration of other Pu isotopes and other nuclides which played important role in the reactor, such as the minor actinide and the poisons. In that case, a starting point is also needed to obtain the parameters in the Eq. (9). In addition, the present calculation constrains the ∆bu to compare with the MC calculations. In principle, the ∆bu can be rather small to obtain an almost continues changes of the reaction rates depending on the burnup and radius. In the case of small ∆bu, the assumption of constant reaction rates in Eq. (7) is more reasonable.
For a more exact study, the present method can be used together with the MC calculation. Such as, starting from bu 1 , the reaction rates and the concentrations can be calculated with rather small ∆bu up to bu 2 . The output is used for the next MC calculation. If the ∆bu is rather small, purely MC calculation is very time consuming from bu 1 to bu 2 . But with the supplement of the present work, more exact and relatively quick results are expected.
In practice, one can use the reaction rates in Eq. (9) with the concentra-tions at bu 1 to calculate the ∆N(r) in Eq. (7) with small ∆bu, such as 0.01 MWd/kgU. The next step is to use the new concentrations to calculate the reaction rates in Eq. (9) and then the next ∆N(r) through Eq. (7). After 285 steps, the concentrations at bu 2 are obtained and can be used for the next MC calculation. It is actually the Euler method with very small steps, which is expected to be more exact. In Eq. (7), the reaction rates (actually the neutron flux and the concentrations) are assumed to be constant during ∆bu. But both the neutron flux and the concentrations change when the burnup changes. For detail, the concentration is in the reaction rate in the numerator of the Eq. (7). The effect of the neutron flux in Eq. (7) is reflected by the integral of the reaction rates (the total fission reaction rate) in the denominator. At the different burnup level, the total fission rate and its distribution in each isotopes are different, which corresponds to the different neutron flux. If the calculation between bu 1 and bu 2 are separated to many small steps, the neutron flux and concentrations are reconsidered for each steps.
The present method is also helpful for solving the properties of the fuel rods at the different positions in an assembly or an reactor core. As discussed before, the relative reaction rates are almost constant in the fuel rod at the different position in the assembly. Only one factor is needed to describe the reaction rates in each fuel rod. One can simulate at an initial level through the MC method to obtain the factors for each fuel rod and investigate the evolution of the local burnup and isotope distributions for all fuel rods at certain average burnup level.

Summary
In conclusion, an analytical and simple formula is suggested to calculate the reaction rates as well as the burnup and the isotope distributions as the function of radius. The parameters of the formula are fitted to the reaction rates of a given burnup level. Starting from the same burnup level, the present formula can give very nice description compared with a MC calculation from the code TRIPOLI-4. The reaction rates depend on the concentrations and the neutron flux, which both vary along the radial direction and at the different burnup levels. The present work finds almost constant relative reaction rates on these two degrees of freedom (except the (n, γ) reaction of 238 U on the degree of freedom of the radius), which provides a solid physical explanation on the simple formula used to calculate the local burnup and the isotope distributions.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this manuscript.