Abstract
Multiphysics coupling of neutronics/thermal-hydraulics models is essential for accurate modeling of nuclear reactor systems with physics feedback. In this work, SCALE/TRACE coupling is used for neutronic analysis and spent fuel validation of BWR assemblies, which have strong coolant feedback. 3D axial power profiles with coolant feedback are captured in these advanced simulations. The methodology is applied to two BWR assemblies (2F2DN23/SF98 and 2F2D1/F6), discharged from the Fukushima Daini-2 unit. Coupling is performed externally, where the SCALE/T5-DEPL module transfers axial power data in all axial nodes to TRACE, which in turn calculates the coolant density and temperature for each of these nodes. Within a burnup step, the data exchange process is repeated until convergence of all coupling parameters (axial power, coolant density, and coolant temperature) is observed. Analysis of axial power, criticality, and coolant properties at the assembly level is used to verify the coupling process. The 2F2D1/F6 benchmark seems to have insignificant void feedback compared to 2F2DN23/SF98 case, which experiences large power changes during operation. Spent fuel isotopic data are used to validate the coupling methodology, which demonstrated good results for uranium isotopes and satisfactory results for other actinides. This work has a major challenge of lack of documented data to build the coupled models (boundary conditions, control rod history, spatial location in the core, etc.), which encourages more advanced methods to approximate such missing data to achieve better modeling and simulation results.
1. Introduction
The ability to accurately and efficiently predict end of life (EOL) fuel composition is essential in ensuring safe and economic management of spent nuclear fuel. In addition, verification, validation, and uncertainty quantification of neutronics/depletion codes with spent fuel data are important in order to evaluate the accuracy of such physics codes. Unfortunately, fuel depletion calculations are expensive to perform, making brute-force uncertainty analysis methods difficult to apply [1]. In single physics or decoupled simulations, it can often be difficult to directly simulate a particular experiment lacking detailed data. For example, time-dependent axial coolant density distributions are needed for light water reactor simulations, especially for boiling water reactors (BWRs) but data on this are rare. Approximations have been used in [2, 3] where an analytical form for axial coolant temperature distribution for a pressurized water reactor (PWR) is used. For BWRs, even with reported experimental data, major assumptions were needed to be made about the behavior of axial coolant density distributions in order to create verification models in [4]. However, such assumptions can introduce uncertainties in the calculations because void fraction feedback is not captured [5, 6]. The current study seeks to avoid this difficulty by coupling the thermal-hydraulics code TRACE to the T5-DEPL depletion sequence in SCALE for validation of spent fuel composition, thereby avoiding the need for detailed history of void distribution data.
Often, verification studies that use spent fuel compositions as a comparative metric rely on the concept of a ratio of computed to experimentally determined isotopic concentration or code bias (C/E) [2, 4, 7]. The results from these verification studies can be further used for uncertainty quantification (UQ) [8]. Many methods of uncertainty quantification based on Monte Carlo parametric sampling [3, 9, 10] exist; data on C/E can be helpful in approximating uncertainty bounds of the isotopic concentrations in spent fuel analysis, similar to what has been done in [11]. An uncertainty analysis of a spent nuclear fuel storage system was performed in [5]. POLARIS was used to calculate the composition of commercial BWR spent fuel for 77 samples. This study [5] presented and compared BWR spent fuel isotopics for both actinides and fission products. Then, the study compared of a spent fuel cask based on code-predicted isotopic compositions with a cask based on experimentally reported isotopic compositions.
Neutronics and thermal-hydraulics coupling involves exchanging data on axial power distribution and coolant characteristics across physics codes [12]. This process can be done with external coupling [13], which is the method used in this study for validation of spent fuel isotopic composition (see Section 3). A variety of neutronics and thermal-hydraulics coupling studies have been conducted previously, which include coupling to subchannel codes (e.g., CTF), coupling to system codes (e.g., TRACE), and coupling to computational fluid dynamics codes (e.g., ANSYS/FLUENT). Internal coupling between Serpent 2-SUBCHANFLOW was carried out in [14], and the coupling scheme was benchmarked against TRIPOLI4-SUBCHANFLOW and MCNP5-SUBCHANFLOW. In addition, MCS neutronics was coupled to the CTF subchannel code to perform 3D pin-by-pin analysis of the BEAVRS benchmark [15]. Applications of coupling involving system codes have been used to simulate BWR accidents such as anticipated transient without scram using Simulate-3K/RELAP5 [16] and BWR instability accidents with TRACE/PARCS [17]. Lastly, nuclear applications that involve coupling to computational fluid dynamics codes were performed in [18, 19].
Overall, the present study seeks to develop advanced depletion models which include coupling to the thermal-hydraulics system code TRACE to aid in verification and validation of neutronic analysis. The spent fuel isotopic data will be used for validation. 3D neutronic models are developed using SCALE/T5-DEPL in order to capture axial effects in the lattice (e.g., axial power and nonuniform fuel composition). Also, TRACE is coupled to these neutronic models to provide a burnup-dependent axial void distribution. The coupling methodology is applied to two BWR benchmarks with reported isotopic measurements at end of cycle. The remaining sections of this paper are organized as follows. Section 2 presents brief descriptions of the BWR benchmarks used in the validation process. Section 3 presents the SCALE/TRACE coupling methodology. The results of this study and their discussion are presented in Section 4, while the conclusions drawn from this work are given in Section 5.
2. Description of Test Cases
2.1. Fukushima Daini-2 2F2D1 Assembly
The first benchmark used in analysis is the 2F2D1 BWR assembly from the Fukushima Daini-2 reactor. The design is 8 × 8-4 with 60 fuel rods and one water rod [20]. The top portion of Figure 1 presents the radial layout of the assembly. In this figure, the location of a special rod is marked, labeled “F6,” where samples have been taken for isotopic analysis [21]. The F6 rod has an active height of 371 cm with 4.5 wt.% enrichment. Isotopics are measured at three different axial elevations shown in Figure 1(a). Within the dataset, an average void fraction and discharged burnup are reported for each sample. Sample power history is determined by scaling the reported total assembly power to enforce the given burnup to each sample.
(a)
(b)
2.2. Fukushima Daini-2 2F2DN23 Assembly
The second benchmark used in this study is taken from the same unit as the previous. This is the 2F2DN23 assembly from the Fukushima Daini-2 unit [22]. Similarly, this is an 8 × 8-2 design but with two smaller water rods. The bottom portion of Figure 1 presents the radial layout with pin-by-pin enrichment of this assembly. Again, the eight locations in which samples have been removed for analysis from a selected rod, the “SF98” rod, are labeled. For each of these samples, time-dependent power, constant and average void history, and discharged burnup are reported. The sample sections taken from this rod are all 0.5 mm in thickness. Sample 1 is taken from the naturally enriched (0.71 wt.%) region of the rod.
3. Methodology
TRACE [23] is a two-phase flow solver and best-estimate reactor system code developed by U.S. Nuclear Regulatory Commission (NRC) to analyze loss of coolant accident (LOCA), operational transients, and other accident scenarios within light water reactors (LWRs) such as PWR and BWR. TRACE is capable of simulating thermal-hydraulic phenomena in LWRs with high accuracy. It is a 1D (z-direction) code that solves the time-averaged and area-averaged two-phase two-fluid model equations. KENO-V.a is a 3D Monte Carlo neutronic solver in the SCALE code system developed by Oak Ridge National Laboratory [24, 25]. KENO-V.a performs 3D Monte Carlo criticality calculations for arbitrary geometries, and it operates in multigroup or continuous energy mode. KENO-V.a is able to determine , neutron lifetime, power and flux densities, and many other important reactor physics quantities.
3.1. Multiphysics Coupling
In this study, the nuclear fuel depletion sequence (T5-DEPL) is coupled to the thermal hydraulic code TRACE. By coupling a neutronics and thermal hydraulics code, there is no longer a need to rely on axial void fraction distributions—as these data can be missing or imprecise. A parametric description of the coupling scheme between these two simulated systems is shown in Figure 2. Internal to the T5-DEPL depletion sequence, KENO-V.a uses Monte Carlo neutron transport to calculate the axial power profile in all nodes, while ORIGEN is used to deplete the fuel in each node based on its respective power level. Externally, TRACE uses the nodal powers to calculate coolant density and temperature and return them to the T5-DEPL sequence for a new neutronics calculation. The coupling procedure is shown in Figure 3 and can be summarized as follows:(1)The TRACE model is executed first with uniform power in all axial nodes. After completion, coolant density and temperature at each axial node are extracted from the TRACE output.(2)The coolant density and temperature in all axial nodes are used to update the SCALE/T5-DEPL model. The updated neutronics input is executed and the axial power profile is extracted and normalized.(3)The updated axial power distribution is used in a new TRACE calculation, which yields updated coolant temperature and density.(4)Steps 2-3 are repeated until the convergence criterion is met. The convergence criterion uses the average relative difference between two successive iterations in all axial nodes. When this quantity falls below a threshold value (e.g., 6%) for axial power, coolant density, and coolant temperature, the convergence criterion is satisfied and the loop is terminated.(5)Steps 2–4 are repeated until all burnup steps are completed. The next burnup step is initiated upon completion of step 4.
3.2. Multiphysics Modeling
SCALE/T5-DEPL is used to model both the 2F2D1 and 2F2DN23 fuel assemblies in 3D with a reflective boundary condition radially and a vacuum boundary condition axially. The assembly is modeled in a square channel without rounded corners. The Monte Carlo parameters used are 1100 cycles with 20,000 neutrons per cycle, with an initial 100 cycles skipped. These parameters in KENO-V.a ensure about ∼10 pcm uncertainty in .
For 2F2D1, twenty three nodes are used in the axial direction, which include three 3.5 cm thick axial slices to represent the samples (TU101, TU102, and TU106) in the F6 rod. For 2F2DN23, twenty one axial nodes are used, and they include eight 0.5 mm thick axial slices to represent the samples in the SF98 rod. For each sample tested, an irradiation/power history must be given to SCALE/T5-DEPL. Fortunately, there is an option in this code to ensure that the sample power can be specified. The option “power basis” in SCALE/T5-DEPL is used for each measured sample in the rods SF98/F6. To explain, the total power for all materials in the problem will be normalized such that the power in the sample material matches the sample power history specified. The burnup of each sample is measured during experiment and is reported in Figure 1. Matching this burnup value is essential in SCALE/T5-DEPL to achieve accurate validation of fuel isotopics. Fortunately, for the SF98 rod, the sample power is provided directly from the benchmark for each of the eight samples, and it is shown in Figures 4(b)–4(d).
(a)
(b)
(c)
(d)
However, only the total assembly power history is given for rod F6. Therefore, the assembly power provided by the benchmark is scaled such that the final burnup matches the sample burnup (i.e., 18.2, 16.1, and 14.0 GWD/MTU in Figure 1). Only the specific power values (kW/kg) are changed, while the total cycle time (418 days for 2F2D1) and time step length are maintained as in the benchmark. The adjusted sample power history is shown in Figure 4(a). The original power history and other neutronic parameters of this benchmark can be found in the SFCOMPO database [21].
TRACE is used to model the two assemblies as a BWR channel (CHAN) component with inlet boundary condition as a FILL component and outlet boundary condition as a BREAK component. The inlet flow rate is set to 17.56 kg/s, while the inlet and outlet coolant temperatures are set to 550 K and 559 K, respectively. The total number and size of the axial nodes in TRACE match those in KENO-V.a. We ensured that mesh sizes in both KENO-V.a and TRACE are appropriate to avoid large numerical uncertainties in the calculations. The extraction of the coolant temperature and density from TRACE is done using the PyPost postprocessor. Notice that since we model the assembly in TRACE as a single component, then the coolant properties calculated by TRACE are generalised for all rods including the SF98/F6 rods in the assembly. To have pin-by-pin resolution of the coolant distribution, subchannel codes (e.g., COBRA-TF) need to be used.
Since this problem involves depletion for multiple cycles, burnup-dependent simulation increases the complexity of the coupling process. Based on preliminary tests, we found that TRACE parameters tend to converge after 300–500 seconds following any power changes from KENO-V.a. Therefore, to reduce TRACE computational time, only a small time is simulated for coolant properties. Moving on, restart capabilities in SCALE are activated to keep track of isotopic composition. The converged values of material isotopic composition of the previous burnup step are used to initiate the calculations of the next time step. Then, on the last burnup step, the converged values of isotopic composition are used for validation. Lastly, it is worth acknowledging that the multiphysics feedback between the two physics is performed only when the power level changes. For 2F2D1, 8 power changes occur, while 2F2DN23 experiences 13 power changes; all are shown in Figure 4. Therefore, the convergence of all coupling parameters is ensured 8 times for 2F2D1 and 13 times for 2F2DN23, over the whole depletion time.
4. Results and Discussion
The results of this study are presented in two sections. First, verification of the coupling process is performed by analyzing the convergence and axial power behavior at the assembly level (i.e., the power basis option is deactivated and all assembly materials are depleted by specifying assembly power). In the second section, validation of the coupling process is performed by analyzing fuel isotopics at the sample level (i.e., the power basis option is activated by specifying sample power).
4.1. Neutronic Analysis (Assembly-Based)
Convergence results at three burnup steps (beginning, middle, and end of irradiation time) are shown in Figure 5 for the 2F2DN23/SF98 benchmark and in Figure 6 for the 2F2D1/F6 benchmark. The convergence of the three coupling parameters: axial power, coolant density, and coolant temperature, is monitored. The convergence criterion is set to ensure that the average relative difference of all axial nodes falls below a threshold of 6% with at least 3 iterations executed. The threshold value is expressed in relative form regardless of the coupling parameter unit and it balances between accuracy and computational costs. If the average relative difference of all axial nodes falls below 6% on the second iteration, another iteration will still be performed to reach 3 in total. The third iteration is used to identify the convergence, i.e., convergence cannot be ensured with two iterations only. However, if three iterations are not sufficient to converge, additional iterations are executed until all average relative differences fall below 6%.
(a)
(b)
(c)
(a)
(b)
(c)
For 2F2D1/F6, the results show fast convergence for all parameters as they all converge after 3 iterations (the minimum limit) at the three selected burnup steps. It is clear that the coupling parameters seem to vary little between iterations, implying that the feedback from coolant is not significant for this benchmark. This can be justified by the almost steady power level in Figure 4(a), which causes the void fraction to be almost constant with burnup. Therefore, for 2F2D1/F6, the evolution in the axial power profile captured by the SCALE/T5-DEPL 3D model is more significant than the evolution of the coolant properties in TRACE when considering reactor behavior over burnup. On the other hand, for 2F2DN23/SF98, the results behave differently, as the system may converge within 3, 4, or 5 iterations, implying strong feedback from the coolant. This can be attributed to the higher variability in the power history associated with 2F2DN23/SF98, as can be observed in Figures 4(b)–4(d).
After confirming the convergence behavior of the coupling process, it is good to verify the change of some neutronic parameters with burnup. Figure 7 shows the variation of assembly with burnup. For 2F2DN23, we can notice that decreases significantly at beginning of cycle due to the xenon effect. Afterward, starts to increase due to the depletion of gadolinium poison, until it reaches the reactivity peak, which is close to 10 GWD/MTU. Next, continues to decrease toward the end of cycle due to the depletion of U-235. This is a typical behavior for BWR assemblies that are loaded with gadolinium rods. This behavior can be observed also for the 2F2D1 assembly, except that the simulated cycle length is shorter. We can see in Figure 7(b) that the reactivity peak is very close to 12 GWD/MTU, which is the end of cycle for 2F2D1. It is worth noticing that for some depletion steps in Figure 7(a), seems not to converge (e.g., between 15 and 20 GWD/MTU), since we did not assign a constraint on converging as a condition to exit the coupling loop, as described in Figure 2. Doing so will increase the computational costs significantly as is a sensitive parameter. Nevertheless, the depletion trend is still clearly seen for this problem in spite of some incompletely converged values.
(a)
(b)
The evolution of the axial power for the two sampled rods is plotted in Figure 8 for different burnup steps. For the SF98 rod, the samples span the axial direction. Therefore, the locations of the eight samples are used in Figure 8(a) for plotting. Since the irradiation time of the SF98 rod is long, the power profile experiences significant changes during the cycle time. On the other hand, for F6, all axial nodes are plotted due to the limited number of measured locations. The power profile for this benchmark did not experience a large change due to the short irradiation time. The bottom-peaked shape of the axial power is caused by the void distribution feedback of TRACE. Only small changes in the power peak at the channel bottom occur during burnup.
(a)
(b)
4.2. Isotopic Analysis (Sample-Based)
To begin, all actinides and fission products are tracked, but six major actinides are plotted as a function of burnup at the sample level. The plots for these results are shown in Figure 9 for 2F2D1/F6 (three samples) and Figure 10 for 2F2DN23/SF98 (eight samples). As expected, the U-235 and U-238 concentrations decrease over burnup, while the remaining 4 actinides tend to build up as burnup increases. The differences in burnup for the last reported concentration in each sample are due to the differences in the sample burnup (see Figure 1). In general, uranium isotopes for all samples tend to change with the same slope. However, it is clear that plutonium isotopic concentrations are more sensitive to depletion conditions, as these isotopes build differently for each sample. The decay time before the rod measurement ( days after rod discharge) causes a sharp decrease in Pu-241 concentration at end of life. This is also the case for the small increase in Pu-239 at end of cycle, due to the decay of other fission products into Pu-239. In addition, Figure 10 shows that sample 1 discharge burnup is much smaller than other samples, due to its location in the natural uranium region (0.71 wt.%).
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
The validation results of fuel isotopics are demonstrated in Figure 11 for selected isotopes in the two benchmarks. In addition, numerical values of the C/E ratio are listed in Table 1 for 2F2DN23/SF98 and in Table 2 for 2F2D1/F6. For completeness, the experimentally determined isotopic data are reported in Tables 3 and 4 of Appendix A, which can be used in conjunction with Tables 1 and 2 to obtain the calculated results. In the scatter plot of experimental versus calculated, the points that are closer to the diagonal line indicate better agreement with experiment. In this plot, the uranium isotopes as well as Nd-148 tend to have good agreement with data as can be inferred from their good C/E values. The good agreement for uranium isotopes can be attributed to their large content in the fuel and their low measurement error. Also, it was seen early that these isotopes are relatively insensitive to depletion conditions. The plutonium isotopes tend to have fair agreement with data, as their relative error can be as low as 2%, but as high as 64%.
(a)
(b)
In terms of uncertainty sources associated with this work, we believe uncertainty comes from three major sources: first, SCALE/T5-DEPL and the nuclear data libraries. Uncertainty in nuclear data, especially for minor actinides and fission products, can be a large factor particularly because the 56 multigroup energy library was used, which could be less accurate as compared to using continuous energy cross section library. The 56-group energy library was used because SCALE/T5-DEPL depletion calculations tend to be very slow with continuous energy cross sections. This coupling process requires multiple SCALE/T5-DEPL runs, so using continuous cross section libraries becomes impractical in terms of computational costs. The second uncertainty source is associated with the TRACE model, which is directly connected to the third source about lack of reported data. As this is a neutronics benchmark, very little information is provided to aid in building an accurate thermal-hydraulics model representing the fuel channel. Therefore, some assumptions about boundary conditions based on commercial BWRs (inlet flow rate, outlet temperature, friction coefficient, etc.) had to be made. These assumptions would affect the quality of coolant density and temperature predictions by TRACE. Additional missing pieces associated with neutronics include configuration of the assemblies that surround 2F2D1/2F2DN23 in the core, control blade history, SF98/F6 rod burnup, and fine-scale power history compared to the approximated profile provided by the benchmark. These factors are very important for realistic 3D modeling. Despite the previous challenges, satisfactory results are achieved at some axial elevations as indicated by the last row of Tables 1 and 2. The average relative difference for all isotopes can be less than 15% for some samples but can also be higher than that for other samples. It is worth mentioning that the current validation results are similar to previous efforts such as [4, 5, 8].
Finally, since lack of documentation was a major challenge for this work, future efforts will focus on incorporating machine learning and data science to aid in identifying the missing input parameters such that they can yield the desired output. Techniques such as Gaussian process optimization and deep reinforcement learning can be efficiently used for this purpose. These efforts can help in completing the current benchmarks to achieve accurate modeling and simulation.
5. Conclusions
Multiphysics coupling of SCALE/TRACE is used in this work for neutronic analysis and spent fuel validation of BWR spent fuel isotopics. The coupling is performed externally where SCALE/T5-DEPL module sends axial power profile to TRACE, which returns coolant density and temperature. This process is repeated until convergence of all quantities at each time step. The methodology is applied to two BWR assemblies, discharged from the Fukushima Daini-2 unit. The neutronic analysis on assembly level shows that the coupling process is verified, where the 2F2D1/F6 benchmark seems to have insignificant void feedback. Validation of the coupling process is done through spent fuel isotopics data, which demonstrated good results for uranium isotopes and satisfactory results for other isotopes, given the significant lack of documented data faced in this work. Future work will include applying machine learning methods to help in identifying the missing data such that more accurate modeling and results can be achieved.
Appendix
A. Experimental Data
In this section, the experimentally determined isotopic data are reported in Tables 3 and 4 as reproduced from [21]. These data, in conjunction with Tables 1 and 2, can be used to calculate the simulated data.
Data Availability
The computer codes are from RSICC and U.S. NRC, and the data are from OECD/NEA; any person who is able to obtain codes from RSICC and U.S. NRC and data from OECD/NEA is able to do the same type of coupling and validation.
Conflicts of Interest
The authors declare that they have no conflicts of interest.