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

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

Abreu, M.P. de

Departamento de Modelagem Computacional - Instituto Politécnico – UERJ CP: 97282 - 28610-974, Nova Friburgo, RJ, Brasil. e-mail: [e-mail]

ABSTRACT

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 source-detector 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.

Keywords: mathematical modeling, discrete ordinates, mixed neutron beams, material void fraction.

1 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.

2 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].

2.1 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 Ω; (2ℓ +1)β_ℓ is the ℓth-order component of the Legendre expansion of the neutron scattering phase function, and P_ℓ denotes the ℓth-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 ψ(z, µ) = ψ^u (z µ) + ψ^d (z, µ) , z ∈ Ω, -1< µ < 1. The quantity

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.

2.2 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

2.3 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

where (Z) ≅ ψ^d (z, µ_m) and (Z) ≅ q^u (z, µ_m) 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

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.

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) to get

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

where A is an NxN real matrix with (known) entries

b is an N-dimensional column matrix with (known) entries

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) to get

Equation (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

4 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¹⁷ cm^-2s^-1sr^-1, and we set γ_L(µ) = 10³µ⁴ cm^-2s^-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–clock Intel 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.

5 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.

6 ACKNOWLEDGMENTS

The author is grateful to two anonymous reviewers for their helpful comments and constructive criticisms.

7 REFERENCES

Data de envio: 11/04/2006 Data de aceite: 09/02/2009

• [1] DUDERSTADT, J.J., MARTIN, W.R., Transport Theory, 1 ed., New York, John Wiley & Sons, 1979.
• [2] BARROS, R.C., ABREU, M.P. DE, ZANI, J.H., DE MELLO, J.A.M., "Recent advances in spectral nodal methods applied to one-speed discrete ordinates inverse transport problems in plane geometry", In: Proceedings of the Joint International Conference on Mathematical Methods and Supercomputing for Nuclear Applications, v. 2, pp. 1235-1244, Saratoga Springs, New York, October 1997.
• [3] ZANI, J.H., BARROS, R.C., ABREU, M.P., "Analytical computational method for estimating material void fraction of slabs using neutron beams transport", In: 4th Brazilian Meeting on Nuclear Applications (CD-ROM Proceedings), Poços de Caldas, Brazil, August 1997.
• [4] BARTON, J.P., Neutron Radiography (4): Including Radioscopy and Complementary Inspection Methods Using Neutrons, 1 ed., Pennsylvania, Gordon & Breach, 1993.
• [5] SMITH, F.A., A Primer in Applied Radiation Physics, 1 ed., Singapore, World Scientific, 2000.
• [6] KNOLL, G.F., Radiation Detection and Measurement, 3 ed., New York, John Wiley & Sons, 1999.
• [7] LEWIS, E.E., MILLER Jr., W.F., Computational Methods of Neutron Transport, 2ed., Lagrange Park, IL, American Nuclear Society, 1993.
• [8] 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, v. 104, n. 3, pp. 199-208, March 1990.
• [9] BADRUZZAMAN, A., "Nodal methods in transport theory", In: Lewins, J., Becker, M. (eds), Advances in nuclear science and technology, 1 ed., v. 21, chapter 3, New York, USA, Plenum Press, 1990.
• [10] MCCORMICK, N.J., "Inverse radiative transfer problems: a review", Nuclear Science and Engineering, v. 112, n. 3, pp. 185-198, February 1992.
• [11] ABREU, M.P., "An improved mathematical model for a void fraction problem in nondestructive testing", Nondestructive Testing and Evaluation, v. 20, n. 3, pp. 175-190, September 2005.
• [12] ABREU, M.P., "Mixed singular-regular boundary conditions in multislab radiation transport", Journal of Computational Physics, v. 197, n. 1, pp. 167-185, June 2004.
• [13] ABREU, M.P., "A two-component method for solving multislab problems in radiative transfer", Journal of Quantitative Spectroscopy and Radiative Transfer, v. 85, n. 3-4, pp. 311-336, May 2004.
• [14] DAVISON, B., Neutron Transport Theory, 1 ed., Oxford, Clarendon Press, 1957.
• [15] CHANDRASEKHAR, S., Radiative Transfer, 1 ed., Oxford, Clarendon Press, 1950.
• [16] CASE, K.M., ZWEIFEL, P.F., Linear Transport Theory, 1 ed., Reading, MA, Addison-Wesley, 1967.
• [17] SIEWERT, C.E., "A concise and accurate solution to Chandrasekhar's basic problem in radiative transfer", Journal of Quantitative Spectroscopy and Radiative Transfer, v. 64, n. 2, pp. 109-130, January 2000.
• [18] ABREU, M.P., "On the spectrum of the one-speed slab-geometry discrete ordinates operator in neutron transport theory", Annals of Nuclear Energy, v. 25, n. 15, pp. 1209-1219, October 1998.
• [19] THOMAS, G.E., STAMNES, K., Radiative Transfer in the Atmosphere and Ocean, 1 ed., New York, Cambridge University Press, 1999.
• [20] PRESS, W.H., TEUKOLSKY, S.A., VETTERLING, W.T., FLANNERY, B.P., Numerical Recipes in FORTRAN 77: The Art of Scientific Computing, 2 ed., New York, Cambridge University Press, 1996.
• [21] INGERSOLL, D.T., SLATER, C.O., WILLIAMS, L.R., "BNCT filter design studies for the ORNL tower shielding facility", Transactions of the American Nuclear Society, v. 75, n. 1, pp. 28-30, November 1996.

Publication Dates

• Publication in this collection
  11 May 2009
• Date of issue
  Apr 2009

History

• Accepted
  09 Feb 2009
• Received
  11 Apr 2006