A Numerical-Analytical Solution of the One-Dimensional and Monoenergetic Neutron Transport Equation for Eigenvalue Problems

Uma Solução Numérico-Analítica da Equação de Transporte de Nêutrons Unidimensional e Monoenergética para Problemas de Autovalor

Abstract

The study of neutron behavior and interactions within the core of a nuclear reactor is essential for ensuring criticality control. The most accurate deterministic framework for modeling these phenomena is the Neutron Transport Equation, which, with the exception of stochastic methods such as Monte Carlo, is typically solved using fine-mesh numerical methods. However, such approaches often entail considerable computational costs to achieve high-fidelity results. In this work, we propose a numerical-analytical solution to the one-dimensional, monoenergetic neutron transport equation with isotropic scattering, formulated for eigenvalue problems and employing the discrete ordinates (SN) method. The proposed formulation is presented through a fully explicit and self-contained step-by-step derivation, emphasizing transparency and reproducibility of the numerical-analytical approach. To assess the accuracy and efficiency of the proposed approach, results were compared with those obtained from a purely numerical method developed in this study, as well as with benchmark solutions reported in the literature. In all cases, the proposed method exhibited excellent agreement with reference data, along with a significant reduction in computational time.

Keywords:
Neutron Transport Equation; Numerical-Analytical Solutions; Numerical Methods; Criticality Calculations

Resumo

O estudo do comportamento e das interações dos nêutrons no núcleo de um reator nuclear é essencial para garantir o controle da sua criticalidade. A formulação determinística mais precisa para descrever esses fenômenos é dada pela Equação de Transporte de Nêutrons, a qual, com execeção de métodos estocásticos como Monte Carlo, é tipicamente resolvida por métodos numéricos de malha fina. No entanto, esses métodos geralmente requerem custos computacionais significativos para se obter resultados de alta precisão. Neste estudo, propomos uma solução numérico-analítica para a equação de transporte de nêutrons unidimensional, monoenergética e com espalhamento isotrópico, formulada para problemas de autovalor e utilizando o método das ordenadas discretas (SN). A formulação proposta é apresentada por meio de uma dedução analítica explícita e autocontida, com ênfase na transparência e na reprodutibilidade da abordagem numérico-analítica. Para validar nossa abordagem, os resultados obtidos foram comparados com os de um método puramente numérico desenvolvido neste trabalho e com soluções de referência disponíveis na literatura. Em todos os casos, os resultados obtidos pela solução proposta demonstraram excelente acurácia se comparados com os da bibliografia, além de uma redução substancial no tempo computacional dos cálculos.

Palavras-chave:
Equação de Transporte de Nêutrons; Soluções Numérico-Analíticas; Métodos Numéricos; Cálculo de Criticalidade

1. INTRODUCTION

Nuclear reactor theory focuses on analyzing the behavior of neutrons within the material regions that compose the reactor core. The Neutron Transport Equation (NTE) is widely recognized as the most accurate deterministic formulation for describing such phenomena. However, due to its inherent complexity, even after several approximations, direct solutions remain computationally demanding.

Traditionally, neutron transport problems have been solved using fine-mesh numerical methods, such as the Diamond Difference (DD) scheme. However, numerical-analytical approaches—such as spectral nodal methods introduced by ^{Larsen in 1986}—have also been proposed and remain higly relevant. Despite their algebraic and computational complexity, these methods have shown the ability to deliver accurate results with significantly reduced computation times compared to purely numerical techniques.

In order to further advance the study of numerical-analytical solutions for the NTE, this work aims to present the formulation of a solution that encompasses eigenvalue problems in one-dimensional slab geometry, with isotropic scattering, and formulated for monoenergetic cases. The term “numerical-analytical” is used to denote approaches in which the NTE is solved analytically with respect to the spatial variable x, while the angular variable μ is discretized using the discrete ordenates (S_N) method. In contrast, the purely numerical methods, such as the DD scheme, are characterized by numerical discretization in both the spatial and angular variables.

In addition to its computational efficiency, a distinguishing feature of the present work is its emphasis on a transparent and fully explicit derivation of the numerical-analytical formulation. Unlike several existing spectral nodal or analytical SN approaches , where key steps are often condensed or presented in implicit form, this study provides a step-by-step construction of the solution. This self-contained presentation is intended to improve clarity, reproducibility, and accessibility of numerical-analytical transport methods, particularly for readers and practitioners seeking to implement such formulations from first principles.

For clarity and simplicity, the mathematical formulation is initially developed for a single homogeneous region with known boundary conditions. This simplified setting allows the fundamental aspects of the numerical-analytical approach to be presented without additional geometric or interface-related complexities. The formulation is subsequently extended to treat multi-region configurations, where region-to-region coupling and continuity conditions are explicitly addressed.

Among the references considered in this work, ^{Lahdour et al. (2019)} stands out as a recent study and is used primarily as a benchmark for monoenergetic and multi-region neutron transport calculations. The problem configurations and nuclear data adopted therein provide a consistent basis for assessing the accuracy and efficiency of the numerical-analytical formulation proposed in this study.

2. MATERIALS AND METHODS

The steady-state, monoenergetic, one-dimensional, slab geometry NTE with isotropic scattering, in a domain with uniform material properties and angular discretization via the S_N method is given by Equation (1) (^{DUDERSTADT et al., 1976}):

Eq. 01 $${\mu }_m\frac{d}{dx}{\phi }_m\left(x\right)+{\Sigma }_t{\phi }_m\left(x\right)=\frac{1}{2}({\Sigma }_{s0}+\frac{\nu {\Sigma }_f}{k_{\mathrm{eff}}})\sum _{m\prime =1}^M{\omega }_{m\prime }{\phi }_{m\prime }\left(x\right),$$

for m=1,…,M, and where M is an even number that denotes the Gauss-Legendre quadrature order (i.e. the number of discretized directions), φ_m(x) denotes the angular neutron flux along the discrete direction m; μ_m is the cosine of the angle between direction m and the x-axis; Σ_t, Σ_s0 and Σ_f, respectively, are the total cross-section, the isotropic scattering cross-section, and the fission cross-section of the region; ν is the average number of neutrons produced by fission; ω_m, the Gauss-Legendre weight functions; and k_eff is the effective multiplication factor.

As boundary conditions, we assume known the angular fluxes entering from the left and right boundaries of the region, which has width a. Hereafter, these boundary conditions shall be written, respectively, as φ_m,l and φ_m,r, and they correspond to:

Eq. 02 $${\phi }_{m,l}\equiv {\phi }_m\left(\frac{-a}{2}\right),\forall {\mu }_m>0;$$

Eq. 03 $${\phi }_{m,r}\equiv {\phi }_m\left(\frac{+a}{2}\right),\forall {\mu }_m<0.$$

These incoming fluxes may be zero, corresponding to vacuum boundary conditions, or may assume other prescribed constant values. In the multi-region formulation, the same expressions are also employed as interface conditions between adjacent regions. In such cases, however, the interface fluxes are not known a priori and must be determined as part of the global iterative solution procedure.

2.1 Numerical-Analytical Solution

The formulation presented in this section follows the standard discrete ordinates (S_N) framework used in purely numerical schemes such as the Diamond Difference method, with the angular variable being discretized into a finite set of directions. The key difference lies in the treatment of the spatial variable, which is handled analytically within each homogeneous region, resulting in a numerical-analytical formulation.

By introducing the substitution $\tilde{\Sigma }\equiv \frac{({\Sigma }_{s0}+\frac{\nu {\Sigma }_f}{k_{\mathrm{eff}}})}{2}$ into equation (1) and expanding it for all values of m, we obtain the following system:

Eq. 04 $$\frac{d}{dx}\left(\begin{array}{c}{\phi }_1\left(x\right)\\ ⋮\\ {\phi }_M\left(x\right)\end{array}\right)=\left(\begin{array}{ccc}\frac{\tilde{\Sigma }{\omega }_1-{\Sigma }_t}{{\mu }_1}& \cdots & \frac{\tilde{\Sigma }{\omega }_M}{{\mu }_1}\\ ⋮& \ddots & ⋮\\ \frac{\tilde{\Sigma }{\omega }_1}{{\mu }_M}& \cdots & \frac{\tilde{\Sigma }{\omega }_M-{\Sigma }_t}{{\mu }_M}\end{array}\right)\left(\begin{array}{c}{\phi }_1\left(x\right)\\ ⋮\\ {\phi }_M\left(x\right)\end{array}\right).$$

By defining the vectors Ψ(x)≡(φ₁(x), …,φ_M(x))^T, μ̅≡(1/μ₁, …,1/μ_M)^T and ω≡(ω₁, …,ω_M)^T the system above can be rewritten compactly as:

Eq. 05 $$\frac{d}{dx}\mathbf{\Psi }\left(\mathbf{x}\right)=\mathbf{M}\mathbf{\Psi }\left(\mathbf{x}\right),$$

with

Eq. 06 $$\mathbf{M}\equiv \tilde{\Sigma }\overline{\mathbf{\mu }}{\mathbf{\omega }}^{\mathrm{T}}-\mathrm{diag}\left({\Sigma }_t\overline{\mathbf{\mu }}\right),$$

where diag(Σ_tμ̅) denotes the diagonal and invertible M×M matrix with entries Σ_t/μ_m.

Since matrix M defines a linear transformation, we seek solutions of Ψ(x) that satisfy the eigenvalue equation:

Eq. 07 $$\mathbf{M}\mathbf{\Psi }\left(\mathbf{x}\right)=\lambda \mathbf{\Psi }\left(\mathbf{x}\right).$$

That is, solutions of Ψ(x) which, when acted upon by M, become scaled versions of themselves (^{BOYCE et al., 2015}).

Therefore, letting ξ be a constant vector, we assume a solution of the form:

Eq. 08 $$\mathbf{\Psi }\left(\mathbf{x}\right)=\mathbf{\xi }{e}^{\lambda x},$$

since substituting equation (8) into equation (5) yields:

Eq. 09 $$\lambda \mathbf{\xi }{e}^{\lambda x}=\mathbf{M}\mathbf{\xi }{e}^{\lambda x},$$

which satisfies equation (7).

By eliminating the nonzero factor e^λx, we can rearrange equation (9) as:

Eq. 10 $$(\mathbf{M}-\lambda \mathbf{I})\mathbf{\xi }=\mathbf{0},$$

where I is the M×M identity matrix and 0 is the null vector.

Equation (10) admits non-trivial solutions if, and only if, det(M-λI)=0, or:

Eq. 11 $$\det [\tilde{\Sigma }\overline{\mathbf{\mu }}{\mathbf{\omega }}^{\mathrm{T}}-\mathrm{diag}({\Sigma }_t\overline{\mathbf{\mu }}+\lambda \mathbf{1}\left)\right]=0,$$

where 1 is the M×1 vector of ones.

The evaluation of this determinant leads to a characteristic polynomial, p(λ), of order M, whose roots λ_l (for l=1,…,L) correspond to the eigenvalues of matrix M.

This polynomial can be conveniently constructed by means of the matrix determinant lemma, which states that, if A is an invertible M×M matrix, γ is a scalar, and uv^T is the outer product of two M×1 vectors, then det(A+γuv^T)=(1+γv^TA^-1u)det(A). Thus, by considering A≡-diag(Σ_tμ̅+λ1), γ≡Σ̃, and uv^T≡μ̅ω^T, one can derive that:

Eq. 12 $$p\left(\lambda \right)\equiv \det (\mathbf{M}-\lambda \mathbf{I})=(1-\sum _{m=1}^M\frac{\tilde{\Sigma }{\omega }_m}{{\Sigma }_t+{\mu }_m\lambda })\prod _{m=1}^M(\frac{{\Sigma }_t}{{\mu }_m}+\lambda ).$$

Although the roots λ_l of this polynomial cannot, in general, be expressed in closed form for M>4 (in accordance with the Abel-Ruffini theorem), they can be readily computed using standard numerical root-finding methods, such as the Newton-Raphson or the secant method. Once obtained, each root/eigenvalue can be substituted into equation (10) to determine their associated eigenvectors, ξ^l. In general, the number L of distinct eigenvalues satisfies L≤M, whether they are real or in pairs of complex conjugates.

Consequently, the general solution Ψ(x) is a linear combination of all solutions Ψ^l(x), each associated with an eigenvalue λ_l. Specifically:

Eq. 13 $$\mathbf{\Psi }\left(\mathbf{x}\right)=c_1{\mathbf{\Psi }}^1\left(\mathbf{x}\right)+\cdots +c_L{\mathbf{\Psi }}^L\left(\mathbf{x}\right)=\sum _{l=1}^L{c}_l{\mathbf{\Psi }}^l\left(\mathbf{x}\right),$$

where c_l are real constants determined by the boundary conditions of the problem.

In the case of a real eigenvalue, Ψ^l(x) takes the form:

Eq. 14 $${\mathbf{\Psi }}^l\left(\mathbf{x}\right)={\mathbf{\xi }}^l{e}^{{\lambda }_{l}x}.$$

To generalize the solution Ψ^l(x) for a complex eigenvalue and its associated complex eigenvector, we first consider the equivalences λ_l≡α_l+iβ_l and ξ^l≡A^l+iB^l and apply them to equation (14), which leads to:

Eq. 15 $${\mathbf{\Psi }}^l\left(\mathbf{x}\right)\equiv ({\mathbf{A}}^l+i{\mathbf{B}}^l)e^{({\alpha }_l+i{\beta }_l)x}.$$

Recalling Euler’s relation e^iθ=cosθ+isinθ and considering that the solution must be real and equal to the sum of both its real and imaginary parts (^{BOYCE et al., 2015}), we get the following expression for the complex case:

Eq. 16 $${\mathbf{\Psi }}^l\left(\mathbf{x}\right)=e^{{\alpha }_{l}x}\left[\right({\mathbf{A}}^l+{\mathbf{B}}^l\left)\cos \right({\beta }_{l}x)+({\mathbf{A}}^l-{\mathbf{B}}^l\left)\sin \right({\beta }_{l}x\left)\right].$$

We now turn to the determination of the constants c_l. Let ψ_m^l(x) denote the m-th component of a solution Ψ^l(x), whether it corresponds to a real eigenvalue, as in equation (14), or to a complex eigenvalue, as in equation (16). In this case, we can write:

Eq. 17 $${\phi }_m\left(x\right)=\sum _{l=1}^L{c}_l{\psi }_m^l\left(x\right).$$

Taking into account the boundary conditions established in equations (2) and (3), it follows that:

Eq. 18 $${\phi }_{m,l}=c_1{\psi }_m^1(-\frac{a}{2})+\cdots +c_L{\psi }_m^L(-\frac{a}{2}),\forall {\mu }_m>0;$$

Eq. 19 $${\phi }_{m,r}=c_1{\psi }_m^1\left(\frac{a}{2}\right)+\cdots +c_L{\psi }_m^L\left(\frac{a}{2}\right),\forall {\mu }_m<0.$$

Assuming L=M — a condition satisfied in most practical cases —, equations (18) and (19), written for their respective values of m, give rise to a linear system of order M×M, whose corresponding solution vector contains the constants c₁,…,c_L. Nonetheless, two additional aspects must be addressed before obtaining definitive values for c_l and, consequently, for φ_m(x).

First of all, this formulation was so far restricted to the case of a single homogeneous region, which allows for a simpler notation and highlights the fundamental aspects of the method. Nevertheless, the procedure can and was meant to be generalized to systems composed of N regions with distinct material properties (e.g., Σ_tⁿ, Σ_s0ⁿ, Σ_fⁿ for each region n). In such cases, local solutions φ_mⁿ(x) must be defined for each region, and several parameters acquire the superscript “n”, including the constants c_lⁿ. The main additional challenge arises from the interface conditions, which are not directly known for internal regions of the system. That is, except for φ_m,l¹ (∀ μ_m>0) and φ_m,r^N (∀ μ_m<0), which are the known boundary conditions, the values of φ_m,lⁿ and φ_m,rⁿ are unknown and must be all initially estimated as unity.

However, to ensure continuity of the angular flux across adjacent regions, the outgoing fluxes from the left boundary of region n, φ_m,lⁿ, must equal the incoming angular fluxes from the right boundary of region n-1, φ_m,r^n-1, for all μ_m>0. Similarly, the outgoing fluxes from the right boundary of region n, φ_m,rⁿ, must equal the incoming fluxes from the left boundary of region n+1, φ_m,lⁿ⁺¹, for all μ_m<0. The relationship between these interface conditions plays a central role in coupling adjacent regions in both the numerical-analytical formulation and purely numerical methods, and is schematically illustrated in Figure 1.

Figure 1
Coupling between distinct regions.

Thus, after an initial calculation of λ_lⁿ, ξ_m^n,l and c_lⁿ on a specific region n, the outgoing fluxes from this region can be expressed as:

Eq. 20 $${\phi }_{m,r}^n=c_1^n{\psi }_m^{n,1}\left(\frac{a^n}{2}\right)+\cdots +c_L^n{\psi }_m^{n,L}\left(\frac{a^n}{2}\right),\forall {\mu }_m>0;$$

Eq. 21 $${\phi }_{m,l}^n=c_1^n{\psi }_m^{n,1}(-\frac{a^n}{2})+\cdots +c_L^n{\psi }_m^{n,L}(-\frac{a^n}{2}),\forall {\mu }_m<0.$$

And, by the continuity of the solutions:

Eq. 22 $${\phi }_{m,l}^{n+1}={\phi }_{m,r}^n,\forall {\mu }_m>0,\text{for}n=1,\dots ,N-1;$$

Eq. 23 $${\phi }_{m,r}^{n-1}={\phi }_{m,l}^n,\forall {\mu }_m<0,\text{for}n=2,\dots ,N.$$

However, after sweeping across all regions, the effective multiplication factor k_eff remains unknown. In pratice, the first evaluation of the angular fluxes φ_mⁿ(x) requires an initial estimate k_eff⁽⁰⁾, which is also commonly set to unity. After each global sweep, which constitutes a single iteration, k_eff is then updated using the Power Iteration Method (^{ALVIM, 2007}), as expressed below:

Eq. 24 $$k_{\mathrm{eff}}^{(k+1)}=k_{\mathrm{eff}}^{\left(k\right)}\frac{\sum _{n=1}^N{a}^n\nu {\Sigma }_f^n{\overline{\Phi }}^{n(k+1)}}{\sum _{n=1}^N{a}^n\nu {\Sigma }_f^n{\overline{\Phi }}^{n\left(k\right)}}.$$

Note that this update requires the values of Φ̅ⁿ, the average scalar flux in region n. The initial estimate Φ̅ⁿ⁽⁰⁾ is also set to unity for all n, but at the end of each iteration, Φ̅ⁿ can be obtained from the following summation:

Eq. 25 $${\overline{\Phi }}^n=\frac{1}{2}\sum _{m=1}^M{\omega }_m{\overline{\phi }}_m^n.$$

Here, φ̅_mⁿ denotes the average angular flux in region n, which can be computed by applying the averaging operator to φ_mⁿ(x) over the domain of region n, [-aⁿ/2, aⁿ/2]:

Eq. 26 $${\overline{\phi }}_m^n\equiv \frac{1}{a^n}{\int }_{\frac{-a^n}{2}}^{\frac{a^n}{2}}{\phi }_m^n\left(x\right)dx.$$

Using the expression for φ_mⁿ(x) given in equation (16), φ̅_mⁿ can be evaluated as a summation of the average solutions, ψ̅_m^n,l:

Eq. 27 $${\overline{\phi }}_m^n=\frac{1}{a^n}{\int }_{\frac{-a^n}{2}}^{\frac{a^n}{2}}\left(\sum _{l=1}^L{c}_l^n{\psi }_m^{n,l}\right(x\left)\right)dx=\sum _{l=1}^L{c}_l^n{\overline{\psi }}_m^{n,l}.$$

For the case of a real eigenvalue, the averaged solution ψ̅_m^n,l takes the form:

Eq. 28 $${\overline{\psi }}_m^{n,l}=\frac{2{\xi }_m^{n,l}}{a^n{\lambda }_l^n}\sinh \left(\frac{{\lambda }_l^n{a}^n}{2}\right),$$

where ξ_m^n,l is the m-th component of the real eigenvector ξ^l,n.

In contrast, for complex eigenvalues and eigenvectors, ψ̅_m^n,l​ is given by:

Eq. 29 $${\overline{\psi }}_m^{n,l}={\Gamma }_m^{n,l}\sinh \left(\frac{{\alpha }_l^n{a}^n}{2}\right)\cos \left(\frac{{\beta }_l^n{a}^n}{2}\right)+{\Delta }_m^{n,l}\cosh \left(\frac{{\alpha }_l^n{a}^n}{2}\right)\sin \left(\frac{{\beta }_l^n{a}^n}{2}\right),$$

where the coefficients Γ_m^n,l e Δ_m^n,l are defined as:

Eq. 30 $${\Gamma }_m^{n,l}=\frac{2}{a^n[{\left({\alpha }_l^n\right)}^2+{\left({\beta }_l^n\right)}^2]}\left(A_m^{n,l}\right({\alpha }_l^n-{\beta }_l^n)+B_m^{n,l}({\alpha }_l^n+{\beta }_l^n\left)\right),$$

Eq. 31 $${\Delta }_m^{n,l}=\frac{2}{a^n[{\left({\alpha }_l^n\right)}^2+{\left({\beta }_l^n\right)}^2]}\left(A_m^{n,l}\right({\beta }_l^n+{\alpha }_l^n)+B_m^{n,l}({\beta }_l^n-{\alpha }_l^n\left)\right).$$

Here, A_m^n,l and B_m^n,l denote the m-th components of the real and imaginary parts, respectively, of the complex eigenvector ξ^l,n≡A^l,n+iB^l,n.

The procedure for computing the angular fluxes φ_mⁿ(x), their averages φ̅_mⁿ, the scalar fluxes Φ̅ⁿ, and updating k_eff is repeated over successive iterations until convergence is satisfactorily achieved. Convergence tests are perfomed at the end of each iteration for both Φ̅ⁿ and k_eff, using prescribed max relative errors ε_Φ̅ and ε_keff, respectively.

For Φ̅ⁿ, the convergence test consists of:

Eq. 31 $$\underset{n=1,\dots ,N}{\max }\left|\frac{{\overline{\Phi }}^{n\left(k\right)}-{\overline{\Phi }}^{n(k-1)}}{{\overline{\Phi }}^{n\left(k\right)}}\right|\le {\epsilon }_{\overline{\Phi }},$$

whereas for k_eff:

Eq. 32 $$\left|\frac{k_{\mathrm{eff}}^{(k+1)}-k_{\mathrm{eff}}^{\left(k\right)}}{k_{\mathrm{eff}}^{(k+1)}}\right|\le {\epsilon }_{k_{\mathrm{eff}}}.$$

When both criteria are satisfied, the iterative process is terminated.

Finally, it should be noted that the calculation of eigenvalues and eigenvectors must be redone after each iteration for every multiplying region (νΣ_fⁿ>0), since the term Σ̃ⁿ, which is part of the matrix M, depends on k_eff and is updated along with it.

2.2 Purely Numerical Solution

In order to validate the analytical solutions, we also present a numerical solution for the neutron transport equation, using the principles of the DD Method (Diamond Difference Method) (^{LEWIS et al., 1984}). First, we apply the same averaging operator used in equation (26) to the entire equation (1) within the domain [-aⁿ/2, aⁿ/2]. However, for the numerical solution, we assume that n corresponds to a subdivision of size aⁿ within a region where nuclear parameters remain homogeneous throughout the domain. We will refer to this subdivision as a “node n”.

So, for the averaging operation of the terms of the total cross-section, isotropic scattering, and fission, we can use the definition of φ̅_mⁿ given by equation (26) to obtain:

Eq. 33 $$\frac{1}{a^n}{\int }_{\frac{-a^n}{2}}^{\frac{a^n}{2}}{\Sigma }_t^n{\phi }_m^n\left(x\right)dx={\Sigma }_t^n{\overline{\phi }}_m^n;$$

Eq. 34 $$\frac{1}{a^n}{\int }_{\frac{-a^n}{2}}^{\frac{a^n}{2}}{\tilde{\Sigma }}^n\sum _{m\prime =1}^M{\omega }_{m\prime }{\phi }_{m\prime }\left(x\right)dx={\tilde{\Sigma }}^n\sum _{m\prime =1}^M{\omega }_{m\prime }{\overline{\phi }}_{m\prime }^n.$$

For the escape term, by the Fundamental Theorem of Calculus, we obtain:

Eq. 35 $$\frac{1}{a^n}{\int }_{\frac{-a^n}{2}}^{\frac{a^n}{2}}{\mu }_m\frac{d}{dx}{\phi }_m^n\left(x\right)dx=\frac{{\mu }_m}{a_n}\left[{\phi }_m^n\right(\frac{a^n}{2})-{\phi }_m^n(-\frac{a^n}{2}\left)\right].$$

So, combining all these integrated terms and remembering, from equations (2) and (3), that φ_m,lⁿ≡φ_mⁿ(-aⁿ/2) e φ_m,rⁿ≡φ_mⁿ(+aⁿ/2), respectively, we are left with:

Eq. 36 $$\frac{{\mu }_m}{a^n}({\phi }_{m,r}^n-{\phi }_{m,l}^n)+{\Sigma }_t^n{\overline{\phi }}_m^n={\tilde{\Sigma }}^n\sum _{m\prime =1}^M{\omega }_{m\prime }{\overline{\phi }}_{m\prime }^n.$$

So, if we define the term Q_mⁿ as:

Eq. 37 $$Q_m^n\equiv {\tilde{\Sigma }}^n\left(\sum _{m\prime =1}^M{\omega }_{m\prime }{\overline{\phi }}_{m\prime }^n\right)-{\Sigma }_t^n{\overline{\phi }}_m^n,$$

we can rewrite equation (36) in two ways:

Eq. 38 $${\phi }_{m,r}^n={\phi }_{m,l}^n+\frac{a^n}{{\mu }_m}{Q}_m^n;\forall {\mu }_m>0;$$

Eq. 39 $${\phi }_{m,l}^n={\phi }_{m,r}^n-\frac{a^n}{{\mu }_m}{Q}_m^n;\forall {\mu }_m<0.$$

It is worth noting that, in this numerical solution, just like in the DD Method, the average angular fluxes φ̅_mⁿ, included in Q_mⁿ, can be approximated using the angular fluxes on each face of node n through the following equation:

Eq. 40 $${\overline{\phi }}_m^n=\frac{{\phi }_{m,r}^n+{\phi }_{m,l}^n}{2}.$$

With initialized values of φ_m,rⁿ and φ_m,lⁿ, we define two sweep directions that allow us to update φ_m,rⁿ and φ_m,lⁿ using equations (38) and (39), respectively. In practice:

• The ‘left → right’ sweep computes the angular fluxes exiting the right face of node n using equation (38).

• The ‘right → left’ sweep computes the angular fluxes exiting the left face of node n using equation (39).

By applying the same principle of continuity used in the analytical solution, the angular fluxes entering from the left face of node n+1 (in the ‘left → right’ sweep) and those entering from the right face of node n-1 (in the ‘right → left’ sweep) are updated with equations (22) and (23), respectively.

Once these incoming fluxes are determined, we move on to the subsequent nodes (depending on the sweep direction), update φ̅_mⁿ and Q_mⁿ, and continue calculating the exiting fluxes to be used as incoming fluxes in adjacent nodes. After completing both sweeps, we evaluate Φ̅ⁿ using equation (25), update k_eff with equation (24), and perform convergence tests on Φ̅ⁿ and on k_eff using equations (31) and (32), respectively.

3. RESULTS AND DISCUSSIONS

For the purpose of validating the solution proposed in this article, two computational algorithms were developed in Python—the “Analytical” (which refers to the Numerical-Analytical solution presented in section 2.1) and the “Numerical” (which refers to the Purely Numerical solution presented in section 2.2)—using the Visual Studio Code environment. These algorithms were applied to simulate the three cases defined in problem model 2 of ^{Lahdour et al. (2019)}. The material configurations for each of the three cases are presented in the figure below.

Figure 2
Core configurations adapted from ^{Lahdour et al. (2019)}.

The nuclear parameters of each region and their properties are shown in Table 1.

Table 1
Nuclear Parameters and Properties of Material Regions.

Label	Material	Σt [cm-1]	Σs0 [cm-1]	νΣf [cm-1]	a* [cm]
F	Fuel	0.4150	0.3340	0.1780	2.40963855
R	Reflector	0.3710	0.3340	0.0000	2.69541779
A	Absorber	0.3710	0.0370	0.0000	2.69541779

For all cases, vacuum boundary conditions were applied. That is, φ_m,l¹ (∀ μ_m>0) and φ_m,r^N (∀ μ_m<0) are all equal to 0. The maximum number of iterations chosen, K_max, was 10⁵. The convergence criteria for the scalar flux and k_eff, ε_Φ̅ and ε_keff, respectively, were both established as 10^-9. The quadrature order used in the calculations was S₃₀₀, as in ^{Lahdour et al. (2019)}.

Theoretically, the analytical solution would require no more than one node per region, but due to overflow issues in the exponential terms during the calculations, the results reported in Tables 2 and 3, specifically for Case 2, were obtained using two nodes per region. For the numerical solution, a total of 500 divisions per region were employed, also in accordance to ^{Lahdour et al. (2019)}.

Table 2
Results obtained for k_eff for the 3 cases

Case	Method
	Analytical	Numerical	CP*	SN*	MOC*	GFM*
1	1.17361	1.17361	1.17360	1.17361	1.17361	1.17361
2	0.94267	0.94267	0.94267	0.94268	0.94268	0.94268
3	1.02265	1.02265	1.02264	1.02265	1.02265	1.02265

Table 3
Comparison of computational times and the number of outer iterations.

Method	Case 1		Case 2		Case 3
	Time (s)	Iterations	Time (s)	Iterations	Time (s)	Iterations
Analytical	0.031838	62	0.053843	58	0.029326	63
Numerical	5.381677	12167	6.804430	15348	7.126371	15533

The results obtained in the calculation of k_eff using the analytical and numerical solutions developed in this work were compared with those obtained by ^{Lahdour et al. (2019)} through the CP (Collision Probability), S_N (DD with S_N Method), MOC (Method of Characteristics) and GFM (Green Function Method) approaches, for the three cases under consideration. These comparisons are summarized in Table 2:

Table 3 presents a comparison between the computational times and the number of outer iterations required by the numerical-analytical and purely-numerical solutions developed in this work. For this comparison, the quadrature order used in the calculations was S₁₆, since extremely high quadrature orders such as S₃₀₀ are not representative of typical reactor calculations and lead to unnecessarily large computational costs.

Figures 3, 4, and 5 present the neutron flux distributions corresponding to Cases 1, 2, and 3, respectively. In the analytical solution, five interpolation points per region were employed to improve the clarity of the plots.

Figure 3
Neutron scalar flux for case 1

Figure 4
Neutron scalar flux for case 2.

Figure 5
Neutron scalar flux for case 3.

Tables 2 and 3, as well as Figures 3, 4, and 5, demonstrate satisfactory results regarding both in the validation against reference calculations and the significant relative improvement between the numerical-analytical and the purely numerical solutions.

4. CONCLUSION

In this work, a numerical-analytical solution was developed for the one-dimensional, monoenergetic neutron transport equation with isotropic scattering, employing the discrete ordinates (S_N​) method for eigenvalue problems, together with a purely numerical approach based on the Diamond Difference (DD) method.

The numerical-analytical and the purely numerical solutions proposed in this study were compared with the reference results reported by ^{Lahdour et al. (2019)}. The comparisons showed that the obtained results are highly accurate. Since the accuracy of the numerical-analytical solution is independent of the size of the homogeneous region, it was possible to employ much coarser meshes, thereby reducing computational time and the number of iterations required, as evidenced in Table 3, when compared to the purely numerical method presented herein.

This reduction in computational effort is particularly relevant for engineering design applications, where large numbers of transport calculations must be performed repeatedly. In such contexts, improvements in computational efficiency may lead to substantial time savings in practical reactor analysis and design workflows.

The perspective for future work is to extend the use of numerical-analytical solutions for the neutron transport equation to the multigroup energy formulation, as well as to two-dimensional and/or three-dimensional cases, especially to geometries of practical relevance, such as cylindrical configurations.

Funding

The authors declare that this work did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

DATA AVAILABILITY STATEMENT

Interview transcripts and qualitative data are not publicly available due to confidentiality commitments made to research participants. Aggregated data supporting the findings are available from the corresponding author upon reasonable request and subject to participant privacy protections.

REFERENCES

• ¹ ALVES FILHO, H. Um método espectro-nodal para problemas da autovalor na teoria de transporte de nêutrons segundo a formulação de ordenadas discretas. Tese (Doutorado) — COPPE/UFRJ, 1999.
• ² ALVIM, A. C. M. Métodos numéricos em engenharia nuclear. Rio de Janeiro: Ed. Certa, 2007.
• ³ BARROS, R. C. A spectral nodal method for the solution of discrete ordinates problems in one and two dimensional cartesian geometry. Tese (Doutorado) — The University of Michigan, 1990.
• ⁴ BARROS, R. C.; LARSEN, E. W. A numerical method for one-group slab-geometry discrete ordinates problems with no spatial truncation error. Nuclear Science and Engineering, Taylor and Francis, v. 104, n. 3, p. 199–208, 1990.
• ⁵ BOYCE, W. E.; DIPRIMA, R. C. Equações diferenciais elementares e problemas de valores de contorno. 10. ed. Rio de Janeiro: LTC, 2015.
• ⁶ DE ABREU, M. P. Métodos determinísticos livres de aproximações espaciais para a solução numérica dominante de problemas de autovalor multiplicativo na formulação de ordenadas discretas da teoria do transporte de nêutrons. Tese (Doutorado) — COPPE/UFRJ, 1996.
• ⁷ DUDERSTADT, J. J.; MARTIN, W. R. Nuclear reactor analysis. New York: John Wiley and Sons, 1976.
• ⁸ LAHDOUR, M.; EL BARDOUNI, T.; CHAKIR, E.; BENAALILOU, K.; MOHAMMED, M.; BOUGUENIZ, H.; EL YAAJOUBI, H. NTP-ERSN: A new package for solving the multigroup neutron transport equation in a slab geometry. Applied Radiation and Isotopes, v. 145, 2019, p. 73–84. DOI: https://doi.org/10.1016/j.apradiso.2018.12.004
  » https://doi.org/10.1016/j.apradiso.2018.12.004
• ⁹ LARSEN, E. W. Spectral analysis of numerical methods for discrete-ordinates problems. Transport Theory and Statistical Physics, Taylor & Francis, v. 15, n. 1-2, p. 93–116, 1986. DOI: https://doi.org/10.1080/00411458608210446
  » https://doi.org/10.1080/00411458608210446
• ¹⁰ LEWIS, E. E.; MILLER, W. F. Computational methods of neutron transport. New York: John Wiley and Sons, 1984.
• ¹¹ RAMÍREZ, S. A. Cálculos de criticalidade usando a equação de transporte de nêutrons multigrupo unidimensional na formulação das ordenadas discretas a partir da solução analítica local. Tese (Doutorado em Modelagem Computacional) — Instituto Politécnico, Universidade do Estado do Rio de Janeiro, Nova Friburgo, 2021.

*

[e-mail]

Conflict of Interest

All authors declare that they have no conflicts of interest.

SCIENTIFIC EDITOR:

Prof. Dr. Bernardo Maranhão Dantas https://orcid.org/0000-0002-2388-6073

SCIENTIFIC EDITOR:

Prof. Dr. Alfredo Lopes Ferreira Filho https://orcid.org/0000-0002-0806-1284

Publication Dates

• Publication in this collection
  30 Oct 2026
• Date of issue
  Jul-Sep 2026

History

• Received
  15 Oct 2025
• Accepted
  16 Mar 2026