A Spectral Approach for Solving the Nonclassical Transport Equation

TL;DR

This paper presents a spectral method using Laguerre polynomials to numerically solve the nonclassical transport equation, enabling deterministic solutions for particle transport in complex systems.

Contribution

It introduces a novel spectral approach that transforms the nonclassical transport equation into a classical form using Laguerre polynomial expansion.

Findings

Validated the spectral method with numerical results in slab geometry.

Demonstrated effectiveness for both classical and nonclassical transport problems.

Provided a new deterministic solution technique for nonclassical particle transport.

Abstract

This paper introduces a mathematical approach that allows one to numerically solve the nonclassical transport equation in a deterministic fashion using classical numerical procedures. The nonclassical transport equation describes particle transport for random statistically homogeneous systems in which the distribution function for free-paths between scattering centers is nonexponential. We use a spectral method to represent the nonclassical flux as a series of Laguerre polynomials in the free-path variable $s$, resulting in a nonclassical equation that has the form of a classical transport equation. We present numerical results that validate the spectral approach, considering transport in slab geometry for both classical and nonclassical problems in the discrete ordinates formulation.

Tables

Table 1. Table 1 : Classical scalar fluxes for the slab-geometry classical transport problem.

$\left|\text{x}\right|$	Classical transport equation	Nonclassical transport equation			Relative error
(cm)	(neutrons/cm2s)	(neutrons/cm2s)			(%)
		$M$=0	$M$=50	$M$=100	$M$=0	$M$=50	$M$=100
$c$ = 0.0
0.0	6.734502E+16	6.734502E+16	6.734502E+16	6.734502E+16	-2.061153E-07	8.305695E-07	-2.316133E-09
1.0	1.2672661E+16	1.2672661E+16	1.2672661E+16	1.2672661E+16	-2.0611535E-07	-6.5724395E-07	-1.4390150E-06
2.0	2.664948E+15	2.664948E+15	2.664948E+15	2.664948E+15	-2.061155E-07	-2.378777E-06	2.347877E-06
3.0	6.997904E+14	6.997904E+14	6.997902E+14	6.997905E+14	-2.061153E-07	-2.423450E-05	6.092642E-06
4.0	2.011571E+14	2.011571E+14	2.011571E+14	2.011571E+14	-2.061156E-07	2.127259E-05	3.275557E-05
5.0	6.084981E+13	6.084981E+13	6.084979E+13	6.084978E+13	-2.061152E-07	-4.467129E-05	-5.366992E-05
6.0	1.902634E+13	1.902634E+13	1.902632E+13	1.902633E+13	-2.061153E-07	-8.692087E-05	-5.532571E-05
7.0	6.089504E+12	6.089504E+12	6.089505E+12	6.089537E+12	-2.061152E-07	1.857005E-05	5.378618E-04
8.0	1.983150E+12	1.983150E+12	1.983146E+12	1.983127E+12	-2.061151E-07	-1.905777E-04	-1.142114E-03
9.0	6.546034E+11	6.546034E+11	6.545958E+11	6.546020E+11	-2.061152E-07	-1.166810E-03	-2.140537E-04
10.0	2.184099E+11	2.184099E+11	2.184053E+11	2.184060E+11	-2.061152E-07	-2.118095E-03	-1.786618E-03
$c=0.9$
0.0	2.719908E+17	2.719907E+17	2.719908E+17	2.719908E+17	-1.916088E-05	2.227138E-06	1.924450E-06
1.0	1.537115E+17	1.537114E+17	1.537115E+17	1.537115E+17	-3.315645E-05	3.321180E-06	3.264306E-06
2.0	8.680613E+16	8.680608E+16	8.680613E+16	8.680613E+16	-5.609675E-05	5.742565E-06	5.833409E-06
3.0	5.065063E+16	5.065058E+16	5.065063E+16	5.065064E+16	-8.992331E-05	8.735536E-06	9.487710E-06
4.0	2.978408E+16	2.978404E+16	2.978408E+16	2.978408E+16	-1.396228E-04	1.485754E-05	1.487591E-05
5.0	1.754675E+16	1.754672E+16	1.754676E+16	1.754676E+16	-2.101284E-04	2.128752E-05	2.129112E-05
6.0	1.032058E+16	1.032055E+16	1.032059E+16	1.032059E+16	-3.049518E-04	3.112824E-05	3.108565E-05
7.0	6.024440E+15	6.024414E+15	6.024442E+15	6.024443E+15	-4.229096E-04	4.217371E-05	4.354501E-05
8.0	3.433063E+15	3.433044E+15	3.433065E+15	3.433065E+15	-5.529738E-04	5.372110E-05	5.202548E-05
9.0	1.805655E+15	1.805643E+15	1.805656E+15	1.805656E+15	-6.683423E-04	6.025153E-05	6.203616E-05
10.0	5.845608E+14	5.845566E+14	5.845611E+14	5.845611E+14	-7.046261E-04	5.564129E-05	5.684071E-05

Table 2. Table 2 : Parameters of nonclassical test problems.

Problem		Space domain limits	${\ell }_1$	${\ell }_2$	$x_1$	$x_2$
		(cm)	(cm)	(cm)	(cm)	(cm)
A	A1	[-9.0, 9.0]	0.5	1.0	-0.5	0.5
	A2	[-10.0, 10.0]	1.0	1.0	-0.5	0.5
	A3	[-9.0, 9.0]	1.0	0.5	-0.5	0.5
B	B1	[-9.0, 9.0]	0.5	1.0	-9.0	9.0
	B2	[-10.0, 10.0]	1.0	1.0	-10.0	10.0
	B3	[-9.0, 9.0]	1.0	0.5	-9.0	9.0

Table 3. Table 3 : Ensemble-Averaged scalar fluxes for Problem Set A.

Problem	$\left|\text{x}\right|$	Benchmark	Nonclassical transport equation		Relative error
	(cm)	(neutrons/cm2s)	(neutrons/cm2s)		(%)
			$M$=50	$M$=200	$M$=50	$M$=200
A1	$c$ = 0.0
	0.0	2.941938E+16	2.948194E+16	2.936095E+16	2.126596E-01	-1.985964E-01
	4.5	1.756486E+15	1.743672E+15	1.745583E+15	-7.295201E-01	-6.207207E-01
	9.0	2.252977E+14	2.187016E+14	2.198255E+14	-2.927735E+00	-2.428882E+00
	$c$ = 0.9
	0.0	9.745114E+16	9.595979E+16	9.586788E+16	-1.530351E+00	-1.624665E+00
	4.5	3.389651E+16	3.322013E+16	3.327105E+16	-1.995425E+00	-1.845201E+00
	9.0	7.610961E+15	7.650189E+15	7.647939E+15	5.154120E-01	4.858457E-01
A2	$c$ = 0.0
	0.0	3.890079E+16	3.891017E+16	3.889305E+16	2.411346E-02	-1.988391E-02
	5.0	6.598115E+14	6.660539E+14	6.591396E+14	9.460837E-01	-1.018436E-01
	10.0	2.999265E+13	3.028603E+13	2.994294E+13	9.781649E-01	-1.657549E-01
	$c$ = 0.9
	0.0	1.427736E+17	1.410333E+17	1.410050E+17	-1.218951E+00	-1.238745E+00
	5.0	3.186643E+16	3.218562E+16	3.214855E+16	1.001670E+00	8.853233E-01
	10.0	4.032572E+15	4.365284E+15	4.367367E+15	8.250617E+00	8.302263E+00
A3	$c$ = 0.0
	0.0	5.186772E+16	5.188694E+16	5.186813E+16	3.706624E-02	7.872449E-04
	4.5	4.585612E+14	4.616722E+14	4.585512E+14	6.784263E-01	-2.173120E-03
	9.0	1.251280E+13	1.249927E+13	1.248988E+13	-1.080800E-01	-1.831402E-01
	$c$ = 0.9
	0.0	1.908438E+17	1.896727E+17	1.896511E+17	-6.136817E-01	-6.249941E-01
	4.5	3.302245E+16	3.324623E+16	3.325332E+16	6.776650E-01	6.991184E-01
	9.0	3.180001E+15	3.312502E+15	3.312802E+15	4.166693E+00	4.176118E+00

Table 4. Table 4 : Ensemble-Averaged scalar fluxes for Problem Set B.

Problem	$\left|\text{x}\right|$	Benchmark	Nonclassical transport equation		Relative error
	(cm)	(neutrons/cm2s)	(neutrons/cm2s)		(%)
			$M$=50	$M$=200	$M$=50	$M$=200
B1	$c$ = 0.0
	0.0	9.893582E+16	9.878834E+16	9.880716E+16	-1.490681E-01	-1.300397E-01
	4.5	9.625624E+16	9.607356E+16	9.612430E+16	-1.897782E-01	-1.370743E-01
	9.0	4.998409E+16	4.986236E+16	4.986238E+16	-2.435295E-01	-2.434841E-01
	$c$ = 0.9
	0.0	7.286388E+17	7.136680E+17	7.138745E+17	-2.054628E+00	-2.026288E+00
	4.5	6.367297E+17	6.234173E+17	6.235945E+17	-2.090757E+00	-2.062916E+00
	9.0	2.273773E+17	2.296963E+17	2.297052E+17	1.019883E+00	1.023785E+00
B2	$c$ = 0.0
	0.0	9.990035E+16	9.989015E+16	9.989337E+16	-1.021726E-02	-6.988034E-03
	5.0	9.893674E+16	9.893792E+16	9.893049E+16	1.191124E-03	-6.325158E-03
	10.0	4.999981E+16	4.999523E+16	4.999523E+16	-9.160875E-03	-9.156552E-03
	$c$ = 0.9
	0.0	9.038937E+17	8.971970E+17	8.972188E+17	-7.408680E-01	-7.384594E-01
	5.0	8.050648E+17	7.988039E+17	7.987906E+17	-7.776922E-01	-7.793447E-01
	10.0	2.386542E+17	2.485283E+17	2.485310E+17	4.137392E+00	4.138536E+00
B3	$c$ = 0.0
	0.0	9.996817E+16	9.996753E+16	9.996782E+16	-6.494109E-04	-3.513593E-04
	4.5	9.946734E+16	9.946217E+16	9.946661E+16	-5.197731E-03	-7.356050E-04
	9.0	4.999998E+16	4.999966E+16	4.999966E+16	-6.286861E-04	-6.282565E-04
	$c$ = 0.9
	0.0	9.431169E+17	9.410914E+17	9.410939E+17	-2.147673E-01	-2.145026E-01
	4.5	8.560157E+17	8.529432E+17	8.529691E+17	-3.589333E-01	-3.558992E-01
	9.0	2.396932E+17	2.445678E+17	2.445680E+17	2.033667E+00	2.033756E+00

Equations

$$dp={\Sigma }_t(\mathbf{x},E)ds$$

$$\displaystyle s=\begin{array}[]{l}\text{distance traveled by the particle since}\\ \text{its last interaction (birth or scattering).}\end{array}$$

$$\frac{\partial }{\partial s}$$

$$\delta (s)\left[c{\int }_{4\pi }{\int }_0^{\infty }P({\mathbf{\Omega }}^{\prime }\cdot \mathbf{\Omega }){\Sigma }_t({\mathbf{\Omega }}^{\prime },s^{\prime })\Psi (\mathbf{x},{\mathbf{\Omega }}^{\prime },s^{\prime })d{\Omega }^{\prime }d{s}^{\prime }+\frac{Q(\mathbf{x})}{4\pi }\right],\mathbf{x}\in V,\mathbf{\Omega }\in 4\pi ,0<s,$$

$$\Psi (\mathbf{x},\mathbf{\Omega },s)={\Psi }^b(\mathbf{x},\mathbf{\Omega })\delta (s),\mathbf{x}\in \partial V,\mathbf{n}\cdot \mathbf{\Omega }<0,0<s.$$

$$p(\mathbf{\Omega },s)={\Sigma }_t(\mathbf{\Omega },s)e^{-{\int }_0^s{\Sigma }_t(\mathbf{\Omega },s^{\prime })d{s}^{\prime }},$$

$$\mathbf{\Omega }\cdot \nabla {\Psi }_c(\mathbf{x},\mathbf{\Omega })+{\Sigma }_t{\Psi }_c(\mathbf{x},\mathbf{\Omega })$$

$$\mathbf{x}\in V,\mathbf{\Omega }\in 4\pi ,$$

$${\Psi }_c(\mathbf{x},\mathbf{\Omega })$$

$${\Psi }_c(\mathbf{x},\mathbf{\Omega })={\int }_0^{\infty }\Psi (\mathbf{x},\mathbf{\Omega },s)ds.$$

$$\frac{\partial }{\partial s}\Psi (\mathbf{x},\mathbf{\Omega },s)+\mathbf{\Omega }\cdot \nabla \Psi (\mathbf{x},\mathbf{\Omega },s)+{\Sigma }_t(\mathbf{\Omega },s)\Psi (\mathbf{x},\mathbf{\Omega },s)=0,$$

$$\Psi (\mathbf{x},\mathbf{\Omega },0)=c{\int }_{4\pi }{\int }_0^{\infty }P({\mathbf{\Omega }}^{\prime }\cdot \mathbf{\Omega }){\Sigma }_t({\mathbf{\Omega }}^{\prime },s^{\prime })\Psi (\mathbf{x},{\mathbf{\Omega }}^{\prime },s^{\prime })d{\Omega }^{\prime }d{s}^{\prime }+\frac{Q(\mathbf{x})}{4\pi },$$

$$\Psi (\mathbf{x},\mathbf{\Omega },s)\equiv \psi (\mathbf{x},\mathbf{\Omega },s)e^{-{\int }_0^s{\Sigma }_t(\mathbf{\Omega },s^{\prime })d{s}^{\prime }}.$$

$$\frac{\partial }{\partial s}\psi (\mathbf{x},\mathbf{\Omega },s)+\mathbf{\Omega }\cdot \nabla \psi (\mathbf{x},\mathbf{\Omega },s)=0,$$

$$\psi (\mathbf{x},\mathbf{\Omega },0)=S(\mathbf{x},\mathbf{\Omega })+\frac{Q(\mathbf{x})}{4\pi },$$

$$\psi (\mathbf{x},\mathbf{\Omega },s)={\Psi }^b(\mathbf{x},\mathbf{\Omega })\delta (s)e^{{\int }_0^s{\Sigma }_t(\mathbf{\Omega },s^{\prime })d{s}^{\prime }},\mathbf{x}\in \partial V,\mathbf{n}\cdot \mathbf{\Omega }<0,$$

$$S(\mathbf{x},\mathbf{\Omega })=c{\int }_{4\pi }{\int }_0^{\infty }P({\mathbf{\Omega }}^{\prime }\cdot \mathbf{\Omega })p({\mathbf{\Omega }}^{\prime },s^{\prime })\psi (\mathbf{x},{\mathbf{\Omega }}^{\prime },s^{\prime })d{\Omega }^{\prime }d{s}^{\prime }$$

$$\psi (\mathbf{x},\mathbf{\Omega },s)=\sum _{m=0}^{\infty }{\psi }_m(\mathbf{x},\mathbf{\Omega })L_m(s),$$

$$\mathbf{\Omega }\cdot \nabla {\psi }_m(\mathbf{x},\mathbf{\Omega })=\sum _{j=m+1}^{\infty }{\psi }_j(\mathbf{x},\mathbf{\Omega }),m=0,1,2,...,$$

$${\psi }_m(\mathbf{x},\mathbf{\Omega })={\Psi }^b(\mathbf{x},\mathbf{\Omega }),\mathbf{x}\in \partial V,\mathbf{n}\cdot \mathbf{\Omega }<0,m=0,1,2,....$$

$$\sum _{j=m+1}^{\infty }{\psi }_j(\mathbf{x},\mathbf{\Omega })$$

$$\mathbf{\Omega }\cdot \nabla {\psi }_m(\mathbf{x},\mathbf{\Omega })+{\psi }_m(\mathbf{x},\mathbf{\Omega })=S(\mathbf{x},\mathbf{\Omega })+\frac{Q(\mathbf{x})}{4\pi }-\sum _{j=0}^{m-1}{\psi }_j(\mathbf{x},\mathbf{\Omega }),m=0,1,2,...,$$

$${\psi }_m(\mathbf{x},\mathbf{\Omega })={\Psi }^b(\mathbf{x},\mathbf{\Omega }),\mathbf{x}\in \partial V,\mathbf{n}\cdot \mathbf{\Omega }<0,m=0,1,2,....$$

$${\Psi }_c(\mathbf{x},\mathbf{\Omega })={\int }_0^{\infty }\Psi (\mathbf{x},\mathbf{\Omega },s)ds={\int }_0^{\infty }{e}^{-{\int }_0^s{\Sigma }_t(\mathbf{\Omega },s^{\prime })d{s}^{\prime }}\sum _{m=0}^{\infty }{\psi }_m(\mathbf{x},\mathbf{\Omega })L_m(s)ds.$$

$$S(\mathbf{x},\mathbf{\Omega })=c{\int }_{4\pi }P({\mathbf{\Omega }}^{\prime }\cdot \mathbf{\Omega })\sum _{k=0}^{\infty }{\psi }_k(\mathbf{x},{\mathbf{\Omega }}^{\prime })\left[{\int }_0^{\infty }p({\mathbf{\Omega }}^{\prime },s^{\prime })L_k(s^{\prime })d{s}^{\prime }\right]d{\Omega }^{\prime },$$

$${\int }_0^{\infty }p({\mathbf{\Omega }}^{\prime },s^{\prime })L_0(s^{\prime })d{s}^{\prime }$$

$${\int }_0^{\infty }p({\mathbf{\Omega }}^{\prime },s^{\prime })L_1(s^{\prime })d{s}^{\prime }$$

$${\int }_0^{\infty }p({\mathbf{\Omega }}^{\prime },s^{\prime })L_2(s^{\prime })d{s}^{\prime }$$

$$\mu \frac{\partial }{\partial x}{\psi }_m(x,\mu )+{\psi }_m(x,\mu )=S(x)+\frac{Q(x)}{2}-\sum _{j=0}^{m-1}{\psi }_j(x,\mu ),m=0,1,\dots M,$$

$${\psi }_m(-X,\mu )=0,\mu >0,m=0,1,\dots ,M,$$

$${\psi }_m(X,\mu )=0,\mu <0,m=0,1,\dots ,M,$$

$$S(x)=\frac{c}{2}{\int }_{-1}^1{\int }_0^{\infty }{\Sigma }_t({\mu }^{\prime },s^{\prime })e^{-{\int }_0^{s^{\prime }}{\Sigma }_t({\mu }^{\prime },s^{\prime \prime })d{s}^{\prime \prime }}\sum _{k=0}^M{\psi }_k(x,{\mu }^{\prime })L_k(s^{\prime })d{s}^{\prime }d{\mu }^{\prime },$$

$${\Psi }_c(x,\mu )={\int }_0^{\infty }{e}^{-{\int }_0^s{\Sigma }_t(\mu ,s^{\prime })d{s}^{\prime }}\sum _{m=0}^M{\psi }_m(x,\mu )L_m(s)ds.$$

$${\mu }_l\frac{d}{dx}{\psi }_{m,l}(x)+{\psi }_{m,l}(x)=S(x)+\frac{Q(x)}{2}-\sum _{j=0}^{m-1}{\psi }_{j,l}(x),m=0,1,\dots M,l=1,2,\dots N,$$

$${\psi }_{m,l}(-X)=0,m=0,1,\dots ,M,l=1,2,\dots \frac{N}{2},$$

$${\psi }_{m,l}(X)=0,m=0,1,\dots ,M,l=\frac{N}{2}+1,\dots N,$$

$$S(x)=\frac{c}{2}\sum _{n=1}^N{\omega }_n\sum _{k=0}^M{\psi }_{k,n}(x){\int }_0^{\infty }{\Sigma }_{t_n}(s^{\prime })L_k(s^{\prime })e^{-{\int }_0^{s^{\prime }}{\Sigma }_{t_n}(s^{\prime \prime })d{s}^{\prime \prime }}d{s}^{\prime },$$

$${\Psi }_{c,l}(x)=\sum _{m=0}^M{\psi }_{m,l}(x){\int }_0^{\infty }{L}_m(s)e^{-{\int }_0^s{\Sigma }_{t_l}(s^{\prime })d{s}^{\prime }}ds.$$

$$\begin{array}{c}\frac{L}{2}{\int }_{-1}^1{\Sigma }_{t_n}\left(\left(z+1\right)\frac{L}{2}\right)L_k\left(\left(z+1\right)\frac{L}{2}\right)e^{-{\int }_0^{\left(z+1\right)\frac{L}{2}}{\Sigma }_{t_n}(s^{\prime \prime })d{s}^{\prime \prime }}dz=\\ \frac{L}{2}\sum _{i=1}^{G_{\ell }}{\Sigma }_{t_n}\left(\left(z_i+1\right)\frac{L}{2}\right)L_k\left(\left(z_i+1\right)\frac{L}{2}\right)e^{-{\int }_0^{\left(z_i+1\right)\frac{L}{2}}{\Sigma }_{t_n}(s^{\prime \prime })d{s}^{\prime \prime }}{\omega }_i,\end{array}$$

$$S(x)=\frac{cL}{4}\sum _{n=1}^N{\omega }_n\sum _{k=0}^M{\psi }_{k,n}(x)\sum _{i=1}^{G_{\ell }}{\Sigma }_{t_n}\left(\left(z_i+1\right)\frac{L}{2}\right)L_k\left(\left(z_i+1\right)\frac{L}{2}\right)e^{-{\int }_0^{\left(z_i+1\right)\frac{L}{2}}{\Sigma }_{t_n}(s^{\prime \prime })d{s}^{\prime \prime }}{\omega }_i$$

$${\psi }_{c,l}(x)=\frac{L}{2}\sum _{m=0}^M{\psi }_{m,l}(x)\sum _{i=1}^{G_{\ell }}{L}_k\left(\left(z_i+1\right)\frac{L}{2}\right)e^{-{\int }_0^{\left(z_i+1\right)\frac{L}{2}}{\Sigma }_{t_n}(s^{\prime \prime })d{s}^{\prime \prime }}{\omega }_i.$$

Full text

††journal: arXiv

1 Introduction

The theory of nonclassical particle transport, which describes processes in which a particle’s distance-to-collision is not exponentially distributed, has received increased attention in the last decade. It was originally proposed by Larsen [1] to describe measurements of photon path-length in the Earth’s cloudy atmosphere that could not be explained by classical radiative transfer (cf. [2]). The theory has been extended over the last few years [3, 4, 5, 6, 7] and has found applications in other areas, including neutron transport in certain types of nuclear reactors [8, 9, 10], computer graphics [11, 12, 13], and problems involving anomalous diffusion (cf. [14]). Moreover, a similar kinetic equation has been independently derived for the periodic Lorentz gas in a series of papers by Golse (cf. [15]) and by Marklof and Strömbergsson [16, 17, 18, 19].

In classical transport, the quantity $dp$ defined as

[TABLE]

represents the incremental probability that a particle located at spatial point ${\bm{x}}=(x,y,z)$ with energy $E$ will collide while traveling an incremental distance $ds$. Here, the cross section $\Sigma_{t}$ is independent of both (i) the particle’s direction of flight ${\bf\Omega}=(\Omega_{x},\Omega_{y},\Omega_{z})$, and (ii) the free-path $s$ defined as

[TABLE]

In other words, the instant immediately following the moment in which a particle scatters or is born, its free-path $s$ is set to 0.

The assumption that $\Sigma_{t}$ is independent of both ${\bf\Omega}$ and $s$ is generally valid when the spatial locations of the scattering centers in the system are uncorrelated (Poisson distributed). This leads to the Beer-Lambert law, with particle flux decreasing as an exponential function of $s$. However, a nonexponential attenuation law for the particle flux occurs in certain heterogeneous media in which the scattering centers in the system are spatially correlated [2, 8, 9, 10, 11, 13, 14, 15, 16, 17, 18, 19]. The theory of nonclassical particle transport [1, 3, 4, 5, 6, 7] was developed to address this class of problems. In particular, it assumes that $\Sigma_{t}$ is a function of both ${\bf\Omega}$ and $s$, which requires an extended phase space that includes the free-path $s$ as an extra independent variable.

The one-speed nonclassical transport equation with angular-dependent free-paths can be written as [5]

[TABLE]

The angular-dependent nonclassical total cross section $\Sigma_{t}({\bf\Omega},s)$ in Eq. 1.5a satisfies

[TABLE]

where $p({\bf\Omega},s)$ is the free-path distribution function in the direction ${\bf\Omega}$.

If classical transport takes place, $\Sigma_{t}$ is independent of both ${\bf\Omega}$ and $s$. In this case, the free-path distribution reduces to the exponential distribution $p(s)=\Sigma_{t}e^{-\Sigma_{t}s}$, and Eq. 1.5 reduce to the classical linear Boltzmann equation

[TABLE]

Numerical results for the nonclassical theory have been provided for diffusion-based approximations and for moment models in the diffusive regime [8, 3, 9, 21, 10, 22]. To our knowledge, numerical results for the nonclassical transport equation given by Eq. 1.5 are only available for problems in rod geometry [23, 24, 25]. This is in part due to the difficult task of estimating the nonclassical free-path distribution. Another reason is that the $s$-dependence of the equation adds to the numerical cost of finding a solution. For instance, if one interprets the free-path $s$ as a pseudo-time variable, a direct discretization will potentially lead to stiffness since the step size for $s$ will need to be smaller than the spatial step size [23], analogous to the Courant-Friedrichs-Lewy (CFL) condition. Moreover, the improper integral on the right-hand side of Eq. 1.5a will need to be approximated numerically.

The original contributions of this paper are as follows. We introduce a mathematical approach that allows one to numerically solve Eq. 1.5 in a deterministic fashion using classical numerical procedures. This approach uses a spectral method to represent the nonclassical flux as a series of Laguerre polynomials [26] in the variable $s$. The resulting equation has the form of a classical transport equation and can be solved by traditional methods. Furthermore, it gives one the option to analytically solve the improper integral on the right-hand side of Eq. 1.5a, depending on how hard it is to obtain explicit expressions for the moments of the free-path distribution. We also present numerical results that validate the spectral approach, considering transport in slab geometry for both classical and nonclassical problems. To our knowledge, this is the first time deterministic numerical results for Eq. 1.5 are provided for problems in slab geometry.

A summary of the remainder of this paper is given next. In Section 2, we present the spectral approach to the nonclassical transport equation. In Section 3, we describe the numerical methodology used to solve the transport problems in this paper. Numerical results are given in Section 4: in Section 4.1, we validate the proposed approach by solving a classical transport problem, and in Section 4.2 we present numerical results for nonclassical transport problems in a random periodic medium. We conclude with a discussion in Section 5.

2 Spectral Approach

Let us consider Eq. 1.5a in an equivalent “initial value” form [5]:

[TABLE]

and define $\psi$ such that

[TABLE]

We can now rewrite the nonclassical problem as

[TABLE]

and $p({\bf\Omega}^{\prime},s^{\prime})$ is the free-path distribution given in Eq. 1.6.

To apply the spectral method, we represent $\psi$ as a series of Laguerre polynomials [26] in $s$:

[TABLE]

and replace this ansatz into Eq. 2.3. The Laguerre polynomials $\{L_{m}(s)\}$ are orthogonal with respect to the weight function $e^{-s}$ and satisfy $\frac{d}{ds}L_{m}(s)=\left(\frac{d}{ds}-1\right)L_{m-1}(s)$ for $m>0$. Therefore, multiplying Eqs. 2.3a and 2.3c by $e^{-s}L_{m}(s)$ and operating on them by $\int_{0}^{\infty}(\cdot)ds$, we obtain

[TABLE]

Moreover, Eq. 2.3b yields

[TABLE]

where $S({\bm{x}},{\bf\Omega})=c\int_{4\pi}\int_{0}^{\infty}P({\bf\Omega}^{\prime}\cdot{\bf\Omega})p({\bf\Omega}^{\prime},s^{\prime})\textstyle{\sum_{k=0}^{\infty}\limits}\psi_{k}({\bm{x}},{\bf\Omega}^{\prime})L_{k}(s^{\prime})ds^{\prime}d{\bf\Omega}^{\prime}$.

Substituting Eq. 2.6 into Eq. 2.5, we have

[TABLE]

The nonclassical angular flux $\Psi({\bm{x}},{\bf\Omega},s)$ is recovered from Eqs. 2.2 and 2.4. The classical angular flux $\Psi_{c}({\bm{x}},{\bf\Omega})$ is obtained using Eq. 1.7c, such that

[TABLE]

Rearranging the terms in the scattering source such that

[TABLE]

we observe that the improper integral in brackets can be rewritten as a linear combination of the $k$ first raw moments of $p({\bf\Omega},s)$. For example, considering the first three Laguerre polynomials, one obtains

[TABLE]

where $\tilde{s}({\bf\Omega})$ and $\tilde{s^{2}}({\bf\Omega})$ are, respectively, the first (mean free-path) and second (mean squared free-path) moments of $p({\bf\Omega},s)$ in the direction ${\bf\Omega}$. In general, it is always preferable to perform the integral in $s$ analytically. However, $p({\bf\Omega},s)$ may be such that obtaining an explicit expression for the moments might be too challenging, or the moments may not even exist.

In order to obtain a numerical solution for Eq. 2.7 it is necessary to truncate the series of Laguerre polynomials. Moreover, in this paper, we have decided to perform the improper integral in Eq. 2.9 numerically, so that our numerical approach is as problem independent as possible. The resulting equations can be solved by traditional numerical methods as discussed in the next section.

3 Numerical Methodology

As the goal of this paper is to provide a first validation of the proposed method, we focus on obtaining numerical solutions to test problems in slab geometry with isotropic scattering and vacuum boundary conditions. Taking $M$ as the truncation order of the Laguerre expansion, we have:

[TABLE]

Next, we write Eq. 3.1 in the conventional discrete ordinates (SN) formulation [27]:

[TABLE]

Here, we have defined $\psi_{m,l}(x)$ as $\psi_{m}(x,\mu_{l})$, where the direction-of-motion variable $\mu$ has been discretized in $N$ discrete values $\mu_{l}$. Similarly, $\Sigma_{t_{l}}(s^{\prime})$ has been defined as $\Sigma_{t}(\mu_{l},s^{\prime})$, and the angular integral in Eq. 3.1d has been approximated by the angular quadrature formula with weights $\omega_{n}$ in Eq. 3.2d.

Furthermore, we take two steps to approximate the improper integrals in Eqs. 3.2d and 3.2e: (i) we truncate the upper limit to a finite number $L$ by neglecting the integral in the complementary range $(L,\infty)$; and (ii) we perform the simple linear transformation $z=\frac{2}{L}s^{\prime}-1$ to translate the interval $[0,L]$ into $[-1,1]$. We use the Gauss-Legendre quadrature [28] in Eq. 3.2d to write

[TABLE]

where $G_{\ell}$ is the order of the Gauss-Legendre quadrature.

A similar procedure is carried out for the improper integral in Eq. 3.2e. Equations 3.2d and 3.2e appear as

[TABLE]

and

[TABLE]

Equation 3.4 is used in the SN equations (3.2a), which are solved numerically using the conventional fine-mesh diamond difference (DD) method with the source iteration (SI) scheme [27]. In addition, the classical angular flux in direction $\mu_{l}$ is approximated by Eq. 3.5.

At this point we consider a spatial discretization grid on the slab, wherein each cell has width $h_{j}$, $j=1,2,\dots J$. Applying the operator $\frac{1}{h_{j}}\int_{x_{j-1/2}}^{x_{j+1/2}}(\cdot)dx$, with $h_{j}=x_{j+1/2}-x_{j-1/2}$, to the SN equations (3.2a), we obtain the conventional discretized spatial balance SN equations, where we define the cell-average quantity

[TABLE]

As with the DD method, we approximate $\psi_{m,l}(x)$ by piecewise continuous linear functions across the spatial grid. In other words, we consider the approximation

[TABLE]

By substituting Eq. 3.7 into the discretized balance equations, we obtain the nonclassical SN sweep equations

[TABLE]

with

[TABLE]

We use Eq. 3.8a to sweep from left to right $\left(\mu_{l}>0\right)$ and Eq. 3.8b to sweep from right to left $\left(\mu_{l}<0\right)$, updating $\bar{S}_{j},\,j=1:J$, and the summation terms, until a prescribed stopping criterion is satisfied. The stopping criterion adopted is that the relative deviations between two consecutive estimates of classical scalar fluxes in each point of the spatial discretization grid need to be smaller than or equal to a prescribed positive constant $\xi$.

4 Numerical Results

In this section we provide numerical results that validate the spectral approach and the proposed numerical methodology. First, we apply this approach to solve a classical transport problem in a homogeneous slab and show that it correctly estimates the scalar flux. Then, we proceed to solve a nonclassical transport problem in a one-dimensional, random periodic slab. For all numerical experiments in this section we adopt $\xi=1\times 10^{-6}$ for the stopping criterion, $N=20$, $G_{\ell}=300$, and truncate the $s$ range in the length $L$ of the slab. The spatial domain is discretized in each problem to yield $h_{j}=0.005$ cm, $j=1,2,...,J$.

4.1 Validation: Classical Transport

In order to validate the proposed methodology, we apply the numerical approach described in the previous section to solve a classical transport problem in slab geometry. The slab is composed of a homogeneous material (material 1), with total cross section $\Sigma_{t_{1}}=1$ cm*-1*, as depicted in Fig. 1. A neutron source $Q(x)=1\times 10^{17}$ neutrons/cm${}^{3}\cdot$s , with $-0.5$ cm $<x<0.5$ cm, is located in a region at the center of the slab. We emphasize that, as $\Sigma_{t_{1}}$ is independent of $\mu$ and $s$ in all spatial domain, the free-path distribution $p(\mu,s)$ reduces to the exponential distribution $p(s)=\Sigma_{t_{1}}e^{-\Sigma_{t_{1}}s}$.

We are interested in how accurately the nonclassical model predicts the classical scalar flux. To this end, we compare the results obtained with the solution of the classical linear Boltzmann equation in slab geometry. Table 1 displays the classical scalar fluxes for two different values of the scattering ratio $c$ and three different values of the truncation order of the Laguerre expansion $M$. One can see that, even for low values of $M$, the results generated by the nonclassical S20 transport equations using the proposed methodology are in close agreement with the results obtained by solving the classical S20 transport equations. This is expected and provides a first validation of the nonclassical spectral approach, since this choice of parameters means the nonclassical transport equation should reduce to the classical linear Boltzmann equation.

4.2 Nonclassical Transport in a Random Periodic Slab

Let us consider a one-dimensional physical system similar to the one introduced in [29], composed of two distinct materials periodically arranged. The period is given by $\ell=\ell_{1}+\ell_{2}$, with $\ell_{1}$ and $\ell_{2}$ representing the width of each material. Material 1 is a solid with $\Sigma_{t_{1}}=1$ cm*-1*, and material 2 is defined as void, i.e., $\Sigma_{t_{2}}=0$ cm*-1*. This periodic system is randomly placed in the infinite line $-\infty<x<\infty$, such that the probability $P_{i}$ of finding material $i\in\{1,2\}$ in a given point is $\ell_{i}/\ell$. Therefore, material parameters (such as the cross sections) are stochastic functions of space. Figure 2 illustrates this periodic system.

The ensemble-averaged free-path distribution for the problem depicted in Fig. 2 has been analytically calculated in [25] to different material widths, with expressions for $p(\mu,s)$ given by

[TABLE]

Case 1: $\ell_{1}<\ell_{2}$

[TABLE]

Case 2: $\ell_{1}=\ell_{2}$

[TABLE]

Case 3: $\ell_{1}>\ell_{2}$

[TABLE]

where $n=0,1,2,...$ .

We perform numerical experiments considering the two sets of problems (A and B) displayed in Table 2. In each problem, we consider the existence of a neutron source defined as [25]

[TABLE]

where $x_{1}$ and $x_{2}$ are spatial points of the domain.

Here, we are interested in how accurately the nonclassical model predicts the ensemble-averaged scalar flux over all physical realizations of the random medium. To this end, we compare the nonclassical results against benchmark results obtained by averaging the solutions of the classical transport equation over a large number of physical realizations of this random system. The benchmarks were produced by solving $1/h_{j}$ classical transport problems. Details of how to obtain the benchmark solution can be found in [25].

Figure 3 depicts the ensemble-averaged scalar fluxes obtained for problem set A, with Fig. 4 showing the percentage relative error of the nonclassical approach with respect to the benchmark solutions. Similarly, Fig. 5 illustrates the ensemble-averaged scalar fluxes for problem set B, with corresponding percentage relative errors presented in Fig. 6. Explicit values of these results are given in Tables 3 and 4 for easy comparison.

The best accuracy in both sets of problems is obtained in purely absorbing systems, with $c=0$. In these cases, increasing the value of $M$ from $50$ to $200$ shows a clear improvement in the accuracy of the solution. This can be seen in the left-hand plots in Figs. 4 and 6 and from the values in Tables 3 and 4.

For this class of test problems, it has been shown [25] that the accuracy of the nonclassical model will deteriorate as the system becomes more diffusive, underestimating the maximum value located at $x=0$. This is confirmed by the numerical results; the estimates obtained for the cases with $c=0.9$ are more inaccurate than the ones generated for purely absorbing cases. Moreover, there is virtually no improvement in the accuracy of the diffusive solutions upon increasing $M$.

Finally, as $|x|$ increases, one can see a clear difference in terms of accuracy between solutions of problem sets A and B. The sinuous shape seen in the relative errors is a consequence of the periodic structure of the random systems; however, the amplitude of the errors is much larger in problem set A. This is due to the neutron sources in problem set B being inserted upon all spatial domain, which has two direct effects in the solution: (i) it increases the contribution of unscattered neutrons in the scalar flux, smoothing the error; and (ii) it decreases the propagation of boundary effects.

5 Discussion

In this paper, we have introduced a spectral approach that allows us to numerically solve the nonclassical transport equation using traditional (classical) methods. By representing the nonclassical flux as a series of Laguerre polynomials in the free-path variable $s$, we obtain a nonclassical equation that has the form of a classical transport equation. We describe a numerical methodology to solve this equation using a discrete ordinates ($S_{N}$) formulation in combination with the diamond difference method and a source iteration scheme. This was used to solve both classical and nonclassical transport problems in slab geometry, thus proving to be a good tool and validating the proposed approach. To our knowledge, this is the first time deterministic numerical results have been given for the nonclassical transport equation (1.5) in slab geometry.

It is important to note that the goal of this paper is not to investigate the accuracy of the nonclassical transport equation; this has been partly addressed elsewhere [25] and will be the subject of future work. Here, we are concerned with the development of an efficient approach to solve the nonclassical transport equation in a deterministic fashion–this work is a first step in that direction. We remark that the choice of using discrete ordinates or diamond differences is not binding; one can choose several different numerical methodologies to solve the nonclassical equations introduced here.

Modifications of the spectral approach presented in Section 2 are possible and may yield interesting results. For instance, we can modify Eq. 2.2 to define

[TABLE]

If we expand $\hat{\psi}$ by a truncated series of Laguerre polynomials,

[TABLE]

and perform the same steps as described in Section 2, we obtain

[TABLE]

The classical angular flux is given by

[TABLE]

Assuming classical transport ($\Sigma_{t}$ independent of $s$ and ${\bf\Omega}$), we can choose $\alpha=1-\Sigma_{t}$ and $M=0$ to rewrite Eq. 5.3 as

[TABLE]

In this case, Eq. 5.4 yields $\Psi_{c}({\bm{x}},{\bf\Omega})=\hat{\psi}_{0}({\bm{x}},{\bf\Omega})$; therefore, Eq. 5.5 represent the classical transport equation for the classical flux as given by Eq. 1.7. This demonstrates that, for classical transport problems, a proper choice of $\alpha$ and $M$ in Eq. 5.3 shall yield the correct solution without truncation errors in the variable $s$.

We remark that, in this paper, we have validated the proposed approach for problems in which the free-path distribution $p({\bf\Omega},s)$ is either an exponential representing classical transport (Section 4.1) or a nonexponential that decays in an exponential fashion as $s\rightarrow\infty$ (Section 4.2). In these cases, the improper integral in Eq. 2.9 can be written as a linear combination of the raw moments of $p({\bf\Omega},s)$. This choice of nonclassical problem for validation was made for two reasons: (i) it represents a transport problem whose benchmark solution can be obtained deterministically, and (ii) the free-path distribution is given analitically by Eq. 4.1.

Nevertheless, there are several important applications [6, 7] in which the free-path distribution has a large tail, decaying algebraically as $s\rightarrow\infty$:

[TABLE]

In these cases, the moments of order $k$ and larger will not exist. Therefore, the approximation

[TABLE]

where $L$ is a constant larger than the maximum chord length of the system, may compromise the efficiency of the method. We acknowledge this possibility and intend to perform a deeper investigation of the proposed method’s capabilities in addressing such problems. However, since benchmark solutions for this type of transport problem would need to be obtained through a Monte Carlo calculation, this task is beyond the scope of this purely deterministic paper and shall be pursued in future work.

Further work will also need to be done to investigate how well this approach performs in multi-dimensional nonclassical systems. In order to address that, the next steps of this work include: (i) implementing an acceleration scheme to improve the efficiency of the method; (ii) investigating the use of coarse-mesh methods, such as Response Matrix [30]; and (iii) performing a full convergence analysis of the method. We also intend to drop the periodic assumption and investigate results in more realistic random media; however, this will require a sophisticated numerical approach to estimate the free-path distribution $p({\bf\Omega},s)$.

Acknowledgements

R. Vasques acknowledges support under award number NRC-HQ-84-15-G-0024 from the Nuclear Regulatory Commission. R. N. Slaybaugh acknowledges support under award number NRC-HQ-84-14-G-0052 from the Nuclear Regulatory Commission. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. L. R. C. Moraes and R. C. Barros also would like to express their gratitude to the support of Conselho Nacional de Desenvolvimento Científico e Tecnológico - Brasil (CNPq) and Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro - Brasil (FAPERJ).