EFFECTS OF FUEL TEMPERATURE-SHAPING FUNCTIONS ON XENON OSCILLATIONS

In coupled multiphysics simulations, single pin-averaged values are typically used to describe the temperature, power, and burnup within a given fuel pin. However, since xenon oscillations are largely driven by fuel temperatures, radially dependent quantities have been implemented in the Virtual Environment for Reactor Applications. These radial-shaping functions are based on Zernike polynomial expansions and allow information to pass effectively between codes with differing spatial meshes. This work examines the effects of radial fuel temperature-shaping functions on the behavior of axial xenon oscillations. A test problem was developed from full-core, multi-cycle depletions using as-built fuel data. The center 25 assemblies of the full-core case were used to test the radial-shaping function by inducing an axial xenon oscillation using an instantaneous control rod movement. The test case was run with and without the radial shapes, and each component was also run individually. Including the shaping functions significantly impacted the xenon oscillations for this problem; the magnitude and period of the oscillations were altered, and large pin power and soluble boron differences were observed. Testing each component individually showed that the radial fuel temperature-shaping function had the largest effect.


INTRODUCTION
Accurately predicting reactor behavior through modeling and simulation depends greatly on the ability to capture the effects of coupled sets of physics. The phenomenon of axial power oscillations due to xenon buildup, particularly near the end of a fuel cycle, is of interest to reactor designers and operators. Advanced modeling and simulation tools struggle to accurately predict xenon oscillations due to the complex multiphysics nature of the problem. The Consortium for Advanced Simulation of Light Water Reactors (CASL) is developing a suite of modeling capabilities called the Virtual Environment for Reactor Applications (VERA) [1]. This paper documents recent developments in VERA to increase the fidelity of coupled multiphysics problems.
To better predict the xenon oscillation behavior in coupled problems, three quantities of interest were examined: fuel temperature, power, and burnup distributions. These three items are used to create continuous radial functions that capture the radial shape of each. Section 2 describes the background information necessary to perform this work. Section 3 describes the various software packages that were used and the problem model that was used. Section 4 presents the results of this work and Section 5 summarizes the conclusions and describes work that is still ongoing.

BACKGROUND
Xenon oscillations can significantly affect operating reactor performance. When large and undamped, axial power oscillations can challenge design and trip limits for operating pressurized water reactors (PWRs). In CASL, it is important to model these axial power oscillations to better quantify the fuel duty during normal end of life operations and load-following scenarios. Avoidance can increase person time in power maneuver planning and, at worst, could result in reduced electrical output. To better capture these effects during multiphysics simulations, Zernike polynomials were used to describe the radial pin temperature, power, and burnup shapes. The following subsections provide the necessary background on these topics.

135
Xe has a very large thermal neutron absorption cross section, making it a significant isotope to track. 135 Xe is produced in an operating reactor as a directly formed fission product and as a daughter of 135 I decay, another fission product. 135 Xe is removed from the core by either absorbing a neutron and converting to 136 Xe or by undergoing beta decay to 135 Ce with a half-life of roughly 9-hours, which is on the same timescale as operations. Due to its large neutron absorption cross section, 135 Xe removal is primarily dictated by the local neutron flux in the reactor.
If a perturbation in the core causes a sudden shift or tilt in the local power distribution (e.g., by a control rod movement), then a spatial xenon oscillation can be induced. This is caused by the higher flux regions burning out the existing built-up 135 Xe faster than the low flux regions, leading to more growth in the power shift as the 135 Xe burns out. Conversely, the 135 I concentration builds up in the high flux region and decreases in the low flux region. However, once the 135 I builds up to sufficient levels, its subsequent decay into 135 Xe leads to a large decrease in the flux of the initially high flux region. This causes a power increase in the low flux region due to the lack of significant 135 Xe levels in this region, and results in an inverse of the initial tilted flux distribution. If left uncorrected, this cycle can continue to oscillate with a period of approximately 15-30 hours [2]. This oscillatory behavior is undesirable for reactors operating at full power because it can lead to increases in localized power-peaking factors that can challenge design and trip limits. Therefore, the ability to accurately model and predict this behavior is important for modern modeling and simulation tools.

Zernike Polynomials
Zernike polynomials were employed to better translate the radial pin shapes between coupled codes. Developed by Frits Zernike [3], Zernike polynomials are a sequence of polynomials that are orthogonal and continuous over the unit disc making them suitable for cylindrical fuel rods. There are even and odd Zernike polynomials. The even polynomials are defined as and the odd are defined as where ݉ is the angular frequency, ݊ is the radial order with ݊ ≥ ݉, ‫ݎ‬ is the radial distance on the unit disc 0 ≤ ‫ݎ‬ ≤ 1, ߠ is the azimuthal angle, and ܴ are the radial polynomials. The radial polynomials, which are dependent only on the radius, ‫,ݎ‬ are given by

METHODOLOGY AND MODELS
This section describes the software packages and the models used in this study. The multiphysics coupling in VERA is performed between the radiation transport solver, MPACT [4], and the thermal hydraulic solver, CTF [5] using a standard Picard iteration scheme. The models are based on the Byron Nuclear Generating Station Unit 1 and were constructed with help from the plant owner, Exelon Corporation, and the plant designer, Westinghouse.

MPACT
MPACT is one of the neutron transport codes being developed by CASL. MPACT uses a 2D/1D approach to solve the neutron transport problem. The problem is divided into a series of 2D axial planes whose axial transverse leakage is solved using a 1D axial calculation. Each axial plane is solved independently using the method of characteristics (MOC) across the explicitly modeled radial geometry.
Each fuel pin in the problem is discretized both radially and azimuthally. Fuel pins are typically discretized into three equal-volume rings radially and eight azimuthal regions. However, when performing coupled multiphysics simulations, MPACT only communicates single value, pin-averaged power and burnup shapes.
To better approximate the actual radial shape, the values for each radial ring were used to calculate the coefficients used in a Zernike polynomial representation of the radial distribution for both power and burnup. An example radial power shape for a fuel pin with three radial rings is shown in Fig. 1. A similar approach is also performed for the radial burnup shape. The order of the Zernike polynomial expansion is the same as the number of radial regions in the fuel pin. Therefore, the example shown in Fig. 1 uses a third order expansion. These continuous radial functions are then passed to CTF to be used for thermal hydraulic and fuel performance calculations.

CTF/CTFFuel
CTF is a subchannel thermal hydraulic simulation code that is being developed by CASL and is designed for light-water reactor (LWR) vessel analysis. CTF uses a two-fluid solution over fluid film, vapor, and liquid droplet flow fields to provide detailed subchannel results. Adding CTFFuel [6] allows higher fidelity multiphysics coupling. CTFFuel is a standalone fuel performance modeling package that uses the CTF fuel rod models and can be run independently. CTFFuel can model steady-state and transient LWR fuel performance. When coupled to MPACT, CTFFuel models fuel rods using a different spatial mesh. Therefore, the continuous Zernike radial power and burnup shapes from MPACT can be easily mapped onto the radial discretization used by CTFFuel.
In coupled simulations, CTFFuel communicated only single value, pin-averaged fuel temperatures to MPACT to update the cross sections on the transport mesh. Rather than volume averaging the fuel temperatures, coefficients from the radial temperatures can be calculated to create a continuous Zernike representation of the fuel temperature shape, which is similar to the approach taken for power and burnup. An example radial temperature shape for a fuel pin with three radial rings is shown in Fig. 2. As seen in Fig. 2., the Zernike shape preserved the fuel centerline temperature and the volume-averaged temperature. This continuous shape was then mapped onto the MPACT radial mesh, which yielded the Ringwise Radial Shape curve seen in Fig. 2

Model
The Byron Nuclear Generating Station Unit 1 is a Westinghouse 4-loop PWR consisting of 193 fuel assemblies. The models used in this work were generated from as-built data and included burnable poisons, control rod movements, and radial support structure details. Starting with cycle 17, a jump-in methodology was used: the once-and twice-burnt assemblies that were shuffled into the start of cycle 17 were approximated using single-assembly burnup calculations with reflective boundary conditions. These single assemblies were depleted until the assembly-average exposure matched that of the real assembly. The burned assemblies were properly shuffled into the correct locations for the start of cycle 17, along with fresh fuel. Cycle 17 was then depleted for a full cycle, after which the resulting assemblies were shuffled to start cycle 18. This process was repeated for cycles 18-21.
After the cycle 21 depletion, an artificial xenon oscillation was induced using a 25-assembly subset of the fullcore problem. The center 25 assemblies were shuffled into a new model with reflective radial boundary conditions. This reduced-size model was used due to constraints in available computing resources. At the center of this new model, a single control rod assembly was added, fully withdrawn. A diagram of this model is outlined in Fig. 3 where the colors represent a typical loading pattern for a 193-assembly core. To induce the xenon oscillation in this new problem, the control rod bank was partially inserted instantaneously at time ‫ݐ‬ = 0. A series of 1-hour depletion steps were carried out over the next 48 hours, immediately following the insertion. This xenon oscillation problem was performed multiple times, both with and without the Zernike radial coupling shapes enabled. The results from these studies are presented in the next section.

RESULTS
Using the methodologies described in the previous section, five different xenon oscillation cases were examined. The first was a control case, the no radial coupling case, with all the radial coupling shapes disabled, and it used only volume-average values throughout the simulation. The second case, the all radial coupling case, enabled all three radial-shaping functions simultaneously. The final three cases had the burnup (burnup only), power (power only), and temperature (temp only) radial coupling shapes enabled independently from each other to determine their individual impacts on the problem.
After the 48-hour depletions had finished, the power axial offset was calculated for each case as a function of time to show the effects on the axial xenon oscillation. The results for all five cases are shown together in Fig.  4 with important parameters highlighted in Table I.  As seen in Fig. 4, the combined effects of the radial-shaping functions are significant on the axial xenon oscillation. When comparing the all radial coupling case with the no radial coupling case, which relied exclusively on volume-averaged quantities, the axial offset increased by roughly 0.73% at the height of the oscillation when the radial shapes were enabled. Table I shows that the magnitude of the oscillation with the shaping functions increased due to the 1.62 decay ratio, whereas the magnitude for the no radial coupling case remained constant over the 48-hour depletion with a 1.00 decay ratio. Additionally, the period of the all radial coupling case was slightly shorter than the no radial coupling case by roughly 30 minutes. Examining each of the three coupling shapes independently shows that when the burnup only shape was enabled, it had a marginal effect on the overall magnitude of the xenon oscillation. The axial offset followed closely with the no radial coupling case, only shifted downward slightly. Also, the burnup only case was the only result that demonstrated a dampening of the oscillation. The power only case results show an increase in the axial offset by roughly 0.3% throughout the depletion. Although the period of the power only case matches the no radial coupling case exactly, its amplitude was undamped, and it began to increase slightly near the end of the run. The temperature only case results had the largest impact on the shape of the oscillation. The axial offset differs from the no radial coupling case by as much as 0.63%, and it also elongated the oscillation period by an hour. The xenon oscillation of the temperature only case was also undamped, and its magnitude grew after the first period.
Enabling the radial shapes increased the critical soluble boron concentration at the start of the transient by 15.6 ppm over the no radial coupling case. The difference in methodologies led to a normalized root mean square (RMS) pin power difference of 0.747% with a maximum difference of 1.932% at the start of the transient, which then propagated as the oscillation proceeded. A visualization of these pin power differences is shown in Fig. 5.  Figure 5 shows the expected axial tilt in the pin powers due to the axial offset shift introduced by the radial coupling shapes. The largest pin power differences occur near the top and bottom of the problem in high power assemblies. Since the xenon oscillations were conducted at the end of a cycle, the fuel in the middle of the core was more depleted, causing the power and the xenon to be more concentrated near the top and bottom. The corner assemblies do not show pin power differences as large as those along the edge because the edge assemblies are higher in power and, therefore, produce more xenon.

CONCLUSIONS
This work examined the effects of radial-shaping functions on xenon oscillations. Instead of using single volume-averaged values to represent quantities of interest, a continuous Zernike polynomial representation was used. Fuel temperatures, power, and burnup were the quantities of interest in this work. In the coupled multiphysics problems that were studied, MPACT translated the power and burnup shapes from its spatial mesh to create Zernike polynomials. CTF/CTFFuel used these continuous functions to map the values onto its computational mesh. The resulting fuel temperatures from CTF/CTFFuel were then translated back to MPACT using a continuous Zernike representation. Using results from models that were created using as-built reactor data, a 25-assembly test problem was created. Despite its reduced size, this model contained the burnup and power histories from its prior depletions. A xenon oscillation was induced in this problem by instantaneously and partially inserting a control rod bank in the center of the problem, and the oscillation was then allowed to propagate for 48 hours following the rod movement. This was performed with and without the radial coupling shapes enabled and with each shape individually.
This work demonstrated that incorporating radial temperature, power, and burnup shapes can have a significant impact on problems that exhibit xenon oscillations. Using these shapes decreased the period and increased the amplitude of the power axial offset caused by the xenon oscillation over the first 48 hours following the insertion. These changes also led to critical soluble boron concentration increases of 15.6 ppm and maximum pin power differences of 1.932%. Observing each radial shape individually showed that the burnup shape had the least impact on the behavior of the xenon oscillation, whereas the temperature shape had the greatest impact.
Work is ongoing to apply these radial shapes to the full-core cases themselves, with the goal of comparing them with measured plant data.