Transport equations, numerical methods

Methods for solving integro-differential equations describing particle or radiation transport. The equations take the following form for stationary problems:

$$ \tag{1 } \Omega \nabla \phi + \Sigma \phi = \int\limits dv ^ \prime \int\limits d \Omega ^ \prime \phi w ( x, \Omega , \Omega ^ \prime , v, v ^ \prime ) + f, $$

where $ x = ( x, y, z) $, $ \Omega = ( \Omega _ {1} , \Omega _ {2} , \Omega _ {3} ) $ is a unit vector, $ \phi = \phi ( x, \Omega , v) $ is the particle flux at a point $ x $ for particles with velocity $ v \Omega $, and the positive functions $ \Sigma $ and $ w $ describe the interaction of the particles with matter, while $ f $ is the source. Two basic aspects are considered: 1) finding the solution to (1) in a (convex) domain $ D( x, y, z) $ such that at its boundary $ \Gamma $,

$$ \tag{2 } \phi ( x, \Omega , v) = 0 \ \textrm{ for } ( \Omega , n) < 0, $$

where $ n $ is the outward normal to $ \Gamma $; and 2) finding the largest eigenvalue $ \lambda _ {1} $, and the corresponding eigenfunction of (1)–(2), in which

$$ \tag{3 } f = \frac{1} \lambda \int\limits dv ^ \prime \int\limits d \Omega ^ \prime \phi w _ {1} ( x,\ \Omega , \Omega ^ \prime , v, v ^ \prime ). $$

Equation (1) contains six independent variables; as it is essentially higher dimensional, it must be reduced to simpler equations. In (1) and (3) one replaces the integral with respect to $ v ^ \prime $ by a quadrature formula containing $ N $ terms and assumes that the scattering is isotropic, which gives a system of so-called multi-velocity equations:

$$ \tag{4 } \Omega \nabla \phi _ {i} + \Sigma _ {i} \phi _ {i} = \sum_{j=1}^ { N } \Sigma _ {s} ^ {ij} \phi _ {i0} + f _ {i} ,\ \ i = 1 \dots N, $$

where

$$ \phi _ {i} = \phi _ {i} ( x, \Omega ),\ \ \phi _ {i0} = \frac{1}{4 \pi } \int\limits \phi _ {i} d \Omega ^ \prime $$

are the zero moments, while the coefficients $ \Sigma _ {i} $, $ \Sigma _ {s} ^ {ij} $ and $ f _ {i} $ are obtained by applying averaging methods and using the solutions to conjugate problems. For the problem (2) one similarly gets

$$ \tag{5 } f _ {i} = \frac{1} \lambda \chi _ {i} Q( \phi ) = \ \frac{1} \lambda \chi _ {i} \sum_{j=1}^ { N } \Sigma _ {f _ {j} } \phi _ {j0} . $$

For $ N= 1 $, one gets the single-velocity equation

$$ \tag{6 } \Omega \nabla \phi + \Sigma \phi = \Sigma _ {s} \phi _ {0} + f $$

for the function $ \phi = \phi ( x, \Omega ) $. Equation (6) takes the following form for a planar layer $ 0 \leq x \leq H $, where the solution depends only on one coordinate $ x $ and one angular variable $ \mu $, $ | \mu | \leq 1 $:

$$ \tag{7 } \mu l \frac{\partial \phi }{\partial x } + \phi = \ c \frac{1}{2} \int\limits _ {-1}^ { 1 } \phi ( x, \mu ^ \prime ) d \mu ^ \prime + f _ {1} ( x, \mu ), $$

where $ l = 1/ \Sigma $, $ c = \Sigma _ {s} / \Sigma $, $ f _ {1} = f/ \Sigma $. The characteristics of the left-hand side in (6) are all the straight lines $ x = x _ {0} + \xi \Omega $, $ x _ {0} \in D $, and along each of them, (6) takes the form

$$ \tag{8 } \frac{\partial \phi }{\partial \xi } + \Sigma \phi = \Sigma _ {s} \phi _ {0} + f. $$

If one makes the substitution $ u = ( \phi ( x, \Omega ) + \phi ( x, - \Omega ))/2 $ in (6), it becomes

$$ \tag{9 } [- l \Omega \nabla ] ^ {2} u + u = cu _ {0} + F. $$

The solution to (9) minimizes the quadratic Vladimirov functional

$$ \tag{10 } G( u) = ( l \Omega \nabla u , l \Omega \nabla u) + ( u, u) - ( cu _ {0} , u) + $$

$$ + \int\limits _ {\Gamma \times \Omega } | ( \Omega , u) | u ^ {2} d \Omega d \Gamma - 2( u, f ), $$

where

$$ ( u, v) = \int\limits _ { D } \int\limits _ \Omega uv dx d \Omega . $$

Let the boundary value problems be written in operator form:

$$ \tag{11 } L \phi = S \phi + f. $$

A characteristic feature of the transport problem, which is used in numerical algorithms, is that the value of $ L ^ {-} 1 \psi $ is found from a given $ \psi $ by a direct method involving the integration of (8) along the characteristics. On this basis, from (11) one obtains the Peierls integral equation

$$ \tag{12 } S \phi = SL ^ {-} 1 ( S \phi + f ) $$

for the zero moment $ S \phi $.

The method of spherical harmonics (a form of Galerkin's method) has been substantially developed for solving transport problems. An approximate solution (approximation $ P _ {n} $) is found in the form

$$ \tag{13 } \phi ^ {(} n) ( x, \Omega ) = \ \sum_{k=0}^ { n } ( 2k+ 1) \sum _ { i=-k} ^ { k } \phi _ {ki} ( x) Y _ {ki} ( \Omega ), $$

where $ \phi _ {ki} ( x) $ are unknown functions while $ Y _ {ki} ( \Omega ) $ are spherical harmonics of order $ k $. Substituting (13) into (6), multiplying the result by $ Y _ {ki} $ and integrating with respect to $ \Omega $ leads to a system of partial differential equations for $ \phi _ {ki} ( x) $( cf. Spherical harmonics, method of). In approximation $ P _ {1} $, the system takes the form

$$ \tag{14 } \mathop{\rm div} \phi _ {1} + \Sigma \phi _ {0} = f _ {0} ,\ \ \frac{1}{3} \nabla \phi _ {0} + \phi _ {0} + \Sigma \phi _ {1} = f _ {1,} $$

$$ \left . 2( \phi _ {1} , n) - \phi _ {0} \right | _ \Gamma = 0, $$

where $ \phi _ {0} = \phi _ {00} $, $ \phi _ {1} = ( \phi _ {11} , \phi _ {12} , \phi _ {13} ) $. For $ f _ {1} \equiv 0 $( 14) implies the diffusion approximation

$$ \tag{15 } - \mathop{\rm div} D \nabla \phi _ {0} + \Sigma \phi _ {0} = f _ {0} ,\ \ \left . 2D \frac{\partial \phi }{\partial n } - \phi _ {0} \right | _ \Gamma = 0, $$

where $ D = 1/( 3 \Sigma ) $, which is an elliptic problem whose solution can be found by variational or grid methods.

To solve one-dimensional cases, analytic methods have been developed based on expanding the solution in terms of generalized eigenfunctions. The Monte-Carlo method is used to find functionals in the solutions to complex multi-dimensional problems.

Finite-difference approximation methods are widely used for transport equations. For example, one can use a quadrature formula for $ D $ and replace (12) by a system of linear equations. One can approximate the integral operator in (4), (5), (6), or (8) by means of quadrature formulas for a sphere. The Gauss quadrature formula for a sphere is known up to the 35th algebraic order of accuracy. In the method of characteristics, a family of characteristics is drawn through each point in the spatial grid along directions corresponding to the nodes of the quadrature for a sphere, and the differential operator $ L $ in (8) is replaced by a difference one. The difference equations of the $ S _ {n} $ method are obtained by integrating (6) over a grid cell in the phase space on the assumption that the solution is linear in the independent variables within the cell. In Galerkin's method, the solution is sought in the form

$$ \tag{16 } \phi = \sum_{n=1}^ { N } g _ {n} ( \Omega ) \phi _ {n} ( x). $$

If the $ \phi _ {n} ( x) $ are given, one obtains a system of degenerate integral equations for the $ g _ {n} ( \Omega ) $; if the $ \phi _ {n} ( x) $ are functions of compact support, one obtains the finite-element method; and if the $ g _ {n} ( \Omega ) $ are given functions of compact support and (16) minimizes (10), one gets the so-called $ P _ {NJ} $ equations.

Iterative methods for solving difference transport problems have the specific feature that the convergence usually becomes slower as $ c \rightarrow 1 $ $ ( c \leq 1) $, and that to derive the next approximation $ \phi ^ {k+1}$ one uses only part of the information on a preceding approximation $ \phi ^ {k} $ of substantially-fewer dimensions — one stores and uses only the values of $ \phi _ {0} ^ {k} $. In iterative methods an intermediate operation (an operation $ K $) is often that of solving the following problem:

$$ \tag{17 } L \Phi ^ {k} = S \phi ^ {k} + f,\ \ \Phi _ {0} ^ {k} = \frac{1}{4 \pi } \int\limits \Phi ^ {k} d \Omega ^ \prime . $$

Then the error $ \phi - \Phi ^ {k} $ satisfies (11) with source $ S( \Phi ^ {k} - \phi ^ {k} ) $ which is, as well as the discrepancy, independent of $ \Omega $. This feature enables one to accelerate the convergence. Consider a periodic problem for (7) with constant coefficients, with a source that is even in $ \mu $, and let $ H = 2 \pi $. In this application, below the following iterative methods are considered. For (7) one constructs a grid with $ N $ nodes in $ x $ and $ M $ angular directions in $ \mu $. Let

$$ r( t) = t ^ {-} 1 { \mathop{\rm arc} \mathop{\rm tan} } t,\ \ 0 \leq r( t) \leq 1, $$

$$ \epsilon _ {0} ^ {k} = \phi _ {0} - \phi _ {0} ^ {k} = \sum _ { n } \epsilon _ {n} ^ {k} e ^ {inx} , $$

$$ \| \epsilon _ {0} ^ {k} \| = \max _ { n } | \epsilon _ {n} ^ {k} | . $$

For convergent iterative methods, $ \| \epsilon _ {0} ^ {k+1}\| \leq q \| \epsilon _ {0} ^ {k} \| $, where $ 0 \leq q < 1 $. Let $ P _ {0} $ be the price (number of operations) in the operation $ K $, while $ P $ is the price of a complete iteration and $ \Delta = P- P _ {0} $. The following relationships apply for the various methods.

1) Simple iteration: $ \phi _ {0} ^ {k+1}= \Phi _ {0} ^ {k} $, where $ \Delta = 0 $ and $ q = c $.

2) Lyusternik's method: For certain $ k $, one uses the simple iteration

$$ \phi _ {0} ^ {k+1} = \Phi _ {0} ^ {k} + ( \lambda _ {1} - 1) ^ {-} 1 ( \Phi _ {0} ^ {k} - \phi _ {0} ^ {k} ), $$

in which $ \lambda _ {1} > 1 $ is the largest eigenvalue of $ L \phi = \lambda S \phi $; then $ \Delta = O( N) $ and $ q = cr( l) $( $ q \rightarrow 1 $ for $ c \rightarrow 1 $, $ l \rightarrow 0 $).

3) The method of estimating iterative deviations: $ \phi _ {0} ^ {k+1} = \Phi _ {0} ^ {k} + W _ {0} ^ {k} $, where $ W ^ {k} $ is the solution of the equation

$$ l \mu \frac{\partial W ^ {k} }{\partial x } + ( 1- c) W ^ {k} = c( \Phi _ {0} ^ {k} - \phi _ {0} ^ {k} ); $$

then

$$ \Delta = O( NM),\ \ q = \max \left ( cr( l), \frac{\pi \sqrt 2 c ^ {2} }{12} \right ) $$

( $ q \rightarrow 1 $ as $ c \rightarrow 1 $, $ l \rightarrow 0 $).

4) The balancing-multiplier method: $ \phi _ {0} ^ {k+1} = \delta ^ {k} \Phi _ {0} ^ {k} $, where

$$ \delta ^ {k} = \frac{\int\limits _ { 0 } ^ { H } f _ {0} dx }{\int\limits _ { 0 } ^ { H } ( 1- c) \Phi _ {0} ^ {k} dx } . $$

Here $ \Delta = O( N) $ and $ q = cr( l) $( $ q \rightarrow 1 $ as $ c \rightarrow 1 $, $ l \rightarrow 0 $).

5) The mean-flux method (rebalance method):

$$ \phi ^ {k+1} ( x, \mu ) = ( 1+ v ^ {k} ( x)) \Phi ^ {k} , $$

where the function $ v ^ {k} $ is selected to minimize the functional (10) or else to minimize it in some finite-dimensional subspace: $ v ^ {k} = \sigma _ {i} a _ {i} \psi _ {i} $, and then the $ a _ {i} $ satisfy a certain system of equations.

6) The quasi-diffusion method:

$$ - l \frac{d}{dx} l \frac{d}{dx} D ^ {k} \phi _ {0} ^ {k+1} + ( 1- c) \phi _ {0} ^ {k+1} = f _ {0} , $$

where

$$ D ^ {k} = \frac{\int\limits _ {-1}^ { 1 } \Phi ^ {k} \mu ^ {2} d \mu }{\Phi _ {0} ^ {k} } ; $$

then $ \Delta = O( NM) $.

7) Splitting methods:

$$ ( I + \tau \Lambda _ {2} )( I + \tau \Lambda _ {1} ) \phi _ {k+1} = \ ( I + \tau \Lambda _ {2} )( I - \tau \Lambda _ {2} ) \phi ^ {k+2} \tau f, $$

where

$$ \Lambda _ {1} = I - \frac{c}{2} \int\limits _ {-1} ^ { 1 } ( \dots ) d \mu ,\ \ \Lambda _ {2} = l \mu \frac{d}{dx} . $$

The methods 4)–6) are non-linear, and their convergence may slow down as $ c \rightarrow 1 $ and $ l \rightarrow 0 $; method 7) requires the storage of $ \phi ^ {k} ( x, \mu ) $( $ q \rightarrow 1 $ as $ c \rightarrow 1 $, $ l \rightarrow 0 $).

8) $ KP $- methods: The correction $ W ^ {k} $ is determined as the solution in $ D $ of the boundary value problem

$$ \tag{18 } Q _ {n} W ^ {k} = P _ {n} c ( W ^ {k} + \Phi _ {0} ^ {k} - \phi _ {0} ^ {k} ), $$

where $ Q _ {n} , P _ {n} $ are second-order linear differential operators, and one puts $ \phi _ {0} ^ {k+1} = \Phi _ {0} ^ {k} W ^ {k} $. In one of the forms of the $ KP $- method, (18) takes the form

$$ \tag{19 } - \frac{g _ {k} }{3} \left ( l \frac{d}{dx} \right ) ^ {2} W ^ {k} + ( 1- c) W ^ {k} = \ c( HI _ {0} ^ {k} - \phi _ {0} ^ {k} ). $$

Then $ \Delta = O( N) $ for (19); $ q = 0.186 c $ for $ g _ {k} = 0.843 $, while for $ g _ {k} = ( 1+ y _ {k} )/2 $, where $ y _ {k} $ are the roots of the Jacobi polynomial $ P ^ {(- 1/2,2( N+ \beta )- 2/3) } ( y) $, $ \beta > 0 $, the geometric mean of $ q $ over $ N $ iterations is close to $ 0.15 c $. In the $ KP $- method, the convergence of the iterations does not slow down as $ c \rightarrow 1 $, $ l \rightarrow 0 $.

Seidel's iterative method is used to solve the multi-velocity problem (4):

$$ \tag{20 } \Omega \Delta \phi _ {i} ^ {k+1} + \Sigma _ {i} \phi _ {i} ^ {k+1} = \ \sum_{j=1}^ { i } \Sigma _ {s} ^ {ij} \phi _ {j0} ^ {k+1} + \sum _ { j=i+ 1} ^ { N } \Sigma _ {s} ^ {ij} \phi _ {j0} + f, $$

$$ i = 1 \dots N , $$

and the solution $ \phi _ {i} ^ {k+1} $ in each equation in (20) is found by an iterative method for the one-velocity equation.

To solve multi-velocity problems for the eigenvalues of (4) and (5), these two iterative loops are supplemented with a further outer one to find the maximal value $ \lambda = \lambda _ {1} $ and the corresponding eigenfunction $ \phi $. If $ x = Q( \phi ) $ and $ A = Q( L- S) ^ {-} 1 \chi $, then problem (4), (5) becomes

$$ \tag{21 } Ax = \lambda x. $$

To find $ \lambda _ {1} $ and $ \phi $, iterative methods are used with the Chebyshev parameters:

$$ \tag{22 } u ^ {k+1} = Ax ^ {k} - \beta _ {k+1} x ^ {k} ,\ \ x ^ {k+1} = \frac{u ^ {k+1} }{Q( u ^ {k+1} ) } , $$

where

$$ \tag{23 } \beta _ {k} = \frac{1}{2} ( M + m + ( M- m) \cos \pi \omega _ {k} ), $$

$ \omega _ {k} \in [ 0, 1] $ and $ M $ and $ m $ are parameters. One assumes that the spectrum is non-negative and first finds $ \lambda _ {1} $ and $ \lambda _ {2} $, which are the largest eigenvalues of (21), on the assumption that $ m = 0 $, $ M = a $, where $ a \geq 0 $ is a lower bound for $ \lambda _ {1} $, and takes $ \omega _ {k} $ as a $ T $- sequence (see below). The values of $ \lambda _ {1} $ and $ \lambda _ {2} $ are determined by the generalized Aitken method, which incorporates the shifts $ \beta _ {k} $. When $ \lambda _ {1} $ and $ \lambda _ {2} $ have been found, $ \phi $ is found from (22) and (23) with $ M = \lambda _ {2} $. The infinite $ T $- sequence is formed, correspondingly, from the specially ordered roots of the Chebyshev polynomials of the first kind $ T _ {2 \cdot 3 ^ {n} } ( \cos \pi \omega ) $, while the initial segment of the $ T $- segment of length $ 2 \cdot 3 ^ {n} $ consists of all numbers of the form $ ( 2j _ {k} - 1)/4 \cdot 3 ^ {n} $, $ 1 \leq j _ {k} \leq 2 \cdot 3 ^ {n} $:

$$ \frac{1}{4} , \frac{3}{4} , \frac{1}{12} , \frac{11}{12} , \frac{5}{12} , \frac{7}{12} , \frac{1}{36} ,\dots. $$

Any segment of the $ T $- sequence of length $ 2 \cdot 3 ^ {n} $ ensures an optimal suppression in a certain sense of the error and stability in the iterative method (22), (23).

The following are used to solve non-stationary problems

$$ \frac{1} \nu \frac{\partial \phi }{\partial t } + L \phi - S \phi = f : $$

the method of characteristics in $ ( x, t) $- space, Galerkin's method, and finite-difference methods amounting to explicit and implicit difference schemes or to operator splitting methods. In the case of implicit schemes, the solution on the upper layer may be found by the $ KP $- method.

References

[1]	V.S. Vladimirov, "Mathematical methods of uniform-velocity transport theory" Trudy Mat. Inst. Steklov. , 61 (1961) pp. 1–158 (In Russian)
[2]	G.I. Marchuk, V.I. Lebedev, "Numerical methods in the theory of neutron transport" , Harwood (1986) (Translated from Russian)
[3]	V.I. Lebedev, S.A. Finogenov, "Utilization of ordered Chebyshev parameters in iterative methods" USSR Comp. Math. Math. Phys. , 16 : 4 (1976) pp. 70–83 Zh. Vychisl. Mat. i Mat. Fiz. , 16 : 4 (1976) pp. 895–907
[4]	V.I. Lebedev, "An iterative method with Chebyshev parameters for finding the maximum eigenvalue and corresponding eigenfunction" USSR Comp. Math. Math. Phys. , 17 : 1 (1977) pp. 92–101 Zh. Vychisl. Mat. i Mat. Fiz. , 17 : 1 (1977) pp. 100–108

References

[a1]	B. Davison, "Neutron transport theory" , Oxford Univ. Press (1957)
[a2]	G.J. Bell, S. Glasstone, "Nuclear reactor theory" , v. Nostrand (1971)