INLINE THERMAL AND XENON FEEDBACK ITERATIONS IN MONTE CARLO REACTOR CALCULATIONS

In this work, we describe a method for converging nonlinear feedback during the convergence of the neutron fission source in a Monte Carlo reactor simulation. This approach involves updating feedback physics during discard batches in the Monte Carlo simulation rather than fully (or partially) converging the neutronics prior to the nonlinear update. This approach is demonstrated for a single PWR pin with thermal feedback and with both thermal and xenon feedback. Converging these feedbacks inline with the fission source is shown to have the benefit of reducing numerical instability by effectively underrelaxing the tallied quantities in the Monte Carlo simulation and improving computational performance by converging feedback within (or near to) the number of discard batches required to converge the fission source even without any feedback.


INTRODUCTION
There has been significant progress in the development of Monte Carlo tools capable of performing large-scale reactor analysis over the past decade. Realistic three-dimensional problems such as the VERA and BEAVRS problems are routinely being solved with Monte Carlo tools. In addition to solving the neutron transport problem, Monte Carlo reactor analysis tools must also account for feedback mechanisms such as thermal feedback, xenon feedback, or core depletion. As the capability of the analysis tools themselves has advanced, work into the iterative methods used to perform feedback calculations has been an active area of research [1]. Research in this area is motivated both by observed numerical instabilities in power and temperature distributions when solving coupled Monte Carlo and thermal-hydraulic (T-H) problems and the desire to reduce the computational expense of these coupled calculations.
Here, we will consider the further development and testing of an iteration method described in a previous article [1]. Specifically, we perform a number of numerical studies considering the simultaneous convergence of the fission source and thermal feedback (temperature and density). The original motivation for this method was to observe the effects on numerical stability of thermal updates that were performed early and often during the iterative process. Results showed that this was an effective approach to alleviating numerically induced oscillations. Beyond preventing numerical instability, the simultaneous convergence of the fission source and thermal feedback has the potential advantage of potentially fully converging thermal feedback during discard cycles, effectively making the addition of the coupled physics "free" given that the fission source would need to be converged even without feedback. In fact, equilibrium xenon has been calculated in this exact manner for a number of years [2] in MC21. In this work we will not only look at performing thermal feedback iterations in this manner, which we will refer to as "inline" from now on, but consider for the first time the simultaneous convergence of thermal feedback, equilibrium xenon, and the fission source distribution.
After providing a brief description on the specifics of the implementation of this iterative strategy the remainder of this work will focus on the presentation of numerical results. For the case of only thermal feedback, the ability of the inline approach to prevent numerical oscillations will be demonstrated and then a sensitivity study will be performed where the convergence of the coupled problem is assessed relative to a numerically calculated reference solution while input parameters to the iterative process are varied. A similar sensitivity study will then be performed with the thermal and xenon distributions being simultaneously converged with additional consideration given to the relative frequencies of the thermal and xenon updates.

INLINE FEEDBACK ITERATIONS
The most common approach to solving coupled Monte Carlo / T-H problems is simple fixed-point iteration, often called Picard iteration in the case that H is nonlinear [1]. In this equation the fixed point mapping H contains all of our physics (neutron transport, thermal-hydraulic equations) and x n the estimate of the neutron fluxes, multiplication factor, coolant temperature / density, and fuel temperature. A natural iterative strategy would then be to solve the transport problem given some initial temperature distribution, use the computed local powers to solve the T-H equations, producing new estimates of the temperature / density distribution. This iteration proceeds until some convergence criteria are satisfied. This scheme has been observed to be prone to numerical instabilities for certain problems leading to the use of an underrelaxation factor as follows, Sometimes underrelaxation is applied to only the powers from the transport solver or the temperature distribution from the T-H solver. In this iterative strategy the individual transport and T-H equations may be fully or partially converged. Partial convergence can be ramped up as the fixed-point iterations proceed using some heuristic in an attempt to reduce the overall cost of the nonlinear solve. The fixed point mapping may also include additional physics such as xenon/iodine rate equations in the case of equilibrium xenon.
The inline iterative approach we are considering in this work is fundamentally different than even a loose partial convergence because we are exchanging information prior to the convergence of the fission source itself. The mechanics of the iteration are nearly identical, information is simply being exchanged more frequently. Assuming the most recent estimate of the fission source is used each batch, and the feedback updates are being done prior to the convergence of the fission source at the current thermal (or thermal/xenon) conditions, then the tallies used in the feedback calculations will be the result of some combination of the current and previous nonlinear iterations' fission sources. This combination can be thought of as analogous to the combination described in Eq. (2) though rather than using the well defined ω parameter the amount of underrelaxation is controlled by how well the new fission source is converged between feedback updates. This can be controlled by the choice of the number of batches between updates, though the effect of adjusting the number of batches will depend heavily on the specific problem being solved via the dominance ratio.
While the ability of the inline iterative strategy to dampen numerical oscillations is easily shown numerically there is no supporting analysis to justify or quantity the effect of this so-called physical underrelaxation. Interestingly, recent work ([3] [4]) in deterministic methods has analytically shown that for coupled systems with nonlinear feedback using the CMFD low-order acceleration method fully converging the fission source of the low-order equations is detrimental to the overall convergence rate. The relationship between partial convergence and relaxation is examined specifically in [4]. This is indeed analogous to the method we are employing as partial convergence in the deterministic world means both partial convergence of the fission source and partial convergence of the neutronic result (e.g. power). These cannot be partially converged separately as with Monte Carlo and thus when these works talk of partial convergence they are talking about the scheme which we are currently investigating.
The implementation of the inline iteration approach is very straightforward. In this work we will be considering thermal feedback and xenon feedback specifically. The focus is not on the underlying physics models themselves but in the convergence of these models. The thermal feedback model [5] is a simplistic model based on conservation of energy and the equilibrium and xenon feedback [2] is a simple two-nuclide system of rate equations. These models and the inline iteration approach are implemented in the MC21 Monte Carlo code [6]. Specification of the iteration strategy for each feedback type accepts two input parameters: 1) the number of batches to perform between feedback iterations, α, and 2) the number of feedback iterations to perform before accumulating tallies, β.
Tallies are accumulated over α batches, used to update the thermal / xenon distribution, and then flushed until this has been done for β iterations at which point the tallies are accumulated for the remainder of the calculation (continuing to update the feedback every α batches).
This strategy is extremely flexible in that it allows the iteration strategy to vary with each feedback mechanism being employed and it provides control over precisely when accumulation should begin. Similar to discard batches during the convergence of the fission source, early feedback iterations should also be discarded to avoid contaminating the final results. A recent study by Kreher [7] considers an implementation of the inline iteration strategy where thermal feedback only is updated each batch and tallies can be accumulated over moving windows of batches. This alternate implementation of the strategy is complementary to the studies performed in this paper, providing strong evidence that there are significant advantages to such an iterative strategy.

NUMERICAL RESULTS
Numerical results are presented for the same single-pin PWR unit cell used in [1]. feedback. Nuclear data was generated at multiple temperatures, to be used for interpolation in MC21, using ENDF-B-VIII with 10 K spacing between fuel nuclides. The system pressure is 2215 psi, the power is 100 kW, the mass flow rate is 0.4 kg/s and the inlet temperature is 515 K. The clad resistivity is 0.81431 (cm 2 ·K)/W and the fuel thermal conductivity is calculated using a temperature dependent UO 2 correlation from Todreas and Kazimi [8]. This problem is identical to the problem used frequently for performing studies in [1].

Reference Solutions
To assess the convergence behavior of the inline method we generated numerical reference solutions for thermal feedback only (TF) and for the thermal and xenon feedback (TF-XE) case. These solutions were generated using 20 Picard iterations for the thermal feedback with an underrelaxation factor of 0.7 applied to local power passed to the thermal solver and the inline approach for xenon feedback with α = 10 and β = 25. Underrelaxation is applied because the TF variation of this problem has previously [1] been shown to be unstable. Each batch of neutrons consisted of 1 million histories with 250 discard and 500 active batches. Each problem was run 10 times using different random number seeds and the results of these independent simulations were then combined to form the reference solution. The final power, fuel temperature, and coolant density for each reference solution is shown in Figure 1 along with the convergence of the thermal feedback in each case. Convergence is plotted by showing the maximum relative iteration-to-iteration error. The uncertainty in these results based on running 10 independent simulations is smaller than the line thicknesses in the plot and the relative uncertainties are on the order of 10 −3 for power and fuel temperature and 10 −4 for coolant density. Observations of the Shannon Entropy for this problem with no feedback indicate that 250 discard batches is sufficient to converge the fission source.

Inline Iteration Convergence and Sensitivity to Running Strategy
First we consider the TF problem (no xenon). We are interested in understanding whether the inline approach converges to our reference solution, how quickly it converges, and whether convergence is sensitive to the values of α and β that we choose. As in the reference we run 250 discard batches and 1 million neutrons per batch. We ran 2000 active batches in order to observe the longer term convergence behavior. We consider α values of 5, 10, 25, and 50 and β values that result in feedback accumulation beginning at batches -200, -100, 0, 100, 200. All of the strategies tested were shown to result in stable iterative schemes. Figure 2 shows results from a subset of the running strategies examined with the legend showing α and the batch at which tally accumulation begins. The top row of plots shows the maximum relative error over all axial bins with respect to the reference solution at each batch. The bottom row shows the maximum relative error over all axial bins between two successive iterations. The comparison to the reference solution shows that all strategies converge towards the reference solution at a similar rate, regardless of the values of α and β. We see that the maximum relative deviation from the reference flattens to approximately the same level as the relative uncertainty in the reference solution itself, showing that the power shapes for the inline strategies indeed converge towards the correct solution. The iteration errors show that a smaller α value (fewer batches between updates) results in smaller changes between iterations, as we expect. Smaller values of α mean smaller steps will be taken (more underrelaxation). We note that the magnitude of the iteration error is not an indication of the convergence with respect to the reference solution and that the error with respect to the reference flattens even before the iteration error would indicate that error has stopped decreasing. In any case, these results show that convergence to the correct solution is independent of the particular running strategy shown and that the error with respect to the reference solution has nearly flattened by the time active batches begin, showing that the thermal solution can be converged simultaneously with the fission source in which case the only additional computational cost over a calculation without feedback is the cost of the thermal solve itself, which is negligible when using the MC21 model.
In Figure 3 we look at the spatial distribution of the error with respect to the reference solution for each of the inline running strategies considered. The shaded regions represent the error including the combined uncertainty in the reference solution and the tallied power. Generally good agreement is observed though the α = 5 results show a clear tilt with respect to the reference. Because of the normalization to power level (100 kW) any over-prediction in flux at some axial location must be accompanied by an under-prediction at a different location and so the fact that the plots appear centered around 0.0 is not a sign of convergence but a result of the normalization. Observing the shape of the error as the iterative process evolves shows that the error "moves" as the calculation progresses but does not oscillate about 0.0 and instead seems to converge to an error with some shape that may be within the expected uncertainty but is often not. More work is necessary to determine if this is due to a systematic difference introduced by either the method used to generate the reference solution, the inline strategy itself, or a symptom of not enough neutron histories / feedback updates.
Xenon feedback was added to the problem using the α = 10 and β = 25 strategy for xenon used to generate the reference solution. The same thermal feedback strategies were tested as in the thermal feedback only problem, effectively creating iteration strategies where xenon was updated  less frequently (thermal α = 25, 50), as frequently as (thermal α = 10), and more frequently (thermal α = 5) than thermal feedback. The iterations were not staggered at all such that for thermal α = 10 both xenon and thermal were being updated after the same batches. Figure 4 again shows the maximum relative difference from the reference and the maximum iteration-toiteration error. Results are very similar to the TF case -the inline iteration approach appears to be insensitive to the particular running strategy used and the expected level of convergence to the reference solution is observed.

CONCLUSIONS
We have shown that both thermal and xenon feedback can be converged in Monte Carlo simulations simultaneously with the fission source providing several key benefits. This approach has effectively underrelaxed the feedback iterations (specifically power and xenon/iodine reaction rates), demonstrably improving numerical stability for a problem with only thermal feedback known to be unstable using standard Picard iteration. This approach also speeds up the overall simulation time relative to Picard iteration because it allows for converging feedback within (or near to) the number of discard batches required to converge the fission source even without any feedback. A sensitivity study was performed considering the number of batches between feedback updates and the tally flushing strategy and convergence was found to be insensitive to the particular parameters chosen, indicating the general robustness of this approach. The iterative strategy employed in this work considered xenon and thermal feedback but any physics that would fit into a steady-stated Picard iteration framework could be handled in this same manner.