A cell-local finite difference discretization of the low-order quasidiffusion equations for neutral particle transport on unstructured quadrilateral meshes

Abstract

We present a quasidiffusion (QD) method for solving neutral particle transport problems in Cartesian XY geometry on unstructured quadrilateral meshes, including local refinement capability. Neutral particle transport problems are central to many applications including nuclear reactor design, radiation safety, astrophysics, medical imaging, radiotherapy, nuclear fuel transport/storage, shielding design, and oil well-logging. The primary development is a new discretization of the low-order QD (LOQD) equations based on cell-local finite differences. The accuracy of the LOQD equations depends on proper calculation of special non-linear QD (Eddington) factors from a transport solution. In order to completely define the new QD method, a proper discretization of the transport problem is also presented. The transport equation is discretized by a conservative method of short characteristics with a novel linear approximation of the scattering source term and monotonic, parabolic representation of the angular flux on incoming faces. Analytic and numerical tests are used to test the accuracy and spatial convergence of the non-linear method. All tests exhibit $O(h^2)$ convergence of the scalar flux on orthogonal, random, and multi-level meshes.

Introduction

The fundamental transport problem is the linear Boltzmann (transport) equation,²

$$\frac{1}{v(E)}\frac{\partial \psi }{\partial t}+\overrightarrow{\Omega }\cdot \overrightarrow{\mathrm{\nabla }}\psi +{\sigma }_t\psi =\underset{4\pi }{\int }d{\Omega }^{\prime }\underset{0}{\overset{\infty }{\int }}d{E}^{\prime }\sigma ({\overrightarrow{\Omega }}^{\prime }\to \overrightarrow{\Omega },E^{\prime }\to E)\psi ({\overrightarrow{\Omega }}^{\prime },E^{\prime })+q(\overrightarrow{r},\overrightarrow{\Omega },E,t),$$

shown in nuclear engineering nomenclature, with the primary unknown, the angular flux, $\psi (\overrightarrow{r},\overrightarrow{\Omega },E,t)=v(E)n(\overrightarrow{r},\overrightarrow{\Omega },E,t)$, a convenient combination of particle speed, v ($\text{cm}/\text{s}$), and particle density, n ($1/{\text{cm}}^3\text{eV}\text{ster}$). The total cross section, ${\sigma }_t={\sigma }_t(\overrightarrow{r},E)$ has units (${\text{cm}}^2$) and the double differential cross section $\sigma =\sigma (\overrightarrow{r},{\overrightarrow{\Omega }}^{\prime }\to \overrightarrow{\Omega },E^{\prime }\to E)$ has units (${\text{cm}}^2/\text{eV}\text{ster}$). The transport equation is seven-dimensional with three space, $\overrightarrow{r}$, two angle, $\overrightarrow{\Omega }$, one energy, E, and one time, t, independent variables. Conservation, positivity, and a proper diffusion limit are properties typically sought after in numerical solutions of the transport equation, with specific applications setting the priorities. Solving transport problems is difficult because of the high dimensionality of the problem, and because the particle transport phenomenon may exhibit different behaviors depending on properties of the physical system. In purely absorbing or void regions of a physical system, the transport equation behaves like a hyperbolic wave equation; in highly scattering regions it behaves like an elliptic (steady-state) or parabolic (time-dependent) diffusion equation. It is extremely difficult to find discretization methods that are accurate over this wide range of behavior.

Deterministic transport methods that discretize Eq. (1), as opposed to Monte Carlo methods, typically assume a discrete ordinates approximation in angle, where the angular flux is evaluated at discrete directions, ${\psi }_m=\psi ({\overrightarrow{\Omega }}_m)$, and angular integration is approximated by quadrature [1]. This directly leads to the “source iteration” (SI) iterative solution method, where a “transport sweep” is performed for each direction ${\overrightarrow{\Omega }}_m$ with iteration-lagged right hand side (RHS) scattering terms. For optically-thick, scattering-dominated (i.e. “diffusive”) problems, the number of iterations required by SI to converge approaches infinity. Thus transport methods intended for diffusive applications must address “acceleration” in some sense, currently achieved through pure acceleration methods such as diffusion synthetic acceleration (DSA), which effectively damps the slow-converging diffusive modes of the solution, and/or Krylov iterative methods.

In many applications, low-order approximations of the transport equation may be used to obtain solutions within acceptable margins of accuracy. Low-order approximations in widespread use are diffusion, $P_N$, simplified $P_N$ (${\mathit{SP}}_N$), and variable Eddington factor (VEF) methods. The approximation in $P_N$ and ${\mathit{SP}}_N$ is related to the order N, which for low-order approximations are typically $N=1$, 3, or 5.

The quasidiffusion (QD) method [2] for solving the transport equation is unique in that it addresses both acceleration of SI of the transport equation and the need for a low-order problem of particle transport. The low-order model, namely the low-order QD (LOQD) equations, contains non-linear tensor factors (similar to VEF [3]) which can capture any transport physics. The system of nonlinear equations that define the QD method is equivalent to the original linear transport equation. This work presents a new discretization of the LOQD equations for arbitrary quadrilateral meshes in Cartesian XY geometry, including multi-level meshes, such as those that arise in adaptive mesh refinement (AMR) applications. Examples of meshes used in this work are shown in Fig. 1. The randomized meshes considered simply serve as difficult tests of the QD discretization. Depending on the application, arbitrarily rotated/skewed quadrilaterals could be used to follow material interfaces not aligned with the coordinate axes or better approximate curved surfaces. Randomized meshes also approximate meshes that could result from radiation-hydrodynamics simulations where Lagrangian meshes follow material motions and skew to preserve mass in cells. Multi-level meshes provide a way to easily refine important regions of the problem without increasing the number of unknowns in unimportant regions.

The QD method consists of two basic parts: (i) the transport equation (high-order problem) and (ii) the low-order QD (LOQD) equations. For steady-state, one-group transport problems with isotropic scattering and source in two spatial dimensions, the system of QD equations are given as follows.

High-order transport problem

$$\overrightarrow{\Omega }\cdot \overrightarrow{\mathrm{\nabla }}\psi +{\sigma }_t\psi =\frac{{\sigma }_s}{4\pi }\varphi +q_0\text{in domain}G,$$

$$\psi ={\psi }_b\text{for incoming}\overrightarrow{\Omega }\cdot \overrightarrow{n}<0\text{on boundary}\partial G,$$

Low-order quasidiffusion (LOQD) problem

$$\overrightarrow{\mathrm{\nabla }}\cdot \overrightarrow{J}+({\sigma }_t-{\sigma }_s)\varphi =q_0\text{in domain}G,$$

$$\sum _{\beta =x,y}\frac{\partial }{\partial \beta }(E_{\alpha \beta }\varphi )+{\sigma }_t{J}_{\alpha }=q_{\alpha }\text{for}\alpha =x,y\text{in domain}G,$$

$$\overrightarrow{n}\cdot \overrightarrow{J}+J_{in}=C(\varphi -{\varphi }_{in})\text{on boundary}\partial G,$$

where G is the domain of the problem, ∂G is the boundary surface of G, $\overrightarrow{n}$ is the outward normal to ∂G, $\varphi ={\int }_{4\pi }\psi d\Omega$ is the scalar flux, $\overrightarrow{J}={\int }_{4\pi }\overrightarrow{\Omega }\psi d\Omega$ is the current, q is the external source (with zeroth and first angular moments $q_0$ and $q_{\alpha }$), $J_{in}={\int }_{\overrightarrow{n}\cdot \overrightarrow{\Omega }<0}|\overrightarrow{n}\cdot \overrightarrow{\Omega }|{\psi }_{b}d\Omega$ is the incoming partial current, ${\varphi }_{in}={\int }_{\overrightarrow{n}\cdot \overrightarrow{\Omega }<0}{\psi }_{b}d\Omega$, is the incoming partial scalar flux, $E_{\alpha \beta }$ is the αβ-th component of the QD (aka “Eddington”) tensor, and C is the boundary QD factor. The nonlinear $E_{\alpha \beta }$ and C factors of the angular flux ψ are given as,

$${\mathrm{E}}_{\alpha \beta }[\psi ]=\frac{{\int }_{4\pi }{\Omega }_{\alpha }{\Omega }_{\beta }\psi d\Omega }{{\int }_{4\pi }\psi d\Omega },$$

$$\mathrm{C}[\psi ]=\frac{{\int }_{\overrightarrow{n}\cdot \overrightarrow{\Omega }>0}\overrightarrow{n}\cdot \overrightarrow{\Omega }\psi d\Omega }{{\int }_{\overrightarrow{n}\cdot \overrightarrow{\Omega }>0}\psi d\Omega }.$$

One of the appealing features of the QD method is that $E_{\alpha \beta }={\mathrm{E}}_{\alpha \beta }[\psi ]$ and $C=\mathrm{C}[\psi ]$ depend only weakly on the angular flux $\psi (\overrightarrow{r},\overrightarrow{\Omega })$, which leads to fast convergence for a wide range of transport problems. Isotropic scattering is only assumed for simplicity—the QD method can be applied to problems with anisotropic scattering as well [9], [10].

The iteration scheme used to solve the QD equations of Eqs. (2a), (2b), (3a), (3b), (3c) is shown in Algorithm 1, using operator notation L for the left side transport operator and A for the left side LOQD operator. This algorithm can be interpreted as a fixed-point iteration method. In steady-state, the low-order problem (LOQD equations) is an elliptic system of first-order partial differential equations with some features of convection–diffusion and tensor diffusion, but containing extra cross-derivative terms. A group of finite volume (FV) methods for discretization of LOQD equations has been developed [5], [6]. An analysis of these methods on randomized meshes showed that these discretization schemes need modifications to improve the performance for these types of meshes [7], [8].

Developing a QD method on unstructured meshes requires developing both a high-order and low-order discretization. Note that the QD method allows independent discretization of these high-order and low-order problems [11].

The subcell balance discretization of Bakirova, Karpov, and Mukhina [12] is used for the high-order (transport) discretization. The method divides the cell into angularly-dependent subcells in which the characteristic transport equation and incoming angular fluxes are used to determine the angular flux within the cell. Integrating the angular flux within a cell produces the face-average and/or cell-average angular fluxes as needed by a particular scheme. A quadrature set is used to approximate integrals over angle, and the transport sweep is used to visit cells in a sequence determined by the dependency graph of a cell on its upwind neighbors. The approximation error of the subcell balance method is due to the interpolation function for the angular flux, constructed on each incoming face from stored unknowns, as well as the approximation of the source within a cell. Other similar conservative short characteristic methods are the slice balance methods of Grove [13] and the split-cell methods of Miller, Matthews, and Brennan [14].

Methods of characteristics have been used to discretize the high-order transport equation for QD as they yield robust and accurate angular flux shapes under many circumstances, e.g. across large, optically thick cells or through void regions. The concern is that unphysical angular flux shapes can cause a chain reaction of producing unphysical QD factors, that produce unphysical currents and scattering sources, that can cause divergence of the non-linear QD iterations. The short characteristics methods are preferred to long characteristics because short characteristics methods are cell-based and also work better with local mesh refinement. Long characteristics methods must overlay a global set of rays across the entire domain. In [8], a non-conservative, vertex-based short characteristics transport method was initially used as the discretization scheme for the high-order transport equation, producing good results on orthogonal Cartesian grids. On arbitrary quadrilateral grids, however, the face-average and cell-average QD factors calculated from vertex angular fluxes had unacceptable numerical fluctuations. Although significantly more expensive, the conservative subcell balance method presented here yielded more stable face-average and cell-average QD factors with much smaller fluctuations on coarse, randomized grids.

The low-order (LOQD) discretization is based on a cell-local finite difference (FD) approximation for the diffusion equation on meshes of arbitrary polyhedrons proposed by Morel [15]. This scheme utilizes (i) representation of the gradient of the scalar flux $\overrightarrow{\mathrm{\nabla }}\varphi$ averaged over a cell and (ii) approximation of $\overrightarrow{\mathrm{\nabla }}\varphi$ by means of face-average and cell-average scalar fluxes using a corresponding reciprocal basis. These two approximations of the gradient are used to derive the gradient of the scalar flux at faces of a cell as well as the representation for the current at the faces. Cells are “linked” by enforcing continuity of the current and scalar flux at each face using interface conditions.

The diffusion discretization of Morel [15] has been chosen due to a number of desirable properties.

• In the limit of a diffusion problem, $E_{xx}=E_{yy}=1/3$ and $E_{xy}=0$, the LOQD discretization should become an appropriate diffusion discretization. By basing it on a diffusion discretization, this can easily be guaranteed.

• The use of cell-average and face-average variables is consistent with previous LOQD discretizations on orthogonal grids. In fact, on orthogonal grids the finite difference LOQD discretization presented here is identical to Gol'din's finite volume LOQD discretization [5].

• The discretization extends to polygons in two spatial dimensions and polyhedra in three spatial dimensions.

• The diffusion discretization converges on randomized spatial grids. The diffusion scheme in [7] did not possess this property.

Based on the expressions relating the gradient within a cell to the average values at faces, a novel linear representation of the scalar flux within a cell has been developed which can be used to construct a linear scattering source for the high-order problem from the solution of the low-order problem, instead of the standard flat approximation. To help ensure the scattering source remains positive everywhere within a cell, thus ensuring positivity of the method, a slope-limiting procedure has been used.

The remainder of this paper is organized as follows. The proposed LOQD discretization is presented in Section 2 and the transport discretization in Section 3. Numerical results are presented in Section 4 and conclusions and discussion are presented in Section 5.