Consistent Algorithms Marching Along Characteristics for the Numerical Solution of the Boltzmann Equation

Abstract

We introduce algorithms marching over a polygonal mesh with elements consistent with the propagation directions of the particle (radiation) flux. The decision for adopting this kind of mesh to solve the one-speed Boltzmann transport equation is due to characteristics of the domain of the transport operator which controls derivatives only in the direction of propagation of the particles (radiation) flux in the absorbing and scattering media. This a priori adaptivity has the advantages that it formulates a consistent scheme which makes appropriate the application of the Lax equivalence theorem framework to the problem. In this work, we present the main functional spaces involved in the formalism and a description of the algorithms for the mesh generation and the transport equation solution. Some numerical examples related to the solution of a transmission problem in a high-contrast model with absorption and scattering are presented. Also, a comparison with benchmarks problems for source and reactor criticality simulations shows the compatibility between calculations with the algorithms proposed here and theoretical results.

1. Introduction

A line propagation of radiation and particles is observed for low-density noncharged particles as neutron in nuclear reactors and X-rays or 𝛾-rays in nuclear medicine and engineering nondestructive test applications. For charged particles, this collimation is found in accelerator and spallation target technologies for ADS nuclear plants. This natural behavior of collimated radiation and particles introduces in the numerical computation a new paradigm of operator, the one-speed Boltzmann transport operator, which has been extensively studied in the nuclear community Vladimirov [1], Duderstadt and Martin [2], Kaper et al. [3], Lewis and Miller Jr. [4], Dautray and Lions [5], Mokhtar-Kharroubi [6] and astrophysical community Chandrasekhar [7]. Outside the nuclear community, the analysis of radiation transfer and particle transport problems has several relevant applications in different areas such as population dynamics, heat transfer, remote sensing, global warming models, natural water radiative properties estimation and transmission, scattering and optical tomography.

The huge number of degree of freedom associated with the seven-dimensional phase-space-time domain introduces the necessity of numerically preprocessing the equation in order to produce a solution computationally feasibly. This may introduce conceptual errors in the proposed approximated models. In this work, we are trying to avoid one of the most common mistakes in these models, that is, the disrespect with the natural characteristic of the solution, and also, develop a procedure for discriminating radiation flux of different orders of magnitudes. In image processing, the lower-order magnitude radiation flux is treated as noise, and if we have collimated radiation flux, angular averaged numerical methods will introduce errors. For planes and axis symmetric problems involving the linear Boltzmann transport equation, the adoption of a polygonal mesh with elements consistent with the directions of particles and radiation flux gives an a-priori adaptability that is necessary in the study of the propagation of absorbing and scattering media. This kind of mesh has been called, by us, a natural pixel partition of the domain of definition and has been used successfully in problems of image reconstruction from projection Reis and Roberty [8] and Montero et al. [9]. The use for the solution of the Boltzmann equation is now introduced by me and my students and collaborators Chaffin and Roberty [10], Montero [11], Montero et al. [12]. In the present work, we introduce a new formalism for the radiative transfer and linear particle transport problem. The procedure has been studied in two dimensional models, but extension for 3-dimensional problems is straightforward. Recent works explore the respect with the characteristic of the linear Boltzmann equation, but, contrary to the presented work, they are based on the ray tracing technique and not in the domain partition, as this work is. We cite Liu et al. [13], Chetaine et al. [14], Jevremovic et al. [15].

This paper will be structured with the following sequence. The physical formulation for the stationary model is introduced in Section 2. Some mathematical results that are important to the understanding of the solution given to this numerically implemented problem are presented in Section 3. The consequences of the necessity of a special basis that rotates to accomplish the regularities of the angular flux are discussed in Section 4. We introduce the implementation of the rotation pixel and the natural partitions elements in Section 5. The basis generated by this partition is the most relevant and original aspect found in the present work. Questions related to the algorithm implementation are introduced. The two-dimensional discrete equation used for algorithm implementation is shown in Section 7. The algorithms for the transillumination, source, and the eigenvalue problem are detailed in Section 7.1. The numerical results for the three simulated problems are presented and discussed in Section 8. Finally, we conclude and point out the advances introduced by the present work.

2. Physical Formulation for the Stationary Problem

The particle Boltzmann transport equation is a balance of particles (or photons) inside an absorbing, scattering, or multiplicative participating medium such as a reactor, a tomography, or the natural environment. To state the stationary equation, consider an 𝑛-dimensional convex region Ω,𝑛=2,3, delimiting an region with internal radiation sources 𝑆(𝜔,𝑥),(𝜔,𝑥)∈𝑆^𝑛−1×Ω, subjected to externally generated parallel beams of particles. Here, 𝑆^𝑛−1 is the set of directions in 𝑅^𝑛. For the steady-state situation, with one-velocity particles or no spectral dependency in the case of radiation, the following formulation is obtained from the linear Boltzmann equation, Vladimirov [1], Kaper et al. [3], Dautray and Lions [5], which is usual for the mathematical modeling of the interaction of the radiation with the participating medium. Based on these initial assumptions, we consider the following boundary problem: for a given function 𝑔⁻, find the constant velocity radiation intensity 𝜙∈𝐿²(𝑆^𝑛−1×Ω) such that𝜔⋅∇𝜙(𝜔,𝑥)+𝜎_𝑡(𝑥)𝜙(𝜔,𝑥)=_𝑆𝑛−1𝑓𝑥,𝜔^′⋅𝜔𝜙𝜔^′,𝑥𝑑𝜔^′+𝑆(𝜔,𝑥),(𝜔,𝑥)∈𝑆^𝑛−1𝜙()×Ω,𝜔,𝜎=𝑔⁻(𝜔,𝜎)for(𝜔,𝜎)∈Σ⁻, (2.1)(2.2) where Σ^±={(𝜔,𝜎)∈𝑆^𝑁−1×Γ;±𝜈(𝜎)⋅𝜔>0} are, respectively, the phase space surfaces of the incident and emergent radiation on the physical surface Γ=𝜕Ω, which is the smooth boundary of the bounded, possibly convex region Ω⊂ℝ^𝑛.

The coefficients 𝜎_𝑡 and 𝑓 in (2.1) are, respectively, the total extinction coefficient (absorption + outscattering) and the scattering and fission coefficient. We assume that 𝜎_𝑡∈𝐿^∞(Ω) and, since 𝜔^′⋅𝜔=cos𝜃₀, with 𝜃₀ being the angle between the direction of incident radiation 𝜔^′ and the emergent scattered radiation 𝜔, we also assume that𝑓(⋅,𝑠)∈𝐿^∞(Ω),a.e.𝑠∈[−1,1],𝑓(𝑥,⋅)∈𝐿¹([−1,1]),a.e.𝑥∈Ω. (2.3) Under these assumptions, 𝐿₂(Σ^±;|𝜈(𝜎)⋅𝜔|𝑑𝜎𝑑𝜔) become the appropriated spaces for traces of the functions in the boundary value problem (2.1) and (2.2). For a more detailed discussion, see Cipolatti [16] and the references therein.

3. Mathematical Preliminaries

The stationary transport (2.1) is composed of the following operators. The operator𝐴∶𝑊𝑆^𝑛−1×Ω⟶𝐿²𝑆^𝑛−1𝜙()𝜙)(),×Ω𝜔,𝑥⟼𝐴(𝜔,𝑥=𝜔⋅∇𝜙𝜔,𝑥 (3.1) where derivatives are in the sense of distributions. The space 𝐻=𝐿²(𝑆^𝑛−1×Ω) is the closure of 𝐶(𝑆^𝑛−1×Ω) with respect to the norm 𝜙_𝐻=_{𝑆𝑛−1Ω}𝜙(𝜔,𝑥)²𝑑𝑥𝑑𝜔^1/2, (3.2) where 𝑑𝜔 denotes the measure on 𝑆^𝑛−1 associated with the Lebesgue measure ℝ^𝑛. The space 𝑊 is defined from 𝐻 as 𝑊=𝜙∈𝐻;𝐴𝜙∈𝐻. (3.3)  𝑊 is a Hilbert space for the norm 𝜙_𝑊=_{𝑆𝑛−1Ω}𝜙(𝑡,𝜔,𝑥)²+𝐴𝜙(𝑡,𝜔,𝑥)²𝑑𝜔𝑑𝑥^1/2. (3.4)

The second operator 𝜙()⟼)𝐵∶𝐻⟶𝐻𝜔,𝑥𝐵𝜙(𝜔,𝑥=𝜎_𝑡(𝑥)𝜙(𝜔,𝑥) (3.5) is a bounded operator which has continuous inverse for 𝜎_𝑡(𝑥)>0,𝑥∈Ω. The total macroscopic cross-section is, in applications, piecewise smooth (piecewise and infinitely differentiable), bounded, and strictly positive on Ω with 0<inf𝜎_𝑡(𝑥);𝑥∈Ω≤𝜎_𝑡(𝑥)≤𝑐₀=sup𝜎_𝑡(𝑥);𝑥∈Ω. (3.6) The streaming operator is the sum𝐿_𝜎=𝐴+𝐵∶𝑊𝑆^𝑛−1×Ω⟶𝐿²𝑆^𝑛−1𝜙()×Ω𝜔,𝑥⟼𝐿_𝜎𝜙(𝜔,𝑥)=𝜔⋅∇𝜙(𝜔,𝑥)+𝜎_𝑡(𝑥)𝜙(𝜔,𝑥). (3.7)

The third operator𝜙()⟼)𝐶∶𝐻⟶𝐻𝜔,𝑥𝐶𝜙(𝜔,𝑥=−_𝑆𝑛−1𝑓𝑥,𝜔^′⋅𝜔𝜙𝜔^′,𝑥𝑑𝜔^′ (3.8) describes the anisotropic scattering and fission process. Kernel 𝑓 can be expanded in an absolutely and uniformly convergent series 𝑓𝑥,𝜔⋅𝜔^′=^∞_𝑘=12𝑘−1𝜎4𝜋^𝑘(𝑥)𝑃_𝑘−1𝜔⋅𝜔^′+1𝜆𝜎_𝑓(𝑥), (3.9) where 𝑃_𝑘−1 is the surface harmonic of degree 𝑘−1 (Legendre polynomials). The coefficients 𝜎_𝑘, 𝑘≥1 are nonnegative, and the first coefficient verifies the inequality 0≤inf𝜎_𝑡(𝑥)−𝜎₁(𝑥);𝑥∈Ω≤𝑐₀−𝑐₁,where𝑐₁=inf𝜎₁(𝑥);𝑥∈Ω, (3.10) and (1/𝜆)𝜎_𝑓 is the nonnegative fission cross section divided by the criticality parameter. With this, the collision operator takes its usual form𝐶𝜙(𝜔,𝑥)=_𝑆𝑛−1𝜎_𝑠𝜔⋅𝜔^′+1,𝑥𝜆𝜎_𝑓(𝑥)𝜙𝜔^′,𝑥𝑑𝜔^′. (3.11)

The operator 𝐶 with this kernel is compact from 𝐿²(𝑆^𝑛−1)↦𝐿²(𝑆^𝑛−1), but not from 𝐿²(𝑆^𝑛−1×Ω)↦𝐿²(𝑆^𝑛−1×Ω), and this is the main mathematical problem in this engineering problem.

The trace operator𝛾⁻∶𝑊=𝜙∈𝑊∶_Σ𝜙(𝜔,𝜎)²|𝜔⋅𝜈(𝜎)|𝑑𝜎𝑑𝜔<+∞⟶𝐿²(Σ⁻,|𝜔⋅𝜈(𝜎)|𝑑𝜎𝑑𝜔) (3.12) is continuous and surjective, where 𝑊 is a Hilbert space with the norm ‖𝜙‖²_𝑊=‖𝜙‖²_𝑊+∫_Σ|𝜙(𝜔,𝜎)|²|𝜔⋅𝜈(𝜎)|𝑑𝜎𝑑𝜔.

Given the incident flux by the boundary condition (2.2),  𝑔⁻∈𝐿²(Σ⁻,|𝜔⋅𝜈(𝜎)|𝑑𝜎𝑑𝜔), the following prolongation of 𝑔⁻ inside the phase space is defined: 𝐽_𝜎𝑔⁻(𝜔,𝑥)=exp−^𝜏−₀^(𝜔,𝑥)𝜎_𝑡(𝑥−𝑠𝜔)𝑑𝑠𝑔⁻(𝜔,𝑥−𝜔𝜏⁻(𝜔,𝑥)), (3.13) where 𝜏^±(𝜔,𝑥)=sup{𝑡≥0∶𝑥±𝑡𝜔∈Ω}are the distance to the emergent and incident boundary, respectively. Note that this definition implies that 𝐽_𝜎 is the solution operator of the problem 𝐿_𝜎[𝜙]=0,𝜙|_Σ⁻=𝑔⁻.

These operators can be used to formulate the first-order boundary non-homogeneous stationary direct problem given (2.1) and (2.2) as the following equivalent homogeneous problem: to find 𝜑=𝜙−𝐽[𝑔⁻]∈𝑊∩dom(𝐴) such that𝐿_𝜎−𝐶𝜑=𝑆−𝐶𝐽_𝜎𝑔⁻=ℎ∈𝐻, (3.14) where dom(𝐴)={𝜑∈𝑊∶𝜑=0onΣ⁻}.

Theorem 3.1. For each ℎ∈𝐻, the boundary value problem (3.14) has a unique solution witch is given by 𝜙=𝐼−𝐿_𝜎⁻¹𝐶⁻¹𝐿_𝜎⁻¹[ℎ], (3.15) and verifies the following estimate: 𝜙≤1−exp−𝑐₀𝑑(𝐼−𝐿_𝜎)⁻¹𝑆⁻¹‖ℎ‖, (3.16) where 𝑑 is the diameter of Ω.

Proof. See Chapter (III) and Lemma 3.1 of Vladimirov [1].

In Sections IV thought VII of Vladimirov [1], operators 𝐿_𝜎∶𝐷(𝐿_𝜎)↦𝐿²(𝑆^𝑛−1×Ω) and 𝐶∶𝐿²(𝑆^𝑛−1×Ω)↦𝐿²(𝑆^𝑛−1×Ω) are investigated for more general Sobolev spaces. We can easily adapt them to our Hilbert space framework.

Lemma 3.2. The operator 𝐿_𝜎⁻¹∘𝐶∶𝐿²(𝑆^𝑛−1×Ω)↦𝐿²(𝑆^𝑛−1×Ω) is compact.

Proof. See Lemma 5.3 of Vladimirov [1].

This result is fundamental for the reformulation of boundary value problem given by (2.1) and (2.2) in terms of an integral equation based on a compact perturbation of the identity operator used in this work: to find 𝜙(𝜔,𝑥) such that𝜙=𝐿_𝜎⁻¹∘𝐶𝜙+𝐽_𝜎𝑔⁻,(𝜔,𝑥)∈𝑆^𝑛−1×Ω. (3.17)

Remark 3.3. The study of the numerical behavior of compact operators is an very well known subject. The spectral properties of compact operators are very similar to that of finite-dimensional operators, fact that make analysis very easy. For this kind of operator, the subspaces associated with a eigenvalue are finite dimensional. It is based on this behavior our interest in model this problem using integral equation (3.17). A situation similar also occurs in models based on the Laplacian operator paradigm, which fortunately, for bounded domains, is an operator with compact resolvent. Since the one-velocity transport equation has a more complex structure than the Laplacian, it is necessary to do the composition of the inverse of the streaming operator, 𝐿_𝜎⁻¹, and the collision operator, 𝐶, to obtain a compact operator. Note that 𝐿_𝜎⁻¹ presupposes and line integrals through the ray paths of particles inside the domain, and is this fact that will have a remarkable influence on the numerical behavior of the problem. We will see that this connects compactness and characteristic basis.

4. Relation between the Integral Operator and the Characteristic Basis

Due to its very nature, the transport equation for the angular flux may be viewed as a system of ordinary differential equations coupled by the collision operator. Each angular field is a different field which has its own direction of propagation as its characteristic directions. This means what characteristic means: the discontinuities of each angular flux field will separately propagate in its own direction. This observation introduces a natural way to solve the transport equation, Askew [17], that will be computationally more expensive but will produce more accurate calculations. For this, consider the mapping in the angular phase space 𝑆^𝑛−1×𝑅^𝑛𝑄(𝜔,𝑥)=𝜔^′,𝑥_𝜔,𝑥_𝜔⟂ (4.1) with𝜔^′=𝜔,𝑥_𝜔=𝜔⋅𝑥,𝑥_𝜔⟂=𝑥−(𝜔⋅𝑥)𝜔. (4.2) This is an isometry that rotates the axis 𝑥 to a new direction 𝑥_𝜔 coincident with the propagation direction. For each fixed direction 𝜔∈𝑆^𝑛−1, it introduces a new coordinate system (𝑥_𝜔,𝑥_𝜔⟂) with the second coordinate in the projection Π_𝜔[Ω] of the physical domain Ω in a subspace of 𝑅^𝑛 made with directions 𝜔^⟂ and evidences the ordinary character of each angular flux in the equation𝜕𝜙𝜕𝑥_𝜔+𝜎_𝑡𝜙=𝐿_𝜎𝜙=𝐶𝜙+𝑆,𝜔,𝑥_𝜔,𝑥_𝜔⟂∈𝑄𝑆^𝑛−1×𝑄[Ω], (4.3) subject to the same incident boundary condition𝜙𝜔,𝑥⁻_𝜔,𝑥⁻_𝜔⟂=𝑔⁻𝜔,𝑥⁻_𝜔,𝑥⁻_𝜔⟂,(𝜔,𝑥)∈Σ⁻. (4.4) Note that this incident boundary condition can be interpreted as an initial value for this system of ordinary differential equations coupled through the collision operator. For simplicity, we have used the same notation 𝜙(𝜔,𝑥_𝜔,𝑥_𝜔⟂) for the rotated angular flux 𝑄[𝜙](𝜔,𝑥_𝜔,𝑥_𝜔⟂).

As pointed out in Section 3., the inversion of the operator 𝐿_𝜎 and its composition with special class of collision operator found in nuclear applications gives a compact operator𝐿_𝜎⁻¹𝐶=𝐿_𝜎⁻¹𝐶^scattering+1𝜆𝐿_𝜎⁻¹𝐶^fission∶𝐿²𝑆^𝑛−1×Ω⟶𝐿²𝑆^𝑛−1×Ω. (4.5) Of course this important property will pose the direct problem involving recovery of the angular flux from appropriate parameter data 𝜎, 𝜎_𝑠, and 𝜎_𝑓. Some more functionals defined in terms of the rotated coordinates systems need to be introduced in order to make comprehensive the approximations that will be done to solve the problem. First of all, the optical distance from interior point 𝑥 to a point 𝑥₀∈𝑅^𝑛𝛼𝑥,𝑥⁰=_𝑥^𝑥⋅𝜔_0𝜔𝜎_𝑡𝑋_𝜔,𝑋_𝜔⟂𝑑𝑋_𝜔, (4.6) where 𝜔=(𝑥−𝑥⁰)/‖𝑥−𝑥⁰‖ is the direction that aligns 𝑥 and 𝑥⁰, and the extension for points in the possibly bare medium outside the physical domain Ω may be done by assuming 𝜎(𝑥)=0,𝑥∈𝑅^𝑛⧵Ω. When 𝑥₀ is a boundary point, it is given by 𝑥^±(𝜔,𝑥)=𝑥±𝜏_±(𝜔,𝑥)𝜔 which is the point in the boundary where the line passing 𝑥 with direction 𝜔 crosses the emergent or the incident boundary, respectively. It may be used to define the attenuation from the incident boundary𝑎⁻𝜔,𝑥_𝜔,𝑥_𝜔⟂=exp−^𝑥𝜔_𝑥−𝜔𝜎_𝑡𝑋_𝜔,𝑋_𝜔⟂𝑑𝑋_𝜔. (4.7) We also define the right inverse of the trace operator, the extension operator as𝐸𝑔⁻(𝜔,𝑥)=𝑔⁻𝜔,𝑥−𝜔𝜏₋(𝜔,𝑥);(𝜔,𝑥)∈𝑆^𝑛−1×Ω, (4.8) and the extension with attenuation from the incident boundary𝐸𝑔⁻(𝜔,𝑥)𝑎⁻(𝜔,𝑥)=𝑔⁻𝜔,𝑥−𝜔𝜏₋(𝜔,𝑥)𝑎⁻(𝜔,𝑥)=𝐽𝑔⁻(𝜔,𝑥);(𝜔,𝑥)∈𝑆^𝑛−1×Ω. (4.9) By noting that the attenuation is the inverse of the integrating factor for the system (4.3), we may integrate from the incident boundary to the interior point 𝑥 to obtain an integral equation which is an explicit expression of (3.17): 𝜙𝜔,𝑥_𝜔,𝑥_𝜔⟂=^𝑥𝜔_𝑥−𝜔𝑎⁻𝜔,𝑥_𝜔,𝑥_𝜔⟂𝑎⁻𝜔,𝑋_𝜔,𝑋_𝜔⟂⋅𝑆𝜔,𝑋_𝜔,𝑋_𝜔⟂𝑑𝑋_𝜔+𝑔⁻𝑥−𝜔𝜏₋(𝜔,𝑥)⋅𝑎⁻𝜔,𝑥_𝜔,𝑥_𝜔⟂⋅^𝑥𝜔_𝑥−𝜔𝑎⁻𝜔,𝑥_𝜔,𝑥_𝜔⟂𝑎⁻𝜔,𝑋_𝜔,𝑋_𝜔⟂⋅_𝑆𝑛−1𝜎_𝑠𝜔⋅𝜔^′,𝑋_𝜔,𝑋_𝜔⟂+1𝜆𝜎_𝑓𝑋_𝜔,𝑋_𝜔⟂𝜙𝜔^′,𝑋_𝜔,𝑋_𝜔⟂𝑑𝜔^′𝑑𝑋_𝜔, (4.10) for all (𝑥_𝜔,𝑥_𝜔⟂)∈[𝑥⁻_𝜔,𝑥⁺_𝜔]×Π_𝜔[Ω],𝜔∈𝑆^𝑛−1.

The integration may also be done between two arbitrary points, inside or outside the domain Ω, 𝑥, and 𝑥⁰. In this case we must extend the cross-sections outside the domain as zero. By defining 𝜔=(𝑥−𝑥⁰)/|𝑥−𝑥⁰| and 𝑎⁰(𝜔,𝑥_𝜔,𝑥_𝜔⟂)=exp(−𝛼(𝑥,𝑥⁰)), we obtain 𝜙𝜔,𝑥_𝜔,𝑥_𝜔⟂=^𝑥𝜔_𝑥0𝜔𝑎⁰𝜔,𝑥_𝜔,𝑥_𝜔⟂𝑎⁰𝜔,𝑋_𝜔,𝑋_𝜔⟂⋅𝑆𝜔,𝑋_𝜔,𝑋_𝜔⟂𝑑𝑋_𝜔+𝑔⁻𝑥−𝜔𝜏₋(𝜔,𝑥)⋅𝑎⁰𝜔,𝑥_𝜔,𝑥_𝜔⟂⋅^𝑥𝜔_𝑥0𝜔𝑎⁰𝜔,𝑥_𝜔,𝑥_𝜔⟂𝑎⁰𝜔,𝑋_𝜔,𝑋_𝜔⟂⋅_𝑆𝑛−1𝜎_𝑠𝜔⋅𝜔^′,𝑋_𝜔,𝑋_𝜔⟂+1𝜆𝜎_𝑓𝑋_𝜔,𝑋_𝜔⟂𝜙𝜔^′,𝑋_𝜔,𝑋_𝜔⟂𝑑𝜔^′𝑑𝑋_𝜔, (4.11) for all (𝑥_𝜔,𝑥_𝜔⟂)∈[𝑥⁰_𝜔,𝑥⁺_𝜔]×Π_𝜔[Ω],𝜔∈𝑆^𝑛−1.

Remark 4.1. Problems given by the equations (1)(2.1) and (2.2), (2)(3.17), (3)(4.3) and (4.4), and(4)(4.10) and (4.11) are alternative formulations to the one-velocity transport equation model and express different aspects of this equation. Note that the initial value character of the problem will be accomplished if we construct the solution by using a method that marches along each one of the characteristics.

Remark 4.2. The consistence that will be appropriated to the application of the Lax equivalence theorem framework that has been formulated for time transient problems will be accomplished by the adoption of a rotated pixel basis in the next section. The compact operator 𝐿_𝜎⁻¹𝐶 is subcritical or critical, which means the inverse of its eigenvalue is not less than one. So this gives stability to the class of problems been investigated.

As a consequence of consistency, the Lax equivalence theorem applies, Valtchev and Roberty [18], and may be stated as follows.

Theorem 4.3. For a well-posed initial value problem, a consistent numerical scheme is convergent if and only if it is stable.

5. The Natural Partition

The integral form of the transport equation (4.11) has been written in a different coordinate system for each direction of propagation 𝜔. In order to attempt the regularities requirements of each angular flux 𝜙, we must introduce a grid with one direction as 𝜔. This is impossible unless we choose a set of discrete ordinates as representative of all angular directions. This must be done a priori. Let 𝑆 be the surface of the unity sphere in 𝑅³ and 𝜔=𝜔₁,𝜔₂,𝜔₃=sin(𝜃)cos𝜓,sin(𝜃)sin𝜓,cos(𝜃),𝜓∈[0,2𝜋),𝜃∈[0,𝜋), (5.1) where 𝜃 and 𝜓 are the usual angles of the spherical coordinates system. To be more specific, let us consider the two-dimensional equation embedded in 𝑅³. So, we choose a set of 𝑆_𝑗𝑖=𝜔_𝑗𝑖=sin𝜃_𝑖cos𝜓_𝑗,sin𝜃_𝑖sin𝜓_𝑗,cos𝜃_𝑖,𝑗=1∶2𝐽,𝑖=1∶2𝐼 (5.2) of equally spaced angular angles 𝜓_𝑗=((𝑗−1)/2𝐽)2𝜋 and at least symmetrical with the angle 𝜃=𝜋/2 of 2𝐼 collocations directions out of the plane generated by the vectors (cos(𝜓),sin(𝜓)), that is, out of the plane 𝑥𝑦. We now may use these directions to introduces grids, one for each direction, or in the case of a two-dimensional problem, only in the plane directions Π_𝑗={cos(𝜓_𝑗),sin(𝜓_𝑗)}. Note that in the two dimensional problem, since the problem is embedded in the three-dimensional space, angular flux is defined for all directions in the 𝑆² sphere but has spatial variations only in the plane 𝑥𝑦. This introduces, for each direction out of the 𝑥𝑦 plane, the necessity of a sin(𝜃_𝑖) factor correction. This gives the possibility of extending this formulation to investigate axis symmetric problems, but, from now on, we will fix in the two-dimensional problem. We have introduced 2J coordinate systems. Without loss of generality, let us consider each one with a rectangular equally spaced grid. These coordinates systems are supposed to be equally rotated with respect to the others. We have chose an even number for J, so, directions of propagations appear in groups of four, to complain with the four quadrants.

We know that it is necessary to combine the different directions in order to calculate the collision operator for groups of pair the directions. Since the cross sections are properties common for all angular flux, the complete intercessions of the 2J rotated grid is the appropriate partition consistent with the numerical solution of the transport equation. To fix ideas, let us introduce the appropriate notation in 𝑅². The direction 𝜔_𝑗 contains one strip𝜒_𝑛𝑗𝑥,𝑦=1,if𝑛_𝑗≤−𝑥sin𝜓_𝑗+𝑦cos𝜓_𝑗<1+𝑛_𝑗,0,else, (5.3) which may be intercept with a strip in one direction in 𝜔^⟂_𝑗𝜒_𝑚𝑗𝑥,𝑦=1,if𝑚_𝑗≤𝑥cos𝜓_𝑗+𝑦sin𝜓_𝑗<1+𝑚_𝑗,0,else, (5.4) to form a rotated pixel 𝑃_{𝑗𝑚𝑗𝑛𝑗}𝑥,𝑦=1,𝑥,𝑦∈𝜒_𝑛𝑗∩𝜒_𝑚𝑗,0,else. (5.5) The regular spacing ℎ_𝑥𝑗,ℎ_𝑦𝑗=1 has been adopted without loss of generality for simplicity. Note that in the direction 𝜔_𝑗=cos(𝜓_𝑗)𝑒_𝑥+sin(𝜓_𝑗)𝑒_𝑦,the integrated (4.11) must be collocated for points (𝜔_𝑗,𝑥_0𝜔𝑗,𝑥^⟂_𝜔𝑗) and (𝜔_𝑗,𝑥_𝜔𝑗,𝑥^⟂_𝜔𝑗) that will be referenced by the indexes (𝑗,𝑚_𝑗,𝑛_𝑗) and (𝑗,𝑚_𝑗+1,𝑛_𝑗), respectively. This is a consequence of the fact that 𝜔_𝑗 aligns the points and 𝑥^⟂_0𝜔𝑗=𝑥^⟂_𝜔𝑗. As the marching direction is 𝜔_𝑗, the equation index 𝑚_𝑗 will grow from its value corresponding to a point 𝑥_𝜔𝑗 in the incident boundary, where the boundary condition is imposed, to the corresponding point in the emergent boundary. This happens for each transversal index 𝑛_𝑗 in each direction 𝜔_𝑗. We now face the problem of computing the collision integral in (4.11) with a consistent discretization. Note that the measure 𝑑𝜔′𝑑𝑋_𝜔^𝑥𝜔_𝑥0𝜔𝑎⁰𝜔,𝑥_𝜔,𝑥_𝜔⟂𝑎⁰𝜔,𝑋_𝜔,𝑋_𝜔⟂⋅_𝑆𝑛−1𝜎_𝑠𝜔⋅𝜔^′,𝑋_𝜔,𝑋_𝜔⟂+1𝜆𝜎_𝑓𝑋_𝜔,𝑋_𝜔⟂𝑑𝜔^′𝑑𝑋_𝜔 (5.6) will corresponds to a volumetric Lebesgue measure when the two integrals interchanges, that is, when the Fubini theorem is applicable. The natural way to obtain this in the discretized approximation is by considering the partition presented by Roberty [19] that makes a complete intersection of the pixels 𝑃_{𝑗𝑚𝑗𝑛𝑗}. This may be done with the complete interception of the strips defining the pixels to form the natural pixels𝐸_{𝑛𝑗1,…,𝑛𝑗2𝐽}𝑥,𝑦=^2𝐽_𝑗=1𝜒_𝑛𝑗𝑥,𝑦. (5.7) This interception produces a domain partition which is consistent in the sense that if we increase the number of directions and reduce the size of the spacing ℎ_𝑗=1 by using a more refined mesh, the error in the approximation to the exact stationary transport equation will approache zero accordingly. A particle or photon collected in one interception can only come from pixels aligned with the current position thought the discrete directions considered in the model. The approximation made can answer the question “from where the particle came?”, in a way consistent with the exact equation.

Calculations of the vertex, area, centroid, and all geometric information are a classical computational geometry problem of domain decomposition, O’Rourke [20], which may present a great complexity in dimensions greater than two.

As an optimization problem, the mesh generation is the determination of the convex polygonal sets limiting the maximum area that follows linear restrictions of the type: 𝐴𝑥≤𝑏, where 𝑥 is a vector in ℝ², 𝑏 is a vector in some ℝ^𝑁𝑏 related with the strips 𝑛_𝑗, and 𝐴 is a 2×𝑁_𝑏 matrix related with the direction angle 𝜃_𝑗. The search process of the elements in the mesh executes the following steps:(1)decomposition of the domain in search subdomains (2)determination of the minimum and maximum strips for the search subdomain (3)search for the true intersections, (4)calculus of elements vertices, edges, area, and pertinence index for all directions (5)assemble subdomains to form the natural pixel partition.

5.1. Decomposition of the Domain in Search Subdomains

The greater majority of the combinations (𝑛₁,…,𝑛_2𝐽) are false intersections and can be a prior avoided by restricting the search for parameters on regions for which we can expect true intersections. We call these regions search subdomains, and they can be a pixel for some prescribed direction 𝑗, an intersection of some given strip with a sector, or any other type of polygonal subdomain. With them, we avoid a search with a large number of false intersections in the step (3). As an example in Figure 1 we show a detail of the mesh 𝐽64𝑀128 for huge numbers of direction (64) and pixels (128×128), which is typical of image processing and reconstruction problems. Since the rotational symmetry allows us to calculate the interception only in a 𝐽/𝜋 sector, we may have less computation than the case of doing only two interceptions, as illustrated by Figure 2. Finally, in Figure 3, we present the complete mesh used to simulate our solution with test and benchmark problems.

Figure 1 

A detail in the J32M64 mesh.

Figure 2 

Sector mesh generated strip by strip in the 𝑗=16 and 𝑀=15 problem.

Figure 3 

The J6M24 mesh.

5.2. Determination of the Minimum and Maximum Strips for the Search Subdomain

This is a complementary procedure for minimizing the determination of false intersection.

5.3. Search for the True Intersections

In each domain of search, every strip combinations must be tested until the area of the domain is totally filled with natural pixels.

5.4. Calculus of Element Vertexes, Edges, Area, and Pertinence Index for all Directions

This calculation is done simultaneously with the search for true intersections. Note that not all intersections of lines limiting the strips are vertexes of the polygons and not all lines across a vertex are a polygonal edge.

5.5. Assemble Subdomains to Form the Natural Pixel Partition

In each subdomain, the pertinence indexes are calculated independently and are assembled to form the matrix of pertinence index for the entire mesh.

For more information, please see Roberty [19] and in the references there. In Figure 4, the dependence of the numbers of polygonal with the numbers of views and slices is presented. Note that the almost exponential growth observed for growing numbers of slice and views introduces memories problems in the computational implementation of the present methodology. In Figure 5, we present a comparison of mesh 𝐽6𝑀3 with the classical Delaunay triangulation made for the same vertex number. Note that it is inconsistent with the natural partition, that is, it is not a triangular decomposition of it.

Figure 4 

Numbers of polygonal dependence with the numbers of views and slices.

(a)

(b)

Figure 5 

A comparison between the natural and the Delaunay mesh.

6. The Discrete Equation for the Two-Dimensional Problem

As we mention before, we may consider the integral form of the transport (4.11) for direction 𝜔_𝑗, 𝑗=1,2𝐽, in the pixel 𝑃_{𝑗𝑚𝑗𝑛𝑗}, 𝑛_𝑗=1,2𝑀 and 𝑚_𝑗=1,2𝑀. Since we are focusing on the two dimensional problem, it necessary introduces the sin(𝜃_𝑖) correction for each out of the 𝑥𝑦 plane direction. Note that the 𝑥𝑦 plane corresponds to the case sin(𝜃_𝑖)=sin(𝜋/2)=1. There is by hypotheses no variation of the angular flux in the 𝑧=𝑥𝑦^⟂ direction, but it still remains dependent of azimuthal angle. Let us rebuild the notation to make this behavior explicit. Since 𝜔₃⋅∇_𝑥𝜙=𝜕𝜙/𝜕𝑧=0 and all other functional variations in this direction are also zero, we rewrite 𝜔_𝑗𝑖=sin(𝜓_𝑖)⋅𝜔_𝑗 and obtain by straightforward calculations, the optical distance for points with projection 𝑥𝑦 along the 𝜔₃ direction as𝛼𝑥𝑦,𝑥𝑦₀,𝜓_𝑖=∫^{𝑥𝑦⋅𝜔𝑗}_{𝑥𝑦0⋅𝜔𝑗}𝜎𝑥𝑦^′_𝜔,𝑥𝑦_𝜔⟂𝑥𝑦,𝑥𝑦₀𝑑𝑥𝑦^′_𝜔𝑗sin𝜓_𝑖. (6.1) Since all other parameters, as the cross sections, and sources are linear measures of the change in the particle population in the direction 𝜔_𝑗𝑖, it must be corrected accordingly.

As a consequence of the divergence theorem, the angular flux in direction 𝜔_𝑗𝑖 has average value on all lateral segments (surfaces in three dimension), that is, the segments that are transversal to the strip 𝜒_𝑚𝑗(𝑥,𝑦), parallel to 𝜔_𝑗 for all 𝑚_𝑗, are equal to zero. We complement this by introducing the following hypotheses: (1)that angular flux is well represented by its average value in each transversal section to the strip 𝜒_𝑛𝑗(𝑥,𝑦), (2)the cross sections are known by its average values in the natural pixel 𝜒_{𝑒𝑛1,…,𝑛2𝐽}=^2𝐽_𝑗=1𝜒_𝑛𝑗, (6.2)(3)the average cross sections inside each rotated rectangular pixel 𝑃_{𝑗𝑚𝑗𝑛𝑗} are the area (volume in three dimensions) average weight of each natural pixel cross section value, (4)to accomplish with high cross section values, we integrate the exponentials inside the rotated pixel 𝑃_{𝑗𝑚𝑗𝑛𝑗}, (5)the mean value of angular flux inside the pixel 𝑃_{𝑗𝑚𝑗𝑛𝑗} is used in the calculation of scattering and fission reaction rate.

7. The Method of Successive Approximations

With respect to the intensity of the extinction cross section 𝜎_𝑡, we observe two important class of problems: (1)there exists a transillumination of the medium Ω, which means that the extinction cross sections not sufficient to produce an optical path that shields the radiation flux. In this class of problems, a flux radiation in the incident boundary propagates through the medium and can be collected in the emergent boundary. Its modification will be used to interpret the optical properties of the medium inside Ω, (2)there is no transillumination, which means that the extinction cross sections are sufficient to produce an optical path that completely shields the radiation flux. In this case, the shielded radiation is must be replaced by an internal source or the fission in the case of reactors.

Since the extinction coefficients may present high values, we considered exponential flux attenuation inside the rotated pixel and calculated the reaction rates using the average flux expressed in terms of the incident and emergent flux in the rotated pixel. Equation (4.11) becomes 𝜙_{𝑚𝑗+1,𝑛𝑗,𝑗}=𝜙_{𝑚𝑗,𝑛𝑗,𝑗}𝐸_{1,𝑚𝑗,𝑛𝑗,𝑗}+𝑞_{𝑚𝑗,𝑛𝑗,𝑗}𝐸_{2,𝑚𝑗,𝑛𝑗,𝑗}for𝑚_𝑗,𝑛_𝑗=1,2𝑀;𝑗=1,2𝐽. (7.1) Here, the source with contribution of the scattering and which forms fission inside the rotated pixel is 𝑞_{𝑚𝑗,𝑛𝑗,𝑗}=∑^2𝐽_𝑘=1∑_{𝑒∈𝑃𝑗𝑛𝑗𝑗𝑚}𝑎_𝑒⋅𝜎_𝑠,𝑒+1𝜆𝜎_𝑓,𝑒⋅𝐸_{3,𝑚(𝑒,𝑘),𝑛(𝑒,𝑘),𝑘)}𝜙_{𝑚(𝑒,𝑘),𝑛(𝑒,𝑘),𝑘}+𝐸_{4,𝑚(𝑒,𝑘)+1,𝑛(𝑒,𝑘),𝑘}𝜙_{𝑚(𝑒,𝑘)+1,𝑛(𝑒,𝑘),𝑘}𝐸;𝑘=1,2𝐽,_{1,𝑚𝑗,𝑛𝑗,𝑗}=exp−𝜎_{𝑚𝑗,𝑛𝑗,𝑗},𝐸_{2,𝑚𝑗,𝑛𝑗,𝑗}=1−exp−𝜎_{𝑚𝑗,𝑛𝑗,𝑗}𝜎_{𝑚𝑗,𝑛𝑗,𝑗},𝐸_{3,𝑚𝑗,𝑛𝑗,𝑗}=1−1+𝜎_{𝑚𝑗,𝑛𝑗,𝑗}exp−𝜎_{𝑚𝑗,𝑛𝑗,𝑗}𝜎_{𝑚𝑗,𝑛𝑗,𝑗}1−exp−𝜎_{𝑚𝑗,𝑛𝑗,𝑗},𝐸_{4,𝑚𝑗,𝑛𝑗,𝑗,𝑖}=−1+𝜎_{𝑚𝑗,𝑛𝑗,𝑗}+exp−𝜎_{𝑚𝑗,𝑛𝑗,𝑗}𝜎_{𝑚𝑗,𝑛𝑗,𝑗}1−exp−𝜎_{𝑚𝑗,𝑛𝑗,𝑗},𝜎_{𝑚𝑗,𝑛𝑗,𝑗}=_{𝑒∈𝑃𝑗𝑛𝑗𝑗𝑚}𝑎_𝑒𝜎_𝑒. (7.2)

Remark 7.1. Note that the correction for extending this methodology for axis symmetric problems is 𝜎_𝑒,𝑖=𝜎_𝑒sin𝜓_𝑖 (7.3) as mentioned in the introduction of the natural partition in Section 5.

The localization of the natural element 𝑒_{𝑛1,…,𝑛2𝐽} inside the rotated pixel 𝑃_{𝑚𝑗,𝑛𝑗,𝑗} is made by the incident matrix. As we mentioned before, given one pair (𝑒,𝑗), one may calculate the rotate integer coordinates (𝑚(𝑒,𝑗),𝑛(𝑒,𝑗),𝑗) by extracting the integer part of the natural element centroid in each rotate coordinate system. This is a procedure similar to the retro projection that is done in image reconstruction algorithms. The function relating the two integer coordinate systems is a discrete version of transformation (4.2). It is one-to-one between (𝑒,𝑗)→(𝑚,𝑛,𝑗). It is not convenient to make it onto since the implementation of an procedure for the localization of the boundary of Ω will be computationally expensive. We will start the march in the algorithm in the position labeled with 𝑚_𝑗=1, which may correspond to a pixel outside the domain, and march in one of the directions 𝑗 until the pixel labeled 2𝑀, which also may be outside the domain, is arrived. This introduces no problem if we take care in extending outside the domain the cross sections value by zero, the flux by its boundary value and calculates the exponentials accordingly.

As mentioned before, these equations are consistent and introduce no spatial error when the properties are constant; in that case only, the angular discretization error will be presented. The simulation problems that we will solve in Section 8.3 have material parameters very close to this situation.

As we pointed out before, the marching scheme equation (7.1) is a discretization of the first-order form of the first order system representing the stationary transport equation that has initial value in the incident boundary. The fact that the total cross-section is bigger than the scattering assures contractility of the compact operator related to the scattering kernel, and from this, we will obtain the stability and as a consequence a convergent scheme.

7.1. Algorithms Based on Successive Approximation Method

Properties of collision operator equation (3.11) and of the boundary conditions allow us to distinguish the basic class of problems for the stationary transport equation:𝜙=𝐸𝑎₋𝑔⁻+𝐿_𝜎⁻¹𝐶^scattering𝜙+1𝜆𝐿_𝜎⁻¹𝐶^fission𝜙+𝐿_𝜎⁻¹𝑆;(𝜔,𝑥)∈𝑆^𝑛−1×Ω. (7.4)(1)The collision operator and the external source localized inside Ω are negligible. In this case, the equation is integrable, and we have the direct problem for the X-ray transform. (2)The fission part of the collision operator 𝐿_𝜎⁻¹𝐶^fission has the first positive eigenvalue greater than one, and the domain Ω is subcritical external source problem. The incident flux may be zero or not. There are two important cases of nonfission domain, that is, when there is no fissile material inside the domain and 𝜎_𝑓 is properly zero. The external source inside the domain problem, with 𝑞≠0 and 𝑔⁻=0 and transillumination problem with 𝑆=0 and 𝑔⁻≠0. The transillumination problem has a greater interest for tomography and nondestructive testing. (3)The fission part of the collision operator 𝐿_𝜎⁻¹𝐶^fission is defined with support inside Ω, and we want to find his first positive eigenvalue and its positive eigenfunction. Once we find the eigenvalue, we give to 𝜅 the same value and obtain a new fission operator with first eigenvalue equal one, which means that we have found the critical material composition for the domain Ω.

7.2. The External Source Inside the Domain Problem

The algorithm for the angular flux calculation in the external source inside the domain problem follows the following steps. (1)Establish an error tolerance criteria. (2)Input material properties and source: 𝜎_𝑡,𝜎_𝑠,𝑆. (3)Input mesh parameters: 𝑀, 𝐽, 𝐼. (4)Generate mesh: elements area, centroid, and vertex and incidence integers matrices for natural elements inside the rotated pixels. (5)Adjust geometric scale between the regular mesh and the spatial distribution of reacting material. (6)Calculate the attenuation, the exponentials, and external source for each rotated pixel. (7)Initialize the flux and reaction rates. (8)For iteration equal 1 to maximum number of iterations do.(9)Calculate scattering reaction rate based with the previous iteration flux estimate. (10)Calculate the new iteration flux with (7.1) with zero fission part and zero incident flux. (11)Verify if the difference between the old flux and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the iteration before the maximum number of iterations. (12)End the iteration do.

7.3. First Positive Eigenvalue and Eigen Angular Flux Determination Problems

The algorithm for first positive eigenvalue and eigen angular flux in the criticality determination problem follows the following steps.(1)Establish an error tolerance criteria. (2)Input material properties and source: 𝜎_𝑡,𝜎_𝑠,𝜎_𝑓. (3)Input mesh parameters: 𝑀, 𝐽, 𝐼. (4)Generate mesh: elements area, centroid and vertex and incidence integers matrices for natural elements inside the rotated pixels. (5)Adjust geometric scale between the regular mesh and the spatial distribution of reacting material. (6)Calculate the attenuation and the exponentials for each rotated pixel. (7)Set the first eigenvalue to one: 𝜆⁰=1.(8)Initializes with the flux 𝜙⁰ normalized by the fission part of the collision operator. (9)For outer iteration equal 1 to maximum number of outer iterations do.(10)Since the reactions rate is calculated by collecting the contribution of each natural pixel to its averaged value inside the rotated pixels, take care to zero the rates before start. (11)Calculate scattering outer reaction rates and the fission outer reaction rate based on the previous outer iteration flux. (12)Calculate the new external outer iteration flux with (7.1) with zero external source and zero incident flux. (13)Initialize the internal loop flux 𝜙⁰ with the outer loop flux 𝜙^iterouter and normalizes it with the fission operator. (14)For inner iterations equal 1 to maximum number of inner iterations do.(15)Take care to zero the scattering reaction rate, but not the reaction fission rate, since it will not change inside the internal loop. (16)Calculate scattering reaction rate based on the previous inner iteration calculated flux. (17)Calculate the new internal iteration flux with (7.1) with zero external source and zero incident flux and external loop value of the fission reaction rate. (18)Verify if the difference between the old flux and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the inner iteration before the maximum number of inner iterations. Note that we may alternatively use the fission reaction rate error as criteria. (19)End inner iteration do.(20)Take the convergent inner loop value as the outer loop flux. (21)Calculates the new eigenvalue value as a ratio between the fission rate based on the actual outer iteration flux and the fission rate based on the previous outer iteration flux. (22)Verify if the difference between the old eigenvalue and the new iteration flux is respecting the established tolerance criteria and take the associated decision to break the outer iteration before the maximum number of outer iterations. (23)End outer iteration do.

Remark 7.2. The first and the inner part of the second algorithm are based on the contracting behavior of the scattering collision operator, that is, we have a compact operator with spectral radius less than one and the rate of convergence will be proportional to this ratio. Perhaps other rates of convergence may be found with other algorithms. The outer loop of the second algorithm is also based on compactness. The problem is to determine the first eigenvalue in the pencil problem 𝐼−𝐿_𝜎⁻¹𝐶^scattering𝜙=𝜆𝐿_𝜎⁻¹𝐶^fission𝜙 (7.5) in 𝐿²(𝑠^𝑛−1×Ω), that is, a perturbation of the identity by a compact operator in a compact weigh operator pencil problem. Its convergence with a decreasing sequence of eigenvalues may be proved by using the fact that the multiplicative region inside Ω has compact support and the composite operator is positive and compact. The fact that the first eigen angular flux is positive almost everywhere in 𝑆^𝑛−1×Ω is a classical result based on Jentzsch’s theorem for compact operators which leaves nonnegative those functions that are non negative. Once again the consistence has a positive consequence, since numerical artifacts have not been introduced in the approximation, and good properties of the continuous equation are preserved.

8. Numerical Simulations

8.1. The Transillumination Test Problem

In order to start the evaluation of the method proposed, we first study a subcritical transillumination test problem with no fission and no external sources. The flux for test is 𝜙^test=exp−0.0125⋅𝑦²+0.05⋅(𝑥−𝑀)for𝑗=24andzerosforotherviews. (8.1) The algorithm is used to solve problem (2.1) with source 𝑆(𝜔,𝑥)=𝜔⋅∇𝜙^test(𝜔,𝑥)+𝜎_𝑡(𝑥)𝜙^test(𝜔,𝑥)−_𝑆𝑛−1𝑓𝑥,𝜔^′⋅𝜔𝜙^test𝜔^′,𝑥𝑑𝜔^′,(𝜔,𝑥)∈𝑆^𝑛−1×Ω, (8.2) and incident boundary condition 𝜙(𝜔,𝜎)=𝜙^test(𝜔,𝜎)for(𝜔,𝜎)∈Σ⁻. (8.3) A high-contrast inclusion problem as shown in Figure 6 in a 24 views with 30×30 rectangular pixels/view model is solved with the algorithm proposed here which is based on (7.4) accordingly. Material properties are given in Table 1. As can be seen from Figure 7 comparing the given and the reconstructed solution showa that they are almost equal. From Figure 8 showing the absolute error for two views 9 and 24, the error is of order 10⁻⁸, bigger in the region with the high-contrast inclusion, but with the same magnitude. The absolute error variation with 𝐽 varying between 6 to 16 and 𝑀 from 10 to 24 is showed in Figure 9.

Table 1 

Absorption and scattering coefficient in the transillumination model.

Figure 6 

High-contrast nonhomogeneous model with 24 views and 30×30 pixels view.

(a)

(b)

Figure 7 

Given and reconstructed view 24 angular test flux in the model.

(a)

(b)

Figure 8 

Absolute error in view 24 and 9 angular test flux for the model.

Figure 9 

Maximum error variation for 𝐽 and 𝑀.

8.2. A Transillumination Problem

The same transillumination problem with no fission and no external sources is now investigated. The incident radiation boundary condition is 𝑔(𝑛_2𝐽,2𝐽)=1, for 𝑛_2𝐽=1,…,2𝑀 and 𝑔(𝑛_𝑗,𝑗)=0 for all others views. Figure 10 shows selected particles (radiation) flux for 𝑗=9 and 𝑗=24. We can see the depleted effect in the main direction of the incident flux 𝑗=24 due to absorption and outscattering and the buildup effect in the other direction 𝑗=9 due to the scattered flux. Also, it is possible to note the differences in magnitudes in the flux for the incident direction, 𝑗=24, with magnitude order of 1 and the other directions where there is no incident radiation, with magnitude order of 0,02, one up to 50-times the other. Another important information that can be determinated by the solver presented here is emergent radiation (flux) map shown in Figure 11. This set of figures graphically defines the flux in the emergent boundary Σ⁺𝜙=𝑔⁺onΣ⁺, (8.4) and the Cauchy data for the problem 𝐶_{𝑞𝑎,𝑓}=𝑔,𝑔⁺onΣ⁻×Σ⁺, (8.5) characterizes the graph for the albedo operator (the incident to emergent flux mapping) 𝒜_{𝑞𝑎,𝑓}∶𝐿²(Σ⁻,|𝜔⋅𝜈(𝜎)|𝑑𝜎𝑑𝜔)⟼𝐿²Σ⁺,|𝜔⋅𝜈(𝜎)|𝑑𝜎𝑑𝜔. (8.6)

(a)

(b)

Figure 10 

Selected radiation (particles) flux for illumination on view 24 in a high-contrast nonhomogeneous problem.

Figure 11 

Emergent flux for illumination on view 24 in the high-contrast model problem.

In the inverse problem, we ask if it is possible to determine the coefficients functions 𝜎_𝑡 and 𝑓 from the a priori knowledge of the albedo operator, that is, data such as that found in Figure 11 for different incident radiation 𝑔⁻. Finally, the convergent sequence {𝜙^iter;iter=1,…,31,…} is tested as a Cauchy sequence in the complete Hilbert space 𝐻. Figure 12 shows that before the eleventh iteration, we achieved convergence of order 0,008. If we claim for 9,8𝑒−008, it is achieved in the model problem in the 31 iteration.

Figure 12 

Cauchy error for 24 views and 30×30 pixels view.

8.3. Reactor with the EIR-2A and EIR-2B Data

The EIR-2A and EIR-2B are two problems found in Mordant [21] with data about values on cross sections, source and geometric domain definition presented in Table 2 and represented in Figure 13. The total and scattering cross-section is the same for the two problems. The source values is for the EIR-2A problem, and in this case the fission cross-section is zero. The fission cross-section is for the EIR-2B problem, and in this case the source term is zero. Results obtained are in accordance with the theory and are showed in figures. In all cases, the number of axial directions is 𝐼=1. Figure 14 shows the angular flux for direction 𝑗=1 in the source problem. Figure 15 shows the maximum Cauchy error flux for the angular flux function, where a high convergence rate is observed. The Figure 16 shows the maximum Cauchy error in the first eigenvalue function, where a not-so-high convergence rate is observed. Finally, Figure 17 shows the first eigen angular flux for direction 𝑗=1. These results showing a monotonically decreasing sequence of eigenvalues with the smallest accumulation point at 1 and a nonnegative angular flux are in complete agreement with theoretical results pointed out in Remark 7.2. We must finally note that the continuity in the direction perpendicular to 𝜔₁ has been artificially introduced by the graphic algorithm used to present the results. In reality, the calculations does not impose this restriction on the function.

Table 2 

Cross sections and sources values for EIR problems.

(a)

(b)

Figure 13 

EIR-2A and EIR-2B benchmark problems with data as in Table 2.

Figure 14 

Angular flux in J6M24 EIR 2A.

Figure 15 

Cauchy error in flux convergence for J6M24 EIR 2A.

Figure 16 

First eigenvalue in J6M24 EIR-2B.

Figure 17 

First eigenfunction in J6M24 EIR-2B.

9. Conclusions

The present methodology utilizes the natural basis concept, originally developed for image reconstruction, in the numerical solution of very classical and technological problems found in reactor physics calculations. The algorithms presented here may be considered as giving almost exact solutions in case where material properties and sources are constants in a small number of distinct regions inside the domain, as observed in the simulations EIR-2A and EIR-2B implemented. It does not array system matrices, and as a consequence few bigger problems signify more processing time. The algorithms are deduced by introducing, as less as possible, artifacts related with preprocessing the transport equation by doing inconsistent averages and simplifications hypotheses. It respects the basic characteristic of the transport equation, that is, no continuity requirements in the directions transverse to the propagation direction. It is also based on a scheme consistent with the exact equation. The algorithm implementation also shows compatibility with the present status of computer development.

Acknowledgment

The authors are partially funded by the Brazilian agencies CNPq, the COPPETEC Foundation, and CNEN.