Double PN Benchmark Solution for the 1D Monoenergetic Neutron Transport Equation in Plane Geometry

Preprint

Article

Double PN Benchmark Solution for the 1D Monoenergetic Neutron Transport Equation in Plane Geometry

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

Barry D. Ganapol^*

This version is not peer-reviewed

Submitted:

17 June 2024

Posted:

18 June 2024

You are already at the latest version

Alerts

Abstract

Abstract: As more and more numerical and analytical solutions to the linear neutron transport equation become available, verification of numerical results is increasingly important. This presentation concerns the development of another benchmark for the linear neutron transport equation. There are numerous ways of solving the transport equation, such as the Wiener-Hopf method based on analyticity, method of singular eigenfunctions, Laplace and Fourier transforms and analytical discrete ordinates, which is arguably one of the most straightforward, to name a few. Another potential method is the PN method, where the solution is expanded in terms of full range orthogonal Legendre polynomials and with orthogonality and truncation, the moments form a set of second order ODEs. Because of the half-range boundary conditions for incoming particles however, full range Legendre expansions are inaccurate near material discontinuities. For this reason, a double PN (DPN) expansion is more appropriate, where the incoming and exit-ing flux distributions are expanded separately to preserve the discontinuity at material interfaces. Here, a new method of solution for the DPN equations is proposed and demonstrated for an iso-tropically scattering medium.

Keywords:

Subject: Physical Sciences - Mathematical Physics

1. Introduction

Boltzmann's equation of particle transport, indeed presents a significant challenge and noteworthy opportunity to solve because of its complexity and wide range of physical phenomena it describes. Originally, the non-linear integro-differential equation, as prescribed by kinetic theory of particle motion, was considered unsolvable. With time however, and advances in mathematics and physical applications, where, in some cases, non-linearity could be relaxed to give a linear equation, the situation changed. In the early to mid twentieth century, a flurry of analytical solutions were constructed for the linear and linearized Boltzmann equation primarily based on solving partial differential equations (PDEs) with distributions admitted, specifically for one-dimension. Alongside the development of analytical solutions were numerical solutions as well such as Monte Carlo, thus enabling the practical use of Boltzmann’s equation in nuclear physics experiments and nuclear reactor physics. With numerical solutions and applications to sensitive physics, there arose a need for numerical method’s verification, which led to the development of benchmarks and benchmarking. This, in turn, led to a host of numerical benchmark solutions to more relevant transport applications requiring sophisticated benchmarking techniques but still generally limited to model problems. Some readers may have the misconception that benchmarking comprehensive transport algorithms is a wasted exercise since only idealized cases can be treated; however, the opposite is true. Even the most advanced numerical method to solve the transport equation can contain unknown errors that a benchmark, even for a simple problem, can easily find. The following is about developing another high precision benchmark in one-dimension.

2. Proposed DPN Algorithm

The PN solution is a well-known solution to determine the angular flux as described by the transport equation. In the method, one expresses the angular flux to the transport equation in a 1D-plane slab as the full-range PN infinite series approximation. However, this approach is problematic for two fundamental reasons. The first concerns how the equations are closed since there is always one more unknown than equation. The second is how one best represents half-range boundary conditions by a continuous full-range expansion, which is not possible in an exact way. There have been several methodologies proposed on how to treat these outstanding issues [See Refs. 1-4]. The simplest approach to specifying boundary conditions is to avoid the difficulty altogether by resorting to the DPN moments approximation, where the expansions are split between forward and backward neutron directions. Now, the moments over the incoming fluxes on the slab near and far boundaries are exact in the limit. Since the DPN approximation also has the advantage of being a coupled set of ODEs for flux moments, one can solve for them following well-established methods. Furthermore, one way to close the system is through convergence acceleration as N approaches infinity.

2.1. DPN Moment and Parity Equations

We consider the simplest case of 1D plane geometry where scattering is isotropic, for which the transport equation is

[μ \frac{\partial}{\partial x} + 1] ϕ (x, μ) = \frac{ω}{2} \int_{- 1}^{1} d μ^{'} ϕ (x, μ^{'}) + Q (x, μ)

(1a)

\begin{array}{l} ϕ (0, μ) = f (μ) \\ ϕ (a, - μ) = 0. \end{array}

(1b)

ϕ (x, μ)

is the neutron angular flux, x is its position, μ is the neutron direction cosine, ω is the probability of scattering and Q is an imposed external source. The DPN approximation assumes the flux representation

ϕ (x, μ) = \{\begin{matrix} \sum_{l = 0}^{\infty} (2 l + 1) ϕ_{l}^{+} (x) P_{l} (2 μ - 1), μ \geq 0 \\ \sum_{l = 0}^{\infty} (2 l + 1) ϕ_{l}^{-} (x) P_{l} (2 μ + 1), μ \leq 0, \end{matrix}

(2a,b)

where

P_{l} (2 μ \pm 1)

is the half-range Legendre polynomial of degree l. When projected over half-range polynomials

P_{j} (2 μ \pm 1)

over intervals

[0, 1]

and

[- 1, 0]

respectively and from orthogonality

\int_{0}^{1} d μ P_{j} (2 μ \pm 1) P_{l} (2 μ \pm 1) = \frac{1}{2 l + 1} δ_{j, l},

(3)

the coefficients in Eqs(2) become the moments over the positive and negative directions,

ϕ_{j}^{\pm} (x) = \int_{(0, - 1)}^{(1, 0)} d μ P_{j} (2 μ \mp 1) ϕ (x, μ) .

(4)

The integral on the RHS of Eq(1a), called the scalar flux, in terms of moments is

ϕ (x) = \int_{- 1}^{1} d μ^{'} ϕ (x, μ^{'}) = \int_{0}^{1} d μ^{'} ϕ (x, μ^{'}) + \int_{- 1}^{0} d μ^{'} ϕ (x, μ^{'}) = ϕ_{0}^{+} (x) + ϕ_{0}^{-} (x)

(5)

and the transport equation to solve becomes

[μ \frac{\partial}{\partial x} + 1] ϕ (x, μ) = \frac{ω}{2} [ϕ_{0}^{+} (x) + ϕ_{0}^{-} (x)] + Q (x, μ) .

(6)

To find a recurrence relation for the moments, we first multiply Eq(6) over

P_{j} (2 μ \mp 1)

\begin{array}{l} [\frac{\partial}{\partial x} μ P_{j} (2 μ \mp 1) + P_{j} (2 μ \mp 1)] ϕ (x, μ) = \\ = \frac{ω}{2} [ϕ_{0}^{+} (x) + ϕ_{0}^{-} (x)] P_{j} (2 μ \mp 1) + Q (x, μ) P_{j} (2 μ \mp 1); \end{array}

(7a)^±

and from the recurrence of half-range Legendre polynomials, substitute

μ P_{j} (2 μ \mp 1) = \frac{1}{2} [\begin{array}{l} \frac{j + 1}{2 j + 1} P_{j + 1} (2 μ \mp 1) + \\ + \frac{j}{2 j + 1} P_{j - 1} (2 μ \mp 1) \pm P_{j} (2 μ \mp 1) \end{array}],

(7b)^±

to find

\begin{array}{l} \frac{1}{2} \frac{\partial}{\partial x} [\frac{j + 1}{2 j + 1} P_{j + 1} (2 μ \mp 1) + \frac{j}{2 j + 1} P_{j - 1} (2 μ \mp 1) \pm P_{j} (2 μ \mp 1)] ϕ (x, μ) + \\ + P_{j} (2 μ \mp 1) ϕ (x, μ) = \frac{ω}{2} [ϕ_{0}^{+} (x) + ϕ_{0}^{-} (x)] P_{j} (2 μ \mp 1) + Q (x, μ) P_{j} (2 μ \mp 1) . \end{array}

(7c)^±

Finally, on projection over intervals

[0, 1]

and

[- 1, 0]

respectively applying orthogonality

\begin{array}{l} \frac{1}{2} [\frac{j + 1}{2 j + 1} \frac{d ϕ_{j + 1}^{\pm} (x)}{d x} + \frac{j}{2 j + 1} \frac{d ϕ_{j - 1}^{\pm} (x)}{d x} \pm \frac{d ϕ_{j}^{\pm} (x)}{d x}] + ϕ_{j}^{\pm} (x) = \\ = \frac{ω}{2} [ϕ_{0}^{+} (x) + ϕ_{0}^{-} (x)] δ_{j, 0} + Q_{j}^{\pm} (x) . \end{array}

(7d)^±

To arrive at the final result, change j to l and multiply by

{(\pm 1)}^{l}

to give for l=0,1,…,N-1

\begin{array}{l} \frac{1}{2} [\frac{l + 1}{2 l + 1} \frac{d {\tilde{ϕ}}_{l + 1}^{\pm} (x)}{d x} + \frac{l}{2 l + 1} \frac{d {\tilde{ϕ}}_{l - 1}^{\pm} (x)}{d x} + \frac{d {\tilde{ϕ}}_{l}^{\pm} (x)}{d x}] \pm ϕ_{l}^{\pm} (x) = \\ = \pm \frac{ω}{2} [{\tilde{ϕ}}_{0}^{+} (x) + {\tilde{ϕ}}_{0}^{-} (x)] δ_{l, 0} \pm {\tilde{Q}}_{l}^{\pm} (x) . \end{array}

(7e)^±

when

{\tilde{ϕ}}_{l}^{\pm} (x) \equiv {(\pm 1)}^{l} ϕ_{l}^{\pm} (x),

(7f)^±

and for closure, we assume

\frac{d {\tilde{ϕ}}_{N + 1}^{\pm} (x)}{d x} \equiv 0.

(7g)^±

By defining the vectors

\begin{array}{l} {\tilde{ϕ}}^{\pm} (x) \equiv {[\begin{matrix} {\tilde{ϕ}}_{0}^{\pm} & {\tilde{ϕ}}_{1}^{\pm} & ... & {\tilde{ϕ}}_{N}^{\pm} \end{matrix}]}^{T} \\ {\tilde{Q}}^{\pm} (x) \equiv {[\begin{matrix} {\tilde{Q}}_{0}^{\pm} & {\tilde{Q}}_{1}^{\pm} & ... & {\tilde{Q}}_{N}^{\pm} \end{matrix}]}^{T}, \end{array}

(8a,b)^±

Eq(7e) becomes

\begin{array}{l} \frac{1}{2} A \frac{d {\tilde{ϕ}}^{\pm} (x)}{d x} \pm {\tilde{ϕ}}^{\pm} (x) = \pm \frac{ω}{2} [{\tilde{ϕ}}_{0}^{+} (x) + {\tilde{ϕ}}_{0}^{-} (x)] 1_{0} \pm {\tilde{Q}}^{\pm} (x), \end{array}

(9a)^±

where

\begin{array}{l} A \equiv \{δ_{l - 1, j} l / (2 l + 1), δ_{l, j}, δ_{l + 1, j} (l + 1) / (2 l + 1); l, j = 0, .., N - 1\} \\ 1_{0} \equiv {\{1 0 0 .... 0\}}^{T} . \end{array}

(9b,c)

Expressing the RHS more conveniently in terms of

ψ (x)

gives

\begin{array}{l} \frac{1}{2} A \frac{d {\tilde{ϕ}}^{\pm} (x)}{d x} \pm {\tilde{ϕ}}^{\pm} (x) = \pm \frac{ω}{2} δ [{\tilde{ϕ}}^{+} (x) + {\tilde{ϕ}}^{-} (x)] \pm {\tilde{Q}}^{\pm} (x), \end{array}

(10a)^±

where the matrix

δ

δ \equiv 1_{0} 1_{0}^{T} = \{δ_{0, 0}\} .

(10b)

The key to solving for the flux vector are the even-odd parity equations derived next.

2.2. The Even/Odd DPN Parity Equations

We form the even and odd parity moment vectors

\begin{array}{l} ψ (x) \equiv {\tilde{ϕ}}^{+} (x) + {\tilde{ϕ}}^{-} (x) \\ χ (x) \equiv {\tilde{ϕ}}^{+} (x) - {\tilde{ϕ}}^{-} (x) \end{array}

(11a,b)

with the flux moments recovered from

{\tilde{ϕ}}^{\pm} (x) = \frac{1}{2} [ψ (x) \pm χ (x)] .

(11c)^±

Similarly, the source parity moments are

q^{\pm} (x) \equiv Q^{+} (x) \pm Q^{-} (x) .

(11d)^±

The parity equations come from adding and subtracting Eqs(10a)^± to give

\begin{array}{l} \frac{1}{2} A \frac{d ψ (x)}{d x} + χ (x) = q^{-} (x) \\ \frac{1}{2} A \frac{d χ (x)}{d x} + ψ (x) = ω δ ψ (x) + q^{+} (x) \end{array}

(12a,b)

On combining Eqs(12) by differentiating Eq(12a)

\frac{1}{2} A \frac{d^{2} ψ (x)}{d x^{2}} + \frac{d χ (x)}{d x} = \frac{d q^{-} (x)}{d x} .

(13)

Then, multiplying by

A / 2

\frac{1}{4} A^{2} \frac{d^{2} ψ (x)}{d x^{2}} + \frac{1}{2} A \frac{d χ (x)}{d x} = \frac{1}{2} A \frac{d q^{-} (x)}{d x}

(14a)

and introducing Eq(12b)

\frac{1}{4} A^{2} \frac{d^{2} ψ (x)}{d x^{2}} - [I - ω δ] ψ (x) = \frac{1}{2} A \frac{d q^{-} (x)}{d x} - q^{+} (x)

(14b)

and simplifying

\frac{d^{2} ψ (x)}{d x^{2}} - Γ^{2} ψ (x) = 2 A^{- 1} \frac{d q^{-} (x)}{d x} - 4 A^{- 2} q^{+} (x),

(14c)

where

Γ^{2} \equiv 4 A^{- 2} [I - ω δ] .

(14d)

Once

ψ (x)

is found from Eq(14b), then from Eq(12a)

χ (x) = - \frac{1}{2} A \frac{d ψ (x)}{d x} + q^{-} (x) .

(15)

The two moments of boundary conditions for the solution of Eq(14b) are found from the incoming flux moments at the boundaries [see Eq(1b)]

\begin{array}{l} ϕ_{l}^{+} (0) = \int_{0}^{1} d μ P_{l} (2 μ - 1) f (μ) \\ ϕ_{l}^{-} (a) = 0. \end{array}

(16a,b)

The exiting moments will be given by the DPN solution designed to accommodate the half-range condition; therefore, the BCs for

ψ (x)

\begin{array}{l} ψ (0) \equiv {\tilde{ϕ}}^{+} (0) + {\tilde{ϕ}}^{-} (0) \\ ψ (a) \equiv {\tilde{ϕ}}^{+} (a), \end{array}

(17a,b)

are not fully determined. This is where the response matrix enters the analysis.

It should be stated that the approach taken thus far is not unique. Instead of a second order ODE for

ψ (x)

, a second order ODE can equally be found for

χ (x)

. However, there are additional technical issues with the second approach, and therefore will not be further considered.

3. Solution to the Parity Equations

Our task is to solve the following second order inhomogeneous ODE:

\frac{d^{2} ψ (x)}{d x^{2}} - Γ^{2} ψ (x) = Q (x),

(18a)

with inhomogeneous term

Q (x) \equiv 2 A^{- 1} \frac{d q^{-} (x)}{d x} - 4 A^{- 2} q^{+} (x) .

(18b)

The solution naturally decomposes into additive solutions to the homogeneous and the particular parts.

ψ (x) = ψ_{h} (x) + ψ_{p} (x) .

(19)

We begin with the homogeneous parity equation for

ψ_{h} (x)

\frac{d^{2} ψ_{h} (x)}{d x^{2}} - Γ^{2} ψ_{h} (x) = 0,

(20a)

and apply diagonalization to

Γ

Γ = T d i a g \{λ_{k}; k = 1, ..., N\} T^{- 1},

(20b)

where T is a matrix of size N whose columns are eigenvectors T_k of

Γ

corresponding to the eigenvalues

λ_{k}

. With diagonalization, Eq(20a) becomes

\frac{d^{2} y (x)}{d x^{2}} - d i a g \{λ_{k}^{2}\} y (x) = 0,

(21a)

with

y (x) \equiv T^{- 1} ψ (x) .

(21b)

Eqs(21) are scalar ODEs along the diagonal with the convenient solution [5]

y (x) \equiv d i a g \{\frac{\sin h (λ_{k} x)}{\sin h (λ_{k} a)}\} y (a) + d i a g \{\frac{\sin h (λ_{k} (a - x))}{\sin h (λ_{k} a)}\} y (0)

(22)

chosen to directly incorporate the (unknown) boundary conditions. By re-inserting y from Eq(21b), the solution to Eq(20a) is

ψ_{h} (x) = H (x) ψ_{h} (a) + H (a - x) ψ_{h} (0),

(23a)

where the matrix function

H (x)

forms the two solutions to the homogeneous part

H (x), H (a - x)

H (x) \equiv T d i a g \{\frac{\sin h (λ_{k} x)}{\sin h (λ_{k} a)}\} T^{- 1}

(23b)

H (a - x) \equiv T d i a g \{\frac{\sin h (λ_{k} (a - x))}{\sin h (λ_{k} a)}\} T^{- 1} .

(23c)

From Eq(19)

\begin{array}{l} ψ_{h} (0) = ψ (0) - ψ_{p} (0) \\ ψ_{h} (a) = ψ (a) - ψ_{p} (a), \end{array}

(24a)

we find for Eq(23a)

ψ (x) - ψ_{p} (x) = H (x) [ψ (a) - ψ_{p} (a)] + H (a - x) [ψ (0) - ψ_{p} (0)]

(24b)

ψ (x) = H (x) ψ (a) + H (a - x) ψ (0) + ψ_{p} (x) - ψ_{p} (a) - ψ_{p} (0) .

(24c)

Note that the solution contains the unknown boundary conditions.

With the chosen two solutions of Eqs(23b,c) to the homogeneous part and from variation of parameters, the particular solution is [5]

ψ_{p} (x) = W^{- 1} [H (x) \int_{x}^{a} d x^{'} H (a - x^{'}) Q (x^{'}) + H (a - x) \int_{0}^{x} d x^{'} H (x^{'}) Q (x^{'})]

(25a)

with Wronskian

W^{- 1} \equiv - T [d i a g \{\frac{s i n h (λ_{k} a)}{λ_{k}}\}] T^{- 1} .

(25b)

Note that the particular solution has been constructed such that

ψ_{p} (0) = ψ_{p} (a) = 0.

(26)

From Eq(24c)

ψ (x) = H (x) ψ (a) + H (a - x) ψ (0) + ψ_{p} (x);

(27a)

and from Eq(15)

χ (x) = - \frac{1}{2} A [H^{'} (x) ψ (a) + H^{'} (a - x) ψ (0)] + U (x),

(27b)

where

\begin{array}{l} H^{'} (x) \equiv T d i a g \{λ_{k} \frac{\cos h (λ_{k} x)}{\sin h (λ_{k} a)}\} T^{- 1} \\ H^{'} (a - x) \equiv - T d i a g \{λ_{k} \frac{\cos h (λ_{k} (a - x))}{\sin h (λ_{k} a)}\} T^{- 1} \\ U (x) \equiv q^{-} (x) - \frac{1}{2} A {ψ^{'}}_{p} (x) . \end{array}

(27c,d,e)

The flux moments in Eqs(2a,b) will be constructed from Eqs(27); but since the parity vectors

ψ (0), ψ (a)

\begin{array}{l} ψ (0) = {\tilde{ϕ}}^{+} (0) + {\tilde{ϕ}}^{-} (0) \\ ψ (a) = {\tilde{ϕ}}^{+} (a) + {\tilde{ϕ}}^{-} (a) \end{array}

(28)

contain the unknown exiting flux vectors

{\tilde{ϕ}}^{-} (0), {\tilde{ϕ}}^{+} (a)

, these must be determined first.

4. Response Matrix

To begin, we define

\hat{A}, \hat{B}

\begin{array}{l} {H^{'} (x)|}_{0} \equiv T d i a g \{\frac{λ_{k}}{\sin h (λ_{k} a)}\} T^{- 1} \equiv \hat{B} \\ {H^{'} (x)|}_{a} \equiv T d i a g \{λ_{k} \frac{\cos h (λ_{k} a)}{\sin h (λ_{k} a)}\} T^{- 1} \equiv - \hat{A}; \end{array}

(29a,b)

and note the following:

\begin{array}{l} {H^{'} (a - x)|}_{0} \equiv - T d i a g \{λ_{k} \frac{\cos h (λ_{k} (a))}{\sin h (λ_{k} a)}\} T^{- 1} = \hat{A} \\ {H^{'} (a - x)|}_{a} \equiv - T d i a g \{\frac{λ_{k}}{\sin h (λ_{k} a)}\} T^{- 1} = - \hat{B} . \end{array}

(29c,d)

To recover the outgoing positive and negative moments at the boundaries, one introduces x = 0 and a into Eq(27b) to give

\begin{array}{l} χ (0) = - \frac{1}{2} A [{H^{'} (x)|}_{0} ψ (a) + {H^{'} (a - x)|}_{0} ψ (0)] + U (0) \\ χ (a) = - \frac{1}{2} A [{H^{'} (x)|}_{a} ψ (a) + {H^{'} (a - x)|}_{a} ψ (0)] + U (a) \end{array}

(30a,b)

or from Eqs(28a,b)

\begin{array}{l} χ (0) = - \frac{1}{2} A [\hat{B} ψ (a) + \hat{A} ψ (0)] + U (0) \\ χ (a) = \frac{1}{2} A [\hat{A} ψ (a) + \hat{B} ψ (0)] + U (a) \end{array}

(31a,b)

Substituting Eqs(11a,b) at x = 0 and a respectively yields

\begin{array}{l} {\tilde{ϕ}}^{+} (0) - {\tilde{ϕ}}^{-} (0) = - \frac{1}{2} A \{\hat{B} [{\tilde{ϕ}}^{+} (a) + {\tilde{ϕ}}^{-} (a)] + \hat{A} [{\tilde{ϕ}}^{+} (0) + {\tilde{ϕ}}^{-} (0)] + U (0)\} \\ {\tilde{ϕ}}^{+} (a) - {\tilde{ϕ}}^{-} (a) = \frac{1}{2} A \{\hat{A} [{\tilde{ϕ}}^{+} (a) + {\tilde{ϕ}}^{-} (a)] + \hat{B} [{\tilde{ϕ}}^{+} (0) + {\tilde{ϕ}}^{-} (0)] + U (a)\} \end{array}

(32a,b)

and if

β \equiv \frac{A \hat{B}}{2} γ \equiv \frac{A \hat{A}}{2},

(32c,d)

then Eqs(32a,b) become

\begin{array}{l} {\tilde{ϕ}}^{+} (0) - {\tilde{ϕ}}^{-} (0) = - β [{\tilde{ϕ}}^{+} (a) + {\tilde{ϕ}}^{-} (a)] - γ [{\tilde{ϕ}}^{+} (0) + {\tilde{ϕ}}^{-} (0)] + U (0) \\ {\tilde{ϕ}}^{+} (a) - {\tilde{ϕ}}^{-} (a) = γ [{\tilde{ϕ}}^{+} (a) + {\tilde{ϕ}}^{-} (a)] + β [{\tilde{ϕ}}^{+} (0) + {\tilde{ϕ}}^{-} (0)] + U (a); \end{array}

(32e,f)

and re-arranging

\begin{array}{l} β {\tilde{ϕ}}^{+} (a) - [Ι - γ] {\tilde{ϕ}}^{-} (0) = - β {\tilde{ϕ}}^{-} (a) - [I + γ] {\tilde{ϕ}}^{+} (0) + U (0) \\ [Ι - γ] {\tilde{ϕ}}^{+} (a) - β {\tilde{ϕ}}^{-} (0) = [I + γ] {\tilde{ϕ}}^{-} (a) + β {\tilde{ϕ}}^{+} (0) + U (a) \end{array}

(32g,h)

into matrix form

[\begin{matrix} β & - x^{-} \\ x^{-} & - β \end{matrix}] [\begin{matrix} {\tilde{ϕ}}^{+} (a) \\ {\tilde{ϕ}}^{-} (0) \end{matrix}] = [\begin{matrix} - β & - x^{+} \\ x^{+} & β \end{matrix}] [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] + [\begin{matrix} U (0) \\ U (a) \end{matrix}],

(32i,j)

where

x^{\pm} \equiv [Ι \pm γ] .

(32k)

When we multiply both sides by the skew symmetric block identity matrix

[\begin{matrix} - I & 0 \\ 0 & I \end{matrix}] [\begin{matrix} β & - x^{-} \\ x^{-} & - β \end{matrix}] [\begin{matrix} {\tilde{ϕ}}^{+} (a) \\ {\tilde{ϕ}}^{-} (0) \end{matrix}] = [\begin{matrix} - I & 0 \\ 0 & I \end{matrix}] [\begin{matrix} - β & - x^{+} \\ x^{+} & β \end{matrix}] [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] + [\begin{matrix} - I & 0 \\ 0 & I \end{matrix}] [\begin{matrix} U (0) \\ U (a) \end{matrix}],

(33a)

we find

[\begin{matrix} - β & x^{-} \\ x^{-} & - β \end{matrix}] [\begin{matrix} {\tilde{ϕ}}^{+} (a) \\ {\tilde{ϕ}}^{-} (0) \end{matrix}] = [\begin{matrix} β & x^{+} \\ x^{+} & β \end{matrix}] [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] + [\begin{matrix} - U (0) \\ U (a) \end{matrix}] .

(33b)

Finally, multiplying by the inverse of the leading matrix gives the exiting flux

[\begin{matrix} {\tilde{ϕ}}^{+} (a) \\ {\tilde{ϕ}}^{-} (0) \end{matrix}] = R [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] + {[\begin{matrix} - β & x^{-} \\ x^{-} & - β \end{matrix}]}^{- 1} [\begin{matrix} - U (0) \\ U (a) \end{matrix}],

(34a)

where the response matrix is

R \equiv {[\begin{matrix} - β & x^{-} \\ x^{-} & - β \end{matrix}]}^{- 1} [\begin{matrix} β & x^{+} \\ x^{+} & β \end{matrix}] .

(34b)

To complete the expression

[\begin{matrix} - U (0) \\ U (a) \end{matrix}] = [\begin{matrix} - q^{-} (0) \\ q^{-} (a) \end{matrix}] + \frac{1}{2} [\begin{matrix} {ψ^{'}}_{p} (0) \\ - {ψ^{'}}_{p} (a) \end{matrix}],

(35a,b)

with

[\begin{matrix} {ψ^{'}}_{p} (0) \\ - {ψ^{'}}_{p} (a) \end{matrix}] = W^{- 1} [\begin{matrix} \hat{B} \int_{0}^{a} d x^{'} H (a - x^{'}) Q (x^{'}) \\ - \hat{A} \int_{0}^{a} d x^{'} H (x^{'}) Q (x^{'}) \end{matrix}] .

(35c,d)

On the left, we have the moments of the outgoing flux distribution and on the right the moments of the incoming flux distribution. Hence, for a known incoming distribution, the outgoing moments are now known.

5. Final Moments Solution

One then recovers the N spatial moments from Eqs(11c)^± as

\begin{array}{l} {\tilde{ϕ}}^{+} (x) = \frac{1}{2} [ψ (x) + χ (x)] \\ {\tilde{ϕ}}^{-} (x) = \frac{1}{2} [ψ (x) - χ (x)] . \end{array}

(36a)^±

and from Eqs(27a,b) to give the general solution

\begin{array}{l} {\tilde{ϕ}}^{+} (x) = \frac{1}{2} [\begin{array}{l} [H (x) - \frac{1}{2} A H^{'} (x)] ψ (a) + [H (a - x) - \frac{1}{2} A H^{'} (a - x)] ψ (0) + \\ + ψ_{p} (x) + U (x) \end{array}] \\ {\tilde{ϕ}}^{-} (x) = \frac{1}{2} [\begin{array}{l} [H (x) + \frac{1}{2} A H^{'} (x)] ψ (a) + [H (a - x) + \frac{1}{2} A H^{'} (a - x)] ψ (0) - \\ + ψ_{p} (x) - U (x) \end{array}] . \end{array}

(36b)^±

Or with

H^{\pm} (x) \equiv H (x) \pm \frac{1}{2} A H^{'} (x)

(37a)^±

and

P^{\pm} (x) \equiv ψ_{p} (x) \pm U (x),

(37b)^±

we have

\begin{array}{l} {\tilde{ϕ}}^{+} (x) = \frac{1}{2} [H^{-} (x) ψ (a) + H^{-} (a - x) ψ (0) + P^{+} (x)] \\ {\tilde{ϕ}}^{-} (x) = \frac{1}{2} [H^{+} (x) ψ (a) - H^{+} (a - x) ψ (0)] + P^{-} (x) . \end{array}

(37c)^±

With additional manipulation

[\begin{matrix} {\tilde{ϕ}}^{+} (x) \\ {\tilde{ϕ}}^{-} (x) \end{matrix}] = \frac{1}{2} [\begin{matrix} H^{-} (x) & H^{-} (a - x) \\ H^{+} (x) & H^{+} (a - x) \end{matrix}] \{[\begin{matrix} {\tilde{ϕ}}^{+} (a) \\ {\tilde{ϕ}}^{-} (0) \end{matrix}] + [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}]\} + \frac{1}{2} \{ψ_{p} (x) [\begin{matrix} 1 \\ 1 \end{matrix}] + U (x) [\begin{matrix} 1 \\ - 1 \end{matrix}]\},

(38a)^±

or introducing Eq(34a)

\begin{array}{l} [\begin{matrix} {\tilde{ϕ}}^{+} (x) \\ {\tilde{ϕ}}^{-} (x) \end{matrix}] = \frac{1}{2} [\begin{matrix} H^{-} (x) & H^{-} (a - x) \\ H^{+} (x) & H^{+} (a - x) \end{matrix}] \{[R + I] [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] + {[\begin{matrix} - β & x^{-} \\ x^{-} & - β \end{matrix}]}^{- 1} [\begin{matrix} - U (0) \\ U (a) \end{matrix}]\} + \\ + \frac{1}{2} \{ψ_{p} (x) [\begin{matrix} 1 \\ 1 \end{matrix}] + U (x) [\begin{matrix} 1 \\ - 1 \end{matrix}]\} \end{array}

(38b)^±

with

[\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] = [\begin{matrix} 0 \\ \int_{0}^{1} d μ P_{N} (2 μ - 1) f (μ) \end{matrix}]

(38c)^±

and

P_{N} (2 μ - 1) \equiv {[\begin{matrix} P_{0} (2 μ - 1) & P_{1} (2 μ - 1) & ... & P_{N} (2 μ - 1) \end{matrix}]}^{T}

(38d,e)^±

1 \equiv {\{11 .... 1\}}^{T} .

and

U (x)

is from Eq(27e).

6. The DPN Approximation

The DPN solution comes from the convergence of the series solution of Eqs(2)

ϕ (x, μ) = \lim_{N \to \infty} ϕ (x, μ; N),

(39a)

where the DPN approximation is the partial sum

ϕ (x, μ; N) = \{\begin{matrix} \sum_{l = 0}^{N} (2 l + 1) ϕ_{l}^{+} (x) P_{l} (2 μ - 1), μ \geq 0 \\ \sum_{l = 0}^{N} (2 l + 1) ϕ_{l}^{-} (x) P_{l} (2 μ + 1), μ \leq 0. \end{matrix}

(39b,c)

To incorporate the moments vector of Eqs(38b), we let μ negative be −|μ| and trivially multiply by (+1)^l in the partial sum for μ positive to give for

μ \geq 0

ϕ (x, \pm μ; N) = \sum_{l = 0}^{N} (2 l + 1) {\tilde{ϕ}}_{l}^{\pm} (x) P_{l} (2 |μ| - 1),

(40a)

and if

L_{N} \equiv d i a g \{2 l + 1, l = 1, ..., N\},

(40b)

Eq(40a) becomes the N-term inner product

ϕ (x, \pm μ; N) = P_{N}^{T} (2 |μ| - 1) L_{N} {\tilde{ϕ}}^{\pm} (x)

(40c)

to be evaluated in the next section.

7. Numerical Implementation and Demonstration

7.a. Numerical Implementation for an Isotropic Source

Numerical implementation of the DPN algorithm for an isotropic source at the near surface and none at the far surface is relatively straightforward by following the algorithm presented. In particular, matrix diagonalization is through an HQR procedure originating from the LAPACK routine [6] and is one of the most reliable of numerical methods in use today. Once eigenvalues are known, the matrix functions come from common matrix multiplication of Eqs(23b,c). LU decomposition gives the matrix inversions. Our final consideration is convergence in DPN order N through a sequence of DPN approximations [Eq(40c)] and convergence acceleration. To develop a sequence of partial sums, we increment the number of moments by stride ΔN (usually ΔN = 5) and monitor convergence of the sequence through sequential convergence and convergence acceleration until the relative error is below a desired relative error or convergence fails. If failure, we increase the number of sequence elements or change expectations. In this way, convergence provides a measure of closure. There are numerous variants of acceleration to choose from. Here, we first apply common sequential convergence similar to a sensitivity study and then Wynn-epsilon (W-e) acceleration. Sequential convergence starts from DPN order N₀ and compares the convergence standard, called the engineering estimate of the relative error between iterations m−1 and m [i.e., orders

N_{0} + (m - 1) Δ N

and

N_{0} + m Δ N

], for m=1,2,…

ε_{S e q} (x, μ; m) \equiv |\frac{ϕ (x, \pm μ; N_{0} + m Δ N) - ϕ (x, \pm μ; N_{0} + (m - 1) Δ N)}{ϕ (x, \pm μ; N_{0} + m Δ N)}|,

(41a)

to the desired relative error 

ε_{S e q} (x, μ; m) < ε .

(41b)

If satisfied, the sequence has converged sequentially at (x, μ). W-e acceleration, in contrast, is non-linear. Similarly, to sequential convergence, the same sequence of DPN approximations

s_{m} = ϕ (x, \pm μ; N_{0} + m Δ N), m = 0, 1, ...

(42a)

is the initial sequence list. W-e convergence essentially extrapolates a known sequence to give the next sequence element (and an estimate of the limit) from previous elements. The algorithm is recursive and for L+1 iterations, in principle, improved DPN partial sums result for k odd

\begin{array}{l} ε_{- 1}^{(m)} = 0 \\ ε_{0}^{(m)} = s_{m}, m = 0, 1, ..., L \\ ε_{k + 1}^{(m)} = ε_{k - 1}^{(m + 1)} + {[ε_{k}^{(m + 1)} - ε_{k}^{(m)}]}^{- 1}; m = 0, 1, ..., L - k - 1; k = 0, 1, ..., L - 1. \end{array}

(42b)

The algorithm is conveniently written as the following tableau:

Every element of the odd columns should give an improved estimate of the original DPN partial sum in column one. The last term in each column, indicated by the arrow, assumed most precise, gives the following relative error standard after L iterations:

e_{W e} (x, μ; L) \equiv |\frac{ε_{L}^{(0)} (x, μ) - ε_{L - 2}^{(2)} (x, μ)}{ε_{L}^{(0)} (x, μ)}|,

(44a)

considered converged when

ε_{W e} (x, μ; L) < ε; L = 1, 2, ....

(44b)

Overall convergence therefore is a competition between the two modes of partial sum convergence, where convergence occurs for the least relative error

m i n [ε_{S e q} (x, μ; m), ε_{W e} (x, μ; L)] < ε .

(45)

As a first demonstration, we consider an isotropic source, f(μ) = 1 entering the near surface and no incoming flux entering the far surface to give surface flux moments vectors

[\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] = [\begin{matrix} 0 \\ 1_{0} \end{matrix}] .

(46)

From Eqs(34a,b), the exiting flux moments vectors are therefore

[\begin{matrix} {\tilde{ϕ}}^{+} (a) \\ {\tilde{ϕ}}^{-} (0) \end{matrix}] = R [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] = {[\begin{matrix} - β & x^{-} \\ x^{-} & - β \end{matrix}]}^{- 1} [\begin{matrix} β & x^{+} \\ x^{+} & β \end{matrix}] [\begin{matrix} 0 \\ 1_{0} \end{matrix}]

(47a)

and the interior flux moment vectors are from Eqs(38b,c)

[\begin{matrix} {\tilde{ϕ}}^{+} (x) \\ {\tilde{ϕ}}^{-} (x) \end{matrix}] = \frac{1}{2} [\begin{matrix} H^{-} (x) & H^{-} (a - x) \\ H^{+} (x) & H^{+} (a - x) \end{matrix}] [R + I] [\begin{matrix} {\tilde{ϕ}}^{-} (a) \\ {\tilde{ϕ}}^{+} (0) \end{matrix}] .

(47b)

The flux vectors enter Eq(40a) re-written as

ϕ (x, μ; N) = \sum_{l = 0}^{N} a_{l}^{\pm} (x) θ_{l} (μ)

(48a)

with

\begin{array}{l} θ_{l} (μ) = P_{l} (2 |μ| - 1) \\ a_{l}^{\pm} (x) \equiv (2 l + 1) {\tilde{ϕ}}_{l}^{\pm} (x) . \end{array}

(48b)

Since

θ_{l} (μ)

obeys the recurrence

θ_{l + 1} (μ) = = (2 |μ| - 1) \frac{2 l + 1}{l} θ_{l} (μ) - \frac{l}{l + 1} θ_{l - 1} (μ),

(49a)

the Clenshaw sum [7] applies, where

ϕ (x, μ; N) = θ_{0} (μ) a_{0}^{\pm} (x) + θ_{1} (μ) b_{1} (x, μ) + β_{1} (μ) θ_{0} (μ) b_{2} (x, μ) .

(49b)

and

\begin{array}{l} α_{l} (μ) = (2 μ - 1) \frac{2 l + 1}{l} \\ β_{l} (μ) = - \frac{l}{l + 1} \end{array}

(49c)

\begin{array}{l} b_{l + 1} (x, μ) = b_{l + 2} (x, μ) = 0 \\ b_{l} (x, μ) = a_{l}^{\pm} (x) + α_{l} (μ) b_{l + 1} (x, μ) + β_{l + 1} (μ) b_{l + 2} (x, μ) . \end{array}

(49d)

We now turn to the issue of solution precision. To address precision, we compare the DPN solution to a well-established fully analytical discrete ordinates response matrix benchmark RM/DOM [5]. We assume no absorption (ω = 1) and a slab of 1 mfp thickness. Table 1 gives what is believed to be precise angular fluxes to better than one unit in the 8^th place using RM/DOM. The DPN angular flux approximation for N = 100 agrees to all 7- digits of the angular flux from the RM/DOM of quadrature order 250 and to all but three entries (last digit emboldened and underlined) in the 8^th place for N = 150 without acceleration. Thus, the method presented indeed does successfully provide an extreme benchmark to nearly 8 places.

The effectiveness of W-e acceleration is shown by Figure 1, which is the ratio of the relative error with and without acceleration over all directions at the seven spatial coordinates. One observes the W-e relative errors at convergence are generally smaller than the errors of the original sequence for this benchmark solution. This is further confirmation of the high quality-- to one unit in the seventh and nearly one unit in the eighth place-- of the proposed DPN algorithm. 40 of the 72 fluxes converged by W-e showing the significance of the Wynn-epsilon algorithm.

7.b. Numerical Implementation for a Beam Source

The final example is for a beam source entering the near surface and none at the far surface

\begin{array}{l} ϕ (0, μ) = f (μ) = δ (μ - μ_{0}) \\ ϕ (a, - μ) = 0. \end{array}

(50a,b)

For this case, it is convenient to separate the uncollided

ϕ_{0} (x, μ)

from the collided

ϕ_{c} (x, μ)

component

ϕ (x, μ) = ϕ_{0} (x, μ) + ϕ_{c} (x, μ),

(51)

where he uncollided (scalar) component should not be confused with the zeroth vector flux moments

{\tilde{ϕ}}_{0}^{\pm} (x)

, a vector. The uncollided flux satisfies

[μ \frac{\partial}{\partial x} + 1] ϕ_{0} (x, μ) = 0

(52a)

\begin{array}{l} ϕ_{0} (0, μ) = δ (μ - μ_{0}) \\ ϕ_{0} (a, - μ) = 0 \end{array}

(52b)

to give

ϕ_{0} (x, μ) = δ (μ - μ_{0}) e^{- x / μ_{0}} Θ (μ)

(52c)

where

Θ (μ)

is the Heaviside step function required to maintain the uncollided flux in the positive direction. When Eq(51) is introduced into Eq(1a) with Eq(52c)

\begin{array}{l} [μ \frac{\partial}{\partial x} + 1] ϕ_{c} (x, μ) = \frac{ω}{2} \int_{- 1}^{1} d μ^{'} ϕ_{c} (x, μ^{'}) + Q (x) \\ Q (x) \equiv \frac{ω}{2} e^{- x / μ_{0}}, \end{array}

(53a,b)

with boundary conditions

\begin{array}{l} ϕ_{c} (0, μ) = 0 \\ ϕ_{c} (a, - μ) = 0. \end{array}

(53c,d)

From Eq(14c) for the collided component derived from Eq(53b)

\frac{d^{2} ψ_{c} (x)}{d x^{2}} - Γ^{2} ψ_{c} (x) = - 4 ω e^{- x / μ_{0}} A^{- 2} 1_{0}

(54a)

with

\begin{array}{l} ψ_{c} (0) = {\tilde{ϕ}}_{c}^{-} (0) \\ ψ_{c} (a) = {\tilde{ϕ}}_{c}^{+} (a) . \end{array}

(54b,c)

Assuming the particular solution

ψ_{p} (x) = e^{- x / μ_{0}} C,

(55a)

we find

C \equiv - 4 μ_{0}^{2} ω [I - μ_{0}^{2} Γ^{2}] A^{- 2} 1_{0} .

(55b)

From Eq(38b)

[\begin{matrix} {\tilde{ϕ}}_{c}^{+} (x) \\ {\tilde{ϕ}}_{c}^{-} (x) \end{matrix}] = \frac{1}{2} [\begin{matrix} H^{-} (x) & H^{-} (a - x) \\ H^{+} (x) & H^{+} (a - x) \end{matrix}] [\begin{matrix} {\tilde{ϕ}}_{c}^{+} (a) \\ {\tilde{ϕ}}_{c}^{-} (0) \end{matrix}] + \frac{1}{2} \{ψ_{p} (x) [\begin{matrix} 1 \\ 1 \end{matrix}] + \frac{e^{- x / μ_{0}}}{2 μ_{0}} A [\begin{matrix} 1 \\ - 1 \end{matrix}]\},

(56a)

where

[\begin{matrix} {\tilde{ϕ}}_{c}^{+} (a) \\ {\tilde{ϕ}}_{c}^{-} (0) \end{matrix}] = {[\begin{matrix} - β & x^{-} \\ x^{-} & - β \end{matrix}]}^{- 1} [\begin{matrix} - U (0) \\ U (a) \end{matrix}],

(56b)

and

[\begin{matrix} - U (0) \\ U (a) \end{matrix}] = \frac{1}{2} [\begin{matrix} {ψ^{'}}_{p} (0) \\ - {ψ^{'}}_{p} (a) \end{matrix}] = - \frac{1}{2 μ_{0}} [\begin{matrix} I & 0 \\ 0 & e^{- a / μ_{0}} I \end{matrix}] [\begin{matrix} C \\ - C \end{matrix}]

(56c)

to be inserted into

ϕ_{c} (x, \pm μ; N) = P_{N}^{T} (2 |μ| - 1) L_{N} {\tilde{ϕ}}_{c}^{\pm} (x)

(56d)

for the angular flux.

The slab for the beam source is the same as for the isotropic source. Table 2 shows the results from RM/DOM [5], which is a 7-place benchmark. Overlaid is the DPN result where it is observed that DPN gives nearly all 7 places except for two entries in one unit in the last digit. The beam source behaves differently from the isotropic source however, as it is more sensitive to round off error. The DPN results of Table 2 required quadruple precision (QP) to achieve nearly 7-places; whereas, the isotropic source required double precision. For the beam source, double precision only gives at best 4-places. A future effort will attempt to determine where the loss of precision occurs; nevertheless, the DPN provides a solid 6 digits and near 7 with QP.

Table 2. Angular Flux for a perpendicularly (μ₀ =1) entering beam source.

μ\x	0	0.05	0.1	0.2	0.5	0.75	1
-1.000E+00	5.3877491E-01	5.1979897E-01	4.9826415E-01	4.5015758E-01	2.8363970E-01	1.3670184E-01	0.0000000E+00
-8.000E-01	6.1358488E-01	5.9454278E-01	5.7227580E-01	5.2122894E-01	3.3675659E-01	1.6617467E-01	0.0000000E+00
-6.000E-01	7.0705901E-01	6.8953074E-01	6.6778458E-01	6.1558120E-01	4.1317580E-01	2.1164480E-01	0.0000000E+00
-4.000E-01	8.1805757E-01	8.0600276E-01	7.8820201E-01	7.4066647E-01	5.2986404E-01	2.9042486E-01	0.0000000E+00
-2.000E-01	9.1606674E-01	9.1832709E-01	9.1231716E-01	8.8438367E-01	7.1013555E-01	4.5375071E-01	0.0000000E+00
0.000E+00	8.7868708E-01	9.2867409E-01	9.4863473E-01	9.5773321E-01	8.6865051E-01	7.1731075E-01	4.8302802E-01
2.000E-01	0.0000000E+00	2.0129990E-01	3.6476856E-01	5.9744092E-01	8.4223675E-01	7.9950036E-01	6.5012804E-01
4.000E-01	0.0000000E+00	1.0687614E-01	2.0478323E-01	3.7093644E-01	6.5966190E-01	7.2022080E-01	6.6508565E-01
6.000E-01	0.0000000E+00	7.2711601E-02	1.4205915E-01	2.6699419E-01	5.2395290E-01	6.1548983E-01	6.1208033E-01
8.000E-01	0.0000000E+00	5.5092891E-02	1.0870690E-01	2.0824783E-01	4.3113588E-01	5.2849359E-01	5.4968216E-01
1.000E+00	0.0000000E+00	4.4345690E-02	8.8026420E-02	1.7060822E-01	3.6524695E-01	4.6023695E-01	4.9305877E-01

Conclusion

By expressing the solution to the 1D monoenergetic neutron transport equation in plane geometry in an infinite series of half-range Legendre polynomials, a first order coupled set of ODEs for half-range moments in the positive and negative directions follows on truncation. By adding and subtracting, the ODEs transform into a second order form for the even and odd parity moments. At this point, one can choose to follow several solution scheme to solve for the parity moments. Here, we choose to establish the even parity solution to include the unknown boundary conditions. By manipulation of the matrix equations, one finds the relationship between the incoming and the exiting fluxes through the response matrix. With the exiting fluxes known, the fluxes interior to the slab also become known. We form the solution in terms of the diagonalization of the matrix associated with the second order ODE, or, in other words, expressible as eigenvectors and eigenvalues. Convergence of the infinite series is through sequential convergence of the partial sums and Wynn-epsilon convergence acceleration. It was shown, via two examples, that benchmarks of 6 and 7 digits can be constructed. It should be noted that the precision quoted depends on the spatial positions and directions sampled but is thought to be representative of a wide variety of samples.

References

G.I. Bell and S. Glasstone, Nuclear Reactor Theory, Van Norstrand Reinhold, NY, 1975.
B. Davison, Neutron Transport Theory, Oxford University Press, UK, 1957.
G.C. Pomraning, NSE 55, 328 1956.
F.D. Federighi, Nukleonik 6, 277, 1965.
B.D. Ganapol, The response matrix discrete ordinates solution to the 1D radiative transfer equation, JQSRT 154, 72–90, 2015.
http://icl.cs.utk.edu/lapack-for-windows/lapack.
Weisstein, Eric W. "Clenshaw Recurrence Formula." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/ClenshawRecurrenceFormula.html.

Figure 1. Ratio of unaccelerated to accelerated relative errors.

Table 1. Angular Flux for an entering isotropic source.

μ\x	0.0	0.05	0.1	0.2	0.5	0.75	1.0
-1.000E+00	3.41328760E-01	3.20920611E-01	3.01041128E-01	2.62366118E-01	1.53240509E-01	7.07430107E-02	0.00000000E+00
-8.000E-01	3.92084430E-01	3.69683964E-01	3.47820410E-01	3.05071361E-01	1.82180798E-01	8.60214436E-02	0.00000000E+00
-6.000E-01	4.58134371E-01	4.33685363E-01	4.09782681E-01	3.62759573E-01	2.24012630E-01	1.09614879E-01	0.00000000E+00
-4.000E-01	5.43854301E-01	5.17792855E-01	4.92356464E-01	4.42065792E-01	2.88493773E-01	1.50567243E-01	0.00000000E+00
-2.000E-01	6.45967494E-01	6.19276078E-01	5.93756446E-01	5.43978308E-01	3.90966211E-01	2.35919292E-01	0.00000000E+00
0.000E+00	7.58146459E-01	7.22978545E-01	6.94563136E-01	6.42872374E-01	5.00000000E-01	3.81715377E-01	2.41853541E-01
2.000E-01	1.00000000E+00	9.42160751E-01	8.90352375E-01	8.02157479E-01	6.09033789E-01	4.80750231E-01	3.54032506E-01
4.000E-01	1.00000000E+00	9.69316928E-01	9.38639544E-01	8.78613253E-01	7.11506227E-01	5.83062169E-01	4.56145699E-01
6.000E-01	1.00000000E+00	9.79131021E-01	9.57479617E-01	9.12982318E-01	7.75987370E-01	6.60525697E-01	5.41865629E-01
8.000E-01	1.00000000E+00	9.84189969E-01	9.67479934E-01	9.32267600E-01	8.17819202E-01	7.15927250E-01	6.07915570E-01
1.000E+00	1.00000000E+00	9.87275159E-01	9.73675082E-01	9.44578242E-01	8.46759491E-01	7.56515359E-01	6.58671240E-01

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Double PN Benchmark Solution for the 1D Monoenergetic Neutron Transport Equation in Plane Geometry

Abstract

1. Introduction

2. Proposed DPN Algorithm

2.1. DPN Moment and Parity Equations

2.2. The Even/Odd DPN Parity Equations

3. Solution to the Parity Equations

4. Response Matrix

5. Final Moments Solution

6. The DPN Approximation

7. Numerical Implementation and Demonstration

7.a. Numerical Implementation for an Isotropic Source

7.b. Numerical Implementation for a Beam Source

Conclusion

References

MDPI Initiatives

Important Links

Subscribe