A ROBUST SECOND - ORDER MULTIPLE BALANCE METHOD FOR TIME - DEPENDENT NEUTRON TRANSPORT SIMULATIONS

A second - order “Time - Dependent Multiple Balance” (TDMB) method for solving neutron transport problems is introduced and investigated. TDMB consists of solving two coupled equations: (i) the original balance equation (the transport equation integrated over a time step) and (ii) the “balance - like” auxiliary equation (an approximate neutron balance equation). Simple analysis shows that TDMB is second - order accurate and robust (unconditionally free from spurious oscillation). A source iteration (SI) method with diffusion synthetic acceleration (DSA) is formulated to solve these equations. A Fourier analysis reveals that the convergence rates of the proposed iteration schemes for TDMB are similar to those of the common (SI + DSA) schemes for Backward Euler (BE); however, TDMB requires about twice the computational effort per iteration. To demonstrate the theory—accuracy, robustness, and convergence rate—and investigate the efficiency of TDMB, we present results from a discrete ordinates (Sn) research code. Results are discussed, and future work is proposed.


INTRODUCTION
In solving the time-dependent neutron transport equation, time dependency is usually discretized, and an auxiliary equation is required to close the system of equations. Varying auxiliary equations are available, and each comes with different accuracy and stability. Widely-used standard methods include Forward Euler (FE), Crank-Nicholson (CN, a.k.a trapezoid rule), and Backward Euler (BE). FE is simple but inappropriate for typically stiff problems. CN is favorable because of its higher (second) order of accuracy. However, BE remains the favorite method in practice, due to its robustness; it is free from producing spuriously oscillating solutions regardless of time step size.
The Multiple Balance (MB) principle is introduced in [1]. A simple analysis (discussed later) shows that MB strictly produces non-negative solutions (thus is free from spurious oscillation) with second-order accuracy. Its applications to spatial ( [1]) and energy discretization of charged particle transport calculations have been studied [2], but the application to time discretization has yet to be investigated. MB is a potential second-order accurate alternative to BE while retaining robustness. Such a robust second-order method is expected to have higher computational complexity, both in memory and number of "solves". However, if it can be formulated efficiently enough, the benefit of higher accuracy may outweigh the additional computational cost per time step.
In this work, MB application to the time derivative of the neutron transport equation-hereafter referred to as TDMB-is investigated. The paper is organized as follows. In Section 2, the MB principle and its formulation for time derivative discretization are presented. In Section 3, a simple analysis revealing accuracy and stability of the method is performed. In Section 4, a discrete ordinate (Sn) source iteration (SI) scheme with TDMB is proposed, and its convergence rate is analyzed with a Fourier analysis (FA ). In Section 5, a diffusion synthetic acceleration (DSA) scheme is formulated to speed up the convergence. In Section 6, test problems are devised to verify the stability and accuracy and to investigate the efficiency of TDMB. Last, Section 7 concludes and discusses future work.

TIME DEPENDENT MULTIPLE BALANCE
Let us consider a time-dependent neutron transport problem: with initial condition (0) = init ; is the linear transport operator and the independent variables in ( ⃗, Ω, , ) are not entirely shown for simplicity. If we discretize time into time bins with an equal time step Δ , and then integrate Eq. (1) over each time bin, we obtain the balance equation for time bin with initial condition 1/2 = . Equation

ACCURACY AND STABILITY
Let us consider an infinite purely absorbing medium problem without internal source. Equation (1) thus becomes In this problem, we seek the solution +1/2 = ( +1/2 ) from a given −1/2 . The exact solution is where = Σ Δ . We recast the solution in term of the amplification factor ( ). Solving Eq. (4) with BE and TDMB give the following: where Exact ( ) − Method ( ) quantifies the error factor introduced per time step of the method. Equation (6) shows that BE and TDMB are respectively first-and second-order accurate. Additionally, the amplification factors indicate that both methods produce strictly positive solutions regardless of Δ , which justifies their robustness-being free from spurious oscillation. The amplification factor of BE may be negative in the case of Σ < 0 (exponentially growing solution, supercritical system), yet this is not the case for TDMB, which makes it robust regardless of . Nevertheless, Δ is typically very small in a supercritical system due to the accuracy and convergence requirement of the transport source iteration.

Source Iteration Scheme
It is not known yet whether there is an efficient strategy for solving the two coupled equations of TDMB. Several strategies have been investigated, but only one of them that is specifically formulated for neutron transport will be discussed in this paper. We consider a mono-energetic non-fission 1D-slab problem with isotropic scattering. The original balance and the balance-like equations are discretized in time and respectively give Initial and boundary conditions are not shown for simplicity. The superscripts () in Eq. (7) indicate the proposed source iteration (SI) scheme to solve the coupled equations. We cannot describe in detail the solution procedure of this unaccelerated SI (hereafter referred to as TDMB-SI), due to lack of space. The basic idea is to integrate Eq. (7) over a spatial cell, −1/2 ≤ ≤ +1/2 , to produce two balance equations; and then by applying discrete ordinates (Sn) method with weighted-diamond (WD) spatial differencing, two auxiliary equations are introduced: where subscript and respectively denote direction and spatial indexes, and we note that different WD parameters , and ̂, may be used; these yield a system of four discrete equations. For > 0, the unknowns are cell-average , , and +1/2, , , and cell-edge angular fluxes , , +1/2 and +1/2, , +1/2 as shown in Figure 1. At each spatial step in a forward sweep, we solve for these unknowns with the given incoming cell-edge angular fluxes , , −1/2 and +1/2, , −1/2 obtained from the previous sweep step or boundary condition. The process is reversed for the backward sweep. One can go over the algebra and effectively ends up with two sweep operations per spatial step, each for the timeaverage and time-edge + 1/2 quantities. We note that only one sweep operation, for the time-edge quantities, is performed in BE-SI. An important take away here is that the amount of work per iteration of TDMB-SI is about twice of BE-SI. Additionally, the amount of memory requirement is doubled since we need to keep track of the time-average quantities as well.

Fourier Analysis of SI
A Fourier analysis is performed to investigate the convergence rate of the proposed TDMB-SI scheme outlined in Eq. (7). This analysis is similar to that performed in [3], where an infinite homogeneous medium with isotropic scattering is considered. It is found that the iteration matrix eigenvalues 1,2 are a pair of complex conjugates (more on this in Section 6), and the resulting spectral radius can be expressed as shown in Eq. (9), where = Σ Δ and = Σ /Σ . We can replace with κ + (1 − ) for a system with fission, where criticality = Σ /Σ ( is used to avoid confusion with time bin index ). As a comparison, the spectral radius of BE-SI is provided as well in Eq. (9). TDMB-SI = | 1,2 | = | 2 + ± 2 + 2 + 2 | = √2 + 2 + 2 ; BE-SI = + 1 .
The spectral radius of an iteration matrix represents the worst-case estimate of error reduction factor per iteration; the smaller , the faster the convergence. The spectral radii, for = 1, of TDMB-SI and BE-SI are shown on the left of Figure 2. It is found that convergence is guaranteed by both methods regardless of Δ , because is always smaller than 1 ( → 1 as → ∞). This verifies the applicability of the proposed TDMB-SI scheme. Should increase, a smaller Δ is needed for convergence, and this generally applies to any method. TDMB-SI and BE-SI have very similar spectral radii, and thus convergence rates. Recalling from the previous subsection that the amount of work per source iteration of TDMB-SI is about twice that of BE-SI, we may conclude that TDMB-SI ultimately requires about twice the work of BE-SI in producing solution at each time step with a given Δ . Another important finding of the Fourier analysis is that, similar to BE-SI, the most slowly converging error mode of TDMB-SI happens to be the flat mode ( = 0, as shown on the right of Figure 2), which is nearly isotropic; this suggests that we can formulate a diffusion acceleration scheme, which works well in eliminating nearly isotropic mode, to speed up the convergence of the proposed TDMB-SI scheme.

Figure 2. (Color online) Spectral radius (left) and eigenvalues (right) of SI and DSA.
Eigenvalues of TDMB are complex conjugate pairs, and the curves shown on the right figure are their modulus.

DIFFUSION SYNTHETIC ACCELERATION
A diffusion synthetic acceleration (DSA) scheme (hereafter is referred to as TDMB-DSA), is formulated. Correction equations for the time-average and time-edge scalar fluxes, respectively and , are obtained by taking the zeroth and first angular moment, with 1 approximation, of the actual error equations based on Eq. (7). These correction equations are shown in Eq. (10) and (11), and the correction steps are shown in Eq. (12). The unknowns and are separately acted on typical diffusion operators but coupled within the same fixed source equations. Comparing to the typical diffusion operator arising in BE-DSA, here we have four blocks of the typical diffusion operator for solving the system of equations with doubled unknowns. We note that there are three kinds of total cross-sections, which are defined in Eq. (13).
[− 1 3Σ The TDMB-DSA scheme is an extension of TDMB-SI. We start by performing the transport sweep outlined in Eq. (7), but we change all the superscript iteration indexes of ( + 1) into ( + 1/2) to mark them as intermediate iteration solutions. After that, acceleration is performed: Eq. Fourier analysis of TDMB-DSA is performed. The iteration matrix eigenvalues are a pair of complex conjugates, similar to TDMB-SI. The resulting spectral radius and eigenvalues are compared in Figure 2.
It is found that TDMB-DSA successfully eliminates the flat mode ( = 0), which is the most persisting error mode in TDMB-SI. Furthermore, similar to the unaccelerated TDMB-SI and BE-SI, TDMB-DSA and BE-DSA have similar spectral radii, indicating that they take similar number of iterations for convergence. To conclude that TDMB-DSA ultimately requires about twice the work of BE-DSA, we need to compare how much more expensive is solving Eqs. (10) and (11). A simple convergence analysis with Gauss-Seidel-based iteration methods (not presented in this paper) shows that the resulting spectral radii of iteratively solving BE and TDMB lower-order diffusion equations are similar; yet again, TDMB requires about twice the work per iteration.

NUMERICAL RESULTS
The proposed TDMB-SI and TDMB-DSA iteration schemes have been implemented into an Sn research code. Diamond difference spatial discretization is used: , =̂, = 0 in Eq. (8). A representative problem of one-group homogeneous non-fission infinite medium with isotropic scattering, whose solution is shown in Eq. (15), is considered: A large slab ( 1000 cm , = 100 ) with vacuum-reflective boundary conditions is used in the Sn simulations to simulate the infinite medium problem. Total cross-section Σ and neutron speed are set to unity. Number of directional quadratures is set to 32. Additionally, to avoid false convergence during iteration, convergence criterion is set as shown in Eq. (16), where tol is the specified relative error tolerance, and the spectral radius is approximated as the ratio of the consecutive "residuals". Two test problems are presented. The first one is to verify the theoretical spectral radius-Eq. (9), and the second is to verify the stability and accuracy, and to investigate the efficiency of TDMB.

Spectral Radius
A zero-solution problem is considered: initial condition (0) and independent source are set to zero. Scattering ratio and are set to unity and 10.0, respectively. We use a step function in space as the first guess to feed enough Fourier error modes. Knowing that the exact solution is zero, we can numerically estimate the spectral radius from the ratio of the two consecutive exact residuals: ‖ ( +1) ‖ 2 /‖ ( ) ‖ 2 , either of or +1/2 . Selected results are shown in Figure 3. It is found that the ratio of residuals of TDMB-SI does not converge to the expected spectral radius obtained in Fourier analysis. It is observed that upon "convergence", the residual ratios keep changing periodically around the expected : they decay down to , rapidly drop below it, suddenly jump up once producing worse solutions, and then repeat. This anomaly is due to the leading eigenvalues 1,2 being a pair of complex conjugates, something that, to the authors' knowledge, has not appeared previously in any standard iteration methods. A numerical spectral radius estimation that better suits for complex conjugate leading eigenvalues needs to be formulated. TDMB-DSA also has a complex conjugate pair as its leading eigenvalues. However, its residual ratio converges close to . More investigation is needed to observe why TDMB-DSA behaves differently despite its complex conjugate leading eigenvalues. We note that it converges to slightly below the expected spectral radius; this is because is the theoretical spectral radius of a real infinite slab. In this second test problem, scattering ratio , independent source , and relative error tolerance tol are set to 0.5, 0.1, and 10 −10 , respectively. The solution "at" the reflective boundary is evaluated for time = 10 and is compared to the analytical solution Eq. (15) to calculate the error. The number of time steps is varied from 1 to 1000; we note that = Δ , since and Σ are unity. The relative error and total number of transport sweeps sweep (one and two per iteration for BE and TDMB, respectively) are recorded. As a figure of merit, efficiency eff is defined as shown in Eq. (17) and is calculated for each value . Results are shown in Figure 4. It is found (plot not presented in this paper) that TDMB is second-order accurate and unconditionally robust. From the left figure of Figure 4, it is shown that to achieve a certain accuracy, TDMB requires a smaller number of transport sweeps compared to BE; this separately applies for both SI and DSA versions. The figure on the right compares their efficiencies, and it is shown that TDMB is more efficient than BE. We note that the efficiency of TDMB grows better (in decreasing ) than that of BE due to TDMB's higher accuracy order.

CONCLUSION AND FUTURE WORK
Using a physics-based auxiliary equation, TDMB offers a robust second-order accurate timediscretization method. Fourier analysis reveals that the convergence rates of the proposed iteration schemes for TDMB are similar to those of the common (SI + DSA) schemes for BE. Despite the doubled amount of work required per transport iteration, TDMB can achieve similar accuracy with a larger time step compared to BE; this potentially leads to less total amount of work and thus increased efficiency. Simple test problems verify the robustness, second-order accuracy, and estimated relative computational cost of the proposed TDMB-SI and TDMB-DSA methods.
Future work includes formulating an appropriate spectral radius estimation for the TDMB methods due to the complex conjugate iteration matrix leading eigenvalues. Also, a Fourier analysis for discretized direction and space is needed to further characterize the convergence rate of TDMB iteration schemes. Investigation on how TDMB features hold in solving more realistic (higher-dimensional, anisotropic scattering, multigroup) problems needs to be done. The ultimate goal of this work is to formulate a robust second-order accurate loosely-coupled multiphysics method for reactor transient simulations. To achieve this, we need to apply TDMB to the other involved physics (e.g., heat equation and thermal-hydraulics), set up the TDMB scheme for solving the nonlinear problem, and formulate a robust second-order accurate operator splitting scheme.