Three-dimensional transport theory: An analytical solution for the internal beam searchlight problem, II

Annals of Nuclear Energy 36 (2009) 1242–1255

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Three-dimensional transport theory: An analytical solution for the internal beam searchlight problem, II Barry D. Ganapol a,*, Drew E. Kornreich b a b

Department of Mechanical and Aerospace Engineering, University of Arizona, Tucson, AZ 85721-0119, USA Applied Engineering Technology 2, Los Alamos National Laboratory, TSA-7, MS F609, Los Alamos, NM 87545, USA

a r t i c l e

i n f o

Article history: Received 4 August 2008 Accepted 29 January 2009 Available online 21 June 2009

a b s t r a c t Multidimensional semi-analytical particle transport benchmarks to provide highly accurate standards of assessment are few and far between. Because of a well-established 1D theory for the analytical solution of the transport equation, it is sometimes possible however, to ‘‘bootstrap” 1D solutions to give more comprehensive solution representations. Here, we propose the internal searchlight problem in a half space, designated ISLP/HS, as a multidimensional benchmark to be constructed from 1D solutions. This is a variation of the usual SLP/HS where a source emits within the half space rather than striking its surface. Our primary interest is in the exiting intensity at the free surface established through a new Fn formulation. The benchmark features true 2/3D particle transport through integration of a point kernel to simulate 2/3D source emission. In this way, we accommodate a solid or hollow cylindrical source and a general line source in addition to the standard point, ring and disk sources featured in previous investigations. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Multidimensional semi-analytical benchmarks are rare. The reason is obvious since, in a rigorous sense, a semi-analytical benchmark requires an analytical solution as its basis, and multidimensional analytical solutions to the transport equation are virtually non-existent. Because of a well-established 1D theory however, it is sometimes possible to ‘‘bootstrap” a 1D analytical solution to establish a more comprehensive solution representation. The searchlight problem (SLP) (Barichello and Siewert, 2000 and Williams, 2007; see the references therein) is one such case where a 1D solution lies at the heart of a multidimensional transport solution. As a result, we can consider the SLP as a potential multidimensional benchmark. Elliott (1955) first considered a variation of the usual SLP. He treated a monoenergetic point source emitting isotropically in the interior of a semi-infinite medium and found the flux distribution at the free surface. In the following, we refer to this case as the interior SLP for the half space (ISLP/HS) and consider its solution using a 1D transport method. Solving for the flux for the SLP and ISLP, however, is easier said than done, as an accurate numerical evaluation involves a Hankel transform inversion of an ‘‘inner” 1D transport solution. This re-

* Corresponding author. E-mail addresses: [email protected] (B.D. Ganapol), drewek@ lanl.gov (D.E. Kornreich). 0306-4549/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2009.01.021

quires a double integration to arrive at the ‘‘outer” multidimensional solution. Thus, even the 2D/SLP presents a significant challenge to achieve the necessary high accuracy of a benchmark. Here, we establish a series of semi-analytical ISLP/HS benchmarks based on a new Fn formulation for the inner 1D solution. The new formulation avoids the complications of a past Fn form (Siewert and Dunn, 1983) in evaluating the matrix elements by expressing them in a closed form. Our investigation begins with a point source located a distance D above the free surface within a half space. Particular 2D cases will include nadir (perpendicularly downward) and zenith (perpendicularly upward) directed point sources. We impose the monoenergetic approximation and assume an isotropically scattering semi-infinite medium of infinite transverse extent. The point source then serves as the basis for several comprehensive sources that are more relevant for benchmarking numerical transport algorithms. By a change of coordinate, the point source is then reconfigured off-center to become a point kernel for general 2/3D sources to give multidimensional flux variations. The cases we consider find importance in radiative transfer experiments to determine cloud absorption and scattering properties. In addition, the ISLP/HS is at the heart of a numerical inversion strategy to determine source depth for a radiation field with possible Department of Homeland Security applications. We begin with the development of the analytical solution representation for a radially symmetric source and proceed to nonsymmetric sources.

1243

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

2. Theory for a radially symmetric source 2.1. The 3D transport equation for a beam With the notation of Williams (2007), the appropriate particle transport equation in 3D cylindrical geometry assuming isotropic scattering is



l

Δ

 Z @ @ c þx þ 1 Iðz; q; XÞ ¼ dX0 Iðz; q; X0 Þ þ Sðz; q; XÞ; @z @q 4p 4p

Vacuum Fig. 2. Embedded point beam source.

ð1aÞ where the coordinates are shown in Fig. 1. The boundary conditions at the free surface z ¼ 0 and infinitely far away are

Ið0; q; XÞ ¼ 0;

l > 0; / 2 ½0; 2p; ð1bÞ

lim Iðz; q; XÞ < 1:

z!1 jqj!1

2.2. Fourier transformed equation If we apply a 2D cylindrical Fourier transform to Eq. (1a), the following transformed transport equation results:



The radially distributed interior beam source is located a distance D above the surface and emits in direction ð/0 ; l0 Þ

Sðz; q; XÞ ¼ SðqÞdðD  zÞdðl  l0 Þdð/  /0 Þ;

ð1cÞ

where

jlj; jl0 j 6 1;

/0 ; / 2 ½0; 2p:

 @ þ uðl;/; kÞ Iðz; l; /;kÞ @z Z 1 Z 2p c ¼ dl0 d/0 Iðz; l0 ;/0 ; kÞ þ SðkÞdðD  zÞdðl  l0 Þdð/  /0 Þ; 4p 1 0 ð3aÞ

l

subject to

In the following, the dependence on the beam coordinates is suppressed. We also normalize the source such that

Ið0; l; /; kÞ ¼ 0;

Z

z!1

dqSðqÞ  1:

l > 0;

ð3bÞ

lim Iðz; l; /; kÞ < 1;

ð3cÞ

with

All Space

For example, Fig. 2 shows an interior point beam source emitting in a specified direction. Note, for later use, that the point source distribution at the origin in this example is SðqÞ  dðqÞ=2pq. The solution to Eq. (1) for a point beam is central to our 3D analysis. The particle longitudinal and radial positions, z, qðq; hÞ, are measured in units of mean free path, and the particle direction unit vector is X. The scalar flux

Iðz; qÞ 

Z

dXIðz; q; XÞ;

uðl; /; kÞ  1  ikð1  l2 Þ1=2 cosð/  wÞ:

The vector kðk; wÞ of the transform is shown in Fig. 1. The Fourier transforms in Eq. (3a) are

Z Iðz; l; /; kÞ  dqeikq Iðz; q; XÞ; Z SðkÞ  dqeikq SðqÞ:

ð2Þ

Iðz; kÞ 

cos−1 ( μ)

ð4bÞ

Z

dXIðz; l; /; kÞ

or

Iðz; kÞ ¼

Z

dqeikq Iðz; qÞ

and on inversion is

Iðz; qÞ ¼

z

ð4aÞ

The transformed scalar flux is

4p

at the surface ðz ¼ 0Þ is the desired physical quantity. If the beam emits perpendicularly to the surface ðl0 ¼ 1Þ and the source SðqÞ is radially symmetric (independent of azimuth hÞ, then the flux will also be independent of azimuth. In this way, with nadir (1) and zenith (+1) emission, we reduce the problem to two dimensions.

ð3dÞ

Z

1 2

ð2pÞ

dkeikq Iðz; kÞ:

ð5Þ

We are now in position to solve Eq. (3) as a 1D transport equation. Unlike Williams (in press) who used the Wiener–Hopf solution representation, however, we follow the Fn procedure of Siewert and Dunn (1989). This procedure is more numerically oriented and has greater applicability. 2.3. Integral equation for Iðz; kÞ Integrating Eq. (3a) along the particle trajectory and including the boundary conditions, we have for l > 0

ψ

θ

φ

Iðz; l; /; kÞ ¼ SðkÞ þ

k Fig. 1. Geometry of the direct and transform spaces.

and

eðzDÞ=Uðl;/;kÞ

c 4pl

Z 0

l z

0

dðl  l0 Þdð/  /0 ÞHðz  DÞ 0

dz eðzz Þ=U Iðz0 ; kÞ

ð6aÞ

1244

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

Iðz; l; /; kÞ ¼ SðkÞ

eðDzÞ=Uðl;/;kÞ Z

c

þ

l

4pl

with

dðl þ l0 Þdð/  /0 ÞHðD  zÞ

nl

1

0 ðz0 zÞ=U

dz e

0

Iðz ; kÞ;

ð6bÞ

z

Then, integrating Eq. (10) over equation for the scalar flux:

where

Uðl; /; kÞ 

l uðl; /; kÞ

aðlÞ : bðlÞ

;

/ðz; l Þ ¼

u0  1  ikð1  l20 Þ1=2 cosð/0  wÞ: Integrating these equations to form the scalar flux Iðz; kÞ generates the following integral equation:

Iðz; kÞ ¼ SðkÞejDzj=U0 ðjl0 j;/0 ;kÞ Hððz  DÞ=l0 Þ=jl0 j Z c 1 0 þ dz Kðjz  z0 j; kÞIðz0 ; kÞ; 2 0 with

Kðz; kÞ 

Z

1

dl

e z=l

l

0

Kðz; kÞ ¼

1

ð7aÞ

0

l

l ð1 þ k l 2

K  ðzÞ 

Z

dl ez=nðlÞ : l aðlÞ

1

0

ð7bÞ

l

dl e

ejðDzÞ=nj ½f ðl; l ÞHððz  DÞ=nÞ laðlÞ 0 Z c 1 0  þ f ðl; l ÞHððD  zÞ=nÞ þ dz K ðjz  z0 jÞ/ðz0 ; l Þ; 2 0 ð11aÞ dl

f ðl; l Þ  jljaðlÞdðl  l Þ;



/n ðzÞ ¼ ejðDzÞ=n j Hððz  DÞ=n Þ þ

l

2 Þ1=2

c 2

Z

1

0

dz Kðjz  z0 j; kÞI ðz0 ; kÞ:

0

We now solve this equation using 1D pseudo transport theory.

n  l

0

dz K  ðjz  z0 jÞ/n ðz0 Þ;

0

aðl Þ : bðl Þ

Consider the following generalized transport equation in the half-space (Siewert and Dunn, 1983):

 Z @ c 1 0 dl /ðz; l0 ; l Þ þ bðlÞ /ðz; l; l Þ ¼ @z 2 1 þ f ðl; l Þdðz  DÞ;

l > 0;

lim ½/ðz; l; l Þ < 1;

z!1

ð9aÞ ð9bÞ ð9cÞ

where a and b are even functions of l. l is a parameter we define below in addition to f ðl; l Þ. Again, applying integral transport theory, as above, yields for l>0

f ðl; l Þ ðzDÞ=n /ðz; l; l Þ ¼ Hððz  DÞ=nÞ e laðlÞ Z z 1 c 0 0 dz eðzz Þ=n /ðz0 ; l Þ þ laðlÞ 2 0

I ðz; kÞ ¼ SðkÞ/ ðz; kÞ;

ð13Þ

where the k dependence is now explicitly included. Hence, this requires

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 þ k l2 ; 2

l2

for Eq. (11b) to match Eq. (7c). Therefore, if we solve for / ðz; l Þ, I ðz; kÞ is determined for the two sources, and we find the final solution through the Fourier transform inversion of Eq. (5). Note that the source radial configuration remains arbitrary. Only the directional component of the source has been prescribed. 2.5. Singular integral equation for exiting pseudo flux When we replace l by l in Eq. (9a), multiply the result by ez=s and integrate over l and z on ½0; 1Þ and [1, 1], respectively, we find (see Appendix A)

KðsÞ/ ðs; l Þ ¼ s

Z

1

dl0

0

n0 /ð0;l0 ; l ÞseD=s n s 0

Z

1

dl0 f ðl0 ; l Þ ; bð l0 Þ n0 s 1 ð14aÞ

where

ð10aÞ

and

f ðl; l Þ ðDzÞ=n Þ¼ HððD  zÞ=nÞ e laðlÞ Z 1 1 c 0 0 þ dz eðz zÞ=n /ðz0 ; l Þ; laðlÞ 2 z

ð12cÞ

Note that /ðz; l Þ has been replaced by /n ðzÞ to indicate the true dependence of the scalar flux on n and for notational convenience. Formally, therefore, Eq. (8±) are recovered, when K   KðjD  zj; kÞ and n  1, giving

bðlÞ ¼ 1 þ k

/ðz; l; l

1

ð12bÞ

aðlÞ ¼

2.4. The pseudo problem



Z

where

ð8Þ

/ð0; l; l Þ ¼ 0;

c 2

ð7cÞ

:

I ðz; kÞ ¼ SðkÞejðDzÞj Hððz  DÞÞ þ

laðlÞ

ð12aÞ

then

With this form, an equivalent pseudo problem follows (Siewert and Dunn, 1983) allowing direct determination of I. At this point, it is prudent to focus on two common source emissions – upward and downward directed. Each case is obtainable from Eq. (7a). The upward and downward directed sources are for l0  1 giving for Eq. (7a)



ð11bÞ

If

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kz 1  l2 ; J0

ð1þk2 2 Þ1=2 z=

1

where

where J0 is the Bessel function of the first kind. After some manipulation (Rybicki, 1965), we find

Z

Z

l gives the following integral

ð10bÞ



/ ðs; lÞ 

Z

1

dze

z=s

/ðz; l Þ;

ð14bÞ

0

KðsÞ  1  cLðsÞ;   czðsÞ zðsÞ þ 1 LðsÞ  ln ; 2 zðsÞ  1 s zðsÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 1  s2 k

ð14cÞ ð14dÞ ð14eÞ

1245

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

When s is set to s  m  ie and allowed to approach m on the cut ½c; c, where

1

giving the following decomposition into uncollided and collided flux components:

/ð0; l; l Þ ¼ /0 ð0; l; l Þ þ /c ð0; l; l Þ

c  pffiffiffiffiffiffiffiffiffiffiffiffiffi2ffi ; 1þk

and subsequently

as shown in Fig. 3, we find [using Eq. (B3) of Appendix B],

 1 K ðmÞ/ ðm; l Þ ¼ m dl0 n0 P 0  ipdðn0  mÞ /ð0;l0 ; l Þ n m 0   Z 1 dl0 1 0 D=m  me  i p dðn  m Þ f ðl0 ; l Þ: P 0 n0  m 1 bðl Þ ð15Þ Z



1

kðnÞbðlÞ/c ð0; l; l Þ 

K ðmÞ/ ðm; l Þ Z 1 ¼m dl0 P  meD=m

Z

1

1

 n0 dl /ð0; l0 ; l Þ  ipm2  /ð0; lðmÞ; l Þ dn lðmÞ n m  0 dl 1 ip dl 0  f ð l ; l Þ  f ðlðmÞ; l Þ; P bðl0 Þ n0  m bðlðmÞÞ dn lðmÞ 0

ð16Þ with

m

dl0 P

0

n0 / ð0; l0 ; l Þ n n c 0

ð21aÞ

  0 dl0 eD=n  eD=n f ðl0 ; l Þ 0 n0  n 1 bðl Þ Z 1 dl0 1 f ðl0 ; l Þ þ eD=n 0 0 0 bðl Þ n þ n

Z

1

and, as shown in Appendix B,

kðmÞ  1 

  1 þ zðmÞ czðmÞ  ; ln  2 1  zðmÞ

m

zðmÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 1  m2 k

q ðnÞ ¼

Using Eq. (B6a) and adding and subtracting Eq. (16±), respectively, results in

Z

n0 /ð0; l0 ; l Þ n m 0 Z 1 dl0 1 f ðl0 ; l Þ;  meD=m P 0 0 1 bðl Þ n  m 1

dl0 P

eD  eD=n eD=n Hð 1Þ þ Hð1Þ: 1n 1þn

ð21bÞ

We now focus our attention on numerically solving the Eq. (21a) as singular integral equations. It is here where we apply the Fn method.

0

ð17aÞ

2.6. The Fn approximation The Fn spectral approximation is

eD=m / ðm; l Þ ¼ m/ð0; lðmÞ; l Þ  f ðlðmÞ; l Þ: 2bðlðmÞÞ bðlðmÞÞ c

/c ð0; l; l Þ ¼ ð17bÞ

Finally, eliminating / ðm; l Þ between Eqs. (17a) and (17b), yields

  eD=n kðnÞbðlÞ /ð0; l; l Þ  f ðl; l Þ nbðlÞ Z Z c 1 0 n0 c 1 dl0 P /ð0; l0 ; l Þ  eD=n dl P 0 ¼ 2 0 2 1 bðl0 Þ n n 1  0 f ðl0 ; l Þ; n n

eD=n f ðl; l Þ nbðlÞ

N 1 c du X aa;n ðkÞP a ð2u  1Þ; 2 dl a¼0

ð22aÞ

where

u

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 þ k n;

0 6 u 6 1;

dl bðlÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; du 2 1 þ k 1  rðkÞ2 u2

ð22bÞ ð22cÞ

k rðkÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 1þk ð18Þ

where we have replaced m by n. Eq. (18) indicates that the uncollided component is (let c be zero)

/0 ð0; l; l Þ ¼

1

For the nadir/zenith source emission given by Eq. (12a), we find

lðmÞaðlðmÞÞ : bðlðmÞÞ

kðmÞ/ ðm; l Þ ¼ m

Z

where

qn ðnÞ ¼

Performing the integration over the delta function gives

c 2

c ¼ qn ðnÞ; 2



0

ð20Þ

ð19Þ

Pa ð2u  1Þ is the shifted Legendre polynomial, which we assume as our basis functions. The derivative included as the second factor makes the evaluation of the matrix elements most convenient. When we introduce the Fn approximation into Eq. (21a), there results N1 X

Ba ðuÞaa;n ðkÞ ¼ qn ðnÞ;

ð22dÞ

a¼0

where

Re(s)

-γ Fig. 3. Cut s-plane.

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 e a ðuÞ; 1 þ k ð1  rðkÞ2 u2 Þ1=2 kðnÞPa ð2u  1Þ  B e a ðuÞ ¼ c ½da;0  2uQ a ð2u  1Þ: B 2

Ba ðuÞ ¼

Im(s)

γ

ð22eÞ ð22fÞ

Q a ð2u  1Þ is the Legendre function of the second and appears because of the leading derivative factor in Fn approximation. As is usual for singular integral equations, an additional constraint applies because at s ¼ m0

1246

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255



 c m0 xðm0 Þ þ 1 Kðm0 Þ  1  ln ¼ 0; 2 xðm0 Þ  1

3.1. Standard sources 3.1.1. Centered and off center point source We first consider a point source centered at the origin. For this case,

with

m0 xðm0 Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1  k m20

SðqÞ ¼

giving for Eq. (14a)

0¼

Z

1

dl0

0

n0 /ð0; l0 ; n Þ  eD=m0 0 n  m0

Z

1

1

Ba ðu0 Þaa;n ðkÞ ¼ qn ðm0 Þ;

I ð0; qÞ ¼ eD

ð22gÞ

a¼0

where

c Ba ðu0 Þ ¼  ½da;0  2u0 Q a ð2u0  1Þ; 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 u0  1 þ k m0 :

ð22hÞ ð22iÞ

This form of the Fn approximation is different from previous Fn forms. In particular, the matrix elements, Ba , are computed through straightforward recurrence rather than a specialized recurrence or numerical integration. The required scalar flux at the half-space surface is the inversion

Iðz; qÞ ¼

Z

1 ð2pÞ

dkSðkÞ/n ð0; kÞeikq ;

2

ð27Þ

and SðkÞ ¼ 1 giving

dl0 f ðl0 ; l Þ : bðl0 Þ n0  m0

When the Fn approximation is introduced into this equation, there results N1 X

dðqÞ 2p q

ð23aÞ

dðqÞ c Hð 1Þ þ 2pq 4p

Z

1

dkkJ0 ðqkÞa0; ðkÞ:

ð28Þ

0

The flux for the point source allows the construction of more complicated and therefore more meaningful sources. We now consider a point source situated off center at z ¼ D as shown in Fig. 4. The appropriate flux at any radial position at the surface is found from Eq. (28) by noting that

q0  ½q2 þ q20  2qq0 cosðh  h0 Þ1=2 : Thus

I ð0; q;q0 Þ ¼

1 2p

Z

1

0

 dkkJ 0 k½q2 þ q20  2qq0 cosðh  h0 Þ1=2 /0; ð0; kÞ: ð29Þ

We obtain an alternative expression with the substitution

 J 0 k½q2 þ q20  2qq0 cosðh  h0 Þ1=2 ¼

1 X

en J0 ðkqÞJ0 ðkq0 Þ cosðnðh  h0 ÞÞ;

ð30aÞ

n¼0

where the 1D transport inner solution (scalar flux) is

/n ð0; kÞ ¼ /0;n ð0; kÞ þ /c;n ð0; kÞ:

where

en  2  dn;0 to give

ð23bÞ I ð0; q; q0 Þ ¼

The uncollided flux for our two sources is

/0; ð0; kÞ ¼ eD Hð 1Þ:

1 X n¼0

en

1 2p

Z 0

1

dkkJn ðkqÞJ n ðkq0 Þ/0; ð0; kÞ cosðnðh  h0 ÞÞ:

ð23cÞ

ð30bÞ

Note that the k dependence of the 1D solution is now explicitly indicated. When Eq. (22a) is introduced into Eq. (23a), we find analytically

Eq. (29) and to some extent Eq. (30a) suggests a procedure to obtain the flux for more complicated sources as we now demonstrate.

c /c;n ð0; kÞ ¼ 2

Z

1

1

" dl

# N1 du X c  aa;n ðkÞPa ð2u  1Þ ¼ a0;n ðkÞ: dl a¼0 2

ð24Þ

As it stands, the evaluation of Eq. (23a) requires a two-fold integration. This along with the Fn solution begins to become more than just a routine computation. However, if the source distribution is radially symmetric

Sðqðq; hÞÞ ¼ SðqÞ;

3.1.2. Centered and off-center ring sources For the centered ring source emitting Q 0 particles/s

SðqÞ ¼

dðq  q0 Þ Q 0; 2pq0

ð31aÞ

for which the Hankel transform is

ð25Þ

then

I ð0; qÞ ¼ eD SðqÞHð 1Þ þ

c 4p

Z

1

dkkSðkÞJ 0 ðqkÞa0; ðkÞ;

ð26Þ

where the source transform depends only on the magnitude of k. Hence, the flux is independent of azimuth. Eq. (26) represents the fundamental solution that we use to find a variety of 2 and 3D source configurations in the next section. 3. Theory for the 3D ISLP/HS To arrive at a 3D flux profile at the surface z ¼ 0 from a 3D embedded source, we first consider standard sources like the point, ring and disk. Then, based on these sources, we consider the more sophisticated sources that follow, including the pill, and the hollow pill.

ρ

ρ

0

ρ0 θ

θ0

Fig. 4. Off-center point source.

1247

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

SðkÞ ¼ Q 0 J 0 ðq0 kÞ:

ð31bÞ

Thus, when introduced into Eq. (26) along with Eq. (23b), we find

dðq  q0 Þ c IRing ð0; q; q0 Þ ¼ Q 0 eD Hð 1Þ þ 2pq0 4p Z 1  dkkJ0 ðqkÞJ 0 ðq0 kÞa0; ðkÞ:

ð32Þ

0

An alternative derivation for the ring source is found by integration of the point source of Eq. (30b) over the ring source and then dividing by the source length 2pq0 to give

Z 2p Q0 dh0 q0 Ið0; q; q0 Þ 2pq0 0 Z 1 Q dkkJ0 ðkqÞJ 0 ðkq0 Þ/0; ð0; kÞ: ¼ 0 2p 0

3.2. Vertical line source In this section, we consider a vertical line source as shown in Fig. 5. The source is of length l and is centered at z ¼ D. For this case, we integrate the centered point source over the line source and divide by the source length

ILine ð0; q; D; lÞ ¼

Dþl=2

0

dz Ið0; q; z0 Þ:

ð38Þ

Dl=2

With Eq. (28), we find

Q 0 dðqÞ ðDl=2Þ ½e  eðDþl=2Þ Hð 1Þ l q  Z 1 c þ dkkJ0 ðkqÞv0;1 ðk; DÞ ; 4p 0

ILine ð0; q; D; lÞ ¼

IRing ð0; q; q0 Þ ¼

ð33Þ

Z

Q0 l

ð39aÞ

where

Z

Dþl=2

Trivially, Eq. (33) is applicable for the off-center ring source with the substitution

va;1 ðk; DÞ 

q ! ½q2 þ q20  2qq0 cosðh  h0 Þ1=2 :

Note that the dependence of a0; ðkÞ on z must now been included. The determination of va comes directly from Eq. (22d), since

3.1.3. Centered and off-center disk sources For a centered solid disk of radius q0 at z ¼ D,

SðqÞ ¼

Hðq0  qÞ

pq20

Q0

N1 X

ð34aÞ

Ba ðuÞaa;1 ðk; zÞ ¼ q1 ðn; zÞ;

a¼0

N1 X

2 J ðq kÞ: kq0 1 0

ð34bÞ

IDisk ð0; qÞ ¼ Q 0 e

D

Z

1

Hðq0  qÞ

p20

Hð 1Þ þ

T  ðk; DÞ ¼

c 2pq0

dkJ 0 ðqkÞJ 1 ðq0 kÞa0; ðkÞ:

0

ð35Þ

The identical expression is found by integrating the ring source of Eq. (32) over a differential ring ð2pq00 dq00 Þ to radius q0 and dividing by the area of the disk

IDisk ð0; q; q0 Þ ¼ ¼

Q0

pq20 Q0

Ba ðuÞva;1 ðk; DÞ ¼ T  ðk; DÞ;

Z

q0

0

Z

pq

2 0

dq00 2pq00 ½IRing ð0; q; q00 Þ=Q 0 

1

dkkJ0 ðkqÞ

0

Z 0

q0

Q0

pq0

Z 0

1 ðDl=2Þ ðDþl=2Þ ½e e n½eðDl=2Þ=n eðDþl=2Þ=n  Hð 1Þ 1n n þ ½eðDl=2Þ=n eðDþl=2Þ=n Hð1Þ: 1þn ð40bÞ

3.3. Pill source The creation of a line source leads directly to 3D sources. To establish the flux from a beam pill source as shown in Fig. 6, all we require is the integration of the disk source over the vertical source extent in the form

dq00 q00 J 0 ðkq00 Þ/0; ð0; kÞ IPill ð0; q; q0 ; D; lÞ ¼

giving

IDisk ð0; q; q0 Þ ¼

ð40aÞ

a¼0

where

Eq. (28) becomes



ð39bÞ

then by integrating

for which the Hankel transform is

SðkÞ ¼ Q 0

0

dz aa;1 ðk; z0 Þ:

Dl=2

Q0 l

Z

Dþl=2

0

dz Dl=2



 IDisk ð0; q; z0 Þ Q0

to give

1

dkJ 0 ðkqÞJ 1 ðkq0 Þ/0; ð0; kÞ

ð36Þ

and eventually Eq. (35) with the substitution of Eq. (23b). It is now a relatively simple matter to find a centered disk-shell source by superposition

IDisk-Shell ð0; q; qi ; q0 Þ  2

 2   pqi Q0 pq0 ¼ I I ð0; q ; q Þ  ð0; q ; q Þ : Disk Disk 0 i Q0 pðq20  q2i Þ Q 0

l

Δ

ð37Þ

Here, we have superposed an appropriately normalized negative disk source of radius qi onto a positive disk source of radius q0 and divided by the source area. Both the disk and the disk-shell sources are translated off center with the replacement

q ! ½q2 þ q20  2qq0 cosðh  h0 Þ1=2 : Fig. 5. Vertical line source.

1248

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

4.1. Implementation of the Fn solution We implement the Fn method through two numerical methods. The first is collocation where Eqs. (22d) and (22g) are satisfied at the following N collocation points determined by the zeros of the shifted Legendre polynomial of degree N  1

l

Δ

PN1 ð2ub  1Þ ¼ 0;

b ¼ 1; . . . ; N  1:

ð43Þ

The first collocation point, u0 , derives from Eq. (22i). In this way, we close the system for the determination of the blending coefficients a0; . The resulting set of equations,

Table 1 Choice of a. Fig. 6. Pill source.

Δ

Fig. 7. Pill/shell source.

ð41Þ

We can accommodate a hollow pill source, as shown in Fig. 7 by superposing

IHollowpill ð0; q; q1 ; q2 ; D; l1 ; l2 Þ

 2  Q0 pq2 l2 ¼ IPill ð0; q; q2 ; D; l2 Þ 2 2 Q0 pq2 l2  pq1 l1  2   pq1 l1 IPill ð0; q; q1 ; D; l1 Þ  Q0

a

Point Ring and Disk

For q < q0 ; q0 and q > q0 ; q

Table 2a Nadir point source: ðc ¼ 0:99Þ.

l

IPill ð0; q; q0 ; D; lÞ

Q 1 ðDl=2Þ ½e  eðDþl=2Þ Hðq0  qÞHð 1Þ ¼ 0 l pq20  Z 1 c dkkJ0 ðkqÞJ 1 ðkq0 Þv0;1 ðk; DÞ : þ 2pq0 0

Source

ð42Þ

a negative pill source within. Note that source self-absorption is neglected for all 3D sources. Hence, the fundamental point source flux enables a true 3D source, which should provide a stringent test of any transport algorithm. 4. Numerical implementation and demonstration Two fundamental numerical methods are required for the evaluation of the flux for the standard sources. The first procedure is the solution of Eq. (21) via the Fn spectral representation of Eq. (22). The second is performing the radially symmetric Hankel inversion numerically.

q

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255 N1 X

Ba ðub Þaa; ðkÞ ¼ q ðnb Þ;

b ¼ 0; 1; . . . ; N  1

ð44Þ

Bab 

a¼0

q;b 

are solved by matrix inversion. Note that

ub nb ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2 1þk and that only the zero coefficient a0; is required in the Hankel inversion. The second, and preferred method of solution, is the Galerkin method, where Eq. (22d) is weighted by the same shifted Legendre polynomials of the spectral representation

Pb ð2u  1Þ ¼ 0;

b ¼ 0; . . . ; N  2;

to give the following N  1 linear equations: N1 X

Bab aa; ðkÞ ¼ qb ;

b ¼ 0; 1; . . . ; N  2;

ð45aÞ

Z Z

1

duP b ð2u  1ÞBa ðuÞ;

0

0

1249

ð45bÞ

1

duP b ð2u  1Þq ðuÞ:

The integrals in Eq. (45b) are evaluated by Gauss/Legendre quadrature. With the system of Eq. (45a) closed by Eq. (22g), we again solve the resulting linear set. The second option forces Eq. (22d) to hold true at a greater number of u-points in a collective sense and therefore gives more consistent results. To achieve a desired accuracy, we iterate the Fn solutions by increasing N by 2 in both the Fn approximations of Eqs. (44) and (45a) until convergence. In addition, since this procedure generates a sequence of solutions that supposedly converge to the desired solution, one can apply convergence acceleration to arrive at the limit more quickly. In particular, the Wynn-epsilon (We) (Baker and Graves-Morris, 1996) acceleration algorithm applies to the se-

a¼0

where Table 2b Nadir ring source: q0 ¼ 3 (c = 0.99, D ¼ 1)

Table 2c Nadir disk source: q0 ¼ 2 (c = 0.99, D ¼ 1).

q

Surface flux

0.00E01 1.00E01 2.00E01 3.00E01 4.00E01 5.00E01 6.00E01 7.00E01 8.00E01 9.00E01 1.00E+00 1.10E+00 1.20E+00 1.30E+00 1.40E+00 1.50E+00 1.60E+00 1.70E+00 1.80E+00 1.90E+00 2.00E+00 2.10E+00 2.20E+00 2.30E+00 2.40E+00 2.50E+00 2.60E+00 2.70E+00 2.80E+00 2.90E+00 3.00E+00 3.10E+00 3.20E+00 3.30E+00 3.40E+00 3.50E+00 3.60E+00 3.70E+00 3.80E+00 3.90E+00 4.00E+00 4.10E+00 4.20E+00 4.30E+00 4.40E+00 4.50E+00 4.60E+00 4.70E+00 4.80E+00 4.90E+00 5.00E+00

4.8490E02 4.8353E02 4.8341E02 4.8153E02 4.7888E02 4.7542E02 4.7111E02 4.6591E02 4.5975E02 4.5257E02 4.4427E02 4.3474E02 4.2385E02 4.1145E02 3.9734E02 3.8127E02 3.6292E02 3.4187E02 3.1743E02 2.8825E02 2.4818E02 2.0884E02 1.8130E02 1.5900E02 1.4037E02 1.2456E02 1.1101E02 9.9325E03 8.9178E03 8.0325E03 7.2566E03 6.5738E03 5.9707E03 5.4359E03 4.9603E03 4.5360E03 4.1563E03 3.8156E03 3.5092E03 3.2328E03 2.9831E03 2.7568E03 2.5515E03 2.3648E03 2.1947E03 2.0394E03 1.8974E03 1.7674E03 1.6482E03 1.5386E03 1.4379E03

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

1250

quences for aa; generated by successively increasing N. The We accelerator takes the form for the sequence X n ; n ¼ 0; 1; . . . ; N

eðnÞ 1  0; eðnÞ n ¼ 0; . . . ; N; 0  Xn; e

ðnÞ kþ1

¼e

ðnþ1Þ k1

h

þ

ðnþ1Þ k

e

ðnÞ k

e

i1

eWe ;

k ¼ 0; . . . ; N; ¼ 0; . . . ; N  k  1: ð46aÞ

The recurrence forms a tableau

eð0Þ 0 eð1Þ 0 eð2Þ 0

eð0Þ 1 eð1Þ 1

eð0Þ 2 eð1Þ 2

...

eð0Þ eð0Þ N N1 ð1Þ . . . eN1 ... ...

ð46bÞ

ðN2Þ 2

e

...

where each element of an even column gives an estimate of the limit. Convergence comes from interrogation of last term of the even ðNiÞ ; i ¼ 0; 2; . . . ; 2½N=2 columns ei

eðN1Þ 1

  eðNiÞ  eðNi2Þ    i i   < e;   eðNiÞ

i ¼ 2; . . . ; 2½N=2:

ð46cÞ

i

Hence, we continue the Fn approximation through successively higher orders to form the sequence to be accelerated by the We procedure. X n is therefore each sequence of a0; from the solutions of Eqs. (44) and (45a), respectively. This gives four error measures, including the original sequences and their We equivalents. Convergence occurs when any one of the sequences converge. The minimal penalty in computational time by including both Fn forms and their We equivalents is well worth the confidence in convergence provided. We now proceed to the outer Hankel transform inversion.

eðNÞ 0 Table 3a Zenith point source: ðc ¼ 0:99Þ.

q

1.00E01 2.00E01 3.00E01 4.00E01 5.00E01 6.00E01 7.00E01 8.00E01 9.00E01 1.00E+00 1.10E+00 1.20E+00 1.30E+00 1.40E+00 1.50E+00 1.60E+00 1.70E+00 1.80E+00 1.90E+00 2.00E+00 2.10E+00 2.20E+00 2.30E+00 2.40E+00 2.50E+00 2.60E+00 2.70E+00 2.80E+00 2.90E+00 3.00E+00 3.10E+00 3.20E+00 3.30E+00 3.40E+00 3.50E+00 3.60E+00 3.70E+00 3.80E+00 3.90E+00 4.00E+00 4.10E+00 4.20E+00 4.30E+00 4.40E+00 4.50E+00 4.60E+00 4.70E+00 4.80E+00 4.90E+00 5.00E+00

Table 3b Zenith ring source: (c ¼ 0:99, D ¼ 1).

D 0.1

1

2

4

5.6660E01 3.7927E01 2.6941E01 2.0235E01 1.5823E01 1.2740E01 1.0485E01 8.7778E02 7.4490E02 6.3921E02 5.5361E02 4.8325E02 4.2469E02 3.7542E02 3.3359E02 2.9777E02 2.6690E02 2.4011E02 2.1673E02 1.9624E02 1.7818E02 1.6220E02 1.4801E02 1.3536E02 1.2405E02 1.1391E02 1.0479E02 9.6570E03 8.9136E03 8.2402E03 7.6288E03 7.0726E03 6.5656E03 6.1027E03 5.6792E03 5.2913E03 4.9352E03 4.6081E03 4.3069E03 4.0294E03 3.7734E03 3.5368E03 3.3180E03 3.1153E03 2.9274E03 2.7530E03 2.5909E03 2.4402E03 2.3000E03 2.1692E03

5.5125E02 5.4220E02 5.2784E02 5.0909E02 4.8705E02 4.6277E02 4.3723E02 4.1124E02 3.8543E02 3.6028E02 3.3609E02 3.1308E02 2.9137E02 2.7100E02 2.5198E02 2.3428E02 2.1785E02 2.0262E02 1.8853E02 1.7550E02 1.6345E02 1.5232E02 1.4203E02 1.3252E02 1.2373E02 1.1560E02 1.0808E02 1.0112E02 9.4667E03 8.8688E03 8.3142E03 7.7995E03 7.3215E03 6.8772E03 6.4641E03 6.0796E03 5.7215E03 5.3878E03 5.0766E03 4.7862E03 4.5150E03 4.2616E03 4.0246E03 3.8029E03 3.5954E03 3.4009E03 3.2187E03 3.0477E03 2.8872E03 2.7365E03

2.1498E02 2.1380E02 2.1186E02 2.0920E02 2.0588E02 2.0194E02 1.9748E02 1.9255E02 1.8723E02 1.8161E02 1.7574E02 1.6971E02 1.6357E02 1.5737E02 1.5118E02 1.4502E02 1.3895E02 1.3298E02 1.2714E02 1.2146E02 1.1595E02 1.1062E02 1.0548E02 1.0053E02 9.5775E03 9.1218E03 8.6856E03 8.2687E03 7.8707E03 7.4910E03 7.1293E03 6.7849E03 6.4571E03 6.1455E03 5.8492E03 5.5676E03 5.3002E03 5.0463E03 4.8052E03 4.5763E03 4.3590E03 4.1528E03 3.9570E03 3.7712E03 3.5949E03 3.4275E03 3.2687E03 3.1178E03 2.9746E03 2.8386E03

6.8456E03 6.8326E03 6.8109E03 6.7807E03 6.7423E03 6.6958E03 6.6415E03 6.5797E03 6.5108E03 6.4352E03 6.3534E03 6.2656E03 6.1724E03 6.0743E03 5.9717E03 5.8651E03 5.7549E03 5.6417E03 5.5258E03 5.4078E03 5.2880E03 5.1669E03 5.0448E03 4.9221E03 4.7992E03 4.6764E03 4.5539E03 4.4320E03 4.3111E03 4.1913E03 4.0728E03 3.9559E03 3.8406E03 3.7272E03 3.6157E03 3.5063E03 3.3990E03 3.2940E03 3.1913E03 3.0909E03 2.9930E03 2.8974E03 2.8043E03 2.7136E03 2.6253E03 2.5395E03 2.4561E03 2.3751E03 2.2965E03 2.2202E03

q

Surface flux

1.00E01 1.00E01 2.00E01 3.00E01 4.00E01 5.00E01 6.00E01 7.00E01 8.00E01 9.00E01 1.00E+00 1.10E+00 1.20E+00 1.30E+00 1.40E+00 1.50E+00 1.60E+00 1.70E+00 1.80E+00 1.90E+00 2.00E+00 2.10E+00 2.20E+00 2.30E+00 2.40E+00 2.50E+00 2.60E+00 2.70E+00 2.80E+00 2.90E+00 3.00E+00 3.10E+00 3.20E+00 3.30E+00 3.40E+00 3.50E+00 3.60E+00 3.70E+00 3.80E+00 3.90E+00 4.00E+00 4.10E+00 4.20E+00 4.30E+00 4.40E+00 4.50E+00 4.60E+00 4.70E+00 4.80E+00 4.90E+00 5.00E+00

8.8748E03 8.8748E03 8.8929E03 8.9231E03 8.9654E03 9.0197E03 9.0861E03 9.1645E03 9.2549E03 9.3572E03 9.4712E03 9.5966E03 9.7332E03 9.8803E03 1.0037E02 1.0204E02 1.0377E02 1.0557E02 1.0741E02 1.0927E02 1.1110E02 1.1288E02 1.1454E02 1.1603E02 1.1728E02 1.1821E02 1.1874E02 1.1880E02 1.1831E02 1.1722E02 1.1552E02 1.1322E02 1.1036E02 1.0702E02 1.0328E02 9.9244E03 9.5013E03 9.0676E03 8.6309E03 8.1977E03 7.7729E03 7.3600E03 6.9619E03 6.5801E03 6.2158E03 5.8694E03 5.5411E03 5.2305E03 4.9374E03 4.6610E03 4.4007E03

1251

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

a is chosen according to Table 1.

4.2. Numerical Hankel transform inversion

The inversion integral can be recast as the infinite series By far the most sensitive and challenging numerical evaluation is the Hankel transform inversion given by Eq. (5). Fortunately, we will not require the full inversion over all k-space, as we will need only radially symmetric sources to generate our 3D flux variations as will be shown. Because of radially symmetry, evaluation of Eq. (26) requires the following general Hankel transform:

f ðqÞ ¼

1 2p

Z

dkkJ0 ðqkÞf ðkÞ:

0

1 2pa2

1 2pa2

1 Z X j¼0

ajþ1

duuJ0 ðqu=aÞf ðu=aÞ;

ð49Þ

aj

where the integration is between the zeros aj of the Bessel function J0 . This is the most straightforward and general approach to this evaluation. One could specialize the integral to be between the zeros of the entire integrand (Williams, in press), entailing a zero

1

ð47Þ

A more numerically convenient form is found by introducing an extraneous variable a into the integrand with a simple substitution to give

f ðqÞ ¼

f ðq Þ ¼

Z

1

duuJ0 ðqu=aÞf ðu=aÞ;

ð48Þ

0

Table 3c Zenith disk source: q0 ¼ 2 ðc ¼ 0:99; D ¼ 1Þ.

q

Surface flux

2.00E01 1.00E01 2.00E01 3.00E01 4.00E01 5.00E01 6.00E01 7.00E01 8.00E01 9.00E01 1.00E+00 1.10E+00 1.20E+00 1.30E+00 1.40E+00 1.50E+00 1.60E+00 1.70E+00 1.80E+00 1.90E+00 2.00E+00 2.10E+00 2.20E+00 2.30E+00 2.40E+00 2.50E+00 2.60E+00 2.70E+00 2.80E+00 2.90E+00 3.00E+00 3.10E+00 3.20E+00 3.30E+00 3.40E+00 3.50E+00 3.60E+00 3.70E+00 3.80E+00 3.90E+00 4.00E+00 4.10E+00 4.20E+00 4.30E+00 4.40E+00 4.50E+00 4.60E+00 4.70E+00 4.80E+00 4.90E+00 5.00E+00

2.9471E02 2.9565E02 2.9471E02 2.9315E02 2.9096E02 2.8815E02 2.8471E02 2.8066E02 2.7600E02 2.7074E02 2.6489E02 2.5847E02 2.5151E02 2.4402E02 2.3606E02 2.2767E02 2.1892E02 2.0986E02 2.0059E02 1.9118E02 1.8173E02 1.7234E02 1.6307E02 1.5403E02 1.4525E02 1.3681E02 1.2873E02 1.2104E02 1.1375E02 1.0686E02 1.0037E02 9.4279E03 8.8560E03 8.3202E03 7.8187E03 7.3497E03 6.9111E03 6.5013E03 6.1183E03 5.7604E03 5.4259E03 5.1132E03 4.8209E03 4.5474E03 4.2916E03 4.0522E03 3.8280E03 3.6179E03 3.4211E03 3.2365E03 3.0633E03

Table 4a Benchmark values of transects for the two-shell source. x

y = 4.5

y¼0

4.5000E+00 4.1400E+00 3.7800E+00 3.4200E+00 3.0600E+00 2.7000E+00 2.3400E+00 1.9800E+00 1.6200E+00 1.2600E+00 9.0000E01 5.4000E01 1.8000E01 1.8000E01 5.4000E01 9.0000E01 1.2600E+00 1.6200E+00 1.9800E+00 2.3400E+00 2.7000E+00 3.0600E+00 3.4200E+00 3.7800E+00 4.1400E+00 4.5000E+00

1.8854E03 2.2324E03 2.6415E03 3.1193E03 3.6700E03 4.2933E03 4.9811E03 5.7140E03 6.4598E03 7.1700E03 7.7854E03 8.2425E03 8.4871E03 8.4871E03 8.2425E03 7.7854E03 7.1700E03 6.4598E03 5.7140E03 4.9811E03 4.2933E03 3.6700E03 3.1193E03 2.6415E03 2.2324E03 1.8854E03

8.5185E03 1.2723E02 1.9563E02 2.2574E02 2.2473E02 2.0810E02 2.4168E02 3.5236E02 4.4644E02 4.6109E02 3.7173E02 3.1103E02 2.8671E02 2.8671E02 3.1103E02 3.7173E02 4.6109E02 4.4644E02 3.5236E02 2.4168E02 2.0810E02 2.2473E02 2.2574E02 1.9563E02 1.2723E02 8.5185E03

Table 4b Benchmark values of transects for stacked disks. x

y 4:5

0

4.5000E+00 4.1400E+00 3.7800E+00 3.4200E+00 3.0600E+00 2.7000E+00 2.3400E+00 1.9800E+00 1.6200E+00 1.2600E+00 9.0000E01 5.4000E01 1.8000E01 1.8000E01 5.4000E01 9.0000E01 1.2600E+00 1.6200E+00 1.9800E+00 2.3400E+00 2.7000E+00 3.0600E+00 3.4200E+00 3.7800E+00 4.1400E+00 4.5000E+00

1.1739E03 1.3558E03 1.5635E03 1.7988E03 2.0623E03 2.3531E03 2.6680E03 3.0009E03 3.3422E03 3.6784E03 3.9925E03 4.2650E03 4.4760E03 4.6077E03 4.6476E03 4.5908E03 4.4415E03 4.2119E03 3.9207E03 3.5892E03 3.2386E03 2.8870E03 2.5483E03 2.2319E03 1.9434E03 1.6851E03

3.4701E03 4.4826E03 5.8872E03 7.8866E03 1.0827E02 1.5343E02 2.2757E02 5.8646E02 7.9843E02 1.4997E01 3.1506E01 3.4573E01 3.5945E01 3.6099E01 3.5063E01 3.2430E01 1.6587E01 1.4645E01 1.3290E01 1.1946E01 6.4271E02 5.4811E02 4.8463E02 4.2619E02 1.2659E02 8.1935E03

1252

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

search for each source configuration; however, as we observe, this is unnecessary when convergence acceleration is included. We perform the integration using an iterative shifted Gauss– Legendre quadrature where the order is incremented by two until convergence. Again, the opportunity to use the We algorithm on the sequences of integrals created presents itself. The integration is at the heart of three convergence accelerations applied to the original partial sums. We first accelerate the series through Euler acceleration (Press et al., 1988) for an alternating series. The next acceleration is via the We accelerator. Finally, we accelerate the Euler sequence by a second We acceleration. The series has converged when any of the four error measures converge to a desired accuracy. From experience, one can conclude that convergence acceleration enables the Hankel numerical inversion in a most efficient and cost effective fashion. 4.3. Demonstration of standard sources Our first demonstration concerns the three standard centered source configurations considered by Williams (in press) for a medium of c ¼ 0:99. For these cases, the source emits perpendicularly downward. Tables 2a–2c give the collided fluxes for the above numerical implementation. For the point source, 5-digits (highlighted) out of 196 are in disagreement by one unit in the last decimal with those of Williams (in press). For the other two sources, only one digit (near the ring source) is in disagreement. All values are converged values requiring N less than 60 and under 100 terms in the series for the Hankel transform inversion. It is well known that any numerical method will have difficulty near a source because of its singular nature. The above method is no different in that the highest order Fn solutions are required near the sources. However, it is this author’s contention that the values shown for the point source are true to all digits reported. We continue the demonstration, with a zenith source (pointed directly upward. Tables 3a–3c provide the same information for this source as in Table 2 for the nadir source. The values in Table 3 are expected to be correct to all digits posted. Of course, this is

Table 4c Benchmark values of transects for disk with a hole. x

2.2500E+00 2.0700E+00 1.8900E+00 1.7100E+00 1.5300E+00 1.3500E+00 1.1700E+00 9.9000E01 8.1000E01 6.3000E01 4.5000E01 2.7000E01 9.0000E02 9.0000E02 2.7000E01 4.5000E01 6.3000E01 8.1000E01 9.9000E01 1.1700E+00 1.3500E+00 1.5300E+00 1.7100E+00 1.8900E+00 2.0700E+00 2.2500E+00

y 2.25

0

2.25

6.3516E03 7.1726E03 8.0854E03 9.0890E03 1.0177E02 1.1333E02 1.2531E02 1.3734E02 1.4890E02 1.5936E02 1.6805E02 1.7427E02 1.7749E02 1.7737E02 1.7392E02 1.6748E02 1.5859E02 1.4795E02 1.3625E02 1.2411E02 1.1204E02 1.0043E02 8.9535E03 7.9502E03 7.0398E03 6.2229E03

1.7785E02 2.2964E02 6.1898E02 6.6961E02 7.0760E02 7.3704E02 7.5989E02 7.7735E02 7.9020E02 7.9897E02 8.0397E02 8.0537E02 8.0324E02 7.9768E02 7.8893E02 7.7762E02 7.6490E02 7.5200E02 7.3916E02 7.2494E02 7.0690E02 6.8252E02 6.4914E02 6.0240E02 2.1623E02 1.6698E02

6.2229E03 7.0192E03 7.9023E03 8.8701E03 9.9144E03 1.1019E02 1.2155E02 1.3284E02 1.4353E02 1.5299E02 1.6052E02 1.6546E02 1.6731E02 1.6582E02 1.6109E02 1.5359E02 1.4401E02 1.3315E02 1.2172E02 1.1032E02 9.9342E03 8.9032E03 7.9515E03 7.0839E03 6.3001E03 5.5969E03

a conjecture since there are no benchmarks in the literature for comparison. 4.4. Demonstration of 2/3D sources In this section, we finally consider some 3D sources all emitting nadir. Graphical flux maps for 3D sources constructed from 1D sources through superposition are provided. In addition, Tables 4a–4d of values along 1D transects (slices) are provided for benchmarking purposes. All values are expected to be accurate to one unit in the last place. The benchmark specifications of all cases considered are given in Table 5. 4.4.1. Shell and disk sources We first consider the shell source of Eq. (37). A shell source of outer radius 2 mfp and inner radius 1 mfp is centered 1 mfp above the surface and emits 1particle/s perpendicular to the surface as shown in Fig. 8a. Fig. 9a shows the flux map over a square region on the surface only slightly larger than the source region. The shell source is clearly outlined in the surface response. One can construct a more complicated source configuration simply by superposition of an additional source (see Fig. 8a). Fig. 9b shows the

Table 4d Benchmark values of transects for hollow pill. x

y

2.5000E+00 2.3000E+00 2.1000E+00 1.9000E+00 1.7000E+00 1.5000E+00 1.3000E+00 1.1000E+00 9.0000E01 7.0000E01 5.0000E01 3.0000E01 1.0000E01 1.0000E01 3.0000E01 5.0000E01 7.0000E01 9.0000E01 1.1000E+00 1.3000E+00 1.5000E+00 1.7000E+00 1.9000E+00 2.1000E+00 2.3000E+00 2.5000E+00

2.5

0

2.5

4.3787E03 4.9463E03 5.5767E03 6.2692E03 7.0189E03 7.8151E03 8.6398E03 9.4664E03 1.0260E02 1.0977E02 1.1570E02 1.1995E02 1.2216E02 1.2211E02 1.1982E02 1.1549E02 1.0948E02 1.0224E02 9.4258E03 8.5950E03 7.7673E03 6.9693E03 6.2189E03 5.5266E03 4.8971E03 4.3312E03

1.2242E02 1.5628E02 2.0606E02 6.0218E02 6.5565E02 6.9416E02 7.2327E02 7.4542E02 7.6208E02 7.7416E02 7.8225E02 7.8671E02 7.8777E02 7.8553E02 7.8008E02 7.7159E02 7.6035E02 7.4655E02 7.2986E02 7.0902E02 6.8191E02 6.4550E02 5.9394E02 1.9944E02 1.5096E02 1.1814E02

4.3312E03 4.8894E03 5.5085E03 6.1873E03 6.9204E03 7.6964E03 8.4970E03 9.2946E03 1.0054E02 1.0730E02 1.1278E02 1.1651E02 1.1816E02 1.1755E02 1.1472E02 1.0994E02 1.0363E02 9.6296E03 8.8415E03 8.0405E03 7.2574E03 6.5129E03 5.8192E03 5.1824E03 4.6046E03 4.0848E03

Table 5 Specifications for benchmarks in Table 4. Fig. 9

a

b

c

d

c

0.99 1

0.99 1 3

0.99 2 0.5

1

1

1

1 1 1 1

0.99 1 2 3 1 0.75 0.50 0.50 1 0.50 1

q1 q2 q3 D1 D2 D3 t1 t2 l1 l2

1

1 1.2728 1 0.5

1253

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

original source with an additional outer source ring in region 3 6 q 6 4. Again, the sources are clearly imprinted on the surface. We next consider stacked disk sources as shown in Fig. 8b. In Fig. 9c, we observe a layered configuration of the flux on the surface as each source adds to the other. The next source (Fig. 8c) in a disk with a off centre hole. The corresponding surface flux is shown in Fig. 9d. These sources suggest the power of superposition to generate relatively comprehensive benchmarks. Note that superposition enables 3D flux maps to a relatively high accuracy suitable for benchmarking. 4.4.2. Hollow pill source The final source considered is the pill source with a hollow section as shown in Fig. 8d. The inner off-center cylinder is filled with identical scattering material as the surrounding the source but no source. Again, the indentation in the flux representing the lack of source in the off-center cylinder is readily apparent in Fig. 9e. This is a relatively complex source, which should be a suitable test of any 3D transport algorithm.

Finally, the flux from a general line source in 3D space can also be obtained from the point kernel method, making this approach truly general. Appendix A A.1. Derivation of the Fn equations To derive the Fn equations, we first let l be – l in Eq. (9a), multiply the result by ez=s and integrate over z on ½0; 1Þ, to give

  laðlÞ / ðs; l Þ bðlÞ  s c ¼ laðlÞ/ð0; l; l Þ þ / ðsÞ 2

In this work, 3D surface flux maps for the interior searchlight problem are presented using the Fn method and superposition of a point source kernel. Both the Fn method and superposition have existed for a relatively long time. However, this is the first time they have been merged to provide what should be considered valuable additions to the benchmarking literature. While some may consider a half-space problem rather simplistic, such a problem still provides a level of diagnostics that will benefit any serious transport method developer. This work does not end with this publication. One can relatively easily extend the benchmark sources presented to the slab ISLP. In addition, the list of sources is not particularly limited. Additional benchmarks are relatively easily evaluated through superposition.

a ρ2

t1

t2

dl0

0

n0 /ð0; l0 ; l Þ n s 0

ðA1aÞ

where

Z

1

dze

z=s

/ðz; l Þ:

ðA1bÞ

0

Next, dividing through by the leading factor and integrating over all l yields

KðsÞ/ ðs; l Þ ¼ eD=s

Z

"

1

dl0

0

þ eD=s

Z

#

laðl0 Þ /ð0; l0 ; l Þ 0 0 bðl0 Þ  l aðs l Þ "

1

dl0

0

#

f ðl0 ; l Þ

bðl0 Þ  l aðs l Þ 0

0

KðsÞ  1  cLðsÞ;   Z s 1 dl 1 LðsÞ  : 2 1 bðlÞ s  n

ðA2bÞ ðA2cÞ

ρ1

Δ1 ρ2

Δ2

Δ1

ρ3

Δ3

t2 t3

c

t1

r1

d Δ1

ðA2aÞ

;

with

b ρ1

1

þ f ðl; l ÞeD=s ;

/ ðs; l Þ 

5. Discussion

Z

l2

r2

l1

t1

Fig. 8. (a) Two shell sources. (b) Three disk sources. (c) Three disk source with a hole. (d) Hollow pill source.

Δ1

1254

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

a

b

0.05

0.05

0.04

0.04

0.03

2

0.02

Flux

Flux

0.03 1

0.01 2

0.00

0

x

0

2

0.01

x

0.00

4

0.02

4 2

-1

1

-2 0

0

y

y

-2

-1

-2

-2

c

-4

-4

d

0.4

0.3

0.2

0.10

0.0

0.08 2

0.06

Flux

4

0.04 0.02 0.00

2 0

1 -2

0

y

-1

-2 -4 -4

x

4

2

0

-2

-1

0 1 2

x

-2

e

0.10 0.08

Flux

0.06 0.04 0.02 0.00

-2 -1

x

y

Flux

0.1

2 1

0 0

1

-1 2

y

-2

Fig. 9. (a) Shell source. (b) Two shell source. (c) Stacked disk sources. (d) Disk source with a hole. (e) Hollow pill source.

1255

B.D. Ganapol, D.E. Kornreich / Annals of Nuclear Energy 36 (2009) 1242–1255

l

If we introduce

n

laðlÞ : bðlÞ

ðA3Þ

2

bðlÞ ¼ 1 þ k

n  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2 1 þ k l2

l2 ;

gives

into Eq. (A2a), one arrives at the desired departure point for the application of the Fn method

KðsÞ/ ðs; l Þ ¼ s

Z

1

0

n0 /ð0; l0 ; l Þ  seD=s dl0 0 n s

Z

1

1

dl0 bðl0 Þ

f ðl0 ; l Þ :  n0  s

Z

cs KðsÞ  1  2

"

1

dl 1

#

1

sð1 þ k

2

l2 Þ  lð1 þ k2 l2 Þ1=2

or

ðA4Þ

  dl 1 1 þ sþn 0 bðlÞ s  n Z 1 cs2 dl ¼1þ : 2 2 1  k s2 0 l2  1ks 2 s2

KðsÞ  1 

Appendix B

Z

cs 2

1

Therefore, if

B.1. Approach to the cut

s

ffi; xðsÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2

The cut, as shown in Fig. 3, is necessary to keep KðsÞ

cs 2

KðsÞ  1 þ

Z

1

1

 dl 1 ; bðlÞ n  s

1  k s2



ðB1Þ

single valued in the s-plane. The function admits the following boundary values when s  m  ie and the limit as e approached zero is taken:

cm K ðmÞ  lim Kðm  ieÞ ¼ 1 þ e!0 2 

Z



1

1 1 dl bðlÞ n  ðm  ieÞ 1

 :

ðB2Þ

lim e!0



1

P

n  ðm  ieÞ

1  ipdðn  mÞ; nm

ðB3Þ

therefore

Z

cm K ðmÞ ¼ 1 þ 2 

1

1

dl 1 ipcm  P bðlÞ n  m 2

Z

1

1

dl dl dðl  lðmÞÞ; bðlÞ dn ðB4Þ

ðB8Þ

  cm 1 þ xðmÞ ln ; 1  xðmÞ 2

ðB9Þ

and

kðmÞ ¼ 1 

m xðmÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 1  k m2 Appendix C C.1. Determination of the matrix elements

e a ðuÞ  c B 2

lðmÞaðlðmÞÞ m : bðlðmÞÞ Z

cm K ðmÞ ¼ 1  2 

1

1

 dl 1 ipcm dl ; P  bðlÞ m  n 2bðlðmÞÞ dn lðmÞ

Z

0

1

dl0 P

0

du n0 P a ð2u0  1Þ; dl0 n0  n

ðC1Þ

we can write

After integration over the delta function

ðB5Þ

e a ðuÞ  c B 2

Z

1

0

du Pa ð2u0  1Þ þ u

0

Z 0

1

0

du P

 Pa ð2u0  1Þ ; u0  n

ðC2Þ

to give

with

e a ðuÞ ¼ c ½da;0  2uQ a ð2u  1Þ: B 2

dl 1 ; ¼ dn ð1  k2 n2 Þ3=2 when we introduce aðlÞ and bðlÞ from below Eq. (13). The penultimate form for the boundary values of K is

ipcm 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; K ðmÞ ¼ kðmÞ  2 2 1  k n2 

ðB6aÞ

with

cm 2

Z

1

1

dl 1 P : bðlÞ m  n

ðB6bÞ

The final form comes when kðmÞ is analytically determined. Letting Eq. (B1) be

KðsÞ  1 

  cs xþ1 ln 2 x1

Since

where

kðmÞ ¼ 1 

KðsÞ  1 

with

From the Plemelj relations (Roos, 1969)



we have

cs 2

Z

1

1

and substituting

  dl 1 bðlÞ s  n

ðB7Þ

ðC3Þ

References Baker, G.A., Graves-Morris, P., 1996. Pade’ Approximants. Cambridge University Press, NY. Barichello, L., Siewert, C.E., 2000. The searchlight problem for radiative transfer in a finite Slab. J. Comp. Phys. 157, 707–726. Elliott, J.P., 1955. Milne’s problem with a point source. Proc. Roy. Soc. London, Ser. A 228 (1174), 424–433. Press, W.H. et al., 1988. Numerical Recipes. Cambridge University Press, NY. Roos, B.W., 1969. Analytic Functions and Distributions in Physics and Engineering. John Wiley & Sons, NY. Rybicki, G.B., 1965. Radiative Transfer in Stochastic Media. Smithsonian Institute Astrological Observatory, No. 180. Siewert, C.E., Dunn, W.L., 1983. Radiation transport in plane-parallel media with non-uniform surface illumination. ZAMP 34, 627. Siewert, C.E., Dunn, W.L., 1989. The searchlight problem in radiative transfer. JQSRT 41, 467. Williams, M.M.R., 2007. The transport and diffusion theory of a line source in an infinite half-space with internal reflection. J. Math. Phys. Gen. 34, 910. Williams, M.M.R., in press. Three-dimensional transport theory: an analytical solution for the internal beam searchlight problem I. Annal. Nucl. Energy.