Solution of the Generalized Linear Boltzmann Equation for Transport in Multidimensional Stochastic Media

TL;DR

This paper extends the generalized linear Boltzmann equation framework to multidimensional stochastic media with finite boundaries, using discrete ordinates and approximations, verified through numerical experiments.

Contribution

It introduces a multidimensional solution approach for the GLBE with finite boundaries, including discretization techniques and approximations, advancing beyond previous one-dimensional models.

Findings

Verified the solution against Gaussian process media

Demonstrated the effectiveness of sweeping methods in multidimensional settings

Introduced practical approximations for non-classical cross-sections

Abstract

The generalized linear Boltzmann equation (GLBE) is a recently developed framework based on non-classical transport theory for modeling the expected value of particle flux in an arbitrary stochastic medium. Provided with a non-classical cross-section for a given statistical description of a medium, any transport problem in that medium may be solved. Previous work has only considered one-dimensional media without finite boundary conditions and discrete binary mixtures of materials. In this work the solution approach for the GLBE in multidimensional media with finite boundaries is outlined. The discrete ordinates method with an implicit discretization of the pathlength variable is used to leverage sweeping methods for the transport operator. In addition, several convenient approximations for non-classical cross-sections are introduced. The solution approach is verified against random realizations of a Gaussian process medium in a square enclosure.

Equations

$$\frac{\partial \psi }{\partial s}+\hat{\Omega }\cdot \nabla \psi +{\Sigma }_t(s)\psi =\delta (s)c{\int }_{4\pi }{\int }_0^{\infty }P({\hat{\Omega }}^{\prime }\cdot \hat{\Omega }){\Sigma }_t(s^{\prime })\psi (s^{\prime })d{s}^{\prime }d{\Omega }^{\prime }+\delta (s)\frac{Q}{4\pi }$$

$$\frac{\partial \psi }{\partial s}+\hat{\Omega }\cdot \nabla \psi +{\Sigma }_t(s)\psi =0.$$

$$\psi (\vec{r},\hat{\Omega },s=0)=c{\int }_{4\pi }{\int }_0^{\infty }P({\hat{\Omega }}^{\prime }\cdot \hat{\Omega }){\Sigma }_t(s^{\prime })\psi (s^{\prime })d{s}^{\prime }d{\Omega }^{\prime }+\frac{Q(\vec{r})}{4\pi }.$$

$$\Psi (\vec{r},\hat{\Omega })={\int }_0^{\infty }\psi (\vec{r},\hat{\Omega },s)ds,$$

$${\Psi }^{(0)}={\int }_{4\pi }\Psi d\Omega .$$

$${\Psi }^{(1)}={\int }_{4\pi }\Psi \hat{\Omega }d\Omega$$

$${\int }_{4\pi }{\int }_0^{\infty }\left(\frac{\partial \psi }{\partial s}+\hat{\Omega }\cdot \nabla \psi +{\Sigma }_t(s)\psi \right)dsd\Omega =0$$

$$\nabla \cdot {\Psi }^{(1)}=Q(\vec{r})-(1-c){\int }_{4\pi }{\int }_0^{\infty }{\Sigma }_t(s)\psi (s,\hat{\Omega })dsd\Omega$$

$$\psi ({\vec{r}}_w,\hat{\Omega },s)=\delta (s){\Psi }_o({\vec{r}}_w,\hat{\Omega })$$

$$\psi ({\vec{r}}_w,\hat{\Omega },s=0)=\delta (s)\frac{\rho }{\pi }{\int }_0^{\infty }{\int }_{\hat{\Omega }\cdot \hat{n}<0}\psi ({\vec{r}}_w,\hat{\Omega },s)|\hat{n}\cdot \hat{\Omega }|d\Omega ds$$

$$\frac{\partial {\psi }_l}{\partial s}+\hat{\Omega }\cdot \nabla {\psi }_l+{\Sigma }_t(s){\psi }_l=0.$$

$${\psi }_{l+1}(\vec{r},\hat{\Omega },s=0)=c{\int }_{4\pi }{\int }_0^{\infty }P({\hat{\Omega }}^{\prime }\cdot \hat{\Omega }){\Sigma }_t(s^{\prime }){\psi }_l(s^{\prime })d{s}^{\prime }d{\Omega }^{\prime }+\frac{Q}{4\pi }.$$

$${ϵ}_{l+1}=\frac{|{\Psi }_{l+1}-{\Psi }_l|}{|{\Psi }_l|}$$

$$\frac{\partial {\psi }_l}{\partial s}+\mu \frac{\partial {\psi }_l}{\partial z}+{\Sigma }_t{\psi }_l=0$$

$${\psi }_{l+1}(z,\mu ,s=0)=\frac{c}{2}{\int }_{-1}^1{\int }_0^{\infty }{\Sigma }_t{\psi }_l(s^{\prime })d{s}^{\prime }d\mu +\frac{Q}{4\pi }$$

$$\frac{\partial f_l}{\partial s}+\mu \frac{\partial f_l}{\partial z}+{\Sigma }_t{f}_l=0$$

$$f_{l+1}(z,\mu ,s=0)=\frac{c}{2}{\int }_{-1}^1{\int }_0^{\infty }{\Sigma }_t{f}_l(s^{\prime })d{s}^{\prime }d\mu$$

$$f_l=A{\rho }^l\exp (i\lambda {\Sigma }_{t}z)g(s,\mu ).$$

$$\frac{\partial g}{\partial s}+(1+i\mu \lambda ){\Sigma }_{t}g(s,\mu )=0$$

$$g=h(\mu )\exp \left(-(1+i\mu \lambda ){\Sigma }_{t}s\right)$$

$$f_l=A{\rho }^l\exp (i\lambda {\Sigma }_{t}z)\exp \left(-(1+i\mu \lambda ){\Sigma }_{t}s\right)h(\mu )$$

$$\rho h(\mu )=\frac{c}{2}{\int }_{-1}^1{\int }_0^{\infty }\exp \left(-(1+i\mu \lambda ){\Sigma }_t{s}^{\prime }\right)h(\mu ){\Sigma }_{t}d{s}^{\prime }d\mu$$

$$\rho h(\mu )=\frac{c}{2}{\int }_{-1}^1\frac{h(\mu )}{1+i\lambda \mu }d\mu$$

$${\rho }_{\lambda }=\frac{c}{\lambda }{\tan }^{-1}\lambda$$

$${\rho }_{max}=c$$

$$\frac{\partial {\psi }_l}{\partial s}+\hat{\Omega }\cdot \nabla {\psi }_l+{\Sigma }_t(s){\psi }_l=0.$$

$$f_l=\psi -{\psi }_l$$

$$\frac{\partial f_l}{\partial s}+\hat{\Omega }\cdot \nabla f_l+{\Sigma }_t(s)f_l=0,$$

$$f_l(s=0)=-{\psi }_l(s=0)+c{\int }_{4\pi }{\int }_0^{\infty }P(\widehat{{\Omega }^{\prime }}\cdot \hat{\Omega }){\Sigma }_t(s^{\prime })({\psi }_l(s^{\prime })+f_l(s^{\prime }))d{s}^{\prime }d{\Omega }^{\prime }.$$

$${\psi }_{l+1}(s=0)=c{\int }_{4\pi }{\int }_0^{\infty }P(\widehat{{\Omega }^{\prime }}\cdot \hat{\Omega }){\Sigma }_t(s^{\prime })({\psi }_l(s^{\prime })+f_l(s^{\prime }))d{s}^{\prime }d{\Omega }^{\prime }$$

Full text

1 Introduction

Stochasticity in the spatial distribution of participating media can have a large impact on the statistics of particle transport. Pebble-bed nuclear reactor neutron fluxes are modified by the distribution of the fuel pebbles compared to predictions from an equivalent homogenization [1]. In high temperature combustion, turbulence can induce highly non-uniform gas densities and blackbody radiation intensities, leading to enhanced thermal emissions from flames [2, 3]. The stochastic distribution and structure of clouds in the atmosphere have various impacts on climate models of solar and terrestrial radiation [4]. The transmission of solar energy through particle curtains subjected to turbulence is also of concern for solar energy devices [5].

Transmission through a spatially correlated medium has been studied extensively in recent years. It was shown in [6, 7, 8, 9] that depending on whether the medium is positively or negatively correlated at the microscopic level, radiation may undergo sub-exponential or super-exponential attenuation respectively relative to the rate predicted from Beer’s law. In [10], the spectra of different stochastic processes were analyzed and related to observed rates of radiation attenuation.

Due to the deviations from the classic transport predictions, as well as the possible variance in results, it is important to develop methods to capture this variability. Monte Carlo sampling of random media is an unbiased way of gathering this information, but this requires many samples of the media which may not be available, as well as a solution of the transport equation in each sample. Such studies can become prohibitively costly. In some cases, it is possible to perform dimension reduction [11, 12] to reduce the number of samples required, but this is not always guaranteed and requires an identification of the spectrum of the stochastic medium. In the limit of short length-scales and mild variability, the so-called atomic mixing approximation [13] may be valid and allow one to neglect stochastic effects, but this limit is prohibitive and inaccurate for many cases of interest. Alternatively, one might seek an effectively reduced cross-section that captures the variability [6, 14], but this approximation is only valid for mean-free-paths much longer than the length scale of variability. A commonly used approach in mixtures of discrete binary media is the Levermore-Pomraning model [15]. With some approximations in the nature of the transport through the medium, a coupled set of transport equations for the average angular flux may be derived and solved using existing transport methods. Though this model is convenient, it suffers from a loss of accuracy at higher scattering ratios and is obviously restricted to the special case of binary Markovian media.

A recently developed model called the Generalized Linear Boltzmann Equation (GLBE) is related to non-classical particle transport, in which the particles are considered to have a probability of extinction that depends on how far they have already traveled [1]. This model unifies all possible forms of stochastic media into a non-classical cross-section that depends on the particle pathlength and may also incorporate spatial and angular dependence [16]. Solution methods for the GLBE have been analyzed using Monte Carlo transport, explicit pseudo-time stepping [17], and implicit pseudo-time stepping [18].

The GLBE is a promising approach for more general problems, but studies thus far have been limited in scope and have yet to address problems involving finite boundaries with sources, and the problem of computing the non-classical cross-section is a potential roadblock for using this approach. In this work, the GLBE will be reviewed and extended to incorporate finite boundaries. The numerical solution using implicit time-stepping and source iteration is considered in order to leverage discrete ordinates methods for solving the resulting equations. Methods for computing the non-classical cross-section are demonstrated, and analytical approximations for binary Markovian media and low-variability Gaussian processes are derived. Finally, the GLBE solution will be shown for two problems in a square enclosure and validated against a Monte Carlo approach for a medium described by a truncated Gaussian process.

2 Generalized Linear Boltzmann Equation

As derived in [1], the generalized linear Boltzmann equation (GLBE) for the angular flux $\psi(\vec{r},\hat{\Omega},s)$ at a location $\vec{r}$ with direction $\hat{\Omega}$ of particles that have traveled a distance $s$ is given by

[TABLE]

where $\Sigma_{t}(s)$ is the non-classical cross-section, $c$ is the single scattering fraction, $P(\hat{\Omega}^{\prime}\cdot\hat{\Omega})$ is the scattering phase function, $Q$ is an arbitrary particle source term, and $\delta(s)$ denotes the Dirac delta function.

The GLBE may be recast into a more convenient formulation by separately considering the cases where $s>0$ and the limit of $s\to 0$. By first considering the case where $s>0$, we find the GLBE simplifies to an initial value problem:

[TABLE]

In the limit of $s\to 0$, an initial condition is derived by integrating the GLBE over the range $[-\varepsilon,\varepsilon]$ in the limit that $\varepsilon\to 0$:

[TABLE]

This initial condition is implicit in the case where the stochastic medium scatters particles. However, in non-scattering cases, the initial condition is fully specified. With a solution for the angular flux, the classical angular flux may be derived by integrating the angular flux over all particle path lengths

[TABLE]

The particle scalar flux $\Phi^{0}$ and current $\Phi^{1}$ may be derived by taking the zeroth angular moment of the ensemble-averaged angular flux

[TABLE]

[TABLE]

The divergence of the particle current $\nabla\cdot\Psi^{(1)}$ may be derived from taking the moment of Equation 2:

[TABLE]

Carrying out the integrals and applying Equation 3 gives a direct formula for the divergence of the current as

[TABLE]

2.1 Boundary Conditions

Previous discussions [1, 18] of the GLBE have assumed the case of infinite domains to simplify the analysis of transport in stochastic media. As we are concerned with applying the GLBE to realistic problems, it is important to allow for finite boundaries and particle sources entering the domain from a boundary. A specified boundary source $\Psi_{o}$ at boundary locations $\vec{r}_{w}$ may be incorporated into the GLBE by applying the initial conditions

[TABLE]

This boundary condition is analogous to the source term in the GLBE, and is applied before any particles have collided. Although this work only considers absorbing walls, a diffusely reflecting wall with reflectance $\rho$ may be modeled by

[TABLE]

where $\hat{n}$ is the wall normal vector.

3 Solution of the GLBE

An efficient solution method for the GLBE will be described in this section. This method is closely related to the standard source iteration solver that is commonly applied to the steady linear Boltzmann equation.

3.1 Fixed Point Iteration

Equations 2 and 3 together may be used to form an iterative solver for the GLBE. With some guess for the pre-collision angular flux $\psi^{l}(\vec{r},\hat{\Omega},s=0)$ at iteration $l$, the angular flux at all path lengths may be computed by solving the equation

[TABLE]

The angular flux initial condition at iteration $l+1$ is computed by evaluating the following equation:

[TABLE]

The solution may be repeated until convergence is achieved. For the purposes of this work, the convergence metric is taken to be the relative norm of the residual of the classic angular flux $\epsilon^{l+1}$:

[TABLE]

3.2 Convergence

The solution approach described here has the same spectral radius as the source iteration solver. To prove this, the case of classical particle transport is assumed, i.e. $\Sigma_{t}(s)$ is a constant for all $s$. Furthermore, we consider one-dimensional transport in the $z$-direction with angular coordinate $\mu=\hat{z}\cdot\hat{\Omega}$ in an infinite domain and isotropic scattering. In that case, the equations reduce to

[TABLE]

[TABLE]

The system of equations may be recast in terms of the residual $f_{l}(z,\mu,s)=\psi(z,\mu,s)-\psi_{l}(z,\mu,s)$:

[TABLE]

[TABLE]

In order to understand the convergence behavior of the iterative system, the residual function is expanded in Fourier modes with wavenumber $\lambda$ as

[TABLE]

This expansion is similar to the method used in [19] to analyze the convergence of the source iteration method. The extra step needed in this case is to determine the Fourier modes associated with the pathlength variable $s$. By applying the Fourier expansion to Equation 16, we can solve for $g(s,\mu)$:

[TABLE]

[TABLE]

where $h(\mu)$ is a function of $\mu$ resulting from the integration. The residual function is now given by

[TABLE]

It remains now to solve for $h(\mu)$ by using Equation 17:

[TABLE]

Performing the integration over $s$ gives the familiar integral equation

[TABLE]

where $h(\mu)$ is clearly a constant which may be removed from the equation. The result gives the spectral radius for each wavenumber

[TABLE]

with the slowest converging mode, determining the convergence of the iterative algorithm, being in the limit of $\lambda\to 0$

[TABLE]

This result is identical to that of the source iteration algorithm. Evidently the iterative approach for solving the GLBE retains the same convergence behavior as the standard source iteration method in the limit of classical particle transport. For non-classical transport, i.e. the case where $\Sigma_{t}(s)$ is variable in $s$, it is expected that one could use the maximum value of $\Sigma_{t}(s)$ over all $s$ as an upper bound for the convergence analysis and show that the spectral radius is still $\rho_{max}=c$.

The fact that the spectral radius of the iterative approach for solving the GLBE is identical to that of source iteration means that the GLBE solver inherits the same advantages and disadvantages as source iteration. The well-defined spectral radius means that one can predict the approximate number of iterations required to reach a desired level of convergence, and in the limit that the medium is purely absorbing ($c=0$), no iterations are required at all. However, in strongly scattering media ($c\to 1$), the spectral radius approaches $1$, and the iterative solver may become an impractical choice.

3.3 Synthetic Acceleration

As an improvement on source iteration, [19] demonstrated that synthetic acceleration methods based on lower order forms of the transport equation may be successfully employed as preconditioners. Similar equations can be derived for the GLBE iterative solver. At iteration $l$, the following equation is solved to derive the angular flux:

[TABLE]

As in the previous section, the residual between the current iterate and the true solution is defined as

[TABLE]

for which a transport equation may be derived

[TABLE]

along with the initial condition

[TABLE]

Having solved for $f_{l}$, the initial condition for the following iterate is determined by

[TABLE]

and the classical angular flux is given by

[TABLE]

If one were able to solve Equations 28, 29, and 30 exactly, the true solution would be achieved in one step. However, those equations are just as difficult to solve as the original problem. As in [19], one could invoke a lower order approximation for the GLBE operator in Equation 28, such as transport synthetic acceleration or using the diffusion approximation derived in [1] for diffusion synthetic acceleration. However, we note that the integration in the pathlength variable allows for solving Equation 28 using a coarser pathlength stepsize, and that the number of iterations required to solve Equation 29 may also be limited for further efficiency gains in the inner loop. This strategy is here referred to as pathlength synthetic acceleration. Provided that the coarsening is not too aggressive, pathlength synthetic acceleration provides a convenient approach to reuse the GLBE solver to improve the rate of convergence, and may even be combined with transport synthetic acceleration.

3.4 Implicit Discretization

The problem of discretizing and solving the GLBE using the iterative solver is now considered. It is necessary to track the evolution of $\psi_{l}$ as a function of $s$ with high accuracy in order to evaluate the initial condition and classical angular flux accurately. At the same time, it is important to reduce the cost of the integration, as many iterations may be required to achieve a converged solution. If $\Sigma_{t}(s)$ is not a strong function of $s$ or $\psi_{l}$ is slowly changing, it is possible to use larger steps in the pathlength integration $\Delta s$ in achieving a steady solution. Thus we consider implicit trapezoidal integration in $s$ as in a pseudo-time integration, allowing for second order accurate convergence both in the evaluation of $\psi_{l}(s)$ as well as $\Psi_{l}$. An implicit discretization also has the advantage of allowing for a step size that is not limited by the spatial grid size for stability. In this work we assume the step size is fixed without loss of generality. The resulting semi-discrete equation for the angular flux is

[TABLE]

where $\psi_{l}^{n}$ and $s^{n}$ are the angular flux and pathlength at step $n$ respectively. By regrouping terms for the angular flux at $n+1$, the discretized equation is

[TABLE]

which has the same form as the quasi-steady Boltzmann transport equation with pseudo-time dependent cross-section and source terms. By applying an angular discretization via the discrete ordinates method and combining it with a spatial discretization such as diamond differencing or discontinuous Galerkin methods, the resulting equation may be efficiently inverted for $\psi_{l}^{n+1}$ by using full domain sweeps.

With a trapezoidal integration in $s$, the classical angular flux and initial condition are given by

[TABLE]

[TABLE]

Rather than storing $\psi_{l}$ for all $s$, the classical angular flux and initial condition may be updated incrementally as the integration proceeds to reduce computational memory requirements.

The final consideration that needs to be made is for the boundary condition. A simple approach is used here to regularize the Dirac delta function in the boundary condition by approximating it as a rectangular pulse of duration $\Delta s$ and height $1/\Delta s$, leading to

[TABLE]

4 Application to Transport in Two-Dimensional Stochastic Media

The solution of the GLBE will be demonstrated in this section for a few different example problems, and the results will be compared against Monte Carlo runs to demonstrate convergence of the GLBE solution. A process for generating the cross-section $\Sigma_{t}(s)$ is described first. In this section we will also make the simplifying assumption of isotropic scattering, i.e. $P(\hat{\Omega}^{\prime}\cdot\hat{\Omega})=1/4\pi$ without loss of generality.

4.1 Generation of Stochastic Cross Sections

In this work the case of stochastic media determined by Gaussian processes with a spatially constant mean $\Sigma_{o}$ and a Gaussian covariance function is considered. Then independent realizations of stochastic media may be constructed by using Karhunen-Loeve expansions as in [12]. Starting from the Gaussian covariance function with variance $\sigma^{2}$ and correlation length $l$

[TABLE]

the covariance matrix may be evaluated on a discrete grid. With a covariance matrix of size $M\times M$, the $M$ eigenvalues $\lambda_{i}$ and eigenvectors $v_{i}$ are computed from

[TABLE]

Defining $Z_{i}$ as independent Gaussian random variables with zero mean and unit variance, realizations of the cross-section are generated as

[TABLE]

Since this process may lead to negative values of $\Sigma_{t}$, negative values are truncated to zero. This results in an increase in the mean and decrease in the variance.

Once many realizations of $\Sigma_{t}$ have been generated, the pathlength distribution and non-classical cross section may be computed. Following the definition of the pathlength probability distribution function in [1]

[TABLE]

with the brackets indicating an ensemble-average, many one-dimensional stochastic realizations of $\Sigma_{t}(s)$ are generated, and the integral is computed numerically using trapezoidal integration. With a well-resolved $p(s)$ in hand, the non-classical cross-section is computed from

[TABLE]

with interpolation to a finer grid as needed for increased accuracy in the GLBE solution. Care is also taken to ensure the length of the domain is long enough that a very small number of particles survive to the end; without this measure, the above equation may yield unphysical results for the cross-section.

Figure 1 shows a few examples of the non-classical cross-section for different truncated Gaussian processes. The resulting curves reflect the underlying statistics of the problem. Initially, all of the particles enter into the medium with the same rate of attenuation on average. Particles in regions of high material density tend to be extinguished faster, so the effective cross-section that the remaining particles see decreases. Past a distance on the order of the correlation length of the medium, the attenuation rate equilibrates as particles see a mixture of high and low density material, and the cross-section asymptotes to a constant cross-section identical to a homogenized constant. The figures also represent differing magnitudes of variation in the cross-section, with lower coefficients-of-variation corresponding to smaller decreases, and longer correlation lengths corresponding to longer decay lengths.

4.2 Analytical approximations

In some cases, it may not be feasible to estimate the non-classical cross-section directly from Monte Carlo sampling; the number of realizations required, along with the spatial resolution and maximum pathlength required for the analysis, may be prohibitively expensive. It is useful to have convenient approximations available in that case. In other cases, the mean line transmission $\langle\Psi(s)\rangle$ may be available already, so $\Sigma_{t}(s)$ may be estimated directly from

[TABLE]

Analytical non-classical cross-sections are derived in the following sections for the special cases of binary Markovian media and Gaussian processes with low coefficient of variation.

4.2.1 Binary Markovian media

For discrete mixtures of two separate components, the approximation of a binary Markovian medium is commonly made. This approximation assumes that the pathlength of a ray passing though any one component of the mixture has an exponential distribution. Using the results of [15], the mean line transmission for a binary Markovian medium is given by

[TABLE]

where the various parameters $r_{\pm}$, $\tilde{\sigma}$, and $\langle\sigma\rangle$ are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\sigma_{1}$ and $\sigma_{2}$ are the cross-sections for the two components in the mixture, $p_{1}$ and $p_{2}$ are the volume fractions of each component, and $\lambda_{1}$ and $\lambda_{2}$ are the average path-lengths within each fluid.

The choice of a Markovian model for the mixture is convenient due to the simplicity of the Levermore-Pomraning model. However, this model is known to suffer from a loss of accuracy at higher scattering ratios due to the need for an additional closure condition. Instead, the above result may be used to derive the non-classical cross section for use in the GLBE as

[TABLE]

This analytical form would allow for convenient incorporation of a binary Markovian medium into the GLBE framework without a loss of accuracy at higher scattering ratios.

4.2.2 Quasi-Gaussian media

For continuously varying media with low coefficient of variation, [20] has demonstrated that another analytical approximation may be made, the derivation of which will be repeated here. It is assumed that the mean, variance, and two-point correlation function of the medium are known.

Define $\tau$ to be the optical depth across some path length in the random medium with cross-section $\sigma(x)$:

[TABLE]

Since $\sigma$ is a random field, $\tau$ is also a random variable. The mean of $\tau$ is given by

[TABLE]

Taking the decomposition of the field into a stationary mean and a fluctuating component $\sigma=\langle\sigma\rangle+\sigma^{\prime}$, the fluctuation of $\tau$ is

[TABLE]

from which the variance may be computed as

[TABLE]

With some manipulation of this integral, we arrive at a convenient expression for the variance as

[TABLE]

which may be further simplified to

[TABLE]

where $\Phi_{\sigma\sigma}(r)=\langle\sigma^{\prime}(z)\sigma^{\prime}(z+r)$ is the stationary two-point correlation function, independent of location $z$.

With the approximation that $\tau$ has a nearly normal distribution (in which we neglect the probability that $\tau$ is negative), $\Psi(s)/\Psi_{o}$ is log-normal distributed. Using the properties of a log-normal distribution, the mean transmission is given by

[TABLE]

from which we compute the relevant non-classical cross section as

[TABLE]

For any sufficiently compact two-point correlation, the non-classical cross-section is initially equal to the atomic-mixing approximation and asymptotes to a reduced constant value (for a positively correlated medium). This reduction accounts for the effects of correlation in the medium; the particles that stream initially through high density regions of the medium tend to be absorbed and the average cross-section seen by the surviving particles tends to decrease. This asymptotic value was derived in [14, 21] for the case where the two-point correlation was approximately a Dirac delta function or $s\to\infty$. This result is far more general for arbitrary functions with the caveat that $\tau$ has a non-zero probability of being negative. This connection suggests that the full non-classical cross-section for various distributions of media may be used to derive homogenization models for lower fidelity modeling. Alternatively, [20] derived an alternative model for transmission statistics that guarantees strict positivity in the non-classical cross-section. Such a model may prove to be more accurate in future work, though it is analytically complex.

5 Example Problems

In this section, the solution of the GLBE is benchmarked against direct sampling of the stochastic medium solved with the discrete ordinates method. The atomic mixing solution is also shown for reference where the cross-section is the mean of the truncated Gaussian.

In the examples below, level symmetric quadrature was used for the angular discretization, and diamond differencing was used for the spatial discretization without a negative-flux fixup. The domain was taken to be a square of length $1$ on each side discretized into a uniform 20x20 grid, and the walls were purely absorbing. The time advancement was stopped when the increments in the integrated angular flux reached below $10^{-7}$, and the outer scattering loop was stopped when the residual in the angular flux reached below $10^{-6}$.

5.1 Uniform source

In this problem, there are no external sources, and a uniformly distributed source of magnitude $1$ is prescribed in the domain. The medium is a truncated Gaussian process with nominal parameters $\mu=5$, $\sigma^{2}=2$, and $l=0.3$ with a stepsize of $\Delta s=0.008$. A set of realizations of this process are shown in Figure 3. The $S_{8}$ level symmetric quadrature set was used.

Figure 4 shows the integrated scalar flux along the centerline of the domain for the medium without any scattering, comparing the GLBE solution with Monte Carlo sampling and the atomic mix approximation. The GLBE solution agrees closely with the Monte Carlo results, with some discrepancy in the middle of the domain. It is likely that higher resolution and more accurate discretizations would narrow the gap in the solutions. Figure 5 supports this, which shows the results for the same problem with a scattering ratio of $c=0.5$. The difference in the results in this case is very small, and the solutions are nearly identical. However, just over $6000$ sweeps were necessary to complete the solution of the GLBE compared to over $17000$ sweeps in the aggregated Monte Carlo samples, making the GLBE solution substantially more attractive for the same degree of accuracy.

5.2 Emitting surface

In this problem, an isotropically emitting source of magnitude $1$ is affixed to one of the surfaces in the domain, and no internal sources are present. The medium is a truncated Gaussian process with nominal parameters $\mu=1$, $\sigma^{2}=2$, and $l=0.5$ with a stepsize of $\Delta s=0.007$. Example fields for this problem are shown in Figure 6. The $S_{16}$ level symmetric quadrature set was used.

Figure 7 shows the integrated scalar flux for the emitting surface along the axis normal to the wall. Although some discretization artifacts are present, the results from the Monte Carlo sampling and GLBE solution are nearly identical, and show substantially higher transmission through the domain than the atomic mixing model predicts. Figure 8 shows the same problem for the case where $c=0.3$. In the scattering case, the aggregated number of GLBE sweeps was over $5000$, whereas over $11000$ were required for the Monte Carlo samples.

5.3 Acceleration

The pathlength synthetic acceleration strategy was tested in a problem involving a scattering ratio of $0.9$ in a truncated Gaussian process medium of mean $10$, variance $5$, and correlation length $0.5$. This setting is expected to be challenging for a source iteration strategy due to the high optical depth and scattering ratio. The high-resolution cross-section was evaluated with a pathlength stepsize of $\Delta s=0.00214$, and the low-resolution was evaluated with a fivefold coarsening of $\Delta s=0.0107$. Only one inner iteration was performed with the low-resolution cross-section. The problem considered was the uniform source problem described above.

The residual history is shown in Figure 9 for the source iteration and accelerated solvers. The initial residuals for the two are nearly identical, but the asymptotic convergence rate of the accelerated solver is nearly half that of the source iteration solver. This is due to the inner iteration adding a substantial error reduction at each step, effectively acting as an additional outer iteration. These results are promising, but likely the pathlength synthetic acceleration would be best combined with a transport synthetic acceleration or geometric multigrid approach to maximize the performance of the solver.

6 Conclusion

This work has demonstrated the applicability of the GLBE to arbitrary stochastic processes in bounded multidimensional problems, and that it may be solved efficiently using existing deterministic methods. The choice of implicit discretization and source iterations lends itself to discrete ordinates methods. The non-classical cross-section may be computed directly from Monte Carlo statistics or from analytical approaches if the mean transmission is known. The resulting numerical method was tested against Monte Carlo sampling of truncated Gaussian process media in a square enclosure and showed excellent agreement in cases with and without scattering. A strategy for accelerating the convergence of the iterative process was also demonstrated to improve the rate of convergence for a particular case with a high scattering ratio.

Although the example problems only considered two-dimensional media, it is easy to extend this approach to three-dimensional media and unstructured grids, lending itself to immediate use in more practical problems of interest. In order to fully open the GLBE to the application space of interest, it is necessary to develop and demonstrate versions that capture charged-particle transport and reactor criticality computations. In addition, the larger scale of these problems requires that more efficient implementations of the GLBE solver be available, including versions that implement adaptive time integration, synthetic accleration or preconditioning methods, and parallel sweeps. Examining the relative effectiveness of these methods is important and should be studied further.

The complexity of the physics scenarios of interest is also important to consider. While this study and others before it have considered homogeneous media and uniform sources, it is important to remember that many processes are inhomogeneous and anisotropic, and that problems of interest have source terms that are dependent on the stochastic process as well (such as blackbody radiation or fission neutrons). The question of how effective the GLBE solver presented here is for such problems needs to be determined for these cases as well.

Acknowledgments

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholle owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy of the United States Government. The author would also wish to thank Dr. Arlyn Antolak for his help and advice throughout the progress of this work, as well as the LDRD office at Sandia National Laboratories for supporting this work under project number 212586. This document is set for unclassified unlimited release under SAND2019-0528J.