Science and Technology of Nuclear Installations

Volume 2018, Article ID 3965309, 10 pages

https://doi.org/10.1155/2018/3965309

## Analysis of Spatial Discretization Error Estimators Implemented in ARES Transport Code for Equations

School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China

Correspondence should be addressed to Bin Zhang; nc.ude.upecn@nibgnahz

Received 23 April 2018; Accepted 24 June 2018; Published 9 July 2018

Academic Editor: Arkady Serikov

Copyright © 2018 Liang Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The discrete ordinates method (S_{N}) is one of the mainstream methods for neutral particle transport calculations. Assessing the quality of the numerical solution and controlling the discrete error are essential parts of large-scale high-fidelity simulations of nuclear systems. Three error estimators, a two-mesh estimator, a residual-based estimator, and a dual-weighted residual estimator, are derived and implemented in the ARES transport code to evaluate the error of zeroth-order spatial discretization for S_{N} equations. The difference in scalar fluxes on coarse and fine meshes is adopted to indicate the error in the two-mesh method. To avoid zero residual in zeroth-order discretization, angular fluxes within one cell are reconstructed by Legendre polynomials. The error is estimated by inverting the discrete transport operator using the estimated directional residual as an anisotropic source. The inner product of the forward directional residual and the adjoint angular flux is employed to quantify the error in quantities of interest which can be denoted by a linear functional of forward angular flux. Method of Manufactured Solutions (MMS) is adopted to generate analytical solutions for S_{N} equation with scattering and the determined true error is used to evaluate the effectivity of these estimators. Promising results are obtained in the numerical results for both homogeneous and heterogeneous cases. The larger error region is well captured and the average effectivity index for the local error estimation is less than unity. For the series test problems, the estimated goal quantity error can be contained within an order of magnitude around the exact error.

#### 1. Introduction

The Boltzmann transport equation [1] is central to the neutral particle transport problems as encountered in the field of nuclear engineering community. However, the numerical solution of this transport equation for realistic three-dimensional radiation shielding problems still remains a challenge because of the high-dimensional nature of the solution space and the nonsmoothness of the exact solution. In the deterministic transport calculation, a large part of errors stem from the numerical discretization of all its variables; therefore it is important to develop an approximate measure of these discretization errors. With the estimated errors, one can get the confidence level of the numerical solution and avoid unreasonable discretization. Also, the estimated error can be used to drive an adaptive mesh refinement (AMR) algorithm.

Over the past 20 years, the concept of adaptivity has been widely employed to obtain accurate solutions to the transport equation with fewer unknowns. As an essential part of adaptivity, error estimation remains active research area. In 1973, Madsen developed a method for calculating rigorous a posteriori error bounds for steady-state neutron transport problems [2]. In 1998, Jessee adopted the normalized gradient of the radiant intensity to measure the local truncation error in the radiative transport equation and developed an AMR algorithm with this measurement [3]. In 2003, the value of the mean free path in a certain cell was adopted to drive mesh refinement in the code Styx [4]. For the discrete ordinates nodal transport methods, Azmy derived a posteriori error estimator as a function of volume residual and surface jumps, which were approximated by solving local problems on a cell by cell basis. Also, the global norm error is split into a local error indicator to account for the error of each element [5, 6]. Based on the arbitrary high-order discontinuous Galerkin method, Fournier presented an error estimator for the S_{N} transport equation and investigated its adoption to nonconforming meshes and a variable polynomial order [7]. To drive the AMR process, Wang used an error estimator based on a two-mesh approach and a jump-based error indicator. Also, he evaluated the error in the quantity of interest (QoI) with the direct and adjoint errors [8]. For the DGFEM-0 spatial discretization, Azmy estimated the residual using Taylor expansion and inverted the discrete transport operator using the residual as a source to estimate the error [9]. Recently, a dual-weighted residual (DWR) method was employed to estimate the error of the QoI. In 2011, Lathouwers combined the forward residuals and adjoint error to give accurate error measures for the effective multiplication factor in two-dimensional reactor problems [10]. For the diamond difference (DD) scheme, Jeffers derived DWR method to estimate the error in the QoI arising from the spatial discretization for fixed source and eigenvalue neutron transport problems [11].

In this paper, we assess the accuracy of the error estimator implemented in the ARES [12] transport code. The reminder of this paper is organized as follows. In Section 2, we briefly present the theory of the error estimator implemented in ARES transport code, including the two-mesh estimator, the residual-based estimator, and DWR method. Test problem specifications are described and numerical results are discussed in Section 3. Also, we provide concluding remarks in Section 4.

#### 2. Theory of Error Estimator

The discrete ordinate method (S_{N}) is employed to deal with the discretization of angular variables. And the monoenergetic steady-state S_{N} transport equation with isotropic scattering iswhere is the angular flux at position along the discrete direction , and are the macroscopic total and scattering cross sections, respectively, and is the fixed-source term. The scalar flux is calculated usingwhere is the quadrature weight corresponding to the discrete direction and is the total number of discrete directions.

Weighting equation (1) by a series of test function and integrating over each mesh give a series of spatial moment balance equations:where denotes the outward normal direction of mesh boundary .

Multiple spatial discretization schemes are available in ARES. The traditional low-order schemes, including DD, weighted diamond difference (WD) [13], and step characteristic (SC) [14, 15], are based on the cell-balance equation, which can be obtained by setting in (3). Additional auxiliary equations, which assume the behavior of the angular fluxes across the spatial cell, are required to close the cell-balance equation. For example, the DD scheme assumes that the cell-averaged flux is a linear average of opposing surface-averaged fluxes, whereas the SC scheme analytically inverts the transport operator along the discrete directions within one cell. The finite-element method [16, 17] approximates the angular fluxes within one cell by a linear combination of trial functions .Using identical test spaces and trial spaces yields the Galerkin approximation. For example, the linear discontinuous (LD) scheme expands the angular flux in the trial function space , whereas the trilinear discontinuous (TLD) scheme is defined in the trial function space . In this paper, we mainly focus on the error estimators implemented to the zeroth-order spatial discretization.

##### 2.1. Two-Mesh Estimator

Considering the justification that spatial discretization errors decrease as the spatial cell is refined in the asymptotic regime, one can use the difference in the solution obtained on the mesh (coarse) and mesh (refined) to estimate spatial discretization errors. Generally, this estimator underestimates the true error but captures the error distribution accurately. Another drawback is that the accuracy of this estimator is also affected by the occurrence of irregularities [18] in angular fluxes, which excludes some spatial cells from the asymptotic regime.

In the transport sweep process, angular fluxes are not stored; therefore the angular integrated quantity, namely, the scalar flux, is chosen as the basis of this estimator. In this paper, two different calculation schemes are investigated:where are the estimated errors in each cell, and are scalar fluxes calculated on the mesh and mesh, and is the scalar flux projected from the mesh to the mesh. The projection process is based on zeroth spatial moment conservation for low-order discretization. In scheme 1, the absolute value is taken based on the coarse mesh, whereas the absolute value is taken on fine mesh in scheme 2.

##### 2.2. Residual-Based Estimator

We rewrite (1) in operator form:where represents the spatially continuous transport operator and denotes the scattering operator. and are the spatially continuous angular fluxes and fixed source, respectively. Discretizing (7) over the mesh yieldswhere represents the discrete transport operator, is the numerical solution obtained on the mesh, and is the fixed source projected onto the discrete domain. The spatial discretization error on the mesh is defined as , with being the exact solution to the spatial continuous equation projected onto the discrete domain. Noting the linearity of the transport operator and the scattering operator , we obtainWith the residual defined as , the error can then be obtained by solving (9) with the same iterative method for the transport equation.

To calculate residual in the zeroth-order spatial discretization schemes, we approximately reconstruct the profiles of angular fluxes on the mesh with normalized Legendre polynomials:Here, , where , are the th-order normalized Legendre polynomials. And is the vector of coefficients and is obtained through the cell-averaged and six face-averaged angular fluxes. With these assumed profiles in hand, we can calculate the residual on the mesh from:whereFrom (9), we can see that inverting the transport operator on the directional residual returns the estimated error. Note that the integral of the directional residual over each mesh is zero. To avoid the zero-residual term, we solve (9) on the mesh:In the zeroth-order scheme, is the average residual value on the mesh and is obtained by locally integrating on the mesh. Considering that the residual may be negative or positive, the DD scheme is employed to solve (13). The following two schemes are investigated to project the error from the mesh to the mesh:In scheme 1, the error on the mesh is integrated over the mesh first and then the absolute value is taken as the estimated error on the mesh. In scheme 2, the error on the mesh is estimated by integrating the absolute value of the error on the mesh.

##### 2.3. Dual-Weighted Residual Estimator

In realistic shielding calculations, an accurate solution may not be required throughout the entire computational domain but only some linear functionals of the angular fluxes, such as the detector reaction rate and the average scalar flux in the subdomain. The linear functionals of the solution can be denoted by :To estimate the error in QoI, the following adjoint transport equation must be solved firstly:where is the adjoint transport operator and denotes adjoint angular fluxes. The adjoint angular flux is commonly referred to as the importance function of some goals, which is dependent on the adjoint source term .

As implied in (16), the adjoint source term is dependent on the goal of our calculation. With the relationship of forward and adjoint transport equation, the error in QoI can be derived as follows:This expression indicates that the error in the quantity of interest can be estimated by the inner product of the adjoint fluxes and the residual of forward transport equation. In our zeroth-order schemes implementation, the adjoint angular fluxes are also solved on the mesh and reconstructed in the same form as that in (10). where . Then, (18) is deduced toFrom the orthogonality of the Legendre polynomials, the inner product of the residual and adjoint fluxes over one cell along one direction are found to be

#### 3. Numerical Results and Discussion

For an assessment of the accuracy of these error estimators, two series test problems are set up and the reference solutions are obtained by the MMS [19, 20] method. The two-mesh estimator and residual-based estimator are used to obtain local errors, while the DWR estimator is employed to estimate the error of region-integrated fluxes. Also, some measurements are introduced to quantify the effectiveness of the estimated errors.

##### 3.1. Test Problem Description and Estimator Rubrics

In the test problem suites, a box-shaped domain including two alternate variant materials is considered and depicted in Figure 1. The geometry configuration has dimension of 20 cm × 20 cm × 20 cm. In the homogeneous problem, the total cross section is adopted for the two materials, whereas in the heterogeneous problem, the two different total cross sections and are employed for the two materials. A vacuum boundary condition is applied to the , , and faces of the domain and a reflecting boundary condition is applied to the , , and faces. Uniform 1 cm × 1 cm × 1 cm meshes and P_{N}T_{N} S16 quadrature sets [21] are used in the calculation performed in this paper.