Accurate spherical harmonics solutions for neutron transport problems in multi-region spherical geometry

Highlights

• A numerically stable $P_N$ method for solving multi-region sphere problems is proposed.

• The method avoids spatial and angular discretizations.

• A post-processing step is used to further improve the $P_N$ solution.

• A new way of computing accurate post-processed currents is devised.

• Highly accurate numerical results are reported for three test cases.

Abstract

A stable spherical harmonics ($P_N$) method for solving the problem of a multi-region sphere with internal sources subject to an externally incident angular flux is developed. The method consists in solving the problem locally (region by region) and using an iterative sweep technique to connect and update the solutions for the different regions until global convergence is attained. A post-processing step which is very effective in improving the $P_N$ solutions is also developed. Since calculating the current from its definition in terms of the angular flux is prone to loss of precision in spherical geometry, an alternative scheme for computing accurate post-processed currents is devised.

Because the developed method does not rely on space and angle discretizations, it can achieve very accurate results even for problems with solutions that display strong flux gradients. This kind of problem is challenging for methods built on discretization schemes. An example of such a problem is solved and the results obtained demonstrate the good performance of the method. Accurate numerical results are also given for two variants of the example problem that involve changes in the driving term and in the scattering law.

Introduction

Transport problems in spherical geometry occur in many areas. In reactor physics, there are subjects of study that require solving the neutron transport equation in spherical geometry. One could mention, for example, critical systems such as Godiva [1], microspherical fuel particles [2], and spherical fuel elements used in pebble-bed reactors [3]. In this context, the interesting problem of a double heterogeneity [4] arises when the last two appear combined [5]. Other examples of nuclear applications that involve dealing with spherical geometry are space radiation shielding [6], neutron transmission experiments through spherical shells [7], radiation transport simulations of inertial confinement fusion experiments [8], and radiative transfer in a reactor core under severe accident conditions [9]. Radiative transfer in spherical geometry is also important to other areas of application, as discussed next.

In studies of atmospheric radiative transfer, many existing codes either ignore the Earth's sphericity completely or use approximations such as the pseudo-spherical model [10], where the uncollided component of the radiance is treated exactly in spherical geometry and the collided one only approximately, using plane geometry. Both of these approaches may lead to substantially large errors when the sun is positioned at a large zenithal angle in the sky [11], more so when a sub-horizon position is involved. In mechanical engineering, modeling of high-temperature industrial processes often require solving the radiative heat transfer equation in spherical cavities and cavities delimited by concentric spherical surfaces [12], usually coupled with other forms of heat transfer (conduction, convection). Moreover, in an astrophysical setting, solving the radiative transfer equation in spherical geometry may be necessary in many situations, such as when studying infrared objects, stellar atmospheres, planetary nebulae, quasars, and so on [13].

In this work, a new method is developed for solving accurately the problem of a multi-region sphere with internal sources subject to an externally incident angular flux. The method makes use of two recently reported $P_N$ solutions in spherical geometry [14], [15] and consists in solving the problem locally (on a region by region basis) and coupling the local solutions iteratively. For this purpose, a series of spatial sweeps is performed back and forth across the regions to connect and update the local solutions, until global convergence is attained. Solving a multi-region problem this way may be more efficient, especially when the number of regions characterizing the problem is large. In addition, it allows the use of different orders of the $P_N$ approximation in the regions. Although not explored in this work, this may be an interesting option for some problems [16].

With the usage of the stable $P_N$ solutions reported in Refs. [14] and [15], the well-known numerical instability of the standard $P_N$ method in spherical geometry [17] is avoided. Our solution does not involve spatial and angular discretizations and, after being post-processed by means of an integral equation derived from the technique of integration along the neutron path, is capable of reproducing accurately the angular flux discontinuities that result from shadowing effects in spherical geometry [18].

High-quality numerical results are reported for a three-region, isotropically scattering case originally proposed and solved by Williams [19], for a second case identical to the first but with a different driving term and for a variant of the first case that includes anisotropic scattering. We hope that these results can be useful as benchmarks for verification of general purpose codes and newly developed solution methods.