IMPROVEMENTS TO THE MODELING OF THE TREAT REACTOR AND EXPERIMENTS

This paper summarizes the latest improvements and lessons learned from the modeling and simulation of the transient test reactor at Idaho National Laboratory using the MAMMOTH reactor physics application. MAMMOTH is a MOOSE-based, Finite Element Method application that specializes in the analysis of the spatial dynamics behavior of nuclear reactors. Since early 2018 several transient tests have been conducted at TREAT, thus providing the opportunity to apply and benchmark modern modeling and simulation tools. MAMMOTH was used to provide predictions of the power coupling factor between the core and the experiment for various experiments. Even though the power coupling factor predictions agree very well with the experimental data, within the bounds of the experimental uncertainty, one shortcoming was the underprediction of the total energy deposited in the core and experiment. Determination of the sources for this discrepancy is ongoing, but several key problems have been identified and resolved, thus providing valuable insights for future research. This paper discusses several of these lessons learned. First, the heat capacity data for the TREAT fuel has some significant problems due to limitations of the measurement techniques used circa 1960s. The sensitivity of the peak power and the total energy deposition to various representations of the heat capacity is approximately 5%. Second, the effects of the biological shield and thermal column on the modeling of the core are non-negligible, since they affect the mean generation time and the effective reflection of neutrons back into the core, which is suspected to be important during the core heat up. Matching the reactor period resolves the fact that the reduced spatial domain used in the MAMMOTH model underpredicts the mean generation time. The neutron reflection from these regions is marginally improved with the use of an albedo boundary condition. Third, modeling of the control rod movement with a multi-scheme method is introduced and its current limitations are exposed. Fourth, we explore the effects of using a homogenized model with Superhomogenization equivalence and how that differs from fully heterogeneous simulations. Finally, the energy condensation effects for this graphite core are significant. Solutions with 10 and 26 energy groups show the benefits of using a finer coarse group structure.


INTRODUCTION
TREAT is an air-cooled, thermal, heterogeneous test facility that resumed operation in 2018. The major, unique features of TREAT include a large flux integral absorption due to its high heat capacity; inherent, instantaneous temperature-dependent shutdown mechanism; rapid transient rod movement; and visual access to the core center [1]. The TREAT core is driven by highly enriched uranium (HEU) dispersed in a graphite matrix. At the center of the core, fuel is removed allowing for the insertion of an experimental test vehicle. TREAT's design provides experimental flexibility and inherent safety during neutron pulsing. This safety stems from the graphite in the driver fuel having a strong negative temperature coefficient of reactivity as well as graphite acting as a temperature sink. A top view of the reactor configuration is shown in Figure 1. The TREAT core consists of a 19x19 core lattice around an experiment test region that is labeled "test hole" in the figure. This active core region is surrounded by a radial reflector and a biological shield. The radial reflector is ∼61 cm thick and does not constitute an asymptotic reflector due to the large neutron migration area. A graphite thermal column is located on the east side of the core. The core lattice is filled with standard fuel, control rod fuel, dummy and slotted elements, as well as other specialized elements with instrumentation. The MAMMOTH reactor physics [2,3] and Rattlesnake radiation transport applications are developed and used at the Idaho National Laboratory to model the TREAT reactor and associated experiments. Since early 2018 several transient tests have been conducted at TREAT, thus providing the opportunity to apply and benchmark modern modeling and simulation (M&S) tools. MAMMOTH was used to provide predictions of the power coupling factor between the core and the experiment for various experiments. The total power predicted with MAMMOTH and the measured power from three separate detector channels are shown in Figure 2 for two transients. The first is a temperature limited transient, in which the transient control rods are kept out of the core until 60 seconds from the initiation of the transient. The second is a clipped transient, in which the control rods are quickly re-inserted into the core in order to control the total energy deposition. The key measurement for TREAT experiments is the energy coupling factor, which is the specific energy deposited in the fuel sample divided by the total energy deposited in the TREAT core, normally in units of J/g − UO 2 /MJ. The MAMMOTH predictions of the energy coupling for all transients that were evaluated agree very well with the experimental data, within the bounds of the experimental uncertainty. One shortcoming was the underprediction of the total energy deposited in the core and in the experiment, 18% for the temperature-limited and 15% for the clipped transient, respectively. This is manifested in the difference behavior of the tail in the power trace. Determination of the sources for this discrepancy are still ongoing, but several key problems have been identified, thus providing valuable insights for future research. This paper summarizes the latest improvements and lessons learned from the M&S of TREAT.

CODES AND NUMERICAL METHODS
The INL is currently evolving the M&S capability that will enable improved core operation as well as design and analysis of TREAT experiments. This M&S capability primarily uses MAMMOTH, a reactor physics application being developed under the Multiphysics Object Oriented Simulation Environment (MOOSE) framework [4]. MAMMOTH allows the coupling of a number of other MOOSE-based applications including: Rattlesnake [5,6] for neutron transport, RELAP-7 [7] for low-resolution thermal-fluids, and BISON [8] for fuel performance analyses. Two key methods in the Rattlesnake application are deployed in this work: 1) the generation of reactor kinetic parameters for a particular model and 2) the multischeme method [9], which is used to couple low and high-order transport discretizations within the same problem.
Serpent [10], a three-dimensional continuous-energy Monte Carlo reactor physics code developed at VTT Technical Research Centre of Finland, was used for generating the tabulated macroscopic cross sections. It was selected because it offers 3-D spatial homogenization and group constant generation for deterministic reactor simulator calculations. The neutron cross sections used in this work are based on ENDF/B-VII.r1. Cross sections and reference fluxes from Serpent are then employed in a superhomogenization equivalence (SPH) in Rattlesnake, thereby allowing the preservation of the reference reaction rates. The PJFNK-SPH method [11,12] has been implemented in MAMMOTH to use MOOSE's non-linear solver (PJFNK) to easily converge problems that contain reflector regions with imposed void boundary conditions. It is noteworthy that the preparation of cross sections and reference fluxes is performed from full core simulations.

Heat Capacity
The specific heat capacity of the TREAT fuel is the principal thermophysical property when modeling the TREAT dynamics, since it is directly correlated to the amount of temperature feedback and, thus, the power level of the core. The measurements for the TREAT fuel enthalpy were conducted in the 1960s [13] and included two sets of measurements for two temperature ranges using different measurement techniques. The specific heat capacities at various temperatures were calculated from the enthalpy in the same report. These data sets are shown in Figure 3(a). An inconsistency in the data sets occurs in the range between 500 and 600 K, which corresponds to typical average fuel temperatures for TREAT transients. This problem can be mitigated in a variety of ways. One approach is adjusting the raw data to a single datum point and fitting a polynomial equation to the enthalphy [14]. Here, the specific heat capacity is just the derivative of the enthalpy as a function of temperature (C p = dh/dT ). This equation is referred to as a "cubic" representation in this paper. A second approach is obtained by averaging the points in the overlapping regions of Figure 3(a) and fitting a polynomial, a.k.a "sextic". We also include a standard graphite empirical fit [15,16]. A simple Rattlesnake point reactor kinetics (PRK) model is used to evaluate the peak power and total core energy deposition for various reactivity insertions using each heat capacity functionalization.

Biological Shield
Through the use of the PJFNK-SPH the steady state solutions from MAMMOTH are consistent with the Serpent results, but the initial transient simulations performed with MAMMOTH showed significant differences from the measured reactor period. These differences in the MAMMOTH model were initially attributed to the mesh, which does not include the biological shield and the thermal column. The contributions to the time-dependent fission source from these regions, via reflection of neutrons back into the active core, are not negligible and affect multiple key parameters for transient simulations. To better understand the effects of omitting these regions for steady state calculations, two Serpent models, one with the biological shield and thermal column and one without are compared. The same kinetic parameters for the MAMMOTH model are computed with Rattlesnake. Furthermore, a tally is added to the Serpent model to compute the albedo matrix on each side of the core and determine the benefits of using an albedo boundary condition in Rattlesnake.

Control Rods
The TREAT transient rods exhibit a high level of heterogeneity both radially and axially, which is visible in Figures 4(a) and 4(b). The model shown is based on a simplified 3x3 configuration with standard and transient control rod fuel elements. A fully homogenized mesh, like the one shown in Figure 4(c), has been used in the past modeling efforts and is henceforth referred to as model 1.
Unfortunately, the interpolation of SPH factors in homogenized control rod regions does not work well when the control rods are in between two tabulated points and can lead to large errors in the absorption and fission rates. This is mainly due to the fact that the number densities and hence the homogenized cross sections are incorrect as the rods change position. In addition, the use of SPH-diffusion with a more heterogeneous configuration also leads to similar problems in partially rodded elements, where the SPH factors are not tabulated. To circumvent these issues, the new modeling of control rods requires a heterogeneous mesh similar to the one shown in Figure 4(d).
Here the central poison and stainless steel regions define the dynamic part of the control rod, which change position during operation of the reactor. This mesh is normally extruded with the rest of the core mesh to define the various standard fuel element material interfaces in the axial direction. A low or high order transport solution with the multischeme method in this part of the domain can produce significantly better solutions than diffusion. The Rattlesnake rodded neutronics material modifies the compositions as the control rod is moved in and out of the core. The control rod cusping treatment must then ensure that the transport solution is smooth and represents the physical behavior of the rod movement well. This test is primarily concerned with the use of the SPH correction in a multischeme calculation with the partly heterogeneous model, referred to as model 2. The eigenvalue and the axial power distribution are evaluated. This work focuses on three fuel cross section preparation approaches in which temperature effects are also considered. In the first case, two fuel cross section regions are included, regular and borated, and where the temperature distribution is constant. In the second case, one cross section region is prepared for each fuel zone also with a constant temperature distribution. The third case,  The cross sections are condensed into two neutron energy group structures. The reference transient calculation is based on a 26 group German HTR structure [17]. The 10 group structure used for TREAT modeling and simulation is derived from this 26 group structure with 5 thermal, 3 epithermal and 2 fast neutron energy groups.

Spatial Homogenization and Energy Condensation
The cross sections obtained from the various spatial homogenization and energy condensation approaches lead to inconsistent kinetic behaviour for the same reactivity insertion. This is a consequence of the cross section sets themselves which produce different values of the effective delayed neutron fractions and mean generation time. To resolved this issue one computes kinetic parameters from each of the cross section sets using the IQS system in Rattlesnake. A simple solution to the Inhour equation provides the necessary reactivity insertion to match an inverse period of 14 s −1 . This way, all models have a consistent dynamic behavior. Solutions from the coarsely homogenized 10 group case are compared to high-resolution heterogeneous 26 group with S 20 .

Heat Capacity
The results from the sensitivity study for a 10 second transient with the Rattlesnake PKE model are presented in Table 1, where the empirical correlation serves as the reference. The empirical and sextic fits exhibit better consistency than the cubic in both cases, with the exception of the energy deposition in larger transient above 2%Δk/k in which the cubic performs slightly better for the energy deposition. The results suggest a 5% uncertainty on both peak power and energy deposition based on the heat capacity correlation alone. The average fuel temperatures for the transients range from 300.0 to 680.0 K.

Biological Shield
The repercussions from omitting the biological shield and thermal column on the kinetic parameters are presented in Table 2. The TREAT data are the values currently used in the Automatic Reactor Control System (ARCS). The MAMMOTH model is SPH corrected, and thus takes into account some of the effects of the biological shield and thermal column from the reference fluxes. Note that the measured value of the TREAT Λ is based on assumed values of the effective delayed neutron fractions and the neutron precursor decay constants [18]. The ratios ρ/Λ and β ef f /Λ dominate the prompt time constant and are useful in determining the consistency of the data, specifically, β ef f /Λ is directly correlated to the measured value from the 1960s [19]. The eigenvalue bias was quantified at +376 pcm. Another, less obvious, consequence appears from the reflection of neutrons from these external regions back into the active core, but with some additional time delay, which impacts the mean generation time.   Figure 5 shows the time-dependent power for the 1.5% Δk/k reactivity insertion with and without SPH correction and albedo boundary conditions. The non-SPH and SPH curves are different because they are using the same reactivity insertion, which leads to different periods. The results demonstrate that including a traditional albedo boundary condition has little effect on the non-SPH corrected model (<1.0% in total energy) and almost no effect on the SPH corrected one. An improved boundary condition would include the time-dependence (delay) of the neutrons in the peripheral regions before they re-enter the core, but the implementation might be cumbersome and with little gain in accuracy.  Table 3 presents the eigenvalues obtained for the two control rod models. The first model is fully homogenized, whereas the second uses multischeme SPH-diffusion with S 8 in the dynamic part of the control rod. In addition, the multischeme solver employs the control rod cusping treatment [20] with a third order projection. Both models reproduce the eigenvalue within the uncertainty of the reference calculation. The axial power distribution for each assembly type in model 2 is shown in Figure 6. The computed assembly powers are very accurate as well in each SPH zone, since the use of SPH-diffusion guarantees recovery of the reaction rates.

Homogenization and Condensation Effects
The reference solutions from the heterogeneous models are shown in Figure 7. The 26 group solution with 10 fuel zones and a temperature distribution is listed as the reference (26G-ref). The reference solution for 10 energy groups is also included (10G-ref). Table 4 lists various key transient parameters for the three cross sections sets with 26 energy groups. The results in the first two rows show that there are negligible spatial effects in the cross sections with a constant temperature distribution. The introduction of the temperature distribution leads to a lower peak power and energy deposition. This is mainly driven by differences in the leakages between the model with constant and distributed temperature. The 10 group reference overestimates the power peak and energy deposition by 2.8 and 0.8 %, respectively. Note that for all of these heterogeneous models the tail behavior is very similar and the main differences are in the peak and, consequently the FWHM of the power trace.   The solutions from the various homogenized SPH models in 10 energy groups are included in Figure 8 and the key transient parameters are shown in Table 5. Here again, the two solutions with a constant temperature distributions show no significant spatial dependence of the cross sections and SPH factors, but they lead to a stronger, but delayed feedback that produces a lower tail compared to the reference. The homogenized 10 group model with the temperature distribution is consistent with the reference solution in the energy deposition and average fuel temperature, but with a peak power that is 1.23% lower. The important observation is that the tail behavior matches that of the reference. Therefore, using the temperature distribution during the calculation of the reference fluxes for SPH is crucial in the recovery of the correct transient behavior. When compared to the 26 group heterogeneous reference this 10 group homogenized model with a temperature distribution slighly overpredicts all values, but in most TREAT transients the tail behavior accounts for > 15% of the total energy deposition. (b) log power

CONCLUSIONS
The results indicate that a 5% uncertainty in the peak power and core energy deposition are directly attributed to the specific heat capacity of the TREAT fuel. Including the biological shield and thermal column in models has a non-negligible effect on the reactor mean generation time. Matching the reactor period instead of the static reactivity insertion with the MAMMOTH reduced domain size resolves this problem. The neutron reflection from these external regions is marginally improved with the use of a traditional albedo boundary condition only on models without SPH equivalence. A time-dependent albedo could improve the mean generation time, but its implementation would be complicated and with little gain. The initial testing of the multi-scheme method with SPH for a transient control rod with cusping treatment shows that the reference eigenvalues and power distribution are reproduced very accurately. Future work will determine if the multi-scheme approach with SPH-diffusion and high order transport can be successfully used in predicting the control rod movement. The TREAT fuel cross sections show no spatial dependence and a small sensitivity to the temperature distribution. The transient behavior with SPH-diffusion suggests a strong dependence on the temperature distribution used in the preparation of the reference fluxes. This work shows, for a 1-D TREAT test, that the coarse, homogenized SPH-diffusion model with 10 energy groups can produce very similar results to a high-resolution, heterogeneous 26 group high-order model with the appropriate preparation of the reference fluxes for the SPH procedure. This insight might resolve the differences observed in the MAMMOTH predictions of the power tail behavior for TREAT transients.