DEVELOPMENT OF DECAY HEAT MODEL FOR RAST-K

Decay heat (DH) is the heat produced through a radioactive decay of fission products during or after a reactor operation. It is known as the second largest source of power in the core after fission. Being such a strong contributor to reactor power, it should be accurately determined at any time of reactor operation. Currently, there are two main approaches for DH estimation used in reactor simulation codes. One approach is based on careful inventorying of all produced target nuclides and their individual contributions to total power. Alternatively, the other popular approach is based on collapsing all target fission products into a small number of groups similar to delayed neutron estimation techniques. However, the last (multigroup) method currently has limitations when used in some transient scenarios such as transients occurred in fresh fuel. In this study, the multigroup method was further developed for reducing limitations while retaining the advantage in computation speed. Then, it was implemented into Reactor Analysis code for Steady state and Transient (RAST-K) and tested against other codes. As a result, the improved method was found capable of determining DH power at all tested stages of reactor operation under any tested operation scenario. In particular, the test simulations using the improved method showed better results in those scenarios that were under accuracy limitations of the original multigroup method. Overall, the quality of transient calculations in RAST-K was improved when using the newly implemented DH module.


INTRODUCTION
During a reactor operation, the dominant source of power is energy of fission released from fuel nuclides. Once a heavy nuclide undergoes fission, it splits into two fission product nuclides plus several neutrons that could cause fission in other fuel nuclides, maintaining a certain power level inside the reactor core. Since nuclear fission heavily relies on having available neutrons in the core, it can be controlled using neutron absorber mechanisms such as control rods (CR) or soluble boron added to moderator (for certain types of reactors). However, nuclear fission is not the only source of power in a reactor. Most of fission products (FP) and heavy metal nuclides (HM) produced during an operation are radioactive. The energy released during the radioactive decay of those nuclides is added to the reactor power as decay heat (DH).
Since radioactive decay is beyond human control, the total amount of power produced through this process should be carefully evaluated. As pointed out by Gauld [1], the value of decay power can reach up to 7% of total reactor power during operation or even exceed given value in some transient scenarios. Currently, the most frequently studied case in this area is reactor shutdown as discussed by Schrock [2]. In this scenario, nuclear fuel could contain large quantities of radioactive nuclides collected during the whole cycle operation, and therefore must be cooled for preventing undesired outcomes such as fuel melting. The same careful evaluation of a DH generation in nuclear fuel is necessary for handling spent nuclear fuel (SNF) outside of the core (i.e. SNF pool storage on site or cask storage for transportation) as stated by Zerovnik et al [3].
Despite a relatively straightforward nature of a radioactivity phenomenon, direct measurement of DH during a normal operation of a reactor seems to be challenging. As presented by Tobias [4], the existing experimental approaches include working with small samples of nuclear fuel in various compositions after being exposed to neutron irradiation for a relatively short period of time. A more recent overview of real-life experiments though using similar methodology is shown in the work of Gauld. Considering given information, it could be concluded that currently available experimental data is based on post-irradiation measurements of DH for various fuel compositions and separate nuclides.
Though the undoubtful importance of DH evaluation at times following reactor shutdown, there are many more transient scenarios where DH could play an important role. Because of that, modern reactor simulation codes are desired to be capable of estimating the DH fraction of total power at any time step during the reactor operation. Johnson et al [5] provide a description of two currently used methods for DH evaluation in computer simulation codes. The first method is based on a rigorous calculation of a detailed fuel inventory for each time step. Alternatively, the second method is based on following the multigroup DH curves suggested by the American Nuclear Society (ANS) Standard [6]. The key difference between the methods is the total number of independent nuclide groups that should be updated. As stated by Johnson et al, for the case of a detailed inventorying method, 1119 FP and 30 HM nuclide concentrations should be determined and updated over time. As for the other method, the total number is reduced to only 23 condensed nuclide group concentrations. Based on such a noticeable difference, Johnson concludes that the detailed inventorying method requires much more computation time and resources compared to the multi-group condensation method. Besides that, the target reactor simulation code should be capable of tracking all required FP and HM nuclides in order to use the detailed method of calculation.
In this study, a DH model is supposed to be implemented into Reactor Analysis code for Steady state and Transient (RAST-K) [7]. RAST-K is an in-house computer code that is using 3-dimensional (3D) nodal diffusion methods for Pressurized Water Reactor (PWR) simulations. Currently, RAST-K is capable of tracking only a few FP that are either strong absorbers of neutrons (such as isotopes of Xenon or Gadolinium) or the short-lived decay precursors of those strong absorbers. Considering given preconditions, and since the main advantage of modern nodal diffusion codes is the speed of computation, the multigroup method seems to be a better match for being used in such a code.

DEVELOPMENT OF METHODOLOGY FOR RAST-K DECAY HEAT MODULE
Following the information provided in the introduction, the model that is chosen for being used in RAST-K code is the condensed multi-group model described by Johnson et al. Therefore, the first part of this section contains an overview of a theoretical and practical development of the method, while the second part provides the further development of a given method, which has been implemented in RAST-K.

Current Status of Condensed Multi-Group Decay Heat Model
As stated by Johnson et al, ANS Standard provides data for 23 groups of precursors that can replace the detailed estimation procedure for each individual radioactive nuclide in the core. The key milestone of a given method was reached by Dunn [8]. In his work, Dunn described the overall idea of the multi-group method and also provided a modified version of the ANS Standard approach, in which the total number of precursor groups was effectively reduced down from 23 to 12 and 6 groups. Since the errors of fitting ANS data with a smaller number of decay precursor (DP) curves were found within the same range as the original data uncertainty, it was concluded that less than 23 DP groups could be used for simulation of DH in computer codes.
The approach shown by Dunn was implemented in some real nodal diffusion codes such as Purdue Advanced Reactor Core Simulator (PARCS) [9]. In PARCS, 6 DP groups were re-obtained from the ANS Standard showing a slight difference in final values compared to those calculated by Dunn. However, both Dunn and PARCS show their primary interest in determination of decay power specifically after reactor shutdown. Therefore, instead of using real-time tracking of DH precursors, PARCS is assuming equilibrium concentrations that are easier to calculate and that could fit the given goal with decent accuracy. At the same time, this assumption could lead to having an increased error in some transient scenarios such as those that could be found in Section 3 of present study.
In order to eliminate stated limitations and potential errors for any possible transient case, the condensed multi-group method for RAST-K was further improved. The key requirement for the modified method was the capability to calculate DH precursor concentrations (DPC) and corresponding decay power at any time step during a reactor operation or after reactor shutdown.

Derivation of Decay Heat Equation for RAST-K
The derivation of the final formula for DH calculation used in RAST-K starts from the differential equation (DE).
In equation (1) and in further equations, C(t) -DPC for a single group at time t, -a decay constant for a given single group, -group DP yield (from fission), F -fission rate at a given moment of time t in a certain volume of a reactor (a function of neutron flux). Specific indexing for separating different DP groups is omitted for brevity. The final solution of equation (1) is shown below.
The detailed derivation of the formula (2) is similar to the one provided by PARCS Theory Manual. In a similar manner, the value of equilibrium DPC for certain fission rate could be obtained.
For estimating the value of DH power at time t, the following boundary conditions could be applied: the minimum value of DH is zero when all DPC are equal zero; the maximum value of DH energy for given fission rate ( ) is around 13 MeV per one act of fission as stated by Ripani [10]. It could be reached when group DPC reach equilibrium values for given fission rate. Based on that, the equilibrium DH power for a certain fission rate could be written in the following form.
Considering that the equilibrium DPC does not represent the entire fuel cycle, a weighting factor for estimating the current available fraction of decay power for time t and fission rate ( ), is added into equation (4) yielding the final formula for determination of DH.
Since the value of weighting factor could be below or above 1.0, the equation (5) is applicable for any transient scenario including scenarios with freshly loaded fuel. A major simplification of equation (5) could be introduced using a substitution of equation (3) into (5). In addition, as stated by Tobias, the ANS Standard includes group-specific energy-per-fission yields instead of number-based yields. As a result, the value of DH energy could be substituted by ( • ) where is a theoretical recoverable energy released from fuel per one act of fission. This yields the equivalent DH formula, which shows the value of DH produced by one group of DH precursors with concentration stated in energy units.
To determine the total DH value, all group-specific values of DH (6) should be combined together as shown in the following equation, in which n -total number of DP groups, and ( ) -DPC (in energy units) for one group i updated using equation (2).
In this study, in order to compare the produced result with the method used in PARCS code, the group constants were chosen equal to the default values used in PARCS (i.e. 6 DP groups). The simulation cases and the results are described in the following section.

RESULTS OF SIMULATIONS
For a preliminary testing of the DH model implemented in RAST-K, two transient scenarios were chosen. The first case is a full-power insertion of one control bank into the core based on NEACRP core geometry [11]. The second case is a problem similar to a well-known NEACRP Rod Ejection Accident (REA) problem. The used reactor type is a PWR with 2775 MWth power. Values of Critical Boron Concentration (CBC) for all cases are determined in RAST-K and manually input in PARCS. Other simulation-related parameters are described in specific subsections of this section below.

Full-Power Rapid Control Rod Insertion Accident
In the first case, the chosen reactor model is running at Hot Full Power (HFP). As given in the original benchmark description, NEACRP problem is using the initial cycle freshly loaded fuel. All CR are ejected before the start of simulation. At time 20 sec after the beginning of simulation, one CR bank is fully inserted into the core. For RAST-K, several pre-calculations of DPC were simulated (0 sec, 1000 sec, 30 days). As for the PARCS model, only the equilibrium value of DPC could be used for the DH 'ON' mode. The result of simulation is shown on Fig. 1. As provided on Fig. 1, the decrease of power is more significant when DH is not activated. At the same time, since PARCS is assuming the DPC values as equilibrium, it shows the highest theoretically possible value of DH even for fresh fuel. Hence, the more appropriate result is expected to be as shown by "RAST-K Decay Heat 'ON' 0s" curve (light blue). After the increase of pre-calculation time (1000 sec and 30 days as shown on Fig. 1), the value of power after the CR bank insertion is also increasing. Since the value of DPC for the pre-calculation time of 30 days is expected to be relatively close to the equilibrium, the corresponding curve on the Fig. 1 is almost equal to the curve obtained using PARCS DH 'ON' mode.
Finally, a CR-induced reactivity impacted by using certain DH models of PARCS and RAST-K was obtained and compared to each other as shown on the bottom right-hand side of Fig. 1. The resulting curve shows that the absolute value of CR-induced reactivity after the transient is smaller for the cases of large stored DPC such as "PARCS Decay Heat 'ON'", "RAST-K Decay Heat 'ON' 30d". This is due to a smaller relative fraction of fission power in those cases, which is directly affected by the movements of CR. The results for RAST-K pre-calculations for 0 sec and 1000 sec show that there are relatively small amounts of DPC at the Beginning of Cycle (BOC) leading to a higher CR-induced reactivity at that stage of fuel cycle.

NEACRP Control Rod Ejection Accident
The next tested scenario is two separate cases of NEACRP REA Problem. The first case is a Hot Zero Power (HZP) ejection of one CR bank. The second scenario is similar to the first one aside from being simulated at HFP and having a transient after 20 sec of full-power operation. As stated in the description of NEACRP Benchmark, the simulation is performed for the freshly loaded fuel of the initial cycle.
For the HZP case, the behavior of power after a rapid ejection of the chosen CR bank is shown on Fig. 2.

Figure 2. Power change for the NEACRP HZP REA scenario (bottom left-hand side -impact of DH on the value of power, bottom right-hand side -impact of DH on the CR reactivity).
Similar to the previous scenario, a difference of behavior between DH 'ON' and DH 'OFF' cases could be observed. Due to having zero power before the transient and consequent zero DPC, the curve of "RAST-K Decay Heat 'ON' 0s" is covering the curve of "RAST-K Decay Heat 'OFF'" on Fig. 2. However, a distinctive difference between the curves of "PARCS Decay Heat 'ON'" and "PARCS Decay Heat 'OFF'" shows that PARCS is using the equilibrium value of DPC for all cases of DH 'ON', which is not appropriate for the given transient scenario. In order to prove given statement, a comparable difference between "RAST-K Decay Heat 'ON' 30d" and "RAST-K Decay Heat 'OFF'" results could be obtained using a large DPC pre-calculation time such as 30 days (red line on Fig. 2). As for the DH impact on CR-induced reactivity, no significant difference was found assumedly due to having zero power before the CR bank ejection. For the last studied case, which is a HFP ejection case, a larger power difference between the studied models compared to the HZP case could be observed as shown on Fig. 3. This case could be considered as the opposite scenario to the one discussed in Section 3.1 of this paper. Therefore, it shows a reverse behavior of power with a smaller power increase for DH 'ON' cases and a larger increase for DH 'OFF' cases during transient. Since the transient occurred after 20 sec of fullpower operation, DPC for all models including "RAST-K Decay Heat 'ON' 0s" is non-zero. Hence, a differentiation of power levels caused by different DPC for each RAST-K pre-calculation case could be found on Fig. 3 (bottom left-hand side). The difference of CR-induced reactivity between the models is also found more noticeable compared to the HZP case.

CONCLUSIONS
In this study, a multigroup method of DH evaluation was improved in order to close known gaps that could be found while modelling certain reactor dynamics problems. The new version of the method was implemented in RAST-K code and tested against PARCS code using well-known benchmark problems. applicable for such case. On the other hand, the power change curves calculated in RAST-K were found in good accordance with the original expectations. Finally, this case showed the highest impact of DH on the value of CR-induced reactivity among all test cases shown in this paper.
The second and third tested cases were BOC HZP and HFP CR ejection problems, respectively. It was found that the initial value of core power influences the impact of DH on power change and CR-induced reactivity change during the transient. In particular, a smaller difference between DH 'ON' and 'OFF' modes was found for HZP, in which the initial value of power was zero. However, both cases showed a clear distinction between the non-equilibrium-based values obtained using RAST-K DH model, and equilibrium-based values obtained using PARCS. Similar to the first scenario, equilibrium values of DPC were found not suitable for being used in fresh fuel.
In future work, we expect to perform a detailed testing of the new RAST-K DH module based on other steady-state and transient scenarios. As for the future work expected from the industry, a new decay heat benchmark that is using other than reactor shutdown case is desired to be developed (using experimental data). That would allow reactor simulation code developers to further improve the quality of their codes.