Material void fraction evaluation with mixed neutron beams – theory and numerical experiment

The aim of this article is to describe a mathematical model and an associated numerical method for the basic problem of evaluating the void fraction of a slab material from the responses of a typical sourcedetector system for backscattered neutrons. The mathematical model and void fraction evaluation method reported herein: i) enable us to treat explicitly and exactly a mixed neutron beam, i.e. a neutron beam consisting of a monodirectional component with normal incidence upon the slab material and a component that is a smooth function of the angular variable and ii) allow for an arbitrary order of anisotropy of the neutron scattering phase function. In addition, the associated void fraction evaluation method is free from computer overflow exceptions, even for the challenging problems where the angular quadrature set is of high order and/or the slab material is exceedingly thick. We perform a numerical experiment to illustrate the positive features of both the mathematical model and evaluation method reported here, and we conclude this article with a discussion and directions for future work.


INTRODUCTION
We have been working on deterministic models in the subject of radiation transport theory [1] for the development of faster yet accurate alternatives to stochastic models such as Monte Carlo for both theoretical investigations and practical applications.As a result of our effort in this field, we have developed and reported years ago [2,3] on a mathematical model for a basic and important problem in nondestructive testing: that of evaluating the void fraction of a slab material by using the responses of a collimated neutron detector due to a neutron source with known activity [4,5].We hereafter refer to this mathematical model as the ZB (ZANI and BARROS) model for simplicity.The evaluation method derived from the ZB model provides the void fraction of a test slab in terms of the ratio between the total thickness of (undesirable) void regions and the thickness of the test slab.The testing equipment upon which the ZB model was based is very simple: an isotropic neutron beam shines the test slab on one side, and a collimated neutron detector is positioned at the incidence side along a number of angular directions to get the angular fluxes of backscattered neutrons in corresponding angular directions [6].
The ZB model was based on the widely used discrete ordinates (S N ) formulation of neutron transport theory [7] and on a numerical method ⎯ the spectral Green's function (SGF) method ⎯ that is free from spatial truncation error for slab-geometry S N neutron transport problems [8].For the class of admissible problems, the ZB model represented at that time a more accurate alternative to current models based on approximate methods for the spatial dependence of the angular neutron flux [9,10].However, the ZB model displays some theoretical/practical flaws due to the oversimplifying assumptions made in the design of its ingredients-the S N formulation and the SGF method.To do away with some of these flaws, we have recently developed and reported on an improved mathematical model [11], which is based on a more accurate S N formulation and on some analytic and numerical schemes adapted from a hybrid analytic/numerical method developed by the present author for more general radiation transport problems [12,13].The testing equipment associated with this improved model differs from the one in BARROS et al. [2] only in the beam shape.In our earlier work [11], the neutron beam is treated more realistically than before: it is monodirectional and normally incident upon one side of the test slab.In contrast to the ZB model, the model described in DE ABREU [11]: i) is able to treat explicitly and exactly a monodirectional neutron beam with normal incidence; ii) allows for an arbitrary (Legendre) order of anisotropy of the neutron scattering phase function; and iii) provides us with a void fraction evaluation method that is free from computer overflow exceptions even for angular quadrature sets of high order and arbitrarily thick slabs.
In this article we take a step further, and while keeping the positive features (ii) and (iii) above, we incorporate into the model described in DE ABREU [11] the practical case of a mixed beam, i.e. a beam consisting of a monodirectional component with normal incidence upon the slab material (as we have thus far considered) and a component that is a smooth function of the angular variable (not necessarily isotropic).This step was in order because it enables us to model beam transport problems associated with even more realistic neutron beams; an actual collimated beam can be thought of as a beam with a weakly divergent component and a monodirectional component.In addition, the resulting model provides us with a companion void fraction evaluation method effective for mixed beams.
An outline of the remainder of this article follows.In Section 2, we briefly describe a mathematical model that treats explicitly the case of mixed beams for a more accurate treatment of neutron beam transport problems of practical interest.In Section 3, we consider the model in Section 2, and develop a companion void fraction evaluation method.In Section 4, we perform a numerical experiment to illustrate the positive features of both the mathematical model and void fraction evaluation method.In Section 5 we close this article with concluding remarks.

A MATHEMATICAL MODEL FOR MIXED BEAMS
In this section, we outline a mathematical model that treats explicitly and exactly the case of mixed neutron beams.For a detailed discussion, the reader may consult a recent work by the author [12].

Transport Problem and Chandrasekhar Decomposition
We start with the governing equation of the transport of a mixed neutron beam through an anisotropically scattering homogeneous slab material.Let us consider the energy-independent integrodifferential transport equation where z is the spatial variable defined on a homogeneous domain Ω representing a test slab with boundaries at z = 0 (left) and at z = Z (right); μ is the cosine of the angle defined by the direction of neutron motion and the positive (from left to right) z-axis, and the remaining notation and nomenclature are standard in neutron transport theory [14].The quantity ψ(z,μ) is the angular flux of neutrons travelling in direction μ at position z; σ t is the macroscopic total cross section everywhere in Ω, and S(z,μ) is the scattering source given by .
The quantity σ s is the macroscopic scattering cross section everywhere in Ω; (2l+1)β l is the lth-order component of the Legendre expansion of the neutron scattering phase function, and P l denotes the lth-degree Legendre polynomial.The transport equation ( 1) is subject to the boundary conditions where ψ L denotes the intensity of the monodirectional (singular) component of the mixed neutron beam, which is normally incident (μ = 1) upon the left side of the slab; γ L (μ), μ > 0, is a nonnegative smooth function of μ representing the angularly continuous (regular) component of the incident beam and δ is to represent the Dirac delta.The energy-independent assumption can be understood here as either a thermal equilibrium condition [14] or a reference to the negligible loss of neutron energy between collisions in a slab material consisting only of heavy nuclides.Following the decomposition technique considered by CHANDRASEKHAR [15] in solving a basic problem in radiative transfer in planetary atmospheres, we decompose the neutron beam transport problem (1)-( 3) into the uncollided transport problem with the boundary conditions and the diffuse transport problem with the boundary conditions so that .The quantity 1 is the first-collided neutron source, which is a function of the uncollided component ψ u (z,μ) of the angular flux ψ(z,μ).Therefore, the diffuse problem ( 6)-( 7) can only be completely stated if the uncollided component is available.

The Uncollided Problem and the First-Collided Source
Since the uncollided problem ( 4)-( 5) represents an auxiliary problem defined in a purely absorbing medium with transparent boundaries, with no interior source, and with an incident particle beam upon the left side only, we must have ψ u (z,μ) = 0, for μ < 0 [16].For μ > 0, ψ u (z,μ) has the closed form We substitute solution (9) into the first-collided source (8), and after some Calculus we obtain ( )

The Diffuse Problem
At this point we are able to develop an S N formulation of our target problem -the diffuse problem ( 6)-( 7) -by considering a standard S N approximation [7] to equation (6) of the form are S N approximations to the diffuse component ψ d (z,μ) and first-collided source q u (z,μ), respectively; {μ m }, m = 1 : N, is a finite set of angular directions on the real interval [-1,1].In this article, we consider even-order Gauss-Legendre quadrature sets [7].Therefore, the discrete directions μ m , m = 1 : N, in the S N equations (11) are the (symmetric) roots of the Nth-degree Legendre polynomial P N .In addition, we have labelled the directions μ m so that μ m > 0 holds for m = 1 : N/2, μ m < 0 holds for m = N/2+1 : N, μ m-1 < μ m , m = 2 : N/2, and μ m+N/2 = -μ m , m = 1 : N/2.We remark that the nonnegative integer L in equations (11) means that the Legendre expansion of the scattering phase function has been truncated after (L+1) terms.Equations (11) are constrained by the discrete boundary conditions where the subscript -m is to denote the angular direction -μ m .The S N equations (11) and the discrete conditions (12) make up our S N formulation of the diffuse transport problem ( 6)- (7).
Solution to the S N problem ( 11)-( 12) can be written in terms of nontrivial elementary solutions to the homogeneous version of equations ( 11) and of particular solutions [12,13,17] in the open form , N : where α i , i = 1:N, are scalars depending upon the discrete conditions (12).The functions are elementary solutions of the homogeneous version of equations (11); ν i and a m (ν i ), i = 1 : N, m = 1:N, are the separation constants and the angular components of the elementary solutions (14), respectively.They can be determined with spectral characteristic methods [18] or by defining associate algebraic eigenvalue problems [17,18]; z i , i =1: N, are real scaling constants chosen to keep our void fraction evaluation method safe from computer overflow exceptions.A detailed discussion on the choice of the scaling constants z i can be found in DE ABREU [12].In the open form (13), the functions ( ) are variants to particular solutions used by CHANDRASEKHAR [15] to express the general solution of the classic albedo problem of atmospheric radiative transfer, wherein f m , m = 1: N, are constants to be determined upon substitution of the particular solution (15) into equations (11).The resulting linear matrix equation is solved by splitting and analytic inversion techniques [11,17].At this point we are done with the mathematical model for the void fraction evaluation method of Section 3.

A COMPANION VOID FRACTION EVALUATION METHOD
We begin by stating the underlying principles of the void fraction evaluation method of this section.Let us assume that the responses of the collimated detector of the testing equipment described in Section 1 are available for a test slab of known thickness, say, Z'.If we consider the principle of invariance of the angular flux of neutral particles along straight line paths in source-free void regions [16,19], then we can draw the following conclusion: if the angular fluxes of backscattered neutrons in a number of angular directions are less than the corresponding ones for a homogeneous slab material of the same thickness Z', then the test slab has void regions.In this regard, a slab with void regions behaves as a thinner homogeneous slab, i.e. a homogeneous slab whose thickness, say, Z is the difference between the thickness Z' of the test slab and the total thickness of void regions.So, we may think of using the mathematical model in Section 2 for the problem of evaluating the thickness Z of such a thinner homogeneous slab associated with the test slab.The void fraction will then be given by the ratio between the total thickness of void regions (Z'-Z) and the thickness Z' of the test slab.These are the general lines of the void fraction evaluation method of this section.Note that the only remaining task in the method is the evaluation of the thickness Z of the thinner homogeneous slab associated with the test slab.Since the Z-evaluation procedure here closely follows the two-step procedure described in DE ABREU [11], presentation will be cursory and driven by the major differences found here and there.
We begin the first step of our Z-evaluation procedure by recalling that the collimated neutron detector can be positioned along an arbitrary number of directions, providing thus the angular fluxes of backscattered neutrons along corresponding angular directions.Therefore, we can assume that the angular fluxes ψ d (z,μ m ) are available quantities at z = 0 for m = N/2+1 : N. From our S N approximation (13), we may write or, using results ( 14) and ( 15), Now we use the left boundary condition (12) together with our S N approximation (13) Equations ( 17) and ( 18) form a real system of N linear algebraic equations in the N unknowns α i , i = 1 : N. It can be written in the matrix form b α A = (19) where A is an NxN real matrix with (known) entries b is an N-dimensional column matrix with (known) entries (21) and α is an N-dimensional column matrix whose entries are the unknown coefficients α i .The coefficients α i , i = 1: N, in approximation ( 13) can be made available after solving the linear system (19).For this task, we have considered a computer version of the well-established Crout's method with implicit pivoting [20].
We begin the second step of our Z-evaluation procedure by noting that no neutron enters the slab through the right side, cf.boundary condition (3) at z = Z.This implies that the norm of the partial neutron current [14] vanishes at z = Z.Let us consider the S N -P N-1 approximation [7] to the norm of the partial neutron current (22), where and the integral on the right-hand side above can be computed with a half-range Gauss-Legendre quadrature.We substitute our S N approximation ( 13) into (23), and after some algebra we obtain the S N approximation where and are known quantities.Since the norm (22) vanishes at z = Z, we use our S N approximation (25 28) is a nonlinear algebraic equation in the unknown parameter Z, the thickness of the thinner homogeneous slab associated with the test slab.To solve equation (28), we have considered a computer version of the Newton-Raphson method [20] with the initial guess Z'.Finally, the void fraction ρ of the test slab is given by .

NUMERICAL EXPERIMENT
In this section we perform a numerical experiment in connection to a relevant application in radiation therapy.We consider a slab model of a 10 cm-thick plate-type test block made of bismuth.This block is part of a neutron-photon beam filtering/shaping assembly (a stack of filtering/shaping blocks) designed for improving the quality of an incident beam from a nuclear reactor beam port for boron neutron capture therapy (BNCT) applications [21].The major role of the bismuth block is to remove photons from the beam, while largely sparing neutrons in the energy range of interest for BNCT.In this experiment, we are interested in nondestructively testing the bismuth block by evaluating its void fraction along the direction defined by the beam channel.We note that bismuth isotopes are heavy nuclides, with mass numbers as high as 216, and so the energy-independent assumption made in Section 2 is very appropriate here.
The experiment is as follows.We consider a beam of 1.2 MeV neutrons with ψ L = 10 17 cm -2 s -1 sr -1 , and we set γ L (μ) = 10 3 μ 4 cm -2 s -1 to represent the regular component of a weakly divergent neutron beam.We assume that the test block has a void fraction of 0.01, i.e. the thinner homogeneous block associated with the test block is 9.9 cm-thick, cf.expression (29).Now, instead of considering a Monte Carlo code, we use the S N method described in DE ABREU [12] for N = 200 (a high order) to solve the neutron beam transport problem for the 9.9 cm-thick block.This is just to mock up the one hundred (N/2) responses of the collimated detector for the backscattered fluxes ψ d (0,μ m ), m = 101: 200.At this point we may look at the void fraction as an unknown quantity, and we use the evaluation method of Section 3 to get it.The macroscopic total and scattering cross sections of bismuth for 1.2 MeV neutrons are equal to 0.147455 cm -1 and 0.143520 cm -1 , respectively.The Legendre scattering phase function data for L = 6 and 1.2 MeV neutrons for bismuth are given in Table 1.We should note that the results from our numerical experiment were generated with FORTRAN [20] versions of the S N method described in DE ABREU [12] for the backscattered angular fluxes, and of the methods (Crout and Newton-Raphson) embedded in the evaluation method of Section 3 for the void fraction of the test block.Our FORTRAN codes were executed on an IBM-compatible PC (1.4GHz-clockIntel Pentium 4 processor and 256 Mbytes of RAM) running on a GNU/Linux platform.The one hundred values of ψ d (0,μ m ), m = 101: 200, generated with our S N method, and particularly the two hundred values of α i , i = 1: 200, from the (Crout) solution of system (19), are bulky to be properly tabulated here.For this reason, these numerical values can be made available electronically upon request to the author.In Table 2 we display the first ten estimates for the thickness Z in the corresponding first ten Newton-Raphson iterations on the nonlinear equation (28) for the initial guess 10 cm.In Figure 1, we plot a number of estimates for Z to illustrate the stable work of our computer version of the Newton-Raphson method to give the value 9.9 cm.As shown in our earlier work [11], this is in deep contrast to the void fraction evaluation method derived from the ZB model, which suffers from uncontrolled overflow exceptions, and for this reason should be restricted to low order N and moderately thick slabs.We finally used expression (29) to get the expected value 0.01 for the void fraction of the test block.

CONCLUDING REMARKS
We have successfully incorporated into our previous mathematical model the important case of mixed beams for a more realistic treatment of the problem of evaluating void fractions of material slabs.The mathematical model presented here enables us to treat neutron scattering with an arbitrary (Legendre) order of anisotropy, to represent explicitly and exactly a neutron beam with a normally incident, monodirectional component and a weakly divergent regular component, and to introduce the necessary modifications into the void fraction evaluation method described in our earlier work [11].
We would like to stress here two important points.One is that the deterministic model outlined in Section 2 is not merely an alternative to more general, computer time-memory-consuming stochastic models such as Monte Carlo for computer simulation of the transport processes relevant to the application considered here.The model in Section 2 is actually the basis of the void fraction evaluation method in Section 3. Without the model in Section 2, the method in Section 3 does not exist.The other point is that we have shown the applicability of the void fraction evaluation method developed in Section 3 with a computer simulation of a sample problem.However, the application of the method of Section 3 to practical problems, with actual backscattering measurements from a real source-detector system, is expected to be subjected to a number of experimental constraints.For example, for fixed slab material and thickness, and given a backscattering angle, the (detector) counting rate is expected to depend inversely on the distance between the detector and the irradiated side of the slab, since a real test slab is finite.Moreover, the counting rate in this situation depends also on the neutron scattering phase function of the slab material.If neutron scattering is highly backwardly peaked for instance, then the counting rate is expected to increase for increasing backscattering angles.These aspects (and so many others) must be taken into account when considering the application of our void fraction evaluation method to practical problems, since for these problems, angular fluxes of backscattered neutrons come from actual measurements (not from discrete ordinates or Monte Carlo).
We are aware that basic and important points remain unexplored.For example, neutron scattering is dependent upon the energy of the incident neutron, and this important fact of nature has not yet been considered for an even more realistic treatment of the void fraction evaluation problem considered here.We intend to extend the model and the method described here to treat energy-dependent problems through a multigroup formalism that is well established in neutron transport theory [14] in near future.

Figure 1 :
Figure 1: Stable work of our computer version of the NR method.

Table 1 :
Scattering phase function data for bismuth.

Table 2 :
Successive Newton-Raphson estimates for the thickness Z.