Predictor-Corrector Quasi-Static Method for Tightly-Coupled Reactor Multiphysics Transient Calculations

Several methods have recently been developed to solve multiphysics transients using the Improved Quasi-Static method and its derivatives. In order to address some perceived drawbacks of these methods, we have developed a new method for solving multiphysics transient calculations using the Predictor-Corrector Quasi-Static method. Our method involves computing reactivity feedback parameters during each transport timestep in order to enable reactivity feedback on the small timescale used by the point kinetics phase of the calculation. The advantage of this approach is that the transport solver does not need to store local flux information between time steps, potentially making it more appropriate for use with a Monte Carlo solver. We have demonstrated the accuracy of our method by solving a simple model problem that exhibits difficult multiphysics behavior. Additionally, we have compared our results against another recently published multiphysics coupling scheme, confirming that our approach does not negatively affect the accuracy of the transient solution.


INTRODUCTION
In recent years there has been a desire to improve transient neutronics calculations. This is partially driven by increases in computing power and an increased interest in complex multiphysics calculations, which are often transient in nature. In the past decade, a number of methods have been developed to quickly solve reactor transients. In particular, the Improved Quasi-Static (IQS) method of Ott and Meneley [1] has seen renewed interest, with the Predictor-Corrector Quasi-Static (PCQS) method of Dulla et al. [2] being a notable example. The IQS method and its derivatives are based around a factorization of neutron flux into an 'amplitude' and 'shape' component. Following this factorization, the amplitude function is determined using point kinetics relations, and the shape component is calculated using space-dependent transport theory in a manner that is consistent physically between the two solution methods. The advantage of this approach is that the amplitude and shape functions can be treated with separate time scales, allowing high accuracy with only limited use of expensive transport calculations.
The manner in which the amplitude and shape functions are related to each other helps to define the difference between the original IQS method and some of its modern derivatives. In the original IQS formulation, shape and amplitude were connected via a nonlinear coupling. In implementation, this requires an iteration over each transport timestep in order to enforce consistency between point kinetics and transport. The PCQS method is a linearization of the IQS method, allowing a marching solution in time at some expense of accuracy. In the past decade, both the linear and nonlinear variants of IQS have been implemented in a number of code systems, and a number of efforts have been undertaken to quantify their performance. In particular, efforts to apply these methods to multiphyics calculations have had a lot of recent interest.
Two recent papers have presented strategies for coupling IQS-based methods with other physics feedback. In 2017, Patricot et al. [3] investigated the use of these methods for a thermo-mechanical feedback problem representing the Godiva core. They identified three potential coupling strategies for IQS that significantly improved calculation results in their simulation. All three of their coupling strategies are centered around methods to estimate reactivity feedback without performing more transport calculations. They found that these corrections were necessary because the point kinetics calculations were otherwise unable to capture feedback effects due to the thermo-mechanical system quickly enough to correctly simulate the core behavior. As part of this study, they observed that their coupling strategies were inappropriate for PCQS (linear) transient methods because their coupling requires some iteration on each time step to inform the point kinetics method of the reactivity feedback behavior of the system.
A 2019 paper by Prince and Ragusa [4] developed a coupling scheme for temperature feedback in a reactor transient. Their feedback scheme involves three timestep sizes (for point kinetics, thermal feedback, and neutron transport) and allows for the point kinetics parameters to be recalculated during the thermal feedback stage by updating cross sections without necessarily updating the flux shape (i.e., without repeating a transport calculation). They provided results and analysis showing the effectiveness of their scheme for both the IQS and PCQS methods for two reactor benchmark problems with Doppler feedback. One potential downside of this approach is that it requires some way to store local flux information in order to recompute the integrals for reactivity feedback during the point kinetics calculation.
In this paper, we develop another scheme for solving transients with temperature feedback by using the PCQS method, and we demonstrate its use on a contrived model problem with extreme temperature feedback behavior. Our method is similar in structure to the method posed by Prince, allowing a third timestep size for thermal feedback and updating point kinetics parameters each time the thermal physics are updated. It is unique from the method proposed by Prince because it does not require re-computing the integrals for point kinetics parameters at each thermal timestep. Instead, it calculates spatially dependent temperature feedback coefficients during each transport timestep and uses these coefficients to perform reactivity feedback during the point kinetics phase. The advantage of this approach is that the storage of local flux data from the transport calculation is not required, which can be advantageous for transport solution methods that do not inherently store local flux data, such as Monte Carlo simulations. Finally, we provide results to show that our method can successfully solve our example problem and that it has similar accuracy to the method of Prince. A realistic implementation with a Monte Carlo solver is a logical next step, but it is outside of the scope of this paper.

METHOD
The method developed in this paper is a coupling scheme for multiphysics calculations with IQS or PCQS. Unfortunately, there is no widely agreed upon notation for multiphysics problems of this sort. As such, in this section we will provide a brief re-derivation of the IQS and PCQS methods and the relevant physics that our coupling scheme will capture, as well as comparisons with the method developed by Prince.

Background on Predictor-Corrector Quasi-Static Method
Our review of the PCQS method follows closely the derivation given by Dulla et al. [2]. The PCQS and IQS methods differ only in their algorithmic implementation (not in their derivation), so this brief review captures the background for both. However, for brevity, the remainder of our paper focuses only on the PCQS method. This sheds light on our expectation that an eventual realistic implementation of the method with a Monte Carlo solver will rely on the linearity of PCQS. The derivation begins with a factorization of the angular flux into an amplitude and shape component, Here, ψ refers to the angular flux, A to the amplitude, and φ to the flux shape. Note that both A and φ are functions of time, meaning that no approximation has been made in this factorization. Next, consider the time-dependent transport and precursor equations without an external source, For brevity, most function parameters have been omitted (i.e., ψ = ψ(x, E, Ω, t)). These equations assume an unspecified number of delayed neutron groups, indexed with i. The operators H, F p and F d i are respectively the transport, prompt fission, and delayed fission operators, defined as We assume that the initial flux distribution ψ 0 and the initial adjoint flux solution ψ * 0 are known and are given by the solution to steady-state systems, In order to make the assumption of Eq. (1) unique, a normalization condition is assumed. For this paper, and as done by Dulla, flux shape is normalized against the initial adjoint distribution, ψ * 0 , to an arbitrary constant z, with ψ * 0 , φ(t) ≡ z .
After some manipulation, Eqs. (1) -(6) can be used to develop a shape equation and point kinetics equations, The parameters of Eq. (8) are defined, and where ·, · represents an inner product over space, energy, and angle. Equation (7) is not used in the PCQS method, but we have included it so that our method can be extended to the IQS method in general. Except for notation, these equations are identical to those of Dulla et al. [2].

PCQS algorithm
The intent of quasi-static methods for transients is to be able to follow a transient while minimizing the number of full transport calculations that must be performed. To this end, we define two timestep sizes, ∆t for transport, and δt for point kinetics, where δt << ∆t. The PCQS method is as follows: 1. Starting from a known steady-state forward and adjoint solution, evaluate Eq. (2) using the timestep ∆t to obtain an approximation for ψ(∆t).
2. Factor this approximation for flux into its shape φ(∆t) and amplitude using Eq. (6) and apply the flux shape to calculate the point kinetics parameters using Eqs. (9) -(13).
3. Solve the point kinetics equations, Eqs. (8) starting from t = 0 and advancing in steps of δt until t = ∆t, ultimately obtaining a new value for amplitude A(∆t) different from that calculated in step 2.

Use this value of A(∆t)
to correct the flux approximation calculated in step 1 using From this, the precursor distributions C i are corrected using the precursor formula in Eq. (2).

Multiphysics Coupling Scheme
Now we allow temperature feedback in the problem. This requires the introduction of some (general) thermal feedback model that provides temperatures as a function of power distribution, and the allowance that macroscopic cross sections are functions of temperature σ x = σ x (x, t, T ), where x is cross section type and P (x, t) is power distribution.
While Eq. (15) is general, it is clear that it too will require a discretization in time. We allow an intermediate timestep size for the thermal model, ∂t such that δt < ∂t < ∆t. Additionally, the thermal model requires a power distribution P (x, t) that can be evaluated on the time scale of the thermal solver. To accomodate this, we follow the lead of Prince [4] by using the flux amplitude from the point kinetics solve and by interpolating the flux shape linearly between φ(0) and φ(∆t).
We have now introduced enough notation to describe our coupling scheme. As observed by Patricot [3], the key to accurate multiphysics simulations with quasi-static methods is allowing feedback on the point kinetics reactivity term defined in Eq. (12) at the time scale of the thermal feedback. Patricot discussed several methods of enabling reactivity feedback on the ∂t time scale by taking advantage of the iteration over the transport time scale that occurs in the IQS method. In the PCQS framework, however, no such iteration exists (each transport timestep occurs only once), and so these powerful schemes can't be used. The scheme used by Prince cleverly sidesteps this problem by re-computing the integral in Eq. (12) at intermediate timesteps by using updated cross sections and an interpolated flux shape, also with success.
Our method instead involves a direct calculation of temperature coefficients of reactivity α T (x, ∆t) during the transport calculation. Then, these reactivity coefficients are used to update the reactivity term at the time scale of the temperature feedback. A flowchart describing this process for a single transport timestep (from t = 0 to t = ∆t) is given in Fig. 1.

RESULTS
The intent of these results is to demonstrate that our coupling scheme can be used to calculate accurate results to a transient calculation with thermal feedback. We do this by solving a simple model using both our coupling method and the method presented by Prince for comparison. We forego a proper analysis of computation time because our test problem is too simple and because a full-core implementation is outside of the scope of this paper. However, it should be noted that direct comparison between our method and that of Prince is not fair. With each transport iteration, our method requires ∆t ∂t fewer evaluations of the integral in Eq. (12) than the method of Prince, but it requires an evaluation of α T (x), which can be very expensive. This difference suggests that ψ(t = 0), T (t = 0) Transport Timestep and α T calculation Thermal Timestep Adjust reactivity with α T our method is inappropriate for use with deterministic transport solvers, but that it may be ideal for use with solvers (such as Monte Carlo solvers) that cannot easily store the information required to re-evaluate the reactivity integral between transport iterations.

Contrived Model
We demonstrate our method using a contrived model with exceptionally difficult feedback behavior. Our model problem is a two-material slab problem in 1D with temperature-dependent cross sections, two temperature regions, one delayed group, one energy group, isotropic scatter, and vacuum boundary conditions. Table I lists the unitless problem parameters. For thermal feedback, each temperature region i is modeled as a well-mixed constant-volume holdup tank such that The steady-state is solved with T in = 1 and then the transient begins by abruptly changing T in to 0.1 at time t = 0.
With the parameters in Table I, the left material has a negative temperature coefficient while the Because we are not interested in the performance of the point kinetics solver nor the barelyphysical thermal model, we fix their time step sizes to small values such that 100δt = ∂t, and 100∂t = ∆t and use an implicit Euler discretization for each. For all of our transport solutions, we use discrete ordinates with 8 angles and 200 spatial meshes. Our reference case is a transport calculation with ∆t = 2 × 10 −5 . Our quantity of interest for accuracy is chosen to be the time-integrated error in the temperature distribution, defined Table II summarizes our results for a variety of time step sizes ∆t when using our method on this problem. Additionally, the results of using the method presented by Prince are presented for comparison. As evidenced by the results of Table II, our method has similar accuracy to the method of Prince.

CONCLUSIONS
We have developed a new scheme for solving multiphysics transients using PCQS. Our scheme involves the calculation of reactivity feedback parameters during each transport calculation of PCQS to enable accurate reactivity feedback during the point kinetics portion of the calculations. We have demonstrated that this scheme can obtain accurate results by testing it on a simple multiphysics model and comparing its accuracy with other existing coupling schemes. Because it requires the calculation of spatially dependent reactivity coefficients during each transport calculation, our coupling method is unlikely to be a desirable way to solve transients using a deterministic transport solver. But because it does not require the storage of local flux information, it may be more appropriate than other coupling methods for Monte Carlo implementations of PCQS. Implementation of the method with a realistic full-core Monte Carlo solver remains as future work.