In this paper, we have used Case's analysis of the neutron transport equation to obtain a new set of quadrature coefficients for the discreteordinates method. We perform the transport calculations by this set of quadratures, dependent on the medium. We also use the orthogonality relations in the discrete case to derive the full-range formulation of the half-range problem. This solution can, indeed, be profitably used in the
1. INTRODUCTION
In this paper, we consider the regular and singular solutions of the monoenergetic neutron transport equation in plane geometry. As it is well known [1], the singular solutions are in the form of Dirac-delta distributions and the Cauchy principal value. Because of the singular nature of cigenfunctions, the exact method of neutron transport theory may be considered to yield only a formal solution and is not effective numerically. The present work will address an important natural question of how one can systematically use information from the exact solution in an approximate method. We have chosen to work with the discrete-ordinates (\left(S_{N}\right)) method [2] because of its suitability in practical transport computations.
In the context of spherical harmonics, various attempts [3-7] have bcen made to improve the tranditional (P_{N}) solutions of the transport equation. The basic improvement was achieved by the use of an exact asymptotic cigenvalue of the transport equation. The replacement of other transport parameters by their exact values (such as the extrapolated end-point and linear extrapolation length) was derived in a systematic manner through the change of boundary conditions. Surprisingly, the corresponding improvement of the discrete-ordinates method does not seem to have received equal attention. It is natural to expect that the regular solutions of the (S_{N}) and (P_{N}) methods will approach the distributional solutions of the transport equation in a suitable sense as the order of approximation (N) tends to infinity. In fact,
Davison [8] showed the convergence of the eigenvalues of the (P_{N}) equations to those of the exact theory. Apart from this work of Davison, the relationship of the regular solutions of the (P_{N}) and (S_{N}) methods with the exact solution of the transport boundary value problem has recently been established by Larsen [9].
We observe that the angular nodes in the (S_{N})-approximation are the zeros of Legendre polynomials (P_{N}(\mu)). However, the use of zeros of (P_{N}(\mu)), which are solutions of a SturmLiouville problem, may not be always optimal as the Legendre polynomials do not form a natural basis in transport theory calculations. Thus, we have used orthogonal polynomial with respect to a weight function other than unity [10] and succeeded in improving the discreteordinates calculations in general. These orthogonal polynomials were generated by using an algorithm of Golub and Welsch [11]. The weight function was derived by demanding certain transport parameters of the approximate scheme be exact. We have, however, experienced some difficuities in obtaining converged solution of the non-linear system of equations involving the parameters of the weight function for the case of a multiplying medium. We had to also deviate slightly from requiring that one of the eigenvalues of the discrete-ordinates equations be the exact asymptotic eigenvalue for highly absorbing cases to ensure positivity of the weight function. In this paper, we employ a different approach which will be described next. In most of our numerical calculations, we observed that the present approach performs better. However, it is to be noted that the former approach yields better results in some cases and we believe that the use of the weight function in another suitable form may improve the calculational results.
In the usual formulation of the Gauss-Legendre quadrature set, the angular nodes are the zeros of the Legendre polynomials (P_{N}(\mu)) and the weights are obtained by demanding the exact evaluation of the integrals for polynomials up to order (2 N-1). In our formulation, the weights and nodes for the quadrature approximation are obtained by applying two constraint equations to be discussed later,
together with the usual equations for the Gauss-Legendre set. Specifically, to obtain the (N)-point Gauss Legendre set, we use the following (2 N) equations:
[
\begin{aligned}
\sum_{j=1}^{N} w_{j} \mu_{j}^{n} & =\int_{-1}^{1} \mu^{n} d \mu, & n=0,2,4, \ldots, 2 N-2 \
w_{j} & =w_{N-j+1}, & j=1,2,3, \ldots, N / 2 \
\mu_{j} & =-\mu_{N-j+1}, & j=1,2,3, \ldots, N / 2
\end{aligned}
]
In our case, the last two equations for obtaining the Gauss-Legendre set (first equation with (n=2 N-4) and (n=2 N-2) ) are discarded and replaced by the two constraint equations. We have also considered the replacement of one constraint equation only. We observe that the Gauss-Legendre set integrates polynomials exactly up to order ((2 N-1)). With the inclusion of the constraint equations, we expect that the transport-theoretic properties will more than compensate for the loss of integrating properties by discarding equations from the usual set to obtain Gauss-Legendre quadrature coefficients. At the same time, we have the same computational advantage enjoyed by the traditional discrete-ordinates within the framework of coupled ordinary differential equations.
It now remains to derive the required constraint equations in a meaningful way. These can be derived from (i) purely physical considerations, (ii) purely mathematical relations derived from the half-range and full-range orthogonality relations, and (iii) a suitable combination of (i) and (ii). In the present work, we derive the constraint equations in Section 2 from the physical considerations. We also use the discrete-ordinates orthogonality relations to derive constraint equations. This has been achieved by a full-range type formulation for the discrete-ordinates scheme. We obtain one of the constraints by the replacement of the discrete-ordinates asymptotic eigenvalue by its exact counterpart. By demanding this, we ensure that the largest eigenvalue corresponding to the asymptotic flux does not lie in the range of ((-1,1)) in agreement with the exact analysis. The second constraint will be derived by demanding that either (a) exact leakage from a half-space with constant source or (b) exact extrapolated end-point of the standard source free Milne problem. For a nonmultiplying medium, either (a) or (b) can be used. However, (b) is the only choice for a multiplying medium because of the fact that the source-free Milne problem can be analytically extended for this case. It is worthwhile to investigate if the physical constraints can be interpreted from the analysis of the transport equation. For this purpose, we also derive in this section suitable orthogonality relations in the framework of the (S_{N}) method. The identities satisfied by the discrete case have exact counterparts from Case's solution. By the use of these relations, we will justify that one of the constraints, namely, the requirement of exact asymptotic eigenvalue in the discrete-ordinates scheme is necessiated by rigorous analysis.
Section 3 of the paper provides the new quadrature nodes and weights based on the present formulation. We then compare the performance of the new weights and nodes with the standard set for various problems involving isotropic and anisotropic scattering. We note that improvement of the (S_{N}) method may also be attempted by imposing only one constraint equation. We have done that with partial success, as the numerical results in Section 3 will indicate. Finally, we hope that the present approach will find applications in neutron transport calculations in geometries other than the slab geometry and in other areas of mathematical physics. In Appendix A, we derive the appropriate interface conditions for more realistic problems involving interfaces. We also present numerical results on a four region interface problem and a fuel moderator cell problem. In Appendix B, the quadrature coefficients are derived with constraints obtained from a full-range type formulation. The numerical results on several benchmark problems are also provided.
2. ANALYSIS
We consider the linear transport equation in plane geometry:
(\mu \frac{\partial \psi(x, \mu)}{\partial x}+\psi(x, \mu)=c \int_{-1}^{1} K\left(x, \mu, \mu^{\prime}\right) \psi\left(x, \mu^{\prime}\right) d \mu^{\prime}),
where (\psi(x, \mu)) is the angular flux in direction (\mu) at position (x) and (c) is the number of secondaries per primary. Here, we consider isotropic scattering such that (K\left(x, \mu, \mu^{\prime}\right)=\frac{1}{2}). We first summarize the exact analysis of the transport equation as our analysis of the discrete-ordinates equation is similar in many respects.
By application of the ansatz [1]
[
\psi(x, \mu)=e^{-x / v} \phi(v, \mu)
]
to the homogeneous equation (4) with the normalization (\int_{-1}^{1} \phi\left(v, \mu^{\prime}\right) d \mu^{\prime}=1), the following eigenvalue equation for (\phi(v, \mu)) is obtained:
[
(v-\mu) \phi(v, \mu)=c v / 2
]
The solution is
[
\phi(v, \mu)=\frac{c v}{2} P \frac{1}{v-\mu}+\lambda(v) \delta(v-\mu)
]
with the eigenvalue spectrum consisting of a continuum
(-1 \leqslant v \leqslant 1) and a discrete spectrum (\pm v_{0}) with the eigenfunctions
[
\phi\left( \pm v_{0}, \mu\right)= \pm \frac{c v_{0}}{2} \frac{1}{ \pm v_{0}-\mu}, \quad v_{0} \notin(-1,1)
]
where
[
\frac{c v_{0}}{2} \ln \frac{v_{0}+1}{v_{0}-1}=1
]
Here, (P) denotes the principal value and (\lambda(v)=) (1-c v \tanh ^{-1} v). The Cauchy principle value prescription is chosen to evaluate the integral in the following form:
[
P \int_{-1}^{1} \frac{v}{v-\mu^{\prime}} d \mu^{\prime}=\lim _{\delta \rightarrow 0}\left[\int_{-1}^{v-\delta} \frac{v}{v-\mu^{\prime}} d \mu^{\prime}+\int_{v+\delta}^{1} \frac{v}{v-\mu^{\prime}} d \mu^{\prime}\right]
]
The solution of (4) is given by
[
\begin{aligned}
\psi(x, \mu)= & a_{0+} e^{-x / v_{0}} \phi\left(v_{0}, \mu\right)+a_{0-} e^{x / v_{0}} \phi\left(-v_{0}, \mu\right) \
& +\int_{-1}^{1} A(v) e^{-x / v} \phi(v, \mu) d v
\end{aligned}
]
and
[
\begin{aligned}
\psi(x, \mu)= & a_{0+} e^{-x / v_{0}} \phi\left(v_{0}, \mu\right) \
& +\int_{0}^{1} A(v) e^{-x / v} \phi(v, \mu) d v
\end{aligned}
]
for the full-range and half-range cases, respectively. The coefficients (a_{0 \pm}) and (A(v)) are obtained from boundary conditions. The corresponding orthogonality relations are
[
\int_{-1}^{1} \mu \phi(v, \mu) \phi\left(v^{\prime}, \mu\right) d \mu=0, \quad v \neq v^{\prime}
]
and,
[
\int_{0}^{1} W(\mu) \phi(v, \mu) \phi\left(v^{\prime}, \mu\right) d \mu=0, \quad v \neq v^{\prime}
]
In Eq. (8) (v, v^{\prime}) can be either (\pm v_{0}) or any value in ((-1,1)), whereas (v, v^{\prime}), corresponding to (9), will take either (v_{0}) or any value in ((0,1)). The weight function (W(\mu)) in the half-range has been found for the neutron transport theory by the method of Case [1]. Analogous orthogonality relations are derived in the present investigation for the discreteordinates method. As in the case of distributional solutions [1], these orthogonality relations are expected to be useful to construct suitable approximate methods.
2.1. Discrete-Ordinates Equations
We consider the discrete-ordinates equations derived directly from (4) by using a quadrature formula to approximate the integral. Thus, we obtain
(\mu_{i} \frac{d \psi_{i}(x)}{d x}+\psi_{i}(x)=\frac{c}{2} \sum_{j=1}^{N} w_{j} \psi_{j}(x), \quad i=1,2, \ldots, N),
where, (\psi_{i}(x) \approx \psi\left(x, \mu_{i}\right)) and (\mu_{i}, w_{i}) are the quadrature weights and points. As usual, we impose the following constraints on the quadratures:
(i) the weights are positive and the quadrature sets integrate constant functions (neutron conservation condition) exactly. That is,
[
\begin{gathered}
w_{j}>0, \quad j=1,2, \ldots, N \
\sum_{j=1}^{N} w_{j}=2
\end{gathered}
]
(ii) the discrete-ordinates are distinct
[
\mu_{i} \neq \mu_{j} ; \quad i \neq j
]
(iii) the quadrature sets are symmetric about (\mu=0) :
[
\begin{array}{ll}
\mu_{i}=-\mu_{N+1-i}, & i=1,2, \ldots, N / 2 \
w_{i}=w_{N+1-i}, & i=1,2, \ldots, N / 2
\end{array}
]
Analogus to the procedure used in the exact method of solving (4), we can solve (10) by expanding the solution in terms of the eigenfunctions of the homogeneous solution and obtaining expansion coefficients from the boundary conditions. This procedure, originally considered by Chandrasekhar [12], was developed with the Gauss-Legendre quadrature set. Sykes [13] proposed a double-Gauss formula to determine the emergent angular distribution more reliably. In this paper, we will, however, consider a different quadrature set dependent on the medium necessitated by the analysis and some physical constraints from transport theory. Our quadrature set will satisfy the constraints mentioned above. We now investigate the analysis more thoroughly to establish relations satisfied by the "discrete" eigenfunctions and eigenvalues of the regular (S_{N}) method similar to the corresponding relations for exact singular solutions.
Equation (10) can be written in the form
[
\frac{d \psi}{d x}=A \psi
]
where (\psi=\left(\psi_{1}, \psi_{2}, \ldots, \psi_{N}\right)^{\mathrm{T}}) is a vector of length (N) and (A) is
an (N \times N) matrix. Through tedious calculation it is possible to show that the eigenvalues of (A) satisfy
[
\prod_{j=1}^{N / 2}\left(\lambda^{2}-\frac{1}{\mu_{j}^{2}}\right)+\sum_{j=1}^{N / 2}\left(\frac{c w_{j}}{\mu_{j}^{2}} \prod_{\substack{i=1 \ i \neq j}}^{N / 2}\left(\lambda^{2}-\frac{1}{\mu_{i}^{2}}\right)\right)=0
]
It is clear, therefore, that the eigenvalues of (A) are (\pm \lambda_{1}) for (l=1,2, \ldots, N / 2) where, for convenience, (\operatorname{Re}\left(\lambda_{l}\right)>0), for (l=1,2, \ldots, N / 2). It is also possible to show that if (\mu_{i}), (1 \leqslant i \leqslant N), are distinct, then the eigenvalues are distinct. It is assumed that (\mu_{i} \neq \mu_{j}) if (i \neq j) and therefore (A) has (N) distinct eigenvalues. Thus, the solution of (11) has the form
[
\psi_{j}(x)=\sum_{l=1}^{N} \phi_{l_{j}} e^{-j_{i} x} \quad \text { for } \quad j=1,2, \ldots, N
]
where (\lambda_{l}, 1 \leqslant l \leqslant N), are the (N) distinct eigenvalues of (A).
It will now be shown that (13) is similar in many respects to (6), the solution of the continuous problem. By dividing (12) by (\prod_{j=1}^{N / 2}\left(\hat{\lambda}^{2}-1 / \mu_{j}^{2}\right)) (it is straightforward to show that (\lambda_{l} \neq 1 / \mu_{j}) for every (l) and (j) and hence this expression is nonzero), one obtains the dispersion relation
[
\sum_{j=1}^{N / 2} \frac{w_{j}}{1-\mu_{j}^{2} \lambda^{2}}=\frac{1}{c}
]
We use the same form of eigenfunctions both in the asymptotic and transient ranges in (14). We observe that mathematically more accurate rational function approximation of the singular eigenfunction may be used in the transient range as in the seminal work of discretized spectral approximation [14]. However, in this work, we prefer to use the same form both in the asymptotic and transient ranges for the computational advantage. Now, substituting (13) into (10) results in
[
\begin{aligned}
& \sum_{l=1}^{N}\left(-\phi_{l_{j}} \lambda_{l} \mu_{j}+\phi_{l_{j}}\right) e^{-\lambda_{i} x} \
& =\frac{c}{2} \sum_{k=1}^{N} w_{k} \sum_{l=1}^{N} \phi_{l_{k}} e^{-\lambda_{i} x} \
& \quad \text { for } j=1,2, \ldots, N
\end{aligned}
]
Since the right-hand side of (15) is the same for any (j),
[
\begin{aligned}
& \sum_{l=1}^{N}\left(-\phi_{l_{j}} \lambda_{l} \mu_{j}+\phi_{l_{j}}\right) e^{-\lambda_{l} x} \
& =\sum_{l=1}^{N}\left(-\phi_{l_{m}} \lambda_{l} \mu_{m}+\phi_{l_{m}}\right) e^{-\lambda_{l} x} \
& \quad \text { for } \quad 1 \leqslant j, \quad m \leqslant N
\end{aligned}
]
Therefore, equating coefficients of (e^{-\lambda_{l} x}) in (16), the following result is obtained:
[
\phi_{l_{j}}\left(1-\lambda_{l} \mu_{j}\right)=\phi_{l_{m}}\left(1-\lambda_{l} \mu_{m}\right) \quad \text { for all } j, m, l
]
Therefore, since (\phi_{l_{j}}\left(1-\lambda_{l} \mu_{j}\right)) does not depend on (j), let (b_{l}=\phi_{l_{j}}\left(1-\lambda_{l} \mu_{j}\right)). Then (\phi_{l_{j}}=b_{l} /\left(1-\lambda_{l} \mu_{j}\right)) and by (13),
[
\psi_{j}(x)=\sum_{l=1}^{N} \frac{b_{l}}{1-\lambda_{l} \mu_{j}} e^{-\lambda_{i} x} \quad \text { for } \quad j=1,2, \ldots, N
]
The similarity between (18) and (7) is clear. We refer the readers to Ref. [9] which provides the relationship of the exact distributional solutions with the regular solutions of the discrete-ordinates equations.
Before continuing, some additional interesting and useful relations are derived which are satisfied by the points, weights, and eigenvalues. Equations (15) and (17) imply that
[
\phi_{l_{j}}\left(1-\lambda_{l} \mu_{j}\right)=\frac{c}{2} \sum_{k=1}^{N} w_{k} \phi_{l_{k}} \quad \text { for } \quad 1 \leqslant l \leqslant N
]
But from (17), (\phi_{l_{k}}=\phi_{l_{j}}\left(1-\lambda_{l} \mu_{j}\right) /\left(1-\lambda_{l} \mu_{k}\right)). Substituting this into (19) results in
[
\left(1-\lambda_{l} \mu_{j}\right) \phi_{l_{j}}=\frac{c}{2} \sum_{k=1}^{N} w_{k} \phi_{l_{j}}\left(1-\lambda_{l} \mu_{j}\right) /\left(1-\lambda_{l} \mu_{k}\right)
]
and, hence,
[
\sum_{j=1}^{N} \frac{w_{j}}{1-\lambda_{i} \mu_{j}}=\frac{2}{c}
]
Note, that if the points and weights are symmetric then (21) reduces to (14). Now consider
[
\begin{aligned}
\sum_{j=1}^{N} & \frac{w_{j} \mu_{j}}{\left(1-\lambda_{l} \mu_{j}\right)\left(1+\lambda_{m} \mu_{j}\right)} \
& =\sum_{j=1}^{N} \frac{w_{j} / \lambda_{m}}{1-\lambda_{l} \mu_{j}}-\sum_{j=1}^{N} \frac{w_{j} / \lambda_{m}}{\left(1-\lambda_{i} \mu_{j}\right)\left(1+\lambda_{m} \mu_{j}\right)} \
& =\sum_{j=1}^{N} \frac{-w_{j} / \lambda_{l}}{1+\lambda_{m} \mu_{j}}+\sum_{j=1}^{N} \frac{w_{j} / \lambda_{l}}{\left(1-\lambda_{l} \mu_{j}\right)\left(1+\hat{\lambda}_{m} \mu_{j}\right)}
\end{aligned}
]
Thus, by equation the right-hand sides of (22), one obtains that
[
\begin{aligned}
& \sum_{j=1}^{N} \frac{\left(1 / \lambda_{l}+1 / \lambda_{m}\right) w_{j}}{\left(1-\lambda_{l} \mu_{j}\right)\left(1+\lambda_{m} \mu_{j}\right)} \
& \quad=\frac{1}{\lambda_{l}} \sum_{j=1}^{N} \frac{w_{j}}{1+\lambda_{m} \mu_{j}}+\frac{1}{\lambda_{m}} \sum_{j=1}^{N} \frac{w_{j}}{1-\lambda_{l} \mu_{j}}=\left(\frac{1}{\lambda_{l}}+\frac{1}{\lambda_{m}}\right) \frac{2}{c}
\end{aligned}
]
assuming that the points and weights are symmetric. Hence,
[
\begin{aligned}
& \sum_{j=1}^{N} \frac{w_{j}}{\left(1-\hat{\lambda}_{l} \mu_{j}\right)\left(1+\lambda_{m} \mu_{j}\right)} \
& \quad=\frac{2}{c} \quad \text { for } \quad l, m \text { such that } \lambda_{l} \neq-\lambda_{m}
\end{aligned}
]
and substituting (23) into (22),
[
\begin{aligned}
& \sum_{j=1}^{N} \frac{w_{j} \mu_{j}}{\left(1-\lambda_{l} \mu_{j}\right)\left(1+\lambda_{m} \mu_{j}\right)} \
& \quad=0 \quad \text { for } \quad l, m \text { such that } \lambda_{l} \neq-\lambda_{m}
\end{aligned}
]
In a similar manner one can obtain
[
\begin{array}{ll}
\sum_{j=1}^{N} \frac{w_{j}}{\left(1-\lambda_{l} \mu_{j}\right)\left(1-\lambda_{m} \mu_{j}\right)}=\frac{2}{c} & \text { for } \quad l \neq m \
\sum_{j=1}^{N} \frac{\mu_{j} w_{j}}{\left(1-\hat{\lambda}_{l} \mu_{j}\right)\left(1-\lambda_{m} \mu_{j}\right)}=0 & \text { for } \quad l \neq m
\end{array}
]
The key relations satisfied by the points, weights, and eigenvalues are (12), (16), (21), (23)-(26). In particular, (26) is the discrete version of the orthogonality relation satisfied by the eigenfunctions (\phi(v, \mu)) :
[
\int_{-1}^{1} \mu \phi(v, \mu) \phi\left(v^{\prime}, \mu\right) d \mu=0
]
Indeed (26) can be written as
[
\sum_{j=1}^{N} w_{j} \mu_{j} \phi_{l_{j}} \phi_{m_{j}}=0 \quad \text { for } \quad l \neq m
]
We also note that
[
\int_{-1}^{1} \phi\left(v_{0}, \mu\right) \phi(v, \mu) d \mu=c / 2
]
where (\phi\left(v_{0}, \mu\right)) and (\phi(v, \mu)) are exact eigenfunctions. Equation (25) is the discrete version of (27). The normalization integral of the full-range problem in the discrete version is given by (\sum_{j=1}^{N}\left(w_{j} \mu_{j} /\left(1-\lambda_{i} \mu_{j}\right)^{2}\right)). By equating this to the exact normalization condition, we may obtain a constraint on quadrature weights and nodes.
The half-range orthogonality relations are relatively difficult to formulate because of their half-range characteristic. In the case of the exact neutron transport equation, the halfrange problem has been completely solved by the singular eigenfunction method. The difficult step in the analysis is determination of the half-range weight function, with respect to which the eigenfunctions are orthogonal. In the discrete case, we have imposed a suitable normalization condition and used the half-range orthogonality relations to derive conditions on the weights and nodes of the quadrature set. Because of the presence of complicated functions, we have not been able to use these relations profitably. Fortunately, Siewert's work [15] on the (F_{N}) method and, then, its generalization by Sengupta [16] tell us that the full-range weight function " (\mu) " will be sufficient to solve the half-range problem. Actually, we do not need to use the complicated half-range weight function " (W(\mu)) " and the corresponding orthogonality relations.
2.2. Quadrature Set
(a) Based on physical considerations. We proceed to derive the two constraints which will be used to obtain the new quadrature coefficients. A detailed discussion on this aspect may be found in Refs. ([6,7]). We require that one of the eigenvalues of the matrix (A) be (1 / v_{0}). By demanding this, we ensure that the largest (v_{j}) corresponding to the asymptotic flux does not lie in the range ((-1,1)) for any order of approximation. Thus the first constraint follows from (12) as
[
\sum_{j=1}^{N / 2} \frac{w_{j} v_{0}^{2}}{\left(v_{0}^{2}-\mu_{j}^{2}\right)}=\frac{1}{c}
]
The second constraint is derived by demanding either (a) the exact leakage for the half-space constant source problem or (b) the exact extrapolated end-point of the standard source-free Milne problem. The exact leakage for the half-space constant source problem is [8]
[
L_{\text {exact }}=\frac{2 a}{c}\left[v_{0}-I(c)\right]
]
where
[
I(c)=\frac{c}{2} \int_{0}^{1}\left(1+\frac{c \mu^{2}}{1-\mu^{2}}\right) \mu g(c, \mu) d \mu
]
and
[
g(c, \mu)=\left[\left(1-c \mu \tanh ^{-1} \mu\right)^{2}+\frac{\pi^{2} c^{2} \mu^{2}}{4}\right]^{-1}
]
The functions (I(c)) and (g(c, \mu)) are tabulated in Ref. [17]. The approximate angular flux is
[
\psi_{i}(x)=\sum_{j=1}^{N / 2} \frac{b_{j}}{1-\lambda_{j} \mu_{i}} e^{-\lambda_{j} x}+\frac{a}{1-c}
]
The boundary conditions are
[
\psi_{i}(0)=0, \quad i=1,2, \ldots, N / 2
]
The coefficients (b_{i}) are computed by use of (32) and (33). The leakage in the discrete-ordinates formulation is computed from
[
L_{S_{N}}=-\sum_{j=N / 2+1}^{N} \mu_{j} w_{j} \psi_{j}(0)
]
The second constraint equation is now obtained by demanding that the leakage from (34) be equal to that given by (29). Therefore,
[
-\sum_{j=N / 2+1}^{N} \mu_{j} w_{j} \psi_{j}(0)=\frac{2 a}{c}\left[v_{0}-I(c)\right]
]
We now derive the constraint equation from the sourcefree Milne problem. The exact solution of this problem is
[
\begin{aligned}
\psi(x, \mu)= & a_{0+} e^{-x / v_{0}} \phi_{0+}(\mu)+e^{x / v_{0}} \phi_{0-}(\mu) \
& +\int_{0}^{1} A(v) e^{-x / v} \phi_{v}(\mu) d v
\end{aligned}
]
where
[
a_{0+}=-\exp \left(-\frac{2 z_{0}}{v_{0}}\right)
]
with (z_{0}) defining the extrapolated endpoint. The discrete solution to this problem is
[
\psi_{i}(x)=\sum_{j=1}^{N / 2} \frac{d_{j}}{1-\lambda_{j} \mu_{i}} e^{-i_{j} x}+\frac{c}{2\left(1+\lambda_{N / 2} \mu_{i}\right)} e^{i_{N / 2} x}
]
The coefficients (d_{i}) in (38) are obtained by use of
[
\psi_{i}(0)=0, \quad i=1,2, \ldots, N / 2
]
It is assumed that (\lambda_{N / 2}) is the smallest eigenvalue in (38). By comparison of (38) with (36), it is readily observed that the discrete formulation yields the exact extrapolated end-point (z_{0}) if
[
\frac{2}{c} d_{N / 2}=-\exp \left(-\frac{2 z_{0}}{v_{0}}\right)
]
For nonmultiplying medium ( (c<1) ), either (35) or (40) may be used as the second constraint. We have chosen (35) for (c<1) which produced better numerical results. Since the source-free Milne problem can be analytically extended for (c>1,(40)) is, however, the only choice for the multiplying medium.
We have considered the following non-linear system of equations to obtain a new quadrature set dependent on the medium: (a) Eqs. (1), with (n=2 N-2) case discarded, (2), (3), and (28)
(b) Eqs. (1), with (n=2 N-2) case discarded, (2), (3), and (35)
(c) Eqs. (1), with (n=2 N-2) and (n=2 N-4) cases discarded, (2), (3), (28), and (35) ((c<1)) or (40) ((c>1)).
The detailed numerical results will be presented in the next section for the case (c), which performs remarkably well for a wide class of transport problems with both isotropic and anisotropic scattering. Prior to this, we will briefly discuss the cases (a) and (b) for (N=4) only, in that section.
(b) Based on orthogonality relations. We believe that the orthogonality relations in the discrete case can be suitably used to obtain new quadrature sets. We already noted the difficulty of using half-range weight function (W(\mu)) and the corresponding orthogonality relations. In the next section, we will show the basic steps of full-range formulation of the half-range problem in the discrete-ordinates case.
The full-range solution of half-range problems is based on a reduction to an equivalent full-range problem by the superposition of two half-range problems. Since our analysis is in the framework of discrete-ordinates, it seems to be suited for numerical computations. Instead of deriving the formulation from Plackzek lemma or the fundamental orthogonality relations, we take advantage of an established full-range formulation.
Consider the half-space problem
[
\begin{aligned}
& \mu \frac{\partial \psi(x, \mu)}{\partial x}+\psi(x, \mu) \
& \quad=\frac{c}{2} \int_{-1}^{1} \psi\left(x, \mu^{\prime}\right) d \mu^{\prime}+a, \quad x>0 \
& \psi(0, \mu)=1-a, \quad \mu>0
\end{aligned}
]
When (a=0), this corresponds to the albedo problem and (41) with (a=1) represents the constant source problem. By use of the relation
[
\begin{gathered}
\int_{-1}^{1}\left[\psi(x, \mu)-\frac{a}{1-c}\right] \mu \phi(-\xi, \mu) d \mu \
=0 \quad \text { for } \quad \xi \in(0,1) \cup\left{v_{0}\right}
\end{gathered}
]
Siewert and Benoist [15] obtained the singular integral equation
[
\begin{aligned}
& \frac{2}{c \xi} \int_{0}^{1} \phi(\xi, \mu) \psi(0,-\mu) \mu d \mu \
& \quad=\frac{2 a}{c}+(1-a)\left[1-\xi \log \left(1+\frac{1}{\xi}\right)\right]
\end{aligned}
]
Here, (\phi(\xi, \mu)) are the regular or singular eigenfunctions depending on the value of (\xi). Thus by projecting Case's halfrange solution onto a subspace of the full-range problem, Siewert indirectly solved the transient integral problem [14] involving singular eigenfunctions. By letting (\psi(0,-\mu)=\sum_{\alpha=0}^{N} a_{\alpha} \mu^{\alpha}) and evaluating (43) at selected values of (\xi), the (F_{N}) equations are obtained,
[
\begin{gathered}
\sum_{\alpha=0}^{N} a_{\alpha} B_{\alpha}\left(\xi_{\beta}\right)=\frac{2 a}{c}+(1-a)\left[1-\xi_{\beta} \log \left(1+\frac{1}{\xi_{\beta}}\right)\right] \
\text { for } \beta=0,1,2, \ldots, N
\end{gathered}
]
where
[
B_{\alpha}(\xi)=\xi B_{\alpha-1}(\xi)-\frac{1}{\alpha-1}, \quad \alpha \geqslant 1
]
with
[
B_{0}(\xi)=\frac{2}{c}-1-\xi \log \left(1+\frac{1}{\xi}\right)
]
In this manner, (\psi(0,-\mu)) can be approximated. The (F_{N}) method has been successfully applied to many neutron transport and radiative transfer problems.
The discrete-ordinates formulation of (41) is
[
\begin{aligned}
\mu_{j} \frac{d \psi_{j}(x)}{d x}+\psi_{j}(x) & =\frac{c}{2} \sum_{k=1}^{N} w_{k} \psi_{k}(x)+a \
\psi_{j}(0) & =1-a, \quad 1 \leqslant j \leqslant N / 2
\end{aligned}
]
where we as usual assume that (\mu_{j} \geqslant 0) for (1 \leqslant j \leqslant N / 2) and (\sum_{k=1}^{N} w_{k}=2). An analogous relation to (43) is now derived that is satisfied by solutions of the discrete problem.
The solution of (46) tends to (a /(1-c)) as (x \rightarrow \infty); i.e., (\psi_{j}(x) \rightarrow a /(1-c)) as (x \rightarrow \infty). Additionally, (46) can be written as
[
\frac{d \vec{\psi}}{d x}=A \vec{\psi}+\vec{s}
]
where (A) is the same matrix that was considered in (11). Therefore, the solution of (47) has the form
(\psi_{j}(x)=\sum_{t=1}^{N / 2} \phi_{l_{j}} e^{-\lambda_{i} x}+\frac{a}{1-c} \quad) for (\quad j=1,2, \ldots, N),
assuming that (\operatorname{Re}\left(\lambda_{l}\right)>0) for (l=1,2, \ldots, N / 2). Of course, all the previous relations derived for the (\lambda_{l}, \mu_{j}), and (w_{j}) still hold as well for the (\phi_{l_{j}}). In particular, (\phi_{l_{j}}=b_{l} /\left(1-\lambda_{l} \mu_{j}\right)) for some constant (b_{l}) for each (l) and (j). Hence,
[
\begin{gathered}
\psi_{j}(x)=\sum_{l=1}^{N / 2} \frac{b_{l}}{1-\lambda_{l} \mu_{j}} e^{-\lambda_{l} x}+\frac{a}{1-c} \
\text { for } j=1,2, \ldots, N
\end{gathered}
]
Note that the constants (b_{l}) can be obtained from the boundary conditions (\psi_{j}(0)=1-a) for (j=1,2, \ldots, N / 2).
Now a discrete form of (43) is derived. Multiplying both sides of (49) by (\mu_{j} w_{j} /\left(1+\hat{\lambda}_{m} \mu_{j}\right)) and summing over (j), one obtains that
[
\begin{aligned}
& \sum_{j=1}^{N}\left(\psi_{j}(x)-\frac{a}{1-c}\right)\left(\frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}}\right) \
& \quad=\sum_{l=1}^{N / 2} b_{l} \sum_{j=1}^{N} \frac{w_{j} \mu_{j}}{\left(1-\lambda_{l} \mu_{j}\right)\left(1+\hat{\lambda}_{m} \mu_{j}\right)} e^{-\lambda_{l} x}
\end{aligned}
]
However, assuming that (\lambda_{m}>0), the right-hand side of (50) is zero by (24). Hence,
[
\sum_{j=1}^{N}\left(\psi_{j}(x)-\frac{a}{1-c}\right) \phi_{l_{j}}^{*} \mu_{j} w_{j}=0
]
is the discrete form of (42), where (\phi_{l_{j}}^{*}=b_{l} /\left(1+\lambda_{l} \mu_{j}\right)).
Now let (x=0) in (51). Then,
[
\sum_{j=1}^{N}\left(\psi_{j}(0)-\frac{a}{1-c}\right)\left(\frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}}\right)=0
]
Since (\psi_{j}(0)=1-a) for (1 \leqslant j \leqslant N / 2), one obtains
(\sum_{j=N / 2+1}^{N} \psi_{j}(0) \frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}})
[
=\sum_{j=1}^{N}\left(\frac{a}{1-c}\right)\left(\frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}}\right)-\sum_{j=1}^{N / 2}(1-a) \frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}}
]
However,
[
\sum_{j=1}^{N} \frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}}=\sum_{j=1}^{N} w_{j} \frac{1}{\lambda_{m}}-\sum_{j=1}^{N} \frac{w_{j} / \lambda_{m}}{1+\lambda_{m} \mu_{j}}=\frac{2}{\lambda_{m}}\left(\frac{c-1}{c}\right)
]
Hence, (53) becomes
[
\begin{aligned}
& \sum_{j=N / 2+1}^{N} \psi_{j}(0) \frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}} \
& \quad=\frac{-2 a}{\lambda_{m} c}-\sum_{j=1}^{N / 2}(1-a) \frac{w_{j} \mu_{j}}{1+\lambda_{m} \mu_{j}}
\end{aligned}
]
which is a discrete form of the singular integral equation
(43). In particular, suppose that (a=1) and (\lambda_{m}=1 / v_{0}) for some (m). Then,
[
-\sum_{j=N / 2+1}^{N} \psi_{j}(0) \frac{w_{j} \mu_{j}}{v_{0}+\mu_{j}}=\frac{2}{c}
]
which is very similar to (43) with (a=1) and (\xi=v_{0}). So, from the constant source Milne problem ((a=1)), the present analysis imposes the requirement of exact asymptotic eigenvalue, namely, the constraint (28). We also observe that (54), for (a=0), reduces to (with (\lambda_{m}=1 / v_{0}) )
[
\sum_{j=N / 2+1}^{N} \psi_{j}(0) \frac{w_{j} \mu_{j}}{1+\left(1 / v_{0}\right) \mu_{j}}=-\sum_{j=1}^{N / 2} \frac{w_{j} \mu_{j}}{1+\left(1 / v_{0}\right) \mu_{j}}
]
By comparison of (56) and its continuous analog (43) with (a=0), one may impose another constraint in the form
[
\sum_{j=1}^{N / 2} \frac{w_{j} \mu_{j}}{v_{0}+\mu_{j}}=1-v_{0} \ln \left(1+\frac{1}{v_{0}}\right)
]
We have suitably used the constraints (28) and (57), together with the usual equations for the Gauss-Legendre set, to obtain new quadrature sets. Excellent numerical results are obtained as expected, being a full-range type formulation. In this paper, we have, however, provided results for (N=4) only in Appendix B as we intend to improve the method further in a more general way. To this end, different schemes [16] other than the (F_{N}) equations may be considered. In particular, we [18] have used the set of orthogonal polynomials with respect to (\mu) in ((0,1)) followed by a reflection of the zeroes on the other half ((-1,0)). This makes the direction cosines independent of the medium. The weights may then be obtained from the usual equations, together with constraints.
3. NUMERICAL RESULTS AND CONCLUSIONS
In Table I, the quadrature points and weights of case (a) (as described in the previous section) are given for (N=4). For (N=2), it is straightforward to show that (w_{1}=w_{2}=1) and (\mu=-\mu_{2}=v_{0} \sqrt{1-c}). It is interesting to note that as (c \rightarrow 0), Radau quadrature points and weights are approached and, as (c \rightarrow 1), Gaussian quadrature points and weights are approached. Since (28) has only (v_{0}^{2}) appearing in the equation, case (a) can also be applied for (c>1), in which case (v_{0}) is purely imaginary. In ail the calculations, discreteordinates calculations are performed for spatial mesh widths sufficiently small so that the error in the calculations is dominated by the error due to the angular discretization.
We first consider the criticality problem. Slabs of widths 2.5786, 1.4732 and 1.0240 mean-free paths are considered
TABLE I
Quadrature Weights and Nodes for Case (a)
\begin{tabular}{ccccc}
\hline(c) & (w_{1}) & (w_{2}) & (\mu_{1}) & (\mu_{2}) \
\hline 0.1 & 0.824918 & 0.175082 & 0.442697 & 0.990197 \
0.2 & 0.812110 & 0.187890 & 0.435806 & 0.976305 \
0.3 & 0.793264 & 0.206736 & 0.425595 & 0.957784 \
0.4 & 0.770410 & 0.229590 & 0.413009 & 0.937806 \
0.5 & 0.746664 & 0.253336 & 0.399586 & 0.919335 \
0.6 & 0.723914 & 0.276086 & 0.386280 & 0.903387 \
0.7 & 0.702978 & 0.297022 & 0.373546 & 0.889946 \
0.8 & 0.684088 & 0.315912 & 0.361564 & 0.878671 \
0.9 & 0.667198 & 0.332802 & 0.350381 & 0.869181 \
1.1 & 0.638730 & 0.361270 & 0.330318 & 0.854262 \
1.2 & 0.626752 & 0.373248 & 0.321335 & 0.848338 \
1.3 & 0.616030 & 0.383970 & 0.312973 & 0.843191 \
1.4 & 0.606398 & 0.393602 & 0.305177 & 0.838686 \
\hline
\end{tabular}
for (c=1.2,1.4), and 1.6 , respectively. The exact multiplication factor for each slab is 1.000 . Multiplication factors for the new quadrature set are (0.996,0.988), and 0.979 for (c=1.2,1.4), and 1.6 , respectively. The corresponding multiplication factors for the Gauss quadrature set are 0.994 , (0.980,0.960). Thus, the new quadrature set perform better for the criticality problem.
In the next example, we consider a slab of five mean-free paths with an entering flux distribution of (\psi(0, \mu)=) (1 /\left(v_{0}-\mu\right)) on the left-hand side. This form of distribution is approached as the distance from the source through a homogeneous medium increases. Using discrete-ordinates calculations with (N=4), leakages from the right-hand side of the slab are calculated for (c=0.3). These values are 0.0103 and 0.0339 for the Gauss quadrature and the new quadrature sets. In comparing these leakages with the exact leakage of 0.0338 , it is clear that the new set performs much better. This is, perhaps, because of the unique form of the entering flux distribution.
Next, the leakages from the half-space of a constant course problem are compared. We observe that the new quadrature points and weights do not fare quite as well as using the Gaussian set. This demonstrates that the use of (1 / v_{0}) as one of the eigenvalues of the discrete-ordinates equations does not guarantee improvements in all cases and additional constraint equations need to be applied.
In Table II we provide the weights and nodes for (N=4) corresponding to case (b), where the only constraint equation is derived by demanding exact leakage from a halfspace with a constant source. In the Table III we compare the leakages for a half-space albedo problem ((a=0)) and find that the new set results (\left(S_{4}^{*}\right)) are very good. However, we have observed that the requirement of the exact asymptotic eigenvalue is an important one for many halfspace and criticality problems.
We now consider the important case (c). Here, two
TABLE II
Quadrature Weights and Nodes for Case (b)
\begin{tabular}{ccccc}
\hline(c) & (w_{1}^{\prime}) & (w_{2}) & (\mu_{1}) & (\mu_{2}) \
\hline 0.1 & 0.529358 & 0.470642 & 0.228498 & 0.805932 \
0.2 & 0.530166 & 0.469834 & 0.229494 & 0.806251 \
0.3 & 0.531326 & 0.468674 & 0.230911 & 0.806708 \
0.4 & 0.533434 & 0.466566 & 0.233460 & 0.807543 \
0.5 & 0.534758 & 0.465242 & 0.235041 & 0.808069 \
0.6 & 0.536102 & 0.463898 & 0.236630 & 0.808603 \
0.7 & 0.538653 & 0.461348 & 0.239610 & 0.809622 \
0.8 & 0.541710 & 0.458290 & 0.243117 & 0.810849 \
0.9 & 0.542602 & 0.457398 & 0.244129 & 0.811209
\end{tabular}
TABLE III
Leakage for the Half-Space Albedo Problem
\begin{tabular}{cccc}
\hline(c) & (S_{4}^{*}) & (S_{4}) & Exact leakage \
\hline 0.1 & 0.011113 & 0.012407 & 0.010850 \
0.2 & 0.023674 & 0.026330 & 0.023130 \
0.3 & 0.038072 & 0.042154 & 0.037225 \
0.4 & 0.054907 & 0.060428 & 0.053650 \
0.5 & 0.074887 & 0.081985 & 0.073250 \
0.6 & 0.099368 & 0.108159 & 0.097350 \
0.7 & 0.130854 & 0.141331 & 0.128300 \
0.8 & 0.174203 & 0.186437 & 0.170950 \
0.9 & 0.242603 & 0.257213 & 0.239000 \
\hline
\end{tabular}
TABLE IV
(S_{4}(N)) Nodes and Weights for (c<1)
\begin{tabular}{ccccc}
\hline(c) & (\mu_{1}) & (w_{1}) & (\mu_{2}) & (w_{2}) \
\hline 0.10 & 0.330375 & 0.740076 & 0.985724 & 0.259924 \
0.20 & 0.325227 & 0.725452 & 0.966760 & 0.274548 \
0.30 & 0.317574 & 0.704988 & 0.942809 & 0.295012 \
0.40 & 0.308673 & 0.681086 & 0.917460 & 0.318914 \
0.50 & 0.298424 & 0.655816 & 0.893746 & 0.344184 \
0.60 & 0.288074 & 0.631214 & 0.872827 & 0.368786 \
0.70 & 0.278866 & 0.608724 & 0.854945 & 0.391276 \
0.80 & 0.270240 & 0.588102 & 0.839637 & 0.411898 \
0.90 & 0.259211 & 0.566686 & 0.825466 & 0.433314 \
\hline
\end{tabular}
TABLE V
(S_{6}(N)) Nodes and Weights for (c<1)
\begin{tabular}{ccccccc}
\hline(c) & (\mu_{1}) & (w_{1}) & (\mu_{2}) & (w_{2}) & (\mu_{3}) & (w_{3}) \
\hline 0.10 & 0.994388 & 0.096946 & 0.684753 & 0.47644 & 0.181644 & 0.426614 \
0.20 & 0.985658 & 0.105738 & 0.678209 & 0.471478 & 0.180278 & 0.422782 \
0.30 & 0.973165 & 0.119710 & 0.668147 & 0.462636 & 0.179330 & 0.417654 \
0.40 & 0.959536 & 0.137454 & 0.654237 & 0.453522 & 0.175927 & 0.409024 \
0.50 & 0.947301 & 0.155588 & 0.639544 & 0.444428 & 0.172727 & 0.399984 \
0.60 & 0.937124 & 0.172482 & 0.624848 & 0.437270 & 0.168904 & 0.390248 \
0.70 & 0.928910 & 0.187308 & 0.611248 & 0.431624 & 0.165553 & 0.381066 \
0.80 & 0.922232 & 0.200224 & 0.598451 & 0.428086 & 0.161744 & 0.371688 \
0.90 & 0.916787 & 0.211298 & 0.586816 & 0.425914 & 0.158145 & 0.362788
\end{tabular}
constraint equations (28) and (35) ((c<1)) or (40) ((c>1)) have been used. The nodes and weights for (S_{4}, S_{6}), and (S_{8}) of case (c) (which are denoted as (S_{4}(N), S_{6}(N)), and (S_{8}(N)) ) for (c<1) are given in Tables IV, V, VI, respectively. The nodes and weights for (S_{4}(N), S_{6}(N)) for (c>1) are given in Tables VII, VIII, respectively. It should be noted that handcalculated eigenvalues were used in the derivation of (S_{4}(N)) for (c<1), whereas a numerical subroutine from the IMSL math library was used for the calculation of the eigenvalues for the remaining quadrature sets. Numerical results are provided in Tables IX through XIX for the following problems:
(1) Leakage for the half-space albedo problem.
(2) Integrated scalar flux for an albedo problem in a slab.
(3) Integrated scalar flux for a constant source problem in a slab.
(4) Scalar flux for an albedo problem.
(5) Scalar flux for a constant source problem.
(6) Leakage for a constant source problem in a slab with linear anisotropic scattering and nonentry of neutrons as the boundary condition on both the left and right sides.
(7) Leakage for a constant source problem in a slab with linear anisotropic scattering, nonentry of neutrons as
TABLE VI
(S_{8}(N)) Nodes and Weights for (c<1)
\begin{tabular}{ccccccccc}
\hline(c) & (\mu_{1}) & (w_{1}) & (\mu_{2}) & (w_{2}) & (\mu_{3}) & (w_{3}) & (\mu_{4}) & (w_{4}) \
\hline 0.10 & 0.997071 & 0.048594 & 0.835109 & 0.262756 & 0.498146 & 0.391364 & 0.125158 & 0.297284 \
0.20 & 0.991968 & 0.054094 & 0.830228 & 0.261342 & 0.495092 & 0.389264 & 0.124156 & 0.295298 \
0.30 & 0.984077 & 0.063764 & 0.821779 & 0.258352 & 0.491005 & 0.383776 & 0.124929 & 0.294108 \
0.40 & 0.975599 & 0.076162 & 0.810052 & 0.256078 & 0.483508 & 0.378312 & 0.122817 & 0.289448 \
0.50 & 0.968430 & 0.088310 & 0.797765 & 0.254876 & 0.475649 & 0.371696 & 0.121449 & 0.285116 \
0.60 & 0.962854 & 0.098900 & 0.785982 & 0.255558 & 0.467092 & 0.365714 & 0.119234 & 0.279828 \
0.70 & 0.958628 & 0.107568 & 0.775527 & 0.257422 & 0.459111 & 0.359814 & 0.117710 & 0.275194 \
0.80 & 0.955364 & 0.114656 & 0.766210 & 0.260410 & 0.450938 & 0.355122 & 0.115366 & 0.269812 \
0.90 & 0.952816 & 0.120402 & 0.758137 & 0.263852 & 0.443267 & 0.351062 & 0.113198 & 0.264682 \
\hline
\end{tabular}
TABLE VII
(S_{4}(N)) Nodes and Weights for (c>1)
\begin{tabular}{ccccc}
\hline(c) & (\mu_{1}) & (w_{1}) & (\mu_{2}) & (w_{2}) \
\hline 1.10 & 0.807019 & 0.459138 & 0.251848 & 0.540862 \
1.20 & 0.796999 & 0.476446 & 0.242117 & 0.523554 \
1.30 & 0.786691 & 0.495270 & 0.230515 & 0.504730 \
1.40 & 0.776626 & 0.514180 & 0.218558 & 0.485820 \
1.50 & 0.767001 & 0.532636 & 0.206806 & 0.467364 \
1.60 & 0.757890 & 0.550406 & 0.195493 & 0.449594 \
1.70 & 0.749470 & 0.567036 & 0.185048 & 0.432964 \
1.80 & 0.741510 & 0.583004 & 0.175040 & 0.416996 \
1.90 & 0.734099 & 0.598066 & 0.165695 & 0.401934 \
\hline
\end{tabular}
TABLE VIII
(S_{6}(N)) Nodes and Weights for (c>1)
\begin{tabular}{ccccccc}
\hline(c) & (\mu_{1}) & (w_{1}) & (\mu_{2}) & (w_{2}) & (\mu_{3}) & (w_{3}) \
\hline 1.10 & 0.908875 & 0.228128 & 0.568489 & 0.422570 & 0.154347 & 0.349302 \
1.20 & 0.905896 & 0.234728 & 0.560797 & 0.421948 & 0.152694 & 0.343324 \
1.30 & 0.903437 & 0.240274 & 0.554136 & 0.421646 & 0.151393 & 0.338080 \
1.40 & 0.901259 & 0.245290 & 0.547713 & 0.422224 & 0.149493 & 0.332488 \
1.50 & 0.899255 & 0.249994 & 0.541271 & 0.423724 & 0.146819 & 0.326282 \
1.60 & 0.897375 & 0.254476 & 0.534783 & 0.425982 & 0.143504 & 0.319542 \
1.70 & 0.895698 & 0.258514 & 0.528777 & 0.428336 & 0.140349 & 0.313148 \
1.80 & 0.894075 & 0.262464 & 0.522639 & 0.431298 & 0.136630 & 0.306236 \
1.90 & 0.892572 & 0.266154 & 0.156753 & 0.434418 & 0.132899 & 0.299428 \
\hline
\end{tabular}
the boundary condition on the left side, and reflecting boundary condition on the right side.
(8) Leakage for an albedo problem in a slab with linear anisotropic scattering and reflecting boundary condition on the right side.
(9) Leakage for a constant source problem in a slab with anisotropic scattering, nonentry of neutrons as the boundary condition on the left side, and reflecting boundary condition on the right side.
(10) Leakage for an albedo problem in a slab with anisotropic scattering and reflecting boundary condition on the right side.
(11) Multiplication factor for a criticality problem with critical thickness (b_{c}).
All of the results for the above problems are compared with the Gauss quadrature and the exact values. When the exact results are not available, (S_{32}) results have been taken as exact. The slab thickness in all the calculations is taken as 20 mean-free paths. The leakage is calculated from (34). The scalar flux is given by (\int_{-1}^{1} \psi(x, \mu) d \mu) and the integrated scalar flux is calculated from the formula (\int_{0}^{L} \int_{-1}^{1} \psi(x, \mu) d \mu d x), where (L) is equal to 20 mean-free paths. In problems (6)-(9), the scattering kernel (K\left(x, \mu, \mu^{\prime}\right)=\frac{1}{2}+\mu \mu^{\prime}) and in problems 10 and 11 ,
TABLE IX
Leakage for the Half-Space Albedo Problem
\begin{tabular}{cccccccc}
\hline(c) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & Exact leakage \
\hline 0.10 & 0.0115661 & 0.0124071 & 0.0109877 & 0.0116510 & 0.0108974 & 0.0113377 & 0.010850 \
0.20 & 0.0245593 & 0.0263301 & 0.0234453 & 0.0247684 & 0.0232181 & 0.0241255 & 0.023130 \
0.30 & 0.0393260 & 0.0421535 & 0.0376645 & 0.0397297 & 0.0373955 & 0.0387396 & 0.037225 \
0.40 & 0.0564232 & 0.0604284 & 0.0542206 & 0.0570772 & 0.0538664 & 0.0557198 & 0.053650 \
0.50 & 0.0765743 & 0.0819845 & 0.0739577 & 0.0776296 & 0.0735283 & 0.0758823 & 0.073250 \
0.60 & 0.1011219 & 0.1081591 & 0.0981382 & 0.1027095 & 0.0976370 & 0.1005466 & 0.097350 \
0.70 & 0.1325400 & 0.1413313 & 0.1297183 & 0.1346761 & 0.1286126 & 0.1320680 & 0.128300 \
0.80 & 0.1756702 & 0.1864365 & 0.1718959 & 0.1784364 & 0.1712688 & 0.1753483 & 0.170950 \
0.90 & 0.2436655 & 0.2572127 & 0.2400296 & 0.2476838 & 0.2393398 & 0.2440739 & 0.239000 \
\hline
\end{tabular}
TABLE X
Integrated Scalar Flux for an Albedo Problem in a Slab
\begin{tabular}{cccccccc}
\hline(c) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & (S_{32}) \
\hline 0.10 & 11.38313 & 11.82113 & 11.38313 & 11.58633 & 11.38307 & 11.50033 & 11.39056 \
0.20 & 12.44780 & 12.90123 & 12.45725 & 12.65723 & 12.44776 & 12.56870 & 12.45690 \
0.30 & 13.76440 & 14.23336 & 13.76840 & 13.97913 & 13.76837 & 13.88788 & 13.77386 \
0.40 & 15.44340 & 15.92464 & 15.44137 & 15.65895 & 15.44136 & 15.56475 & 15.44836 \
0.50 & 17.65553 & 18.15677 & 17.65551 & 17.87808 & 17.65537 & 17.78060 & 17.66162 \
0.60 & 20.74179 & 21.26709 & 20.74374 & 20.97349 & 20.74368 & 20.87240 & 20.75064 \
0.70 & 25.42651 & 25.97333 & 25.53009 & 25.66215 & 25.42625 & 25.55715 & 25.43207 \
0.80 & 33.59562 & 34.16811 & 33.59119 & 33.83578 & 33.59141 & 33.72620 & 33.59667 \
0.90 & 52.90842 & 53.53645 & 52.92432 & 53.17741 & 52.92370 & 53.06272 & 52.92695 \
\hline
\end{tabular}
TABLE XI
Integrated Scalar Flux for a Constant Source Problem in a Slab
\begin{tabular}{rrrrrrrr}
\hline(c) & (S_{4}(N)) & \multicolumn{1}{c}{(S_{4})} & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & (S_{32}) \
\hline 0.10 & 876.250 & 875.746 & 876.235 & 876.005 & 876.252 & 876.106 & 876.185 \
0.20 & 984.437 & 983.872 & 984.459 & 984.176 & 984.437 & 984.284 & 984.390 \
0.30 & 1123.191 & 1122.527 & 1123.172 & 1122.891 & 1123.163 & 1122.991 & 1123.123 \
0.40 & 1307.575 & 1306.787 & 1307.599 & 1307.234 & 1307.593 & 1307.368 & 1307.549 \
0.50 & 1564.672 & 1563.668 & 1564.638 & 1564.224 & 1564.674 & 1564.386 & 1564.606 \
0.60 & 1948.124 & 1946.811 & 1948.078 & 1947.541 & 1956.021 & 1947.789 & 1948.006 \
0.70 & 2581.809 & 2579.982 & 2598.853 & 2580.971 & 2581.763 & 2581.365 & 2581.641 \
0.80 & 3831.667 & 3828.736 & 3831.603 & 3830.412 & 3831.689 & 3830.995 & 3831.291 \
0.90 & 7469.030 & 7462.700 & 7468.870 & 7466.208 & 7468.734 & 7467.301 & 7467.402 \
\hline
\end{tabular}
TABLE XIIa
Scalar Flux for an Albedo Problem, (c=0.2) and (c=0.4)
\begin{tabular}{rccccccc}
\hline & \multicolumn{3}{c}{(c=0.2)} & & \multicolumn{3}{c}{(c=0.4)} \
\cline { 2 - 7 } Position & (S_{4}) & (S_{4}(N)) & (S_{32}) & (S_{4}) & (S_{4}(N)) & (S_{32}) \
\hline 1 & 1.0557 & 1.0557 & 1.0557 & 1.1270 & 1.1270 & 1.1270 \
2 & 0.9527 & 0.9447 & 0.8954 & 1.0316 & 1.0229 & 0.9805 \
5 & 0.7070 & 0.6848 & 0.6435 & 0.7979 & 0.7742 & 0.7376 \
10 & 0.4447 & 0.4177 & 0.4183 & 0.5345 & 0.5062 & 0.5056 \
15 & 0.2911 & 0.2686 & 0.2857 & 0.3694 & 0.3463 & 0.3605 \
\hline
\end{tabular}
TABLE XIIb
Scalar Flux for an Albedo Problem, (c=0.6) and (c=0.8)
\begin{tabular}{rccccccc}
\hline & \multicolumn{4}{c}{(c=0.6)} & & \multicolumn{3}{c}{(c=0.8)} \
\cline { 2 - 4 } \cline { 6 - 8 } Position & (S_{4}) & (S_{4}(N)) & (S_{32}) & & (S_{4}) & (S_{4}(N)) & (S_{32}) \
\hline 1 & 1.2251 & 1.2251 & 1.2251 & & 1.3819 & 1.3819 & 1.3819 \
2 & 1.1405 & 1.1313 & 1.0973 & & 1.3146 & 1.3066 & 1.2829 \
5 & 0.9262 & 0.9012 & 0.8715 & & 1.1369 & 1.1147 & 1.0932 \
10 & 0.6686 & 0.6391 & 0.6382 & & 0.9047 & 0.8778 & 0.8759 \
15 & 0.4937 & 0.4698 & 0.4810 & & 0.7299 & 0.7072 & 0.7141 \
\hline
\end{tabular}
TABLE XIIIa
Scalar Flux for a Constant Source Problem, (c=0.2) and (c=0.4)
\begin{tabular}{rccccccc}
\hline & \multicolumn{4}{c}{(c=0.2)} & & \multicolumn{3}{c}{(c=0.4)} \
\cline { 3 - 4 } \cline { 6 - 8 } Position & (S_{4}) & (S_{4}(N)) & (S_{32}) & & (S_{4}) & (S_{4}(N)) & (S_{32}) \
\hline 1 & 1.1803 & 1.1803 & 1.1803 & & 1.4550 & 1.4550 & 1.4550 \
2 & 1.3092 & 1.3192 & 1.3807 & & 1.6140 & 1.6284 & 1.6991 \
5 & 1.6163 & 1.6439 & 1.6956 & & 2.0034 & 2.0429 & 2.1040 \
10 & 1.9442 & 1.9778 & 1.9771 & & 2.4426 & 2.4896 & 2.4906 \
15 & 2.1361 & 2.1642 & 2.1428 & & 2.7177 & 2.7562 & 2.7324 \
\hline
\end{tabular}
TABLE XIIIb
Scalar Flux for a Constant Source Problem, (c=0.6) and (c=0.8)
\begin{tabular}{rccccccc}
\hline & \multicolumn{4}{c}{(c=0.6)} & & \multicolumn{3}{c}{(c=0.8)} \
\cline { 2 - 4 } \cline { 6 - 8 } Position & (S_{4}) & (S_{4}(N)) & (S_{32}) & & (S_{4}) & (S_{4}(N)) & (S_{32}) \
\hline 1 & 1.9371 & 1.9371 & 1.9371 & & 3.0901 & 3.0901 & 3.0900 \
2 & 2.1488 & 2.1718 & 2.2565 & & 3.4269 & 3.4670 & 3.5848 \
5 & 2.6844 & 2.7476 & 2.8210 & & 4.3152 & 4.4260 & 4.5331 \
10 & 3.3284 & 3.4023 & 3.4042 & & 5.4758 & 5.6107 & 5.6195 \
15 & 3.7656 & 3.8254 & 3.7971 & & 6.3500 & 6.4635 & 6.4284 \
\hline
\end{tabular}
TABLE XIV
Leakage for a Constant Source Problem in a Slab with Anisotropic Scattering and Nonentry of Neutrons as the Boundary Condition on Both Left and Right Sides
\begin{tabular}{cccccccc}
\hline(c) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & (S_{32}) \
\hline 0.10 & 0.551932 & 0.574550 & 0.552338 & 0.562613 & 0.552384 & 0.558270 & 0.552790 \
0.20 & 0.615846 & 0.640068 & 0.616694 & 0.627441 & 0.616787 & 0.622916 & 0.617305 \
0.30 & 0.696676 & 0.722632 & 0.698187 & 0.709193 & 0.698321 & 0.704456 & 0.698685 \
0.40 & 0.802326 & 0.829984 & 0.803959 & 0.815568 & 0.804129 & 0.810574 & 0.804603 \
0.50 & 0.945443 & 0.975450 & 0.947617 & 0.959825 & 0.947812 & 0.954511 & 0.948281 \
0.60 & 1.184133 & 1.184147 & 1.153893 & 1.166964 & 1.154108 & 1.161231 & 1.154653 \
0.70 & 1.473440 & 1.509830 & 1.476379 & 1.490522 & 1.476591 & 1.484210 & 1.477134 \
0.80 & 2.050953 & 2.092121 & 2.053865 & 2.069708 & 2.054072 & 2.062527 & 2.054698 \
0.90 & 3.391191 & 3.440495 & 3.395425 & 3.413518 & 3.395517 & 3.405003 & 3.396097 \
\hline
\end{tabular}
TABLE XV
Leakage for a Constant Source Problem in a Slab with Anisotropic Scattering and Nonentry of Neutrons as the Boundary Condition on the Left Side, and Reflecting Boundary Condition on the Right Side
\begin{tabular}{cccccccc}
\hline(c) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & (S_{32}) \
\hline 0.10 & 0.551950 & 0.574556 & 0.552345 & 0.562619 & 0.552391 & 0.558276 & 0.552797 \
0.20 & 0.615877 & 0.640080 & 0.616709 & 0.627454 & 0.616800 & 0.622930 & 0.617319 \
0.30 & 0.696732 & 0.722662 & 0.698219 & 0.709224 & 0.698352 & 0.704487 & 0.698716 \
0.40 & 0.802442 & 0.830061 & 0.804039 & 0.815647 & 0.804207 & 0.810653 & 0.804682 \
0.50 & 0.945732 & 0.975677 & 0.947846 & 0.960054 & 0.948039 & 0.954739 & 0.948509 \
0.60 & 1.184910 & 1.184925 & 1.154663 & 1.167737 & 1.154877 & 1.162003 & 1.155422 \
0.70 & 1.476776 & 1.513026 & 1.479522 & 1.493687 & 1.479733 & 1.487364 & 1.480279 \
0.80 & 2.068188 & 2.109335 & 2.070754 & 2.086726 & 2.070962 & 2.079482 & 2.071597 \
0.90 & 3.541750 & 3.593372 & 3.545581 & 3.564725 & 3.545666 & 3.555672 & 3.546287 \
\hline
\end{tabular}
TABLE XVI
Leakage for an Albedo Problem in a Slab with Anisotropic Scattering and Reflecting Boundary Condition on the Right Side
\begin{tabular}{cccccccc}
\hline(c) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & (S_{32}) \
\hline 0.10 & 0.003961 & 0.004171 & 0.003027 & 0.003574 & 0.002895 & 0.003296 & 0.002859 \
0.20 & 0.008657 & 0.009206 & 0.006835 & 0.007967 & 0.006577 & 0.007399 & 0.006520 \
0.30 & 0.014313 & 0.015408 & 0.011751 & 0.013472 & 0.011389 & 0.012601 & 0.011273 \
0.40 & 0.021358 & 0.023234 & 0.018138 & 0.020540 & 0.017683 & 0.019347 & 0.017563 \
0.50 & 0.030458 & 0.033432 & 0.026785 & 0.029900 & 0.026259 & 0.028367 & 0.026114 \
0.60 & 0.042862 & 0.047300 & 0.038915 & 0.042829 & 0.038329 & 0.040931 & 0.038195 \
0.70 & 0.061239 & 0.067362 & 0.057051 & 0.061815 & 0.056423 & 0.059515 & 0.056273 \
0.80 & 0.09135 & 0.099403 & 0.086805 & 0.092569 & 0.086139 & 0.089815 & 0.086025 \
0.90 & 0.150341 & 0.161873 & 0.146402 & 0.153366 & 0.145698 & 0.150054 & 0.145625 \
\hline
\end{tabular}
TABLE XVII
Leakage for a Constant Source Problem in a Slab with Anisotropic Scattering and Nonentry of Neutrons as the Boundary Condition on the Left Side, and Reflecting Boundary on the Right Side
\begin{tabular}{cccccccc}
\hline(c) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & (S_{32}) \
\hline 0.1 & 0.538861 & 0.562792 & 0.540037 & 0.550572 & 0.540020 & 0.546110 & 0.540398 \
0.2 & 0.585487 & 0.612655 & 0.588052 & 0.599374 & 0.588018 & 0.594577 & 0.588480 \
0.3 & 0.643038 & 0.673877 & 0.647348 & 0.659306 & 0.647287 & 0.654106 & 0.647544 \
0.4 & 0.716408 & 0.751215 & 0.721956 & 0.735036 & 0.721855 & 0.729336 & 0.722199 \
0.5 & 0.812607 & 0.852700 & 0.820075 & 0.834439 & 0.819924 & 0.828100 & 0.820224 \
0.6 & 0.993157 & 0.993181 & 0.955868 & 0.972088 & 0.955658 & 0.964889 & 0.956010 \
0.7 & 1.147261 & 1.204057 & 1.160014 & 1.178800 & 1.159703 & 1.170349 & 1.160002 \
0.8 & 1.494551 & 1.567729 & 1.512212 & 1.535468 & 1.511753 & 1.524945 & 1.512123 \
0.9 & 2.299457 & 2.415617 & 2.334346 & 2.367481 & 2.333514 & 2.352399 & 2.333905 \
\hline
\end{tabular}
TABLE XVIII
Leakage for an Albedo Problem in a Slab with Anisotropic Scattering and Reflecting Boundary Condition on the Right Side
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline(c) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) & (S_{8}(N)) & (S_{8}) & (S_{32}) \
\hline 0.1 & 0.014164 & 0.014759 & 0.014103 & 0.014426 & 0.014027 & 0.014265 & 0.014018 \
\hline 0.2 & 0.029886 & 0.031150 & 0.029757 & 0.030441 & 0.029603 & 0.030103 & 0.029591 \
\hline 0.3 & 0.047564 & 0.049563 & 0.047358 & 0.048425 & 0.047134 & 0.047891 & 0.047092 \
\hline 0.4 & 0.067761 & 0.070550 & 0.067386 & 0.068916 & 0.067092 & 0.068163 & 0.067051 \
\hline 0.5 & 0.091208 & 0.094933 & 0.090668 & 0.092717 & 0.090312 & 0.091715 & 0.090256 \
\hline 0.6 & 0.124012 & 0.124016 & 0.118427 & 0.121099 & 0.118013 & 0.119810 & 0.117960 \
\hline 0.7 & 0.153892 & 0.160078 & 0.152899 & 0.156291 & 0.152428 & 0.154662 & 0.152357 \
\hline 0.8 & 0.199852 & 0.207762 & 0.198513 & 0.202833 & 0.197980 & 0.200779 & 0.197923 \
\hline 0.9 & 0.268913 & 0.279739 & 0.267551 & 0.273132 & 0.266943 & 0.270496 & 0.266900 \
\hline
\end{tabular}
(K\left(x, \mu, \mu^{\prime}\right)=\sin ^{2}\left(\mu \mu^{\prime}\right) /(1-\sin 2 \mu / 2 \mu)). Both these kernels satisfy the normalizing condition (\int_{-1}^{1} K\left(x, \mu, \mu^{\prime}\right) d \mu^{\prime}=1).
In these test problems, we observe that the new quadrature set performs significantly better than the Gauss quadrature set. In many problems, we find that the (S_{4}(N)) is comparable to (S_{32}) with the traditional Gauss-Legendre quadrature. Also, the results compare favorably with those of Refs. ([6,7]). By replacing the last two equations for (n=2 N-2) and (n=2 N-4) of (1) by the constraints ((28), (35), or (40)), the numerical results are significantly improved as expected from transport-theoretic considerations. However, it has not yet been justified from the numerical analysis point of view.
BΓΌckner [19] observes that, if the (m) th derivative of the integrand is discontinuous, there is no advantage in using a formula which is exact for polynomials of order higher than (m-1). This may partially explain why our nodes and weights perform well in transport calculations. It is also interesting to note that our approach yields significantly improved results for anisotropic scattering cases.
Our final conclusion is that the new quadrature set may
TABLE XIX
Multiplication Factor for a Criticality Problem with Critical Thickness (b_{c})
\begin{tabular}{cccccc}
\hline(c) & (b_{c}) & (S_{4}(N)) & (S_{4}) & (S_{6}(N)) & (S_{6}) \
\hline 1.1 & 4.2266 & 0.9999955 & 0.9985735 & 0.9999996 & 0.9995073 \
1.2 & 2.5786 & 0.9999129 & 0.9948241 & 0.9999955 & 0.9983894 \
1.3 & 1.8755 & 0.9997404 & 0.9886890 & 1.0000320 & 0.9965505 \
1.4 & 1.4732 & 0.9994403 & 0.9804660 & 0.9999900 & 0.9937381 \
1.5 & 1.2101 & 0.9991258 & 0.9707364 & 0.9998697 & 0.9899453 \
1.6 & 1.0240 & 0.9988502 & 0.9600388 & 0.9997062 & 0.9852737 \
1.7 & 0.8851 & 0.9986270 & 0.9486828 & 0.9993250 & 0.9797643 \
1.8 & 0.7774 & 0.9984134 & 0.9369345 & 0.9988465 & 0.9735298 \
1.9 & 0.6919 & 0.9984993 & 0.9253452 & 0.9985121 & 0.9670290 \
\hline
\end{tabular}
enable numerical transport calculations for some problems with fewer points and weights than quadrature sets currently used. However, it seems that the new method has certain disadvantages in more realistic problems than the ones considered in this paper. In Appendix A, we have considered the interface problem with the medium dependent quadrature set. The interface condition was achieved by equating the moments of the angular flux at the interface and the overall performance of the new method is better. The implementation of the method to other geometries and multigroup transport would be relatively difficult. However, it may be easier if the weights are only medium dependent. We have started looking into these aspects and the conclusion awaits further work.
APPENDIX A
In this appendix, we consider the suitability of the new quadrature sets for more realistic problems involving interfaces. With standard discrete-ordinates coefficients, the discrete directions do not change at interfaces between different media. Thus, the continuity of the angular flux provides sufficient interface conditions to determine the entering discrete angular fluxes at interfaces. However, the new quadrature coefficients are dependent on the medium and the discrete directions change abruptly at interfaces between different media. In the next section, we derive the appropriate interface conditions required to determine the entering discrete angular fluxes. Finally, we present numerical results on a four region interface problem and a fuel moderator cell problem for (S_{4}) approximation.
The present nodes and weights are dependent on the medium and hence, interface conditions need to be formulated in the discrete ordinate scheme to take advantage of the new set of points. We consider a slab of thickness (L) with an interface at (x_{0}). Also, we assume that the constant (c),
which is the average number of secondaries produced per collision, is different in each of these regions, namely, (c_{1}) and (c_{2}) such that (00,1 \leqslant i \leqslant n), has an inverse; hence the solution to this system can be written as
[
\begin{aligned}
& =\left[\begin{array}{cccc}
\frac{1}{\bar{w}_{1}} & & & \
& \frac{1}{\bar{w}_{2}} & & \
& & \ddots & \
& & & 1 \
& & & \frac{1}{\bar{w}_{n}}
\end{array}\right] \cdot\left[\begin{array}{c}
1 \cdots 1 \
\bar{\mu}_{1} \cdots \bar{\mu}_{n} \
\vdots \
\bar{\mu}_{1}^{n-1} \cdots \tilde{\mu}_{n}^{n-1}
\end{array}\right] \
& \cdot\left[\begin{array}{c}
w_{1} \psi\left(x_{0}, \mu_{1}\right)+\cdots+w_{n} \psi\left(x_{0}, \mu_{n}\right) \
w_{1} \mu_{1} \psi\left(x_{0}, \mu_{1}\right)+\cdots+w_{n} \mu_{n} \psi\left(x_{0}, \mu_{n}\right) \
\vdots \
w_{n} \mu_{1}^{n-1} \psi\left(x_{0}, \mu_{1}\right)+\cdots+w_{n} \mu_{n}^{n-1} \psi\left(x_{0}, \mu_{n}\right)
\end{array}\right]
\end{aligned}
]
The sweep to the left has a similar interface condition, where the summation in (A1) is from (N / 2+1) to (N).
We will report the transport calculations by the new quadrature sets obtained by the approach (A) of this paper and approach (B) of Ref. [10]. We have, however, experienced some difficulties in approach (B) in obtaining a converged solution of the nonlinear system of equations involving the parameters of the weight function for the case of a multiplying medium. We had to also deviate slightly from requiring that one of the eigenvalues of the discreteordinates equations by the exact asymptotic eigenvalue for highly absorbing cases to ensure positivity of the weight function.
In this section, we present the numerical results corresponding to (N=4). For notational convenience, we call the quadrature sets obtained by approaches (A) and (B) by set I and set II, respectively. Correspondingly, we denote our results by (S_{4}\left(N_{1}\right)) and (S_{4}\left(N_{2}\right)). We now consider several problems involving interfaces to test the performance of the new quadrature sets. We compare our results with those of Gauss-Legendre set for (N=4) and (N=16). Since exact results are not available in most cases, we choose the conventional (S_{N}) results with (N=16) for the purpose of comparison.
Our test problem is a four region source problem in a nonmultiplying medium with both isotropic scattering (\left(K\left(x, \mu, \mu^{\prime}\right)=\frac{1}{2}\right)) and linearly anisotropic scattering (\left(K\left(x, \mu, \mu^{\prime}=\frac{1}{2}+\mu \mu^{\prime}\right)\right.). The scattering kernel (K\left(x, \mu, \mu^{\prime}\right)) satisfies the normalization condition (\int_{-1}^{1} K\left(x, \mu, \mu^{\prime}\right) d \mu=1).
FIG. AI. The model problem.
In Fig. AI, we provide the appropriate data for each of the four regions with spatial dimensions indicated in mean-free paths. Vacuum boundary conditions are imposed at both the left and right boundaries of the slab of thickness (L).
The source combination given here will be referred to as Q0101 which means that the constant source of magnitude unity is placed in the second and fourth regions and there are no sources in first and third regions. We have also considered various other source combinations and these are represented in a similar way. We compute the scalar flux in different regions for different source combinations, with both isotropic and linearly anisotropic scattering. The results for a typical case with linear anisotropic scattering and source combination Q0101 are shown graphically in Fig. AII for (S_{4}(N), S_{4}), and (S_{16}) cases. As noted, the difference of our results with the usual Gauss quadrature results is hardly seen on these figures. Therefore, we compare the integrated flux
[
\int_{0}^{L} \int_{-1}^{1} \psi(x, \mu) d \mu d x
]
for different sets with different source combinations.
FIG. AII. Scalar flux for anisotropic case with source combination Q0101.
TABLE AI
Integrated Scalar Flux for the Isotropic Case
\begin{tabular}{crrrr}
\hline & \multicolumn{4}{c}{ Integrated scalar flux } \
\cline { 2 - 5 } Cource & \multicolumn{1}{c}{(S_{4})} & \multicolumn{1}{c}{(S_{4}\left(N_{1}\right))} & \multicolumn{1}{c}{(S_{4}\left(N_{2}\right))} & \multicolumn{1}{c}{(S_{16})} \
\hline Q0101 & 1572.05 & 1581.98 & 1581.77 & 1580.18 \
Q1011 & 2486.59 & 2497.30 & 2495.45 & 2495.86 \
Q1110 & 1075.78 & 1075.97 & 1074.82 & 1075.27 \
Q1010 & 995.16 & 995.64 & 994.25 & 995.47 \
Q0100 & 80.63 & 80.33 & 80.57 & 79.80 \
Q0001 & 1491.42 & 1501.65 & 1501.20 & 1500.38 \
Q1111 & 2567.21 & 2577.63 & 2576.02 & 2575.65 \
\hline
\end{tabular}
In Tables AI and AII we give the values of integrated scalar flux for isotropic and anisotropic cases for the (S_{4}) and (S_{16}) methods. Finally, we compute the leakage from the right and left sides for both isotropic and anisotropic scattering cases. These are
[
L_{\text {right }}=\int_{0}^{1} \mu \psi(L, \mu) d \mu
]
and
[
L_{\text {left }}=\int_{-1}^{0} \mu \psi(0, \mu) d \mu
]
The right and left leakage values are given in Tables AIII(a) and (b) and AIV(a) and (b) for isotropic and anisotropic scattering cases, respectively.
We observe that the results on integrated scalar flux are, in general, improved with the new quadrature coefficients. Also, the set I performs better overall. However, there are a few cases where the set II is superior. The improvement for the anisotropic scattering case is more striking. The overall results on the leakages from the left and right sides are also improved with the new quadrature coefficients. However, the leakage in some cases for the set I worsens slightly, compared to that of conventional Gauss-quadrature. In a few cases, set II performs better. We now consider the
TABLE AII
Integrated Scalar Flux for the Anisotropic Case
\begin{tabular}{crrrr}
\hline & \multicolumn{4}{c}{ Integrated scalar flux } \
\cline { 2 - 5 } Source & \multicolumn{1}{c}{(S_{4})} & (S_{4}\left(N_{1}\right)) & (S_{4}\left(N_{2}\right)) & \multicolumn{1}{c}{(S_{16})} \
\hline Q0101 & 1366.14 & 1380.70 & 1383.21 & 1376.93 \
Q1011 & 2291.95 & 2307.13 & 2307.40 & 2303.95 \
Q1110 & 1086.05 & 1086.07 & 1084.43 & 1085.58 \
Q1010 & 1005.93 & 1006.25 & 1004.31 & 1006.30 \
Q0100 & 80.12 & 79.82 & 80.12 & 79.28 \
Q0001 & 1286.02 & 1300.88 & 1303.09 & 1297.65 \
Q1111 & 2372.07 & 2386.95 & 2387.52 & 2383.23 \
\hline
\end{tabular}
TABLE AIII(a)
Leakage from the Right Side for Isotropic Scattering
\begin{tabular}{|c|c|c|c|c|}
\hline \multirow{2}{*}{\begin{tabular}{c}
Source \
Combinations
\end{tabular}} & \multicolumn{4}{|c|}{ Leakage from the right side } \
\hline & (S_{4}) & (S_{4}\left(N_{1}\right)) & (S_{4}\left(N_{2}\right)) & (S_{16}) \
\hline Q0101 & 2.486 & 2.456 & 2.458 & 2.459 \
\hline Q1011 & 2.507 & 2.477 & 2.459 & 2.480 \
\hline Q1110 & (2.183 \times 10^{-2}) & (2.147 \times 10^{-2}) & (2.071 \times 10^{-2}) & (2.157 \times 10^{-2}) \
\hline
\end{tabular}
TABLE AIII(b)
Leakage from the Left Side for Isotropic Scattering
\begin{tabular}{|c|c|c|c|c|}
\hline \multirow{2}{*}{\begin{tabular}{c}
Source \
Combinations
\end{tabular}} & \multicolumn{4}{|c|}{ Leakage from the left side } \
\hline & (S_{4}) & (S_{4}\left(N_{1}\right)) & (S_{4}\left(N_{2}\right)) & (S_{16}) \
\hline Q0101 & (9.046 \times 10^{-2}) & (1.061 \times 10^{-1}) & (9.982 \times 10^{-2}) & (9.122 \times 10^{-2}) \
\hline Q1011 & (7.292 \times 10^{-1}) & (7.719 \times 10^{-1}) & (7.217 \times 10^{-1}) & (7.040 \times 10^{-1}) \
\hline Q1110 & (8.197 \times 10^{-1}) & (8.780 \times 10^{-1}) & (8.115 \times 10^{-1}) & (7.953 \times 10^{-1}) \
\hline
\end{tabular}
TABLE AIV(a)
Leakage from the Right Side for Anisotropic Scattering
\begin{tabular}{|c|c|c|c|c|}
\hline \multirow{2}{*}{\begin{tabular}{c}
Source \
Combinations
\end{tabular}} & \multicolumn{4}{|c|}{ Leakage from the right side } \
\hline & (S_{4}) & (S_{4}\left(N_{1}\right)) & (S_{4}\left(N_{2}\right)) & (S_{16}) \
\hline Q0101 & 3.021 & 2.978 & 2.978 & 2.984 \
\hline Q1011 & 3.106 & 3.062 & 3.057 & 3.068 \
\hline Q1110 & (8.512 \times 10^{-2}) & (8.321 \times 10^{-2}) & (7.876 \times 10^{-2}) & (8.393 \times 10^{-2}) \
\hline
\end{tabular}
TABLE AIV(b)
Leakage from the Left Side for Anisotropic Scattering
\begin{tabular}{ccccc}
\hline & \multicolumn{4}{c}{ Leakage from the left side } \
\cline { 2 - 5 } Source \
\cline { 2 - 5 } Combinations & (S_{4}) & (S_{4}\left(N_{1}\right)) & (S_{4}\left(N_{2}\right)) & (S_{16}) \
\hline Q0101 & (1.243 \times 10^{-1}) & (1.425 \times 10^{-1}) & (1.217 \times 10^{-1}) & (1.249 \times 10^{-1}) \
Q1011 & (7.657 \times 10^{-1}) & (8.112 \times 10^{-1}) & (7.561 \times 10^{-1}) & (7.381 \times 10^{-1}) \
Q1110 & (8.900 \times 10^{-1}) & (9.537 \times 10^{-1}) & (8.778 \times 10^{-1}) & (8.631 \times 10^{-1}) \
\hline
\end{tabular}
TABLE AV
Disadvantage Factor for Different Methods
\begin{tabular}{cc}
\hline Method & Disadvantage factor \
\hline(S_{4}) & 2.28 \
(S_{4}\left(N_{1}\right)) & 2.33 \
(S_{4}\left(N_{2}\right)) & 2.29 \
(S_{16}) & 2.38 \
\hline
\end{tabular}
problem of evaluating the disadvantage factor for a fuel moderator cell.
An infinite homogeneous array of (2 w)-wide fuel plates separated by a moderator with a centerline spacing of (W) is now considered. The disadvantage factor (d) in such an array is defined as the ratio of average flux in the moderator to the average flux in the fuel, i.e.,
[
d=\frac{(1 /(W-w)) \int_{w}^{w} \overline{\phi(x)} d x}{(1 / w) \int_{0}^{w} \phi(x) d x}
]
where (\phi(x)) is the scalar flux for the fuel and (\overline{\phi(x)}) is that for the moderator. Using the following data, we compare the disadvantage factor for (S_{4}\left(N_{1}\right)) and (S_{4}\left(N_{2}\right)) with that for the Gauss-Legendre set with (N=4) and (N=16). The later two sets will be represented by (S_{4}) and (S_{16}). We now consider the following data:
[
\begin{aligned}
c_{1} & =0.55 \
c_{2} & =0.99 \
w & =0.635 \
W & =15.240
\end{aligned}
]
The results, shown in Table AV, indicate that the disadvantage factor is improved significantly by set I. Set II also performs slightly better.
The comparison of numerical results on the interface problems reveals that the present approach yields results improved over the ones of conventional Gauss quadrature. As for the comparison between set I and set II, it is found that set I performs better in many situations. The reason for this is because, to guarantee the positivity of weight function for low values of (c), we had to replace (v) by (0.9 v), and this caused some inconsistencies in the numerical results. Also the weight function in Ref. [10] is quite arbitrary as the full-range weight function is the simple odd polynomial " (\mu)." However, it appears worthwhile to consider the weight function in approach (B) of the form (W(\mu)=) (\left(1-\alpha \mu^{2}-\beta \mu^{4}\right)^{\gamma}). The weight function in this proposed form may make approach (B) a viable competing candidate with (A). However, the calculation involving the determination of weights and nodes will be tedius. The half-range and fullrange orthogonality relations may also be suitably used to derive appropriate constraint equations for obtaining different quadrature sets. We are presently studying the convergence analysis of the new discrete-ordinates method for interface problems.
APPENDIX B
In this appendix, we will describe how the constraints (28) and (57) may be profitably used in transport calcula-
TABLE BI
Quadrature Weights and Nodes
\begin{tabular}{ccccc}
\hline(c) & (w_{1}) & (w_{2}) & (\mu_{1}) & (\mu_{2}) \
\hline 0.1 & 0.7282 & 0.2717 & 0.3092 & 0.9851 \
0.2 & 0.7124 & 0.2876 & 0.3028 & 0.9654 \
0.3 & 0.6906 & 0.3094 & 0.2941 & 0.9405 \
0.4 & 0.6652 & 0.3348 & 0.2837 & 0.9141 \
0.5 & 0.6388 & 0.3612 & 0.2730 & 0.8894 \
0.6 & 0.6132 & 0.3868 & 0.2626 & 0.8675 \
0.7 & 0.5894 & 0.4106 & 0.2528 & 0.8485 \
0.8 & 0.5672 & 0.4328 & 0.2436 & 0.8321 \
0.9 & 0.5470 & 0.4530 & 0.2351 & 0.8179 \
\hline
\end{tabular}
tions. We consider the case (N=4). The quadrature coefficients are obtained from
[
\begin{aligned}
& \sum_{i=1}^{N} w_{i} \mu^{n}=\int_{-1}^{1} \mu^{n} d \mu, \quad n=0,2 \
& w_{i}=w_{N-i+1}, \quad i=1,2 \
& \mu_{i}=-\mu_{N-i+1}, \quad i=1,2
\end{aligned}
]
together with (28) and (57).
TABLE BII(a)
Leakage for the Half-Space Albedo Problem
\begin{tabular}{cccc}
\hline(c) & (S_{4}) & (S_{4}^{*}) & Exact leakage \
\hline 0.1 & 0.0124 & 0.0112 & 0.0109 \
0.2 & 0.0263 & 0.0237 & 0.0231 \
0.3 & 0.0421 & 0.0380 & 0.0372 \
0.4 & 0.0604 & 0.0545 & 0.0537 \
0.5 & 0.0820 & 0.0741 & 0.0733 \
0.6 & 0.1082 & 0.0981 & 0.0974 \
0.7 & 0.1413 & 0.1288 & 0.1283 \
0.8 & 0.1864 & 0.1714 & 0.1710 \
0.9 & 0.2572 & 0.2392 & 0.2390 \
\hline
\end{tabular}
TABLE BII(b)
Leakage for Constant Source Half-Space Problem
\begin{tabular}{cccc}
\hline(c) & (S_{4}) & (S_{4}^{*}) & Exact leakage \
\hline 0.1 & 0.5654 & 0.5351 & 0.5435 \
0.2 & 0.6187 & 0.5870 & 0.5957 \
0.3 & 0.6844 & 0.6516 & 0.6608 \
0.4 & 0.7681 & 0.7338 & 0.7436 \
0.5 & 0.8786 & 0.8431 & 0.8528 \
0.6 & 1.0327 & 0.9961 & 1.0047 \
0.7 & 1.2664 & 1.2281 & 1.2389 \
0.8 & 1.6741 & 1.6331 & 1.6459 \
0.9 & 2.6401 & 2.5794 & 2.6103 \
\hline
\end{tabular}
TABLE BIII
Disadvantage Factor for Different Methods
\begin{tabular}{cc}
\hline Method & Disadvantage factor \
\hline(S_{4}) & 2.28 \
(S_{4}^{*}) & 2.36 \
(S_{16}) & 2.38 \
\hline
\end{tabular}
Table BI gives the weights and nodes for different values of (c). We now present our results in Tables BII(a) and BII(b) on leakage for a half-space albedo and constant source problems. The improved results are denoted by (S_{4}^{*}) and the conventional results by (S_{4}). The numerical calculations reveal that the inclusion of (28) and (57) in the derivation of the quadrature set does improve the results for several benchmark problems. The results compare favorably with the (F_{N}) method, considering the numerical advantage of the discrete-ordinates scheme.
Finally, we test the new quadrature set on a more realistic problem involving interface. We consider an infinite homogeneous array of (2 w)-wide fuel plates separated by a moderator with a centerline spacing of (W). The following data- (c_{1}=0.55, c_{2}=0.99, w=0.635), and (W=15.24) are used to calculate the fuel-moderator disadvantage factor. Table BIII shows the disadvantage factor is improved by employing (28) and (57).
ACKNOWLEDGMENT
\begin{abstract}
The authors gratefully acknowledge suggestions by the referees to improve the presentation of the paper.