Three-dimensional transport theory: An analytical solution of an internal beam searchlight problem-I

Annals of Nuclear Energy 36 (2009) 767–783 Contents lists available at ScienceDirect Annals of Nuclear Energy journal homepage: www.elsevier.com/loc...

Download PDF

458KB Sizes 1 Downloads 117 Views

Report

PDF Reader
Full Text

Annals of Nuclear Energy 36 (2009) 767–783

Contents lists available at ScienceDirect

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

Three-dimensional transport theory: An analytical solution of an internal beam searchlight problem-I M.M.R. Williams * Computational Physics and Geophysics, Department of Earth Science and Engineering, Imperial College of Science, Technology and Medicine, Prince Consort Road, London SW7 2BP, UK

a r t i c l e

i n f o

Article history: Received 12 September 2008 Received in revised form 4 January 2009 Accepted 1 February 2009 Available online 3 March 2009

a b s t r a c t We describe a number of methods for obtaining analytical solutions and numerical results for threedimensional one-speed neutron transport problems in a half-space containing a variety of source shapes which emit neutrons mono-directionally. For example, we consider an off-centre point source, a ring source and a disk source, or any combination of these, and calculate the surface scalar ﬂux as a function of the radial and angular co-ordinates. Fourier transforms in the transverse directions are used and a Laplace transform in the axial direction. This enables the Wiener–Hopf method to be employed, followed by an inverse Fourier–Hankel transform. Some additional transformations are introduced which enable the inverse Hankel transforms involving Bessel functions to be evaluated numerically more efﬁciently. A hybrid diffusion theory method is also described which is shown to be a useful guide to the general behaviour of the solutions of the transport equation. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction One of the most challenging remaining problems in analytical solutions of the transport equation must be that concerning multi-dimensional geometry. When the transport equation involves more than one space variable there is an increase in complexity due to interaction of the effects caused by these additional variables. The earliest exact solution known to the author on threedimensional transport theory is due to Elliott (1952, 1955). This work describes the problem of a point source on the surface of a half-space; i.e. an early form of the so-called searchlight problem. Later work involved a variety of source conﬁgurations and even the quarter space problem, although admittedly this did not lead to a closed form solution (Lam and Leonard, 1971). Multidimensional solutions are notoriously difﬁcult to obtain and the usual procedure is to take Fourier transforms in the two, inﬁnite transverse directions, and a Laplace transform in the semi-inﬁnite direction. This leads to a solution either by the Wiener–Hopf technique or as the solution of a singular integral equation using the FN, or other methods (Ganapol, 2008a,b; Naz and Loyalka, 2008). The main practical problem is the inversion of these transforms. Whilst the Wiener–Hopf method generally enables one of the inversions to be done analytically, the algebraic complexity of the Fourier transform usually prevents any neat inversion, except in one speciﬁc case (Williams, 1982). The quest for such solutions is well worth * Address: 2a Lytchgate Close, South Croydon, Surrey CR2 0DX, UK. Tel.: +44 208 680 4450. E-mail address: [email protected] 0306-4549/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2009.02.001

while however since they form the basis of essential benchmark problems for validating large transport theory computer codes. The techniques and philosophy of such techniques and associated numerical methods has been discussed by Ganapol (2008a). We should also mention that problems involving spherical and cylindrical bodies have been studied, details of which may be found in Cassell and Williams (2007a,b). In pursuit of the above goals, Ganapol and Kornreich (2008) have recently developed a technique based on the FN method for solving one-speed transport problems involving interior monodirectional beam sources. Such sources cover point, multi-ring and disk geometries with the particles directed towards the surface. The co-ordinate system used was cylindrical and so the problem is essentially a two-dimensional one. The radial co-ordinate is transformed out via a Hankel transform and the resulting pseudoone-dimensional equation in the axial direction is treated by the FN method. In the present work we wish to extend the scope of the Ganapol–Kornreich work to three-dimensional problems and also to introduce some useful ideas for reducing the numerical problems which arise due to the presence of Bessel functions in the inverse Hankel transform. Some numerical results are given but the emphasis is on technique. A useful hybrid-diffusion theory method is devised which enables the directional effects of the source to be included as a ﬁrst collision source in the classic diffusion equation. We also note that Siewert and co-workers and Crosbie and co-workers have made signiﬁcant contributions to the searchlight problem and refer the reader to Williams (2007a) for references.

768

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

In arriving at Eq. (8), we have noted that

2. General theory The transport equation to be studied may be written in the form Davison (1957)

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @Iðx; y; z; l; uÞ @ @ l þ 1 l2 cos u þ sin u Iðx; y; z; l; uÞ @z @x @y c I0 ðx; y; zÞ þ Sðx; y; zÞdðl l0 Þdðu u0 Þ ð1Þ þ Iðx; y; z; l; uÞ ¼ 4p where I(x, y, z, l, /) is the angular ﬂux of particles (neutrons or photons), the scattering is isotropic and we have scaled distances in units of the mean free path. The scalar ﬂux in the above equation is deﬁned as

I0 ðx; y; zÞ ¼

Z 2p

du

Z

0

1

dlIðx; y; z; l; uÞ

1

The source term is an internal beam of spatial distribution S(x, y, z) emitting particles from each point in the source in the direction deﬁned by the unit vector X0(l0, /0). We further note that

Sðx; y; zÞ ¼ Sðx; yÞdðz z0 Þ

ð2Þ

This denotes a plane source at z = z0 with lateral distribution S(x, y). In addition the boundary condition of no incident particles can be written as I(x, y, 0, l, /) = 0 for l > 0 and /(0, 2p). We further deﬁne the Fourier transform

Z

Iðk; z; l; uÞ ¼

1

dx

Z

1

1

dye

ikx xiky y

Iðx; y; z; l; uÞ

ð3Þ

1

@Iðk; z; l; uÞ þ uðk; l; uÞIðk; z; l; uÞ @z c I0 ðk; zÞ þ SðkÞdðz z0 Þdðl l0 Þdðu u0 Þ ¼ 4p

¼

ð4Þ

Iðk; z; l; uÞ ¼

4pl

0

dz I0 ðk; z0 Þe

0 u0 ðzz0 Þ=l0

dðu u0 Þe

l

þ

ð5Þ

SðkÞ

l0

dðl l0 Þ

Hðz z0 ÞHðl0 Þ

ð6Þ

Iðk; z; l; uÞ ¼ c 4pl

1

0

0

dz I0 ðk; z0 Þeuðz zÞ=l

z

l

0

dl

l

dueuz=l

Z qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ez=l J 0 ðkz 1 l2 =lÞ ¼

1

0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1þk2 w2 z=w dwe pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 w 1 þ k w2

Z

1

ð9Þ

dz0 epz0Iðk; z; l; ujz0 Þ

ð10Þ

0

where we have introduced z0 as a conditional variable in I. Thus Eq. (4) becomes

l

@~Iðk; z; l; ujpÞ þ uðk; l; uÞ~Iðk; z; l; ujpÞ @z c ~ ¼ I0 ðk; zjpÞ þ SðkÞepz dðl l0 Þdðu u0 Þ 4p

ð11Þ

Now introduce the Laplace transform

^~ Iðk; s; l; ujpÞ ¼

Z

1

dze

sz~

Iðk; z; l; ujpÞ

ð12Þ

0

and Eq. (11) becomes

^ l~Iðk; 0; l; ujpÞ þ ðsl þ uÞ~Iðk; s; l; ujpÞ SðkÞ c ~^ I0 ðsjpÞ þ ¼ dðl l0 Þdðu u0 Þ 4p pþs

Z

Z 2p

^ l~Iðk; 0; l; ujpÞ þ Vðs; kÞ~I0 ðk; sjpÞ sl þ u 0 1 Z 2p SðkÞ Z 0 Þdðu u0 Þ dðl þ l dl du ¼ sl þ u p þ s 1 0 0

dl

ð13Þ

du

SðkÞ

l0

ð14Þ

c 4p

du

Z

1

0

1

ð15Þ

The term on the right hand side of Eq. (14) may now be simpliﬁed to

ð16Þ

Let us now deﬁne

dðl l0 Þ

dðu u0 Þeu0 ðz0 zÞ=l0 Hðz0 zÞHðl0 Þ

Z 2p

dl sl þ uðk; l; uÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ# " 2 c 1 þ s2 k pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ log 2 2 2 s2 k 1 s2 k

Vðs; kÞ ¼ 1

SðkÞ 1 p þ s u0 sl

where we have used the boundary condition noted above. For l < 0

Z

1

~Iðk; z; l; ujpÞ ¼

We may now integrate Eq. (4) with respect to z to get for

uðzz0 Þ=

Z 2p

We shall return to Eq. (8) later but for now, in order to solve Eq. (4), we use the Wiener–Hopf technique (Williams, 1971) and deﬁne the Laplace transform

l>0 z

dl

where

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ uðk; l; uÞ ¼ 1 þ i 1 l2 ðkx cos u þ ky sin uÞ

Z

Z

1

0

0

with

c

Z

1 2p

Dividing by (sl + u) and integrating over l(1, 1) and /(0, 2p), we then have

where k = (kx, ky). Applying this transform to Eq. (1) we ﬁnd

l

Kðz; kÞ ¼

gðk; sjpÞ ¼

Z

0

1

ð7Þ

dl

Z 2p 0

du

l~Iðk; 0; l; ujpÞ sl þ u

ð17Þ

and write Eq. (14) as

In theqﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ above,ﬃ H(x) is the Heaviside unit function and u0 ¼ 1 þ i 1 l20 ðkx cos u0 þ ky sin u0 Þ. We may integrate Eqs.

SðkÞ ^ ðp þ sÞ~I0 ðk; sjpÞVðs; kÞ ¼ ðp þ sÞgðk; sjpÞ þ u0 sl

(6) and (7) over l and / to get the following integral equation for the scalar ﬂux I0 ðk; zÞ,

Now using the Wiener–Hopf factorisation (Williams, 2007a) we may write

I0 ðk; zÞ ¼ c 2

Z 0

1

0

dz I0 ðk; z0 ÞKðjz z0 j; kÞ þ

SðkÞ eu0 ðz0 zÞ=l Hðz0 zÞ l ð8Þ

In the above equation we have assumed that all particles from the source are travelling towards the surface, i.e. in the negative . l direction, and have accordingly set l0 ¼ l

sðs; kÞ ¼

s2 1 k

2 2

s2 m 2 k

Vðs; kÞ ¼

ð18Þ

sþ ðs; kÞ s ðs; kÞ

where s±(s; k) are deﬁned in Williams (2007a) and m is the root of the equation

1¼

c 1þm log 2m 1m

ð19Þ

769

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

Rearranging Eq. (18) we ﬁnd

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 ^~ s þ m2 þ k 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ðp þ sÞI0 ðk; sjpÞ 2 s þ 1 þ k s ðs; kÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 SðkÞ s 1þk 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ ðp þ sÞgðk; sjpÞ þ 2 u0 sl s m2 þ k sþ ðs; kÞ

Thus we may write the Laplace inverse of Eq. (27) as

SðkÞ 1 ðp þ sÞgðk; sjpÞ þ ¼ C0 sþ ðs; kÞ u0 sl

L

epz0 1 SðkÞ u0 z0 =l H ;k 1 þ e p u0 þ pl l

More concisely,

Z z0 SðkÞ 0 0 0 I0 ðk; 0jz0 Þ ¼ SðkÞ dz0 Gðkjz Þeðz0 z0 Þu0 =l þ ez0 u0 =l 0 l 0 l

1 2p

ð21Þ

2 Z

1

dkx

Z

1

1

dky eikx xþiky yI0 ðk; 0jz0 Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ SðkÞ p þ 1 þ k2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ s ðp; kÞ p þ m2 þ k2 u0 þ pl

The second term in Eq. (32) is written

ð22Þ

Iu ðx; y; 0jz0 Þ 2 Z 1 Z 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 2 dkx dky SðkÞeikx xþiky yi 1l z0 ðkx cos u0 þky sin u0 Þ=l ¼ ez0 =l 2p 1 1 ð34Þ

ð23Þ

where we have used the relation s+(p; k) = 1/s(p; k). Inserting C0 into Eq. (22) we ﬁnd

SðkÞ ^~ 1 1 1 I0 ðk; sjpÞ ¼ H ;k H ;k s p p þ s u0 þ pl

ð24Þ

ð25Þ

is a generalised Chandrasekhar H-function (Williams, 1982). We are interested in the value of the ﬂux at the surface z = 0. This is readily obtained from Eqs. (12) and (24) by noting that

~I ðk; 0jpÞ ¼ lim s^~I ðk; sjpÞ ¼ 0 0 s!1

SðkÞ 1 H ;k p u0 þ pl

Iu ðx; y; 0jz0 Þ ¼ ez0 =l

1

; u0 Þ dy0 Sðx0 ; y0 ÞGu ðx x0 ; y y0 jz0 ; l

1

Physically, Eq. (35) denotes the un-collided ﬂux and describes the fact that the particle travels in a straight line ﬃdeﬁned by the pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 co-ordinates pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ on the surface viz: x z0 1 l cos u0 =l and . If S(x0, y0) = d(x0)d(y0), then these lines 2 sin u0 =l y z0 1 l emanate from the point (0, 0, z0). 2.1. Radial source

ð27Þ

There p is ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ an alternative representation of Eq. (35) if 2 þ y2 Þ, i.e. depends on the radial co-ordinate Sðx; yÞ ¼ Sð x pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ q ¼ x2 þ y2 only. We may then write

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 cos u0 =l R cos r 1l qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 sin u0 =l R sin r ay ¼ y z0 1 l kx ¼ k cos g;

ð28Þ

In fact this function, when divided by c, corresponds to Elliott’s point source problem (Elliott, 1952, 1955) and is the solution of the pseudo-problem, obtained by integration of Eqs. (6) and (7), viz:

ð29Þ

ky ¼ k sin g

1

ð30Þ

1

In fact, the solution of Eq. (29) obeys the reciprocity relation W0(k, z|z0) = W0(k, z0|z) .

ð37Þ ð38Þ ð39Þ

Hence

kx ax þ ky ay ¼ kR cosðg rÞ 2

where R ¼

a2x

þ

a2y

2

and k ¼

ð40Þ 2 kx

þ

2 ky .

Moreover, if we set

x ¼ q cos #

ð41Þ

y ¼ q sin #

ð42Þ

it is readily seen that

; # u0 Þ R2 ¼ R2 ðq; z0 ; l qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 cosð# u0 Þ=l þ ð1 l 2 Þz20 =l 2 ¼ q2 2qz0 1 l

where W(k, 0, w|z0) = 0 for w > 0 and

dwWðk; 0; wjz0 Þ

Z

where

1 ~ GðkjpÞ ¼ H ;k 1 p

Z

dx0

ð35Þ

ax ¼ x z0

Gðkjz 0 Þ ¼ W0 ðk; 0jz0 Þ ¼

1

1

As we will see, the second term above corresponds to the uncollided ﬂux from the source. The ‘1’ in the ﬁrst term ensures convergence as |p| ? 1. We shall also deﬁne a Green’s function as

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @ Wðk; z; wjz0 Þ 2 2 w 1 þ k w2 þ ð1 þ k w2 ÞWðk; z; wjz0 Þ @z c ¼ ðW0 ðk; zjz0 Þ þ dðz z0 ÞÞ 2

Z

ð26Þ

It is convenient to re-write Eq. (26) as

~I ðk; 0jpÞ ¼ SðkÞ H 1 ; k 1 þ SðkÞ 0 p u0 þ pl u0 þ pl

We may now use the convolution theorem for Fourier transforms and get

; u0 Þ Gu ðx; yjz0 ; l qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ d y z0 1 l ð36Þ 2 cos u0 =l 2 sin u0 =l ¼ d x z0 1 l

where

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 sþ 1þk pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ s ðs; kÞ H ;k ¼ 2 s s þ m2 þ k

ð33Þ

1

To ﬁnd C0 we set s = p in Eq. (21), leading to

C0 ¼

ð32Þ

In order to regain the dependence on x and y, we now invert the Fourier transform, viz:

I0 ðx; y; 0jz0 Þ ¼

and

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 ^~ C0 s þ 1 þ k pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ s ðs; kÞ I0 ðk; sjpÞ ¼ p þ s s þ m2 þ k2

dp

ð31Þ ð20Þ

Each side of Eq. (20) is analytic in an overlapping half-plane, which is a necessary condition for application of the Wiener–Hopf method (Noble, 1958). Now using Liouville’s theorem, we note that as |s| ? 1, each side of Eq. (20) tends to a constant C0. We have therefore

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 s 1þk pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 s m2 þ k

Z

I0 ðk; 0jz0 Þ ¼ SðkÞ 1 2p i

ð43Þ

We then write Eq. (34) as

ez0 =l

1 2p

2 Z

1

1

dkx

Z

1 1

dky SðkÞeikR cosðgrÞ

ð44Þ

770

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

or more concisely

ez0 =l

1 2p

Z

1

0

(c) A Gaussian distribution

dkkJ0 ðkRÞSðkÞ

SðqÞ ¼

ð45Þ

and

and the angular difwhich shows that the ﬂux depends on q; z0 ; l ference # u0. The Fourier–Hankel transform of S(q) is

SðkÞ ¼ 2p

Z

dqqJ 0 ðkqÞSðqÞ

ð46Þ

and we shall normalise the source such that

Sð0Þ ¼ 2p

Z

1

dqqSðqÞ ¼ 1

ð47Þ

0

If we have a point source such that

SðqÞ ¼

SðkÞ ¼ ek2 =4a2

1

0

dðqÞ 2pq

ð48Þ

then SðkÞ ¼ 1, i.e. a point source, emitting particles in the direction ; u0 Þ, which leads to an un-collided ﬂux of X0 ðl

; u0 Þ ¼ ez0 =l Gu ðqjz0 ; l

dðRÞ 2pR

ð49Þ

We may now use Eq. (45) to write the ﬁrst term in Eq. (32) as

I0 ðq; #; 0jz0 Þ ¼

1 l

Z

z0

0

0

0

dz0 eðz0 z0 Þ=l

1 2p

Z

1 0

0 SðkÞGðkjz dkkJ0 ðkRÞ 0Þ

a2 a2 q2 e p ð56Þ

We note that should S(x, y) = S(q, #), then the above theory must be modiﬁed. An example of this would be a disk source with a segment cut out, or a ‘broken’ ring. It should be re-emphasised that the particles from these source shapes are emitted in the direction ; u0 Þ. If we had an off-centre point source then S(q, #) = X0 ðl ; u0 Þ is much d(q q0)d(# #0) and the intensity I0 ðq; #jq0 ; #0 ; z0 ; l more complicated; its derivation is given in Appendix. However, ¼ 1, the ﬂux assumes the same form as that for the simple case of l qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ in Eq. (51) but with q replaced by

q2 þ q20 2qq0 cosð# #0 Þ;

this is obvious from geometric considerations. The other limiting ¼ 0. Both cases are given in Appendix case of interest is when l A and we shall examine the behaviour of these functions below. Finally, we observe that the ﬂux arising from an array of point sources around a circle of radius q0 can be obtained by summing the solutions for each desired value of #0. Similarly, one can calculate the ﬂux due to point sources at different radii, q0, for a given #0, or indeed a combination of both. It is also possible to consider sources at several values of z0 and sum the results.

ð50Þ z00 .

where R is given by Eq. (43) with z0 ! z0 It is convenient to set u0 = 0, as only the difference between this and the surface angular co-ordinate is of signiﬁcance, except if there is a varying emission distribution in the azimuthal direction; a situation that will be discussed below. In the special case of a beam normal to the surface and along the ¼ 1, we have z-axis, in which l

I0 ðq; 0jz0 Þ ¼

Z

z0

0

0

0

dz0 eðz0 z0 Þ

1 2p

Z 0

1

0 dkkJ0 ðkqÞGðkjz 0Þ

ð51Þ

which is independent of # as expected. The corresponding uncollided ﬂux at the surface is

Gu ðqjz0 Þ ¼ ez0

dðqÞ 2pq

ð52Þ

It should be noted that this is a line source of variable strength depending on the distance of the source from the surface. We shall see below that this line source leads to a singularity of the form O(1/q) as q ? 0 in the collided ﬂux. Details of this calculation may be found in Appendix C. The spatial form of the source is of some interest and Ganapol and Kornreich (2008) have considered a disk source, a set of concentric annular rings and a Gaussian distribution. These may be written mathematically as (a) N-ring sources distances qi from the axis, i.e.

SðqÞ ¼

N X i¼1

ai

dðq qi Þ 2p q

N X

ai J 0 ðkqi Þ

; u0 Þ ¼ ez0 =l Gu ðqjz0 ; l

dðR q0 Þ 2pq0

ð57Þ

For the disk source of Eq. (55), we must evaluate

; u0 Þ ¼ ez0 =l Gu ðqjz0 ; l

1 pa

Z

1

0

dkJ 1 ðkaÞJ 0 ðkRÞ

ð58Þ

This is a standard integral and leads to

; u0 Þ ¼ ez0 =l Gu ðqjz0 ; l ¼ 0;

1

pa2

R
;

ð59Þ

R>a

Finally, the Gaussian source leads to

; u0 Þ ¼ Gu ðqjz0 ; l

a2 z0 =l a2 R2 e e p

ð60Þ

where in each case R is given by Eq. (43). 2.3. Radial source emission One situation of particular interest is if the source emits parti ¼ 0. To deal cles in the radial direction only, i.e. the case where l with this, we change variables in Eq. (50) such that , whence w ¼ ðz0 z00 Þ=l

ð54Þ

i¼1

I0 ðq; #; 0jz0 Þ ¼

1 2p

Z

z0 =l

dwe 0

w

Z

1

^ SðkÞGðkjz dkkJ 0 ðkRÞ 0 lwÞ

0

ð61Þ

ai being speciﬁed weight functions. (b) Uniform disk of radius q0 and strength 1=pq20

SðkÞ ¼ 2 J ðkq Þ 0 kq0 1

In Eq. (45), we have obtained an explicit expression for the ﬂux of particles that have undergone no collisions. It is of some interest to obtain the associated analytical expressions for this ﬂux corresponding to the point, ring and disk sources described above. In Eq. (49) we have the point source case and it is clear that for a ring source for which SðkÞ is given by Eq. (54), we have

ð53Þ

for which

SðkÞ ¼

2.2. The un-collided ﬂux distribution

where

ð55Þ

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ^ 2 ¼ q2 2qw 1 l 2 cosð# u0 Þ þ ð1 l 2 Þw2 R

ð62Þ

771

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

and we have suppressed the dependence on /0 in Eq. (61). For l ¼ 1, we recover Eq. (51) and in the limit l ¼ 0, Eq. (62) becomes

1 I0 ðq; #; 0jz0 Þ ¼ 2p

Z

1

dwe

w

Z

0

1

~ SðkÞGðkjz dkkJ 0 ðkRÞ 0Þ

ð63Þ

0

with

~ 2 ¼ q2 2qw cosð# u Þ þ w2 R 0

ð64Þ

It is tempting to try to perform the integral

Z

1

dwe

w

~ J 0 ðkRðwÞÞ

0

in Eq. (63) but in practice this does not help. Interestingly, one could also simulate a broken ring by deﬁning the radial emission source Q(/0) as

1 ; 0 < u0 < 2p w 2p ¼ 0; 2p w < u0 < 2p

Qðu0 Þ ¼

ð65Þ

where w is the angular range of the cut segment. Thus we would have to calculate

Z 2p 0

du0 Q ðu0 ÞI0 ðq; #; 0ju0 ; z0 Þ

ð66Þ

where I0 is given by Eq. (61) and we have explicitly shown the ~ from Eq. (64). dependence on /0 via R 3. The Green’s function As seen above, the crucial quantity required to determine the surface ﬂux is the Green’s function Gðkjz0 Þ as deﬁned by Eq. (28). Thus we require

1 Gðkjz 0Þ ¼ 2p i

Z

pz0

dpe

L

1 H ;k 1 p

Z

1 2p

0

1

0 ðkjz0 Þ dkkJ0 ðkqÞSðkÞG

Eq. (68) may also be used in Eq. (63) to obtain the radial emission solution. We shall discuss numerical values of I0(q, 0|z0) below, however it is important to note that, for q ? 0, the integrand of Eq. (70) for SðkÞ ¼ 1 leads to a singularity of the form 1/q. This feature is implicit in the work of Davison (1957) in his study of the inﬁnite line source and re-appears in the present half-space problem. It is due to the effective line source resulting from the ﬁrst collision particle ﬂux along the z-axis as in Eq. (32). Mathematically, the singularity arises from the integral term in Eq. (71) so does not occur for diffusion theory. Full details are in Appendix C where it is also shown that the behaviour of the ﬂux for large q is independent of the source shape S(q). 4. The non-absorbing case The non-absorbing case of c = 1, m = 0, has a number of interesting features. An important one is the fact that it leads to a useful analytic solution for large values of q. This type of solution was ﬁrst pointed out by Elliott (1952, 1955). In order to arrive at the appropriate equation we must take the limit carefully as c ? 1 by using the approximation from Eq. (19), viz: m2 3(1 c). Then from (74) we ﬁnd

1 0 ðkjz0 Þ ¼ 3 G ðekz0 ez0 Þ 2k Hð1=k; kÞð1 kÞ Z 1 1 dtgð1; tÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ þ 2 2 2 2 0 1 þ k t2 H t= 1 þ k t 2 ; k t 1 þ k t2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 e 1þk t z0 =t ez0 ð74Þ We now seek an expansion of this transform in powers of k. Using the properties of the H(1/k; k) function, we ﬁnd that the ﬁrst term on the right hand side of (74) becomes

asy ðkjz0 Þ ¼ G

ð68Þ Z0 ¼ 1 a2 ¼ ð69Þ

The function H(1/p; k) is deﬁned in the Appendix of Williams (2007a). We may now insert Eq. (68) into Eq. (51), to get

I0 ðq; 0jz0 Þ ¼

ð70Þ

where

pﬃﬃﬃ 3ðekðz0 þZ0 Þ 1 2 3 4 ez0 kZ0 Þ 1 þ k þ k þ ð7 6a2 Þk þ Oðk Þ 6

1

p

Z

1 2p

Z

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 0 ðkjz0 Þ ¼ A0 e m2 þk2 z0 ez0 G þ

c 2

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 2 2 dtAðtÞ e 1þk t z0 =t ez0

ð75Þ

ð71Þ

0

m2 ð1 m2 Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð72Þ A0 ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 m2 þ k H 1= m2 þ k2 ; k ðc þ m2 1Þ 1 m2 þ k2

1

dw tan1 KðwÞ ¼ 0:71044609 . . .

ð76Þ

0 1

2

dww tan1 KðwÞ ¼ 0:0769847 . . .

ð77Þ

0

with

KðwÞ ¼

pw=2

1 w2 log 1þw 1w

ð78Þ

In order to evaluate the quadratures in Eqs. (76) and (77) we ﬁnd it useful to write

tan1 KðwÞ ¼

Z

ð73Þ

where

where

2 2 1 ct 1þt c pt þ ¼ 1 log gðc; tÞ 2 1t 2

gðc; tÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ AðtÞ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 2 1 þ k t 2 H t= 1 þ k t 2 ; k t 1 þ k t2

ð67Þ

As noted above, the ‘1’ is required to ensure convergence of the integrand as |p| ? 1. It also conﬁrms the result G ¼ 0 when c = 0. Following the methods discussed in Williams (1982, 2007a,b) we ﬁnd

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 m2 ð1 m2 Þ 2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e m þk z0 Gðkjz 0 Þ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 2 2 2 m þ k H 1= m þ k ; k ðc þ m 1Þ Z pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 c 1 dtgðc; tÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e 1þk t z0 =t þ 2 0 t 1 þ k2 t 2 H t= 1 þ k2 t 2 ; k

and

p 2

tan1

1 KðwÞ

Inserting the series (75) into Eq. (70) with SðkÞ ¼ 1, and using the standard relations for integrals over the Bessel function we ﬁnd using the relation

Z 0

1

p

dkk J 0 ðkRÞekz ¼ p!

Pp ðcos #Þ Y pþ1

and tan # = R/z, Y2 = R2 + z2, the following expression.

772

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

Iasy ðq; 0jz0 Þ pﬃﬃﬃ P2 ðsÞ P3 ðsÞ P4 ðsÞ 3 P1 ðsÞ ¼ þ 2 þ 6 þ 4ð7 6a Þ þ 2 2p X 2 X5 X3 X4 # pﬃﬃﬃ " P2 ðs0 Þ P 3 ðs0 Þ P4 ðs0 Þ 3 P1 ðs0 Þ þ 2 þ 6 þ 4ð7 6a Þ þ ez0 2 2p X 50 X 20 X 30 X 40

X 20 ¼ q2 þ Z 20

It remains to invert the Laplace transform, which is straightforward if tedious. Then the Hankel transform is needed. We shall not pursue this matter any further here but note it to illustrate that the emergent distribution is available.

ð80Þ

6. Diffusion theory

ð81Þ

It is always of value to examine the usefulness of diffusion theory in problems of the type discussed above because of its potential simplicity. Thus if we write the equation for a point source with cylindrical symmetry we have

and

z0 þ Z 0 s¼ ; X

Z0 s0 ¼ X0

For the second term involving the integral transient in (74) we ﬁnd the expansion in powers of k to be

Z

1

dtgð1; tÞ z =t e 0 ez0 þ tðez0 =t ez0 Þk HðtÞðt 1Þ 0 i 0 ez0 Þk2 þ ðB1 ez0 =t B 1 ez0 Þk3 þ Oðk4 Þ þðB0 ez0 =t B

trans ðkjz0 Þ ¼ 1 G 2

ð82Þ Inserting this into Eq. (70), we note that only terms in even powers of k survive leading to

Itrans ðq; 0jz0 Þ ¼

Z

1

1

dttgð1; tÞ z0 =t 9 ez0 Þ þ ðe 4pq 0 HðtÞðt 1Þ 8pq5 Z 1 dtgð1; tÞ 1 ez0 Þ ðB1 ez0 =t B 0 HðtÞðt 1Þ 3

ð83Þ

where H(t) is the classic Chandrasekhar H-function and

t2 ð1 ðt 1Þðt 1 þ z0 þ Z 0 ÞÞ t1 2 1 ¼ t ð1 ðt 1Þðt 1 þ Z 0 ÞÞ B t1

@ 2 Iðq; zÞ 1 @ @Iðq; zÞ 3dðqÞdðz z0 Þ m2 Iðq; zÞ þ þ ¼0 q @z2 @q 2pq q @q

ð89Þ

where in this case m2 = 3(1 c). This is subject to the boundary condition of zero return current, viz:

Iðq; 0Þ ¼

2 @Iðq; zÞ 3 @z z¼0

ð90Þ

There is, however, a defect in the physics of the source term in Eq. (89). It describes a source which emits particles isotropically whereas in the transport problem of interest the particles are emitted uni-directionally. To account for this feature and still use diffusion theory it is necessary to employ a hybrid form of diffusion theory in which the source term is replaced by the ﬁrst collision term from the transport equation. To see how this is done we return to Eq. (1) and write the ﬂux as

Iðx; y; z; l; uÞ ¼ Wðx; y; z; l; uÞ þ Aðx; y; zÞdðl l0 Þdðu u0 Þ ð91Þ

B1 ¼

ð84Þ

The total asymptotic ﬂux is given by the sum of Eqs. (79) and (83) and is useful for large values of q. For absorbing materials the asymptotic behaviour is more difﬁcult to obtain, however for 1 c small, we ﬁnd

Iasy ðq; 0jz0 Þ Iasy ðq; 0jz0 Þc¼1 ð1 þ mXÞe

ð88Þ

ð79Þ where

X 2 ¼ q2 þ ðz0 þ Z 0 Þ2 ;

Speciﬁcally,

Iðk; 0; l; ujl ; u 0 ; z0 Þ l 1 Z c Hð1=p; kÞepz0 ¼ dp SðkÞH ; k Þðu þ plÞ 4p 2pi L u ðu0 þ pl

mX

Inserting this into Eq. (1) and collecting up the delta, and nondelta function terms, we ﬁnd

l

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @ @ Wðx; y; z; l; uÞ 1 l2 cos u þ sin u @x @y c c þ Wðx; y; z; l; uÞ ¼ W0 ðx; y; zÞ þ Aðx; y; zÞ ð92Þ 4p 4p

@ Wðx; y; z; l; uÞ þ @z

with A(x, y, z) given by

l0

@Aðx; y; zÞ þ @z

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @Aðx; y; zÞ @Aðx; y; zÞ 1 l20 cos u0 þ sin u0 @x @y

¼ Sðx; y; zÞ

5. Emergent angular distribution

ð93Þ

Taking the Fourier transform in x and y leads to In addition to the surface ﬂux, an important quantity as far as actual measurement is concerned is the emergent angular distribu ; u0 ; z0 Þ. From Eq. (7), we tion from the surface, i.e. Iðq; #; 0; l; ujl can write this as

Iðk; 0; l; ujl ; u 0 ; z0 Þ ¼

Z

1 c 0 0 dz I0 ðk; z0 jz0 Þeuz =l 4pl 0 SðkÞ Þdðu u0 Þeu0 z=l þ dðl l l

c ^ I0 ðk; u=ljz0 Þ

4pl

ð85Þ

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ zÞ ¼ Sðk; zÞ 1 l20 ðkx cos u0 þ ky sin u0 ÞAðk;

ð94Þ

, the solution of this equation for Sðk; zÞ ¼ For l0 ¼ l SðkÞdðz z0 Þ is

ð95Þ

We now return to Eq. (92) and reduce it to the equivalent diffusion approximation. This leads to

r2 W0 ðx; y; zÞ m2 W0 ðx; y; zÞ þ 3cAðx; y; zÞ ¼ 0

ð96Þ

Taking Fourier transforms in the x and y directions, we ﬁnd

ð86Þ

2

0 ðk; zÞ d W dz

which from Eq. (24) is

l 1 c SðkÞ 1 Iðk; 0; l; ujl ; u0 ; pÞ ¼ H ;k H ;k 4p u þ pl u0 þ pl p u

zÞ @ Aðk; þ @z

zÞ ¼ SðkÞ Hðz0 zÞeu0 ðz0 zÞ=l Aðk; l

The ﬁrst term on the right hand side of Eq. (85) is related to the Laplace transform (12), viz:

Iðk; 0; l; ujl ; u 0 ; z0 Þ ¼

l0

2

3c 2 ðm2 þ k ÞW SðkÞHðz0 zÞeu0 ðz0 zÞ=l ¼ 0 0 ðk; zÞ þ l ð97Þ

ð87Þ

and /0 is where the implicit dependence of W0(k, z) on z0 ; l suppressed.

773

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

or

We may write the solution of Eq. (97) as

0 ðk; zjz0 Þ W ¼

3c SðkÞ l

Z

z0

D ðk; zjzÞeðz0 zÞ=li dzG

Gðkjz 0Þ ¼

pﬃﬃﬃﬃﬃﬃﬃﬃﬃ

0

where we have restored the dependence on z0. The Greens function GD satisﬁes the equation 2

dz

ð106Þ

2 ðkx cos u0 þky sin u0 Þðz0 zÞ=l 1l

ð98Þ

2 d G zÞ D ðk; zj

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 2 2 k þ k21 k þ m21 e k þm1 z0

2 ðm2 þ k ÞG zÞ þ dðz zÞ ¼ 0 D ðk; zj

where a possible choice is k21 ¼ 3 and that Eq. (50) or Eq. (102) become

m21 ¼ 3ð1 cÞ. It is then clear

; u 0 ; z0 Þ W0 ðq; #jl ¼

ð99Þ

Z 1 0 0 ; # u0 Þ dz0 eðz0 z0 Þ=l dkkJ 0 ðkRðq; z0 z00 ; l 0 2pl 0 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 0 2 2 ð107Þ SðkÞ k þ k21 k þ m21 e k þm1 z0 1

Z

z0

with the boundary condition or more conveniently

D ðk; 0jzÞ ¼ 2 dGD ðk; zjzÞ G 3 dz z¼0

ð100Þ

; u 0 ; z0 Þ W0 ðq; #jl ¼

0 ðk; 0jz0 Þ, we need As we are interested only in the surface ﬂux W only calculate GD ðk; 0jz0 Þ. This is found by elementary means and leads to

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 2 e m þk z0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ GD ðk; 0jz0 Þ ¼ 3 1 þ 2 m2 þ k2 3

ð101Þ

which is related to G in Eq. (28) by G ! 3cGD . For c near unity this is a reasonable approximation for many problems. Then after some manipulation following the techniques leading to Eq. (45), we ﬁnd

; u0 ; z0 Þ ¼ W0 ðq; #jl

c pl

Z 0

z0

0

0

dz0 eðz0 z0 Þ=l

Z

1

Z

This result is similar in form to Eq. (102). There are some advantages in using this technique for the H-function as shown in Loyalka and Williams (in press). 7. Diffusion theory: non-absorbing case

ð102Þ

W0 ðqjz0 Þ ¼

¼ 1, Eq. (102) simpliﬁes to For the special case of l

c

Z

p

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 dkkSðkÞJ ðkqÞ 2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e m þk z0 ez0 2 2 1 m2 þ k 1 þ 23 m2 þ k

1

0

¼

W0 ðq; #ju0 ; z0 Þ Z

p

1 w

0

e

Z

we pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

0

1

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dkkJ 0 k q2 2qw cosð# u0 Þ þ w2 SðkÞ

0

2 3

ð103bÞ

which does depend on #. It is worth noting at this stage that there is a useful approximation to the H-function appearing in Eq. (25). This is discussed fully in Williams (2007a) and the result is due to Rybicki (1971). Thus we may write

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 Y n pþ k þ k2j 1 1 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ H ; k Hn ;k ¼ 2 p p j¼1 p þ k þ m2j

ð104Þ

where kj and mj are constants deﬁned in Williams (2007a,b). Thus a ~ ﬁrst approximation to GðkjpÞ deﬁned in Eq. (28) will be

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 k þ k21 k þ m21 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 p þ k þ m21

p X2 "

þ2

1 P1 ðs0 Þ

p

X 20

P2 ðsÞ X3

þ2

ð105Þ

þ 22

P2 ðs0 Þ X 30

ð109Þ

P3 ðsÞ X4

þ 22

þ

728 P4 ðsÞ þ 27 X 5

P3 ðs0 Þ X 40

þ

# 728 P4 ðs0 Þ þ ez0 27 X 50 ð110Þ

where Pn(s) are the Legendre polynomials and

X 2 ¼ q2 þ ðz0 þ 2=3Þ2 ;

m2 þk2 z0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1þ m2 þ k2

~ GðkjpÞ ¼

1

dkkJ0 ðkqÞðekðz0 þ2=3Þ 11 2 91 3 4 ez0 2k=3 Þ 1 þ k þ k þ k þ Oðk Þ 9 81

p

1 P1 ðsÞ

¼ 0, we ﬁnd which as expected is independent of #. For the case of l

¼

Z

W0 ðq; 0jz0 Þ

ð103aÞ

c

1

Carrying out the integrals over k leads to

W0 ðqjz0 Þ ¼

z0

In the special case of c = 1 and m = 0, we may expand the integrand in Eq. (103a) in powers of k to get, for SðkÞ ¼ 1,

dkkJ0 ðkRðq; z0 0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 0 2 e m þk z0 0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ z0 ; l; # u0 ÞÞSðkÞ 2 1 þ 23 m2 þ k

Z 1 0 0 ; # u0 Þ dz0 eðz0 z0 Þ=l dkkJ 0 ðkRðq; z0 z00 ; l 0 0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ 2 2 0 e k þm1 z0 ð108Þ SðkÞ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 k þ k21 þ k þ m21 3c 2pl

X 20

2

¼ q þ 4=9;

s ¼ ðz0 þ 2=3Þ=X

s0 ¼ ð2=3Þ=X 0

ð111Þ

These results are analogous to those in Eq. (79) and it is useful to note that for large pﬃﬃﬃ q and z0 >> 1, the ratio of surface ﬂuxes ‘transport/diffusion’ is 3=2 ¼ 0:8660 . . .. We shall discuss this matter again in Section 10. 8. An alternative approach to the Hankel inversion I Evaluating the Hankel inversion integral as deﬁned in Eqs. (70) and (103a) requires quadrature involving the Bessel function J0(x). This is oscillatory and while the quadrature can be carried out, it involves a slowly converging series which must be summed by using the Euler, Shanks or other methods (Davis and Rabinowitz, 1967; Ganapol, 2008a). There is, however, an alternative procedure which can sometimes have some beneﬁts, especially when the source S(q) is no longer a point. We shall illustrate below the general procedure and then apply it to some of the problems discussed above.

774

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

Consider the Hankel inversion integral

GðqÞ ¼

Z

1 2p

1 0

dkkJ0 ðkqÞSðkÞFðkÞ

ð112Þ

FðxÞ ¼

1 2p

1

FðxÞ ¼

ikx dke FðkÞ

ð113Þ

1

Thus

ð126Þ

Z

b

ð127Þ

q

sjxj

dsgðsÞe

ð114Þ

Now if we invert (113) by writing

Z

dq0 q0 Sðq0 ÞGðqjq0 Þ

Z q dssgðsÞ K 0 ðsqÞ dq0 q0 Sðq0 ÞI0 ðsq0 Þ 0 a Z 1 þI0 ðsqÞ dq0 q0 Sðq0 ÞK 0 ðsq0 Þ

GðqÞ ¼ 2

b

a

FðkÞ ¼

1

0

It is often possible by considering appropriate contours, to evaluate Eq. (113) in the form

Z

Z

GðqÞ ¼ 2p

Now deﬁne the Fourier integral

Z

The case of a uniform disk is also of interest. To obtain the appropriate form, we note from Eq. (125) that we may write G(q) for any arbitrary S(q) as

1

ikx

dxe

We must now evaluate

2pK 0 ðsqÞ

FðxÞ

Z

q

dq0 q0 Sðq0 ÞI0 ðsq0 Þ þ 2pI0 ðsqÞ

0

ð115Þ

Z

1

dq0 q0 Sðq0 ÞK 0 ðsq0 Þ

q

ð128Þ

1

For the case of a disk source of radius a such that

we ﬁnd with (114),

FðkÞ ¼2

Z

b

ds

a

sgðsÞ

ð116Þ

2

k þ s2

SðqÞ ¼

1

1

p

Z

b

dssgðsÞ

Z

a

1

dk 0

ð129Þ

q>a

¼ 0;

Inserting Eq. (116) into (112), leads to

GðqÞ ¼

q
;

pa2

where the source has been normalised to Eq. (47), we ﬁnd

kSðkÞJ 0 ðkqÞ

ð117Þ

2

k þ s2

For the special case of a point source where

2 K 0 ðsqÞ a2

Z

a

dq0 q0 Iðsq0 Þ ¼

0

2 I1 ðasÞK 0 ðsqÞ; as

q>a

ð130Þ

and

dðqÞ SðqÞ ¼ 2pq

ð118Þ

we ﬁnd SðkÞ ¼ 1 and hence using the relation

Z

1

dk

0

kJ 0 ðkqÞ 2

k þ s2

¼

¼ K 0 ðsqÞ

ð119Þ

G(q) becomes

GðqÞ ¼

1

p

Z

b

dssgðsÞK 0 ðsqÞ

ð120Þ

a

This integral is usually far more convenient numerically than that of Eq. (112) which contains the oscillating function J0(kq). This procedure was ﬁrst used by Williams (1982) to solve a transport problem. In addition, suppose we have a ring source of the form

SðqÞ ¼

2 K 0 ðsqÞ a2

dðq q0 Þ 2pq

ð121Þ

for which

0

dk

dq0 q0 I0 ðsq0 Þ þ

0

2 I0 ðsqÞ a2

1 ð1 saK 1 ðasÞI0 ðsqÞÞ; a2 s2

Z q

a

dq0 q0 K 0 ðsq0 Þ

q
ð131Þ

None of these functions involving In or Kn is oscillatory and the integrals over s in Eq. (127) may be carried out efﬁciently by standard methods. One other source of some interest is that of the Gaussian, viz:

SðqÞ ¼

a2 a2 q2 e p

ð132Þ

Eq. (128) then becomes

2a2 K 0 ðsqÞ

Z

q

0

dq0 q0 ea

2 q2

I0 ðsq0 Þ þ2a2 I0 ðsqÞ

Z

1

q

dq0 q0 ea

2 q2

K 0 ðsq0 Þ ð133Þ

but it is not possible to reduce this expression to anything simpler.

ð122Þ

We note in Eq. (102) that the essential factor is (with z00 ¼ w)

and we ﬁnd the integral 1

q

8.1. Application to a diffusion theory problem

SðkÞ ¼ J ðkq Þ 0 0 Z

Z

kJ 0 ðkq0 ÞJ 0 ðkqÞ 2

k þ

ð123Þ

s2

But this integral is given by Gradshteyn and Ryzhik (1994) as

I0 ðsq0 ÞK 0 ðsqÞ; I0 ðsqÞK 0 ðsq0 Þ;

q > q0 q < q0

ð124Þ

Thus we may write

GðqÞ ¼ Gðqjq0 Þ ¼

1

p

Z a

b

I0 ðsq0 ÞK 0 ðsqÞ; dssgðsÞ I0 ðsqÞK 0 ðsq0 Þ;

q > q0 q < q0

ð125Þ

For a multi-ring source with different weights and radii, we simply add the appropriate solutions of (103a) with q0 = qi (i = 1, 2, . . . , N).

FðkÞ ¼

e

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

m2 þk2 w

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 þ 23 m2 þ k

ð134Þ

If this can be put into a suitable form then the integrations over w and k follow easily. Thus consider

Z

1 FðxÞ ¼ 2p

1

ikx

dke 1

e

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

m2 þk2 w

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 þ 23 m2 þ k

ð135Þ

We note that the integrand has branch points at k = ±im. Thus deforming the contour along the real axis to the imaginary axis each side of the cut from im to i1, it is easily found that

FðxÞ ¼

1

p

Z 0

1

ﬃ sin uw þ 23 u cos uw pﬃﬃﬃﬃﬃﬃﬃﬃﬃ duu m2 þu2 jxj pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e 1 þ 49 u2 m2 þ u2

ð136Þ

775

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

By inverting this transform as described in Section 8, we may write

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2 Z 2 sin uw þ 23 u cos uw e m þk w 2 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ duu

2 2 p 0 1 þ 49 u2 ðk þ m2 þ u2 Þ 1 þ 23 m2 þ k

p

Z

1

duu

0

ð137Þ

Z 1 sin uw þ 23 u cos uw dkkSðkÞJ 0 ðkRÞ

4 2 2 1 þ 9u 0 ðk þ m2 þ u2 Þ

ð138Þ

We may now use any of the sources described above. For the simple case of a point source, for which SðkÞ ¼ 1, we ﬁnd

Z pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sin uw þ 23 u cos uw 2 1

K0 duu u2 þ m2 R 4 2 p 0 1 þ 9u

ð139Þ

where in general

; # u0 Þ R2 ¼ R2 ðq; z0 ; l qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 cosð# u0 Þ=l þ ð1 l 2 Þw2 =l 2 ¼ q2 2qw 1 l

ð140Þ

The integral over w (z00 in Eq. (102)), can then be carried out by ¼ 1, i.e. a source quadrature. However, in the special case of l pointing normally towards the surface, we can do this integral analytically to get

2c W0 ðqjz0 Þ ¼ 2 3p

Z

1

0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ duuK 0 u2 þ m2 q

fuez0 þ ð3 þ 2u2 Þ 1 þ 49 u2 ð1 þ u2 Þ

sin uz0 u cos uz0 g

ð141Þ

This is to be compared with the original form of

W0 ðqjz0 Þ ¼

c

p

Z 0

1

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 e k þm z0 ez0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dkkJ0 ðkqÞ 2 2 1 þ 23 k þ m2 1 k þ m2

Even in Eq. (141), it is still necessary to evaluate oscillatory functions (in this case over z0) but the standard library routines for such trigonometric integrals are well-developed (IMSL libraries) and because of the presence of the K0 function the integrand decays rapidly and only a few oscillations have any signiﬁcant effect. For the case of the ring source at q = q0, Eq. (141) remains the same but with the replacement

9

= q > q0 > ; q < q0 > ð143Þ

Although the Gaussian source leads to a difﬁcult expression, as we saw in Eq. (133), this is one case where the direct Hankel transform quadrature may be preferable, because in this case 2 SðkÞ ¼ expðk =4a2 Þ and the original equation becomes

Z

k2 =4a2

1

J 0 ðkqÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ e 2 1 þ 23 k þ m2

dkke

0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

k2 þm2 z0

ð147Þ

which will converge rapidly for small a. Finally we note the case of an annular ring with internal and external radii R0 and R1, respectively. In that case we have

2

SðkÞ ¼

kðR21 R20 Þ

½R1 J 1 ðkR1 Þ R0 J 1 ðkR0 Þ

ð148Þ

which is related to the problem of the difference between two single disk problems. Also, as expected, as R1 ? R0, then SðkÞ ! J 0 ðkR0 Þ, which is the ring source as derived above in Eq. (122). Applications of this technique to the transport from of G, as deﬁned by Eq. (68), presents some difﬁculties and we have yet to complete the task fully. However the appropriate form for the asymptotic term in Eq. (68) has been obtained. The method may also be applied with little difﬁculty in principle to the form of the H-function deﬁned by Eq. (104) but at the expense of considerable algebra.

It is of interest to examine the form taken by the ﬂux when the source leads to an angular variation of the ﬂux at a given value of ¼ 1 leads the radial co-ordinate q. For example, whilst the case l to a ﬂux which depends on q only, the point, ring source and the < 1, lead to ﬂuxes which are functions of disk sources, for l (q, #). From Eqs. (139) and (140), we may write the associated ﬂux on the surface using Eq. (102). This expression explicitly shows the dependence on the azimuthal angle /0. Normally, it is only necessary to consider the difference # /0 and so setting /0 = 0 leads to no loss of generality. However, suppose the source emits particles non-uniformly in /0, e.g. with a distribution Q(/0). In that case the ﬂux would take the form

Z 2p

; z0 Þ ¼ W0 ðq; #jl

0

; u0 ; z0 Þ du0 Q ðu0 ÞW0 ðq; #jl

ð149Þ

In the special case where particles are emitted radially and

l ¼ 0, we have

Similarly, Eq. (142) is the same but with

J 0 ðkqÞ ! J 0 ðkq0 ÞJ 0 ðkqÞ

ð146Þ

8.2. q # Distributions for diffusion theory

ð142Þ

8 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ > < I0 u2 þ m2 q0 K 0 u2 þ m2 q ; pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ K0 u2 þ m2 q ! > : I0 u2 þ m2 q K 0 u2 þ m2 q0 ;

2 J ðkaÞJ 0 ðkqÞ ak 1

J 0 ðkqÞ !

Now consider the integral over k in Eq. (102) to get

2

and for Eq. (142),

ð144Þ

One sees immediately the advantage of the new approach because if Eq. (142) were to be evaluated numerically with (144), then the zeros of the integrand, which are needed for the series summation, are now no longer those of J0(x) = 0 but of the product J0(x)J0(ax) or, for the disk problem below, the roots of J0(x)J1(ax). This is a much more tedious affair. For the disk source, we have (141) with pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ K0 u2 þ m2 q 9 8 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 > > ﬃ I1 u2 þ m2 a K 0 u2 þ m2 q ; q > a = < pﬃﬃﬃﬃﬃﬃﬃﬃﬃ a u2 þm2 ! pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ p ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ > ; : 2 22 2 1 a u2 þ m2 I0 u2 þ m2 q K 1 u2 þ m2 a ; q < a > a ðu þm Þ

ð145Þ

W0 ðq; #jz0 Þ ¼

Z

c

p

1

w

dwe

0

SðkÞ

e

Z

1

dkk

0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ

Z 2p 0

^ q; # u ; wÞÞ du0 Q ðu0 ÞJ 0 ðkRð 0

m2 þk2 z0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 þ 23 m2 þ k

ð150Þ

^ 2 ¼ q2 2qw cosð# u Þ þ w2 . where R 0 In terms of the K0 representation,

W0 ðq; #jz0 Þ ¼

Z

c

p2

Z 0

Z 2p sinðuwÞ þ 23 u cosðuwÞ

dwe duu du0 Q ðu0 Þ 1 þ 49 u2 0 0 ^ dkkSðkÞJ 0 ðkRÞ ð151Þ 2 2 k þ m þ u2

1

0 1

w

Z

1

776

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

We have discussed above the various forms taken by the integral over k for different source shapes, but for the point source, we have

Z 2p 0

du0 Q ðu0 ÞK 0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ^ q; # u ; wÞ m2 þ u2 Rð 0

ð152Þ

It is readily seen how one can introduce a source which emits in a given range of /0 only by suitable deﬁnition of Q. For the general case of a uni-directional source d(/ /0), this approach is also applicable to the transport equation for radial emission, Eq. (63), which we may re-write in the form

1 2p

I0 ðq; #; 0ju0 ; z0 Þ ¼

Z

1

dkkSðkÞGðkjz 0Þ

0

Z

Inserting this into Eq. (157) leads to

I0 ðq; 0jz0 Þ ﬃﬃﬃﬃﬃﬃﬃﬃ2 Z p1 m 2c udu pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ 2 p 0 ð1 þ 4u2 =9Þ 1 m2 u2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 u2 þ m2 ; q sinðuz0 Þ þ u cosðuz0 Þ F 3 Z 2c 1 udu 2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ cosðuz0 Þ u sinðuz0 Þ þ 2 pﬃﬃﬃﬃﬃﬃﬃﬃ 3 p 1m2 ð1 þ 4u2 =9Þ m2 þ u2 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 F u þm ;q ð160Þ where

1

w

dwe 0

~ q; #; u ; wÞÞ J 0 ðkRð 0

Fðs; qÞ ¼

dk

kSðkÞJ 0 ðkqÞ

ð161Þ

2

s2 þ k

Then following the procedure described in Section 8.1, we can write F(s, q) in terms of the Bessel functions K0 and I0.

where

~ 2 ¼ q2 2qw cosð# u Þ þ w2 R 0

ð154Þ

we shall see that the integration over w has some numerical advantages as it damps out the oscillations of the Bessel function. As a simple example of Eq. (149), let us assume that Q(/0) = 1/2p over the range /0(0, 2p), i.e. the emission is isotropic in the radial direction. In that case we must evaluate the term in Eq. (50) as

Z 2p

^ q; # u ; wÞÞ du0 J 0 ðkRð 0 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 ðz0 z00 Þ=l Þ ¼ J 0 ðkqÞJ 0 ðk 1 l 0

ð155Þ

9. An alternative approach to the Hankel inversion II The method described in the previous section is not convenient for dealing with transport problems except in very special cases (Williams, 1982, 2007a) as it is very difﬁcult to obtain a numerically useful Fourier inversion of Eq. (68). We therefore describe another method for transforming the Hankel transform into a more convenient form for numerical evaluation by noting that the Bessel function can be written as (Gradshteyn and Ryzhik, 1994)

J n ðzÞ ¼

1 I0 ðq; 0jz0 Þ ¼ 2p

Z

z 0 =l

w

Z

Z

0

1

kSðkÞ dk pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Gðkjz 0 ÞJ 0 ðkqÞ 2 1þk

ð156Þ

2nþ1

sin

ðtÞ

e2z cotðtÞ

J 0 ðzÞ ¼ gðzÞ sinðzÞ þ hðzÞ cosðzÞ

ð162Þ

ð163Þ

where

ð157Þ

ð158Þ

Although it is not difﬁcult to evaluate the quadrature over J0 in Eq. (157), there are some advantages in using the Fourier transform technique described in Section 8. This may be done by noting that Eq. (157) has branch points at k = ±im, ±i, thus by deforming the contour onto the imaginary axis, we ﬁnd

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 e k þm z0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 1þk 1 þ 2 k þ m2 =3 Z 2 1 sds pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ p m ð1 þ 4ðs2 m2 Þ=9Þ 1 s2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 s2 m2 cos sin s2 m 2 z0 þ s2 m2 z0 2 3 s2 þ k Z h pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 sds pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ cos þ s2 m2 z0 p 1 ð1 þ 4ðs2 m2 Þ=9Þ s2 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ s2 m2 sin s2 m 2 z0 ð159Þ 2 3 s2 þ k

2

Z p=2

cosðt=2Þ dt pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e2z cotðtÞ cosðtÞ sinðtÞ rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Z ﬃ 2 1 e2t 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ dt qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ t2 þ z2 þ t 2 p 0 tðt 2 þ z2 Þ r ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ Z 2 1 e2zt 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dt qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ t2 þ 1 þ t ¼ 2 p 0 tðt 2 þ 1Þ

gðzÞ ¼

We may illustrate this with diffusion theory, using

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ k2 þm2 z0 e pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Gðkjz 0 Þ ¼ 2c 2 1 þ 2 k þ m2 =3

dt

0

From which we may write

0

As expected, Eq. (156) is independent of # due to the radial sym ¼ 1, Eq. (156) reduces to Eq. (70). On the other hand metry. For l ¼ 0 and we ﬁnd for radial emission, l

1 2p

Z p=2

cosn1=2 ðtÞ sinðz nt þ t=2Þ

1

dkkSðkÞGðkjz 0 qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ wÞJ 0 ðkqÞJ 0 ðk 1 l 2 wÞ l dwe

0

2nþ1 zn Cð1=2ÞCðn þ 1=2Þ

whence, after a simple transformation

I0 ðq; 0jz0 Þ ¼

1 0

ð153Þ

1 2p

Z

p

and

ð164Þ

Z p=2

sinðt=2Þ dt pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e2z cotðtÞ cosðtÞ sinðtÞ 0 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Z ﬃ 2 1 e2t 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ dt qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð t2 þ z2 tÞ 2 p 0 tðt 2 þ z2 Þ r ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ Z 2 1 e2zt 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ q dt ð t þ 1 tÞ ¼ p 0 tðt 2 þ 1Þ 2

hðzÞ ¼

2

0

p

ð165Þ

It is useful to note that

lim

z!1

pﬃﬃﬃﬃﬃﬃ pzgðzÞ ¼ 1

and

lim

z!1

pﬃﬃﬃﬃﬃﬃ pzhðzÞ ¼ 1

Also as z ? 0, h(0) = 1 and g(z) ? 2 log (z)/p. A further useful property is that the functions h(z) and g(z) are monotonically decreasing functions of z. Thus in Eq. (163) we have extracted

777

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

the oscillatory behaviour of the Bessel function and written it in terms of simple trigonometric functions. There is yet another representation of the Bessel function which arises from the relation

J 0 ðzÞ ¼

2

Z

p

1

1

The disk source deﬁned by Eq. (55), requires the representation of J1(z) via the use of Eq. (162). We deﬁne

J 1 ðzÞ ¼ pðzÞ sin z qðzÞ cosðzÞ where

sinðzuÞ du pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ u2 1

hðzÞ ¼

2

p

1

0

sinðztÞ dt pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ tðt þ 2Þ

gðzÞ ¼

2

p

Z

0

1

cosðztÞ dt pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ tðt þ 2Þ

Although these integrals have oscillatory integrands they are in a form that is particularly useful for expansion about z = 0 (Grosjean, 1965). These expansions are used in our numerical analysis discussed in Section 10.

I0 ðq; 0jz0 Þ ¼

Z

1 2p Z

1

0 1

0 ðkjz0 Þ þ dkk sinðkqÞgðkqÞSðkÞG

ð166Þ

This result has the very useful property that there exist a number of very efﬁcient quadrature library routines for solving the Fourier integrals; more of which below. In the case of a point source, SðkÞ ¼ 1 and the integral is a straightforward sum of the two Fourier terms. In the case of a ring of radius q0, however, SðkÞ ¼ J 0 ðkqÞ and we must re-write the term J 0 ðkqÞSðkÞ ¼ J 0 ðkqÞJ 0 ðkq0 Þ as

J 0 ðkqÞJ 0 ðkq0 Þ ¼ ½gðkqÞ sinðkqÞ þ hðkqÞ cosðkqÞ½gðkq0 Þ sinðkq0 Þ ð167Þ

Multiplying this out and re-arranging the trigonometric functions we have

2J 0 ðkqÞJ 0 ðkq0 Þ ¼ ðgg 0 þ hh0 Þ cos kðq q0 Þ þ ðgh0 hg 0 Þ

Z 1 1 0 ðkjz0 Þ dkkðgg 0 þ hh0 Þ cos kðq q0 ÞG 4p 0 Z 1 1 0 ðkjz0 Þ þ dkkðgh0 hg 0 Þ sin kðq q0 ÞG 4p 0 Z 1 1 0 ðkjz0 Þ þ dkkðhh0 gg 0 Þ cos kðq þ q0 ÞG 4p 0 Z 1 1 0 ðkjz0 Þ ð169aÞ þ dkkðgh0 þ hg 0 Þ sin kðq þ q0 ÞG 4p 0

This formulation is very useful because it avoids having to evaluate the highly oscillating and irregular integrand J0(kq)J0(kq0), and reduces the problem to evaluating trigonometric Fourier integrals which have a long numerical history. The case of q = 0 needs special attention and reduces to

1 I0 ð0; 0jz0 Þ ¼ 2p

Z 0

1

0

8

Z p=2

ð171Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ cosðtÞ sinðt=2Þ

e2z cotðtÞ 3 sin ðtÞ vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! u Z 1 ut 8 t 1=4 ¼ dtðz2 þ t 2 Þ t 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e2t 2 pz 0 t 2 þ z2

qðzÞ ¼

p

dt

0

ð172Þ

We ﬁnd that

lim zpðzÞ ¼

2

p

2

lim qðzÞ ¼

p

z!0

pﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃ Also as z ? 1, we ﬁnd pzpðzÞ ! 1 and pzqðzÞ ! 1. Further useful representations of p(z) and q(z) are

pðzÞ ¼

2

p

Z

1

dt

0

ðt þ 1Þ sinðztÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ; tðt þ 2Þ

qðzÞ ¼

2

p

Z

1

dt 0

ðt þ 1Þ cosðztÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ tðt þ 2Þ

from which we may obtain series solutions in powers of z. In order to use this form to represent the disk source, Eq. (55), we write

J 0 ðkqÞJ 1 ðkq0 Þ ¼ ½gðkqÞ sinðkqÞ þ hðkqÞ cosðkqÞ½pðkq0 Þ sinðkq0 Þ qðkq0 Þ cosðkq0 Þ and after some algebra this leads to

cos kðq þ q0 Þ þ ðhp0 gq0 Þ sin kðq þ q0 Þ

ð168Þ

where we have used an obvious notation. We see therefore that the quadrature for the ﬂux is the sum of four Fourier integrals, which are each readily evaluated by the sub-routines mentioned above. Explicitly, the ﬂux takes the form

I0 ðq; 0jz0 Þ ¼

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ cosðtÞ cosðt=2Þ

2J 0 ðkqÞJ 1 ðkq0 Þ ¼ ðgp0 hq0 Þ cos kðq q0 Þ ðgp0 þ hq0 Þ

sin kðq q0 Þ þ ðhh0 gg 0 Þ cos kðq þ q0 Þ þ ðgh0 þ hg 0 Þ sin kðq þ q0 Þ

dt

and

0

þ hðkq0 Þ cosðkq0 Þ

p

z!0

1 2p

0 ðkjz0 Þ dkk cosðkqÞhðkqÞSðkÞG

Z p=2

and

9.1. Application to transport and diffusion theory The Hankel transform of the surface ﬂux is given by Eq. (70) with the appropriate form of SðkÞ and G0 ðkjz0 Þ. Thus we may use Eq. (163) and re-write Eq. (70) as

8

e2z cotðtÞ 3 sin ðtÞ vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! u Z 1 u 8 t 2 2 1=4 t t ¼ dtðz þ t Þ 1 þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e2t 2 pz 0 t 2 þ z2

pðzÞ ¼

Now setting u = t + 1, we ﬁnd after expanding the sin function that

Z

ð170Þ

0 ðkjz0 Þ dkkðh0 cosðkq0 Þ þ g 0 sinðkq0 ÞÞG ð169bÞ

ðgq0 þ hp0 Þ sin kðq q0 Þ

ð173Þ

again in an obvious notation. As for the ring source in Eq. (169a) the disk source is also the sum of four Fourier terms and again reduces the problem to a trigonometric one as shown below

I0 ðq; 0jz0 Þ ¼

Z

1

1

0 ðkjz0 Þ dkðgp0 hq0 Þ cos kðq q0 ÞG 2pq0 0 Z 1 1 0 ðkjz0 Þ dkðgp0 þ hq0 Þ cos kðq þ q0 ÞG 2pq0 0 Z 1 1 0 ðkjz0 Þ dkðhp0 gq0 Þ sin kðq þ q0 ÞG þ 2pq0 0 Z 1 1 0 ðkjz0 Þ ð174aÞ dkðgq0 þ hp0 Þ sin kðq q0 ÞG 2pq0 0

Again, the case of q = 0 needs special attention and we ﬁnd

I0 ð0; 0jz0 Þ ¼

1

pq0

Z

1 0

0 ðkjz0 Þ dkðp0 sinðkq0 Þ q0 cosðkq0 ÞÞG ð174bÞ

The efﬁciency of this technique may be greatly improved if accurate rational approximations for the universal functions, g(z), h(z), p(z) and q(z) can be obtained. In this respect, it is useful to

778

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

note the work of Grosjean (1965) which describes a method for developing the functions g(z), h(z), p(z) and q(z) in a series about z = 0 which converges well up to z = 1. We will use this series expansion for the numerical results below, but it should be stressed that it was necessary to use Mathematica to get the required accuracy in the coefﬁcients.

In order to illustrate the above theory we have considered a few problems of general interest numerically and, where possible, have compared them with similar calculations by Ganapol (2008b). In all cases studied, we have found excellent agreement with the work of Ganapol and also with some independent results of Loyalka (2008). The only case in which deviations of one or two digits in the fourth signiﬁcant ﬁgure arise are for Ganapol’s ring source near the singularity. This is a very difﬁcult region to

treat numerically and the analytic solution derived in Appendix C is useful for checking such methods. However, we stress that in this work the emphasis is not on numerical results per se but rather on technique and method. Thus only a few illustrative numerical results are given. Comparisons with diffusion theory are also made where appropriate. In all cases we give the surface ﬂux. To illustrate the point source problem deﬁned by Eq. (70) with SðkÞ ¼ 1, we present Table 1a for the transport theory ﬂux as a function of q in the range (0.1, 5) and z0 = 0.1, 1.0, 2.0, 4.0 with c = 0.99. These values will be useful for benchmarking problems and, subject to the comments made below, agree to the ﬁgures given with the results of Ganapol (2008b). In order to compare the results with diffusion theory, or to be precise modiﬁed diffusion theory according to Eq. (103a), we also present Table 1b. Finally, to compare transport and diffusion theory, Table 1c shows the ratio of transport to diffusion ﬂuxes. These results deviate from unity the most near the source as expected. What is rather surprising,

Table 1a Point source c = 0.99, transport theory.

Table 1b Point source c = 0.99, diffusion theory.

10. Numerical results and discussion

q/z0

0.1

1

2

4

q/z0

0.1

1

2

4

0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50 3.60 3.70 3.80 3.90 4.00 4.10 4.20 4.30 4.40 4.50 4.60 4.70 4.80 4.90 5.00

6.1978E01 1.8717E01 8.7346E02 4.9919E02 3.2030E02 2.2138E02 1.6110E02 1.2175E02 9.4708E03 7.5365E03 6.1085E03 5.0268E03 4.1897E03 3.5301E03 3.0024E03 2.5746E03 2.2237E03 1.9330E03 1.6900E03 1.4851E03 1.3114E03 1.1629E03 1.0353E03 9.2498E04 8.2922E04 7.4569E04 6.7251E04 6.0816E04 5.5137E04 5.0107E04 4.5638E04 4.1656E04 3.8098E04 3.4910E04 3.2047E04 2.9469E04 2.7143E04 2.5040E04 2.3134E04 2.1404E04 1.9831E04 1.8397E04 1.7088E04 1.5892E04 1.4797E04 1.3793E04 1.2871E04 1.2023E04 1.1242E04 1.0522E04

5.8696E01 3.1852E01 2.1818E01 1.6309E01 1.2753E01 1.0254E01 8.4052E02 6.9898E02 5.8800E02 4.9941E02 4.2764E02 3.6882E02 3.2011E02 2.7941E02 2.4515E02 2.1609E02 1.9130E02 1.7001E02 1.5165E02 1.3572E02 1.2185E02 1.0971E02 9.9057E03 8.9665E03 8.1358E03 7.3989E03 6.7431E03 6.1580E03 5.6344E03 5.1647E03 4.7424E03 4.3619E03 4.0182E03 3.7072E03 3.4252E03 3.1690E03 2.9358E03 2.7233E03 2.5292E03 2.3518E03 2.1892E03 2.0401E03 1.9032E03 1.7773E03 1.6613E03 1.5543E03 1.4556E03 1.3644E03 1.2799E03 1.2017E03

2.5478E01 1.5545E01 1.1760E01 9.6096E02 8.1552E02 7.0728E02 6.2190E02 5.5196E02 4.9319E02 4.4292E02 3.9934E02 3.6123E02 3.2763E02 2.9786E02 2.7136E02 2.4768E02 2.2645E02 2.0738E02 1.9020E02 1.7469E02 1.6066E02 1.4795E02 1.3642E02 1.2594E02 1.1640E02 1.0771E02 9.9775E03 9.2523E03 8.5887E03 7.9807E03 7.4230E03 6.9109E03 6.4402E03 6.0069E03 5.6077E03 5.2395E03 4.8996E03 4.5855E03 4.2949E03 4.0259E03 3.7765E03 3.5453E03 3.3305E03 3.1310E03 2.9455E03 2.7728E03 2.6119E03 2.4620E03 2.3221E03 2.1915E03

4.8190E02 3.4691E02 2.9478E02 2.6443E02 2.4317E02 2.2665E02 2.1294E02 2.0107E02 1.9050E02 1.8088E02 1.7201E02 1.6374E02 1.5599E02 1.4867E02 1.4174E02 1.3516E02 1.2890E02 1.2293E02 1.1725E02 1.1182E02 1.0665E02 1.0171E02 9.6996E03 9.2497E03 8.8205E03 8.4110E03 8.0203E03 7.6479E03 7.2927E03 6.9542E03 6.6316E03 6.3242E03 6.0314E03 5.7525E03 5.4869E03 5.2341E03 4.9933E03 4.7641E03 4.5460E03 4.3383E03 4.1407E03 3.9526E03 3.7735E03 3.6031E03 3.4410E03 3.2866E03 3.1397E03 2.9998E03 2.8667E03 2.7400E03

0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50 3.60 3.70 3.80 3.90 4.00 4.10 4.20 4.30 4.40 4.50 4.60 4.70 4.80 4.90 5.00

2.7223E01 1.1954E01 6.7077E02 4.2653E02 2.9267E02 2.1144E02 1.5859E02 1.2240E02 9.6641E03 7.7725E03 6.3483E03 5.2532E03 4.3963E03 3.7154E03 3.1672E03 2.7206E03 2.3531E03 2.0479E03 1.7923E03 1.5766E03 1.3934E03 1.2368E03 1.1021E03 9.8577E04 8.8470E04 7.9652E04 7.1926E04 6.5131E04 5.9131E04 5.3816E04 4.9092E04 4.4881E04 4.1115E04 3.7739E04 3.4705E04 3.1970E04 2.9501E04 2.7265E04 2.5237E04 2.3394E04 2.1715E04 2.0184E04 1.8784E04 1.7503E04 1.6328E04 1.5249E04 1.4256E04 1.3342E04 1.2499E04 1.1720E04

3.8766E01 2.7359E01 2.0962E01 1.6650E01 1.3505E01 1.1113E01 9.2448E02 7.7605E02 6.5654E02 5.5931E02 4.7949E02 4.1345E02 3.5843E02 3.1229E02 2.7336E02 2.4033E02 2.1216E02 1.8802E02 1.6722E02 1.4924E02 1.3361E02 1.1999E02 1.0806E02 9.7581E03 8.8345E03 8.0178E03 7.2933E03 6.6488E03 6.0739E03 5.5597E03 5.0986E03 4.6842E03 4.3108E03 3.9737E03 3.6687E03 3.3922E03 3.1410E03 2.9124E03 2.7040E03 2.5137E03 2.3397E03 2.1802E03 2.0339E03 1.8994E03 1.7757E03 1.6617E03 1.5566E03 1.4594E03 1.3695E03 1.2863E03

1.9075E01 1.4805E01 1.2333E01 1.0590E01 9.2469E02 8.1592E02 7.2508E02 6.4768E02 5.8082E02 5.2248E02 4.7123E02 4.2596E02 3.8580E02 3.5007E02 3.1819E02 2.8968E02 2.6413E02 2.4119E02 2.2056E02 2.0198E02 1.8521E02 1.7007E02 1.5637E02 1.4396E02 1.3270E02 1.2247E02 1.1316E02 1.0469E02 9.6958E03 8.9900E03 8.3447E03 7.7539E03 7.2125E03 6.7156E03 6.2591E03 5.8393E03 5.4528E03 5.0964E03 4.7676E03 4.4639E03 4.1830E03 3.9230E03 3.6821E03 3.4586E03 3.2512E03 3.0585E03 2.8792E03 2.7124E03 2.5570E03 2.4120E03

4.2848E02 3.6996E02 3.3530E02 3.1006E02 2.8981E02 2.7261E02 2.5748E02 2.4385E02 2.3136E02 2.1977E02 2.0895E02 1.9877E02 1.8915E02 1.8004E02 1.7140E02 1.6318E02 1.5535E02 1.4790E02 1.4081E02 1.3405E02 1.2760E02 1.2147E02 1.1562E02 1.1005E02 1.0475E02 9.9709E03 9.4908E03 9.0340E03 8.5996E03 8.1865E03 7.7938E03 7.4205E03 7.0657E03 6.7286E03 6.4083E03 6.1041E03 5.8150E03 5.5405E03 5.2797E03 5.0320E03 4.7968E03 4.5733E03 4.3611E03 4.1594E03 3.9679E03 3.7859E03 3.6130E03 3.4487E03 3.2926E03 3.1442E03

779

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

however, is that for larger values of z0, diffusion theory tends to become less accurate overall. For example, at q = 2.5, the values of the ratio of transport/diffusion are equal to 0.9372. . ., 0.9209. . .,

Table 1c Point source c = 0.99, ratio of transport to diffusion theory ﬂuxes.

q/z0

0.1

1

2

4

0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50 3.60 3.70 3.80 3.90 4.00 4.10 4.20 4.30 4.40 4.50 4.60 4.70 4.80 4.90 5.00

2.2767E+00 1.5658E+00 1.3022E+00 1.1703E+00 1.0944E+00 1.0470E+00 1.0158E+00 9.9469E01 9.8000E01 9.6963E01 9.6223E01 9.5689E01 9.5301E01 9.5014E01 9.4799E01 9.4634E01 9.4502E01 9.4391E01 9.4294E01 9.4194E01 9.4114E01 9.4025E01 9.3931E01 9.3833E01 9.3729E01 9.3618E01 9.3500E01 9.3375E01 9.3244E01 9.3107E01 9.2964E01 9.2815E01 9.2661E01 9.2503E01 9.2342E01 9.2177E01 9.2009E01 9.1839E01 9.1667E01 9.1494E01 9.1321E01 9.1146E01 9.0972E01 9.0798E01 9.0624E01 9.0452E01 9.0281E01 9.0111E01 8.9943E01 8.9776E01

1.5141E+00 1.1642E+00 1.0408E+00 9.7951E01 9.4435E01 9.2277E01 9.0918E01 9.0069E01 8.9561E01 8.9290E01 8.9188E01 8.9205E01 8.9309E01 8.9474E01 8.9681E01 8.9915E01 9.0166E01 9.0425E01 9.0686E01 9.0944E01 9.1196E01 9.1438E01 9.1669E01 9.1887E01 9.2091E01 9.2281E01 9.2449E01 9.2617E01 9.2764E01 9.2896E01 9.3014E01 9.3120E01 9.3212E01 9.3293E01 9.3362E01 9.3420E01 9.3468E01 9.3506E01 9.3535E01 9.3556E01 9.3569E01 9.3575E01 9.3574E01 9.3567E01 9.3555E01 9.3537E01 9.3515E01 9.3488E01 9.3458E01 9.3424E01

1.3357E+00 1.0500E+00 9.5354E01 9.0743E01 8.8194E01 8.6685E01 8.5769E01 8.5221E01 8.4914E01 8.4772E01 8.4746E01 8.4804E01 8.4922E01 8.5085E01 8.5281E01 8.5501E01 8.5736E01 8.5982E01 8.6234E01 8.6489E01 8.6744E01 8.6996E01 8.7244E01 8.7486E01 8.7722E01 8.7950E01 8.8169E01 8.8380E01 8.8581E01 8.8773E01 8.8956E01 8.9129E01 8.9292E01 8.9446E01 8.9591E01 8.9728E01 8.9855E01 8.9974E01 9.0085E01 9.0188E01 9.0283E01 9.0369E01 9.0453E01 9.0528E01 9.0597E01 9.0659E01 9.0716E01 9.0768E01 9.0815E01 9.0856E01

1.1247E+00 9.3770E01 8.7915E01 8.5283E01 8.3908E01 8.3139E01 8.2700E01 8.2458E01 8.2339E01 8.2302E01 8.2321E01 8.2379E01 8.2465E01 8.2573E01 8.2695E01 8.2828E01 8.2969E01 8.3117E01 8.3268E01 8.3422E01 8.3578E01 8.3735E01 8.3891E01 8.4047E01 8.4202E01 8.4355E01 8.4507E01 8.4656E01 8.4803E01 8.4947E01 8.5088E01 8.5227E01 8.5362E01 8.5494E01 8.5622E01 8.5747E01 8.5869E01 8.5987E01 8.6102E01 8.6214E01 8.6322E01 8.6426E01 8.6528E01 8.6626E01 8.6720E01 8.6811E01 8.6899E01 8.6984E01 8.7066E01 8.7145E01

Flux

transport diffusion

0.02

0.01

0.005

0

1

2

ρ

3

4

Fig. 1. Flux for a ring source: c = 0.99, q0 = 3.0, and z0 = 1.0.

5

Table 2 Ring source: c = 0.99, z0 = 1.0, and q0 = 3.0.

q

Transport

Diffusion

Trans./Diff

0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 2.91 2.92 2.93 2.94 2.95 2.96 2.97 2.98 2.99 3.00 3.01 3.02 3.03 3.04 3.05 3.06 3.07 3.08 3.09 3.10 3.20 3.30 3.40 3.50 3.60 3.70 3.80 3.90 4.00 4.10 4.20 4.30 4.40 4.50 4.60 4.70 4.80 4.90 5.00

5.1647E03 5.1728E03 5.1973E03 5.2383E03 5.2964E03 5.3720E03 5.4661E03 5.5796E03 5.7139E03 5.8704E03 6.0509E03 6.2577E03 6.4934E03 6.7611E03 7.0643E03 7.4074E03 7.7955E03 8.2348E03 8.7326E03 9.2977E03 9.9407E03 1.0675E02 1.1516E02 1.2486E02 1.3611E02 1.4932E02 1.6508E02 1.8449E02 2.0996E02 2.4938E02 2.5504E02 2.6130E02 2.6830E02 2.7629E02 2.8562E02 2.9690E02 3.1128E02 3.3132E02 3.6520E02 inﬁnity 3.6396E02 3.2907E02 3.0812E02 2.9289E02 2.8080E02 2.7070E02 2.6197E02 2.5426E02 2.4733E02 2.4101E02 1.9602E02 1.6633E02 1.4366E02 1.2537E02 1.1019E02 9.7401E03 8.6505E03 7.7152E03 6.9071E03 6.2052E03 5.5926E03 5.0555E03 4.5828E03 4.1651E03 3.7948E03 3.4653E03 3.1712E03 2.9081E03 2.6719E03

5.5597E03 5.5689E03 5.5966E03 5.6432E03 5.7092E03 5.7952E03 5.9022E03 6.0316E03 6.1847E03 6.3634E03 6.5698E03 6.8066E03 7.0767E03 7.3837E03 7.7316E03 8.1254E03 8.5705E03 9.0734E03 9.6415E03 1.0283E02 1.1009E02 1.1828E02 1.2754E02 1.3800E02 1.4979E02 1.6305E02 1.7794E02 1.9457E02 2.1307E02 2.3360E02 2.3577E02 2.3797E02 2.4019E02 2.4244E02 2.4471E02 2.4701E02 2.4934E02 2.5170E02 2.5409E02 2.5652E02 2.5322E02 2.4998E02 2.4679E02 2.4365E02 2.4056E02 2.3751E02 2.3450E02 2.3154E02 2.2862E02 2.2573E02 1.9890E02 1.7540E02 1.5485E02 1.3691E02 1.2128E02 1.0768E02 9.5839E03 8.5526E03 7.6532E03 6.8671E03 6.1786E03 5.5741E03 5.0419E03 4.5721E03 4.1561E03 3.7869E03 3.4582E03 3.1649E03 2.9023E03

9.2895E01 9.2888E01 9.2864E01 9.2825E01 9.2770E01 9.2698E01 9.2611E01 9.2507E01 9.2388E01 9.2253E01 9.2102E01 9.1937E01 9.1758E01 9.1568E01 9.1369E01 9.1164E01 9.0958E01 9.0758E01 9.0573E01 9.0414E01 9.0298E01 9.0248E01 9.0293E01 9.0479E01 9.0872E01 9.1577E01 9.2775E01 9.4820E01 9.8541E01 1.0676E+00 1.0817E+00 1.0980E+00 1.1170E+00 1.1396E+00 1.1672E+00 1.2020E+00 1.2484E+00 1.3163E+00 1.4373E+00 inﬁnity 1.4373E+00 1.3164E+00 1.2485E+00 1.2021E+00 1.1673E+00 1.1397E+00 1.1171E+00 1.0981E+00 1.0818E+00 1.0677E+00 9.8553E01 9.4827E01 9.2776E01 9.1570E01 9.0857E01 9.0455E01 9.0261E01 9.0208E01 9.0252E01 9.0361E01 9.0515E01 9.0697E01 9.0894E01 9.1099E01 9.1305E01 9.1507E01 9.1702E01 9.1886E01 9.2060E01

780

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783 Table 3 Disk source: c = 0.99, z0 = 1.0, and q0 = 2.0.

0.1 Flux

Diffusion transport

0.01

0.001 0

1

2

ρ

3

4

5

Fig. 2. Flux for a disk source: c = 0.99, z0 = 1.0, and q0 = 2.0.

0.8773. . . and 0.8420. . . for z0 equal to 0.1, 1.0, 2.0 and 4.0, respectively. However, for c = 1 and large q and z0, the ratio does pﬃﬃﬃ tend to the exact value deduced from Eqs. (79) and (110), viz: 3=2. The explanation for this anomalous effect must lie in the fact that there is a long range, ‘deep-penetration’ effect which diffusion theory is unable to describe well. This is not unexpected because we are not considering the conventional ‘deep-penetration’ problem at inﬁnite depth. In this case, the ﬂux is that at the surface of the half-space and, no matter how far away the source, this will always differ from diffusion theory. The hybrid diffusion Eq. (96) is a reasonable attempt at describing this situation but clearly it has limitations due, probably, to the highly directional nature of the source. The numerical results have been obtained by using the formalism of Eq. (166) and also by direct use of Eq. (70) using a Shanks summation method. The results from the two methods agree to the places given in Tables. In all cases the quadratures are carried out with the IMSL library routines which have been thoroughly validated. We consider now the case of a ring source of radius q0 in which all neutrons are moving towards the surface. This is described by setting SðkÞ ¼ J 0 ðkq0 Þ in Eq. (70). Fig. 1 shows the surface ﬂux in the range q(0, 5) when the plane of the ring is at z0 = 1.0, q0 = 3.0 and c = 0.99. Both diffusion and transport theory results are shown using Eqs. (141), (169a) and (169b). The most noticeable effect is the increase in ﬂux near the source. This is not unexpected because, as shown in Appendix C, the integral transient term in Eq. (70) leads to a singularity as q ? q0 of the form log (|q q0|) whereas the diffusion theory case is ﬁnite. We also note that as q0 ? 0, the result tends to that of the point source and proves the existence of the singularity at q = 0 for that case. Table 2 shows the values of the ﬂux for transport and diffusion theory and also the ratio of transport/diffusion. We have noted, although no results are given, that as for the point case larger values of z0 lead to less accurate values in diffusion theory. We also note that as the singularity is approached the numerical effort required for a given accuracy becomes much greater. The other case of some interest is the disk source, sometimes called a ’ﬂashlight’ source. In Fig. 2 we show the surface ﬂux for transport and diffusion theory for c = 0.99, z0 = 1.0 and q0 = 2.0. Eqs. (141), (145), (174a), and (174b) are used. There is no singularity at q = q0 for this case and diffusion theory appears to be reasonably accurate. Table 3 shows numerical results for transport and diffusion theory.

q

Transport

Diffusion

Trans/Diff

0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50 3.60 3.70 3.80 3.90 4.00 4.10 4.20 4.30 4.40 4.50 4.60 4.70 4.80 4.90 5.00

4.8490E02 4.8453E02 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

4.9696E02 4.9654E02 4.9528E02 4.9317E02 4.9017E02 4.8626E02 4.8140E02 4.7553E02 4.6858E02 4.6047E02 4.5111E02 4.4039E02 4.2817E02 4.1430E02 3.9862E02 3.8094E02 3.6103E02 3.3868E02 3.1363E02 2.8561E02 2.5431E02 2.2360E02 1.9710E02 1.7420E02 1.5439E02 1.3723E02 1.2233E02 1.0937E02 9.8064E03 8.8178E03 7.9507E03 7.1880E03 6.5151E03 5.9196E03 5.3912E03 4.9209E03 4.5013E03 4.1257E03 3.7888E03 3.4859E03 3.2128E03 2.9661E03 2.7426E03 2.5399E03 2.3556E03 2.1878E03 2.0345E03 1.8945E03 1.7662E03 1.6486E03 1.5405E03

9.7573E01 9.7580E01 9.7603E01 9.7641E01 9.7696E01 9.7770E01 9.7862E01 9.7977E01 9.8116E01 9.8284E01 9.8483E01 9.8718E01 9.8993E01 9.9312E01 9.9678E01 1.0009E+00 1.0052E+00 1.0094E+00 1.0121E+00 1.0092E+00 9.7587E01 9.3399E01 9.1984E01 9.1274E01 9.0917E01 9.0769E01 9.0751E01 9.0817E01 9.0939E01 9.1094E01 9.1270E01 9.1456E01 9.1643E01 9.1829E01 9.2007E01 9.2177E01 9.2336E01 9.2483E01 9.2618E01 9.2740E01 9.2850E01 9.2947E01 9.3031E01 9.3105E01 9.3167E01 9.3218E01 9.3260E01 9.3292E01 9.3316E01 9.3331E01 9.3340E01

Acknowledgements The author wishes to thank Professors Barry Ganapol and Sudarshan Loyalka for many stimulating discussions about this work and other related matters. We also wish to note that this paper is part of a series of such publications on this topic and future contributions will be made by Barry Ganapol and Sudarshan Loyalka. Appendix A. Derivation of ﬂux for off-centre point The Fourier transform of the collided ﬂux for a general source may be written from Eq. (32) as

Z z0 0 0 ðz0 z00 Þu0 =l I0 ðk; 0jz0 Þ ¼ SðkÞ dz0 Gðkjz 0 Þe l 0

ðA1Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 ðkx cos u0 þ ky sin u0 Þ. The Fourier transwhere u0 ¼ 1 þ i 1 l form of the source is given by

781

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

SðkÞ ¼ Sðkx ; ky Þ ¼

Z

1

dx

1

Z

1

dySðx; yÞeikx xiky y

ðA2Þ

1

In terms of the radial co-ordinates deﬁned in Section 2.1, we have

Sðk; gÞ ¼

Z

1

dqq 0

Z 2p

d#Sðq; #Þeikq cosðg#Þ

ðA3Þ

0

Thus we may write the inverse of Eq. (A1) as

; u0 Þ I0 ðq; #jq0 ; #0 ; z0 ; l 2 Z z0 Z 1 Z 1 1 1 0 0 0 dz0 eðz0 z0 Þ=l dkx dky Sðk; gÞGðkjz ¼ 0Þ l 2p 0 1 1 p ﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ p ﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ 0 0 2 2 eikx ðxðz0 z0 Þ 1l cos u0 =l Þþiky ðyðz0 z0 Þ 1l sin u0 =l Þ ðA4Þ Using the deﬁnitions in Section 2.1 this reduces to

; u0 Þ ¼ I0 ðq; #jq0 ; #0 ; z0 ; l

1

1

2 Z

l 2p

Z 2p 0

z0

0

0

0

dz0 eðz0 z0 Þ=l

Z

1

dkk

0

0 ikR cosðgrÞ dgSðk; gÞGðkjz 0 Þe

ðA5Þ

where R is deﬁned by Eq. (43), viz:

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 cosð# u0 Þ=l þ ð1 l 2 Þðz0 z00 Þ2 =l 2 R2 ¼ q2 2qðz0 z00 Þ 1 l Now for an off-centre point source at (q0, #0), we have from Eqs. (A3) and (A4),

; u0 Þ ¼ I0 ðq; #jq0 ; #0 ; z0 ; l

2 Z z0 Z 1 1 0 0 dz0 eðz0 z0 Þ=l dkk l 2p 0 0 Z 2p 0 ikðRcosð grÞq0 cosðg#0 ÞÞ dgGðkjz 0 Þe 1

I0 ðq; #jq0 ; #0 ; z0 ; u0 Þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Z 1 Z 1 1 w dkkGðkjz Þ dwe J q2 þ q20 2qq0 cosð#0 rÞ ¼ 0 0 k 2p 0 0 ðA13Þ where

q cos # w cos u0 cos r ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ q2 2wq cosð# u0 Þ þ w2 q sin # w sin u0 sin r ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 q 2wq cosð# u0 Þ þ w2

ðA14Þ

Appendix B. An integral equation for the Green’s function It is useful to reduce the basic Green’s function deﬁned by Eq. (28) to the solution of a singular integral equation which is suitable for treatment by the FN method. We do this by introducing the Laplace transform of the spatial ﬂux and showing how it may be related to the emergent angular distribution. The details follow. The surface ﬂux due to any source can be related to some integral over the function W0(k, 0|z0) where W(k, z, w|z0) obeys the equation (see Eq. (29) of the text)

wQ

@W c þ Q 2 W ¼ ðW0 þ dðz z0 ÞÞ 2 @z

where Q ¼

ðB1Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 þ k w2 and

W0 ðk; zjz0 Þ ¼

Z

1

dwWðk; z; wjz0 Þ

ðB2Þ

1

and the boundary condition is

0

ðA6Þ But

R cosðg rÞ q0 cosðg #0 Þ ¼ A cosðg wÞ

ðA7Þ

Wðk; 0; wjz0 Þ ¼ 0;

where

A2 ¼ R2 þ q20 2Rq0 cosð#0 rÞ tan w ¼

ðA8Þ

R sin r q0 sin #0 R cos r q0 cos #0

ðB3Þ

Wðk; z; wjz0 Þ ¼

c 2wQ

Z

z

0

0

dz ðW0 ðk; z0 jz0 Þ þ dðz0 z0 ÞÞeQðzz Þ=w

ðB4Þ

0

and for w < 0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ q cos # ðz0 z00 Þ 1 l 2 cos u0 =l cos r ¼ R

Wðk; z; wjz0 Þ ¼

Z

c 2wQ

1

0

0

dz ðW0 ðk; z0 jz0 Þ þ dðz0 z0 ÞÞeQðz zÞ=w

z

ðA9Þ

ðB5Þ Now for w < 0, we may set z = 0 in (B5) to get

and

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ q sin # ðz0 z00 Þ 1 l 2 sin u0 =l sin r ¼ R

ðA10Þ

Wðk; 0; wjz0 Þ ¼

Thus Eq. (A6) becomes

; u0 Þ ¼ I0 ðq; #jq0 ; #0 ; z0 ; l

w>0

Let us cast Eq. (B1)) into a form suitable for the application of the FN method. By integration of (B1) we may write For w > 0

Z

c 2wQ

1

0

0

dz ðW0 ðk; z0 jz0 Þ þ dðz0 z0 ÞÞeQz =w

ðB6Þ

0

Let us now deﬁne the Laplace transform

Z

1

z0

0

0

ðk; s; wjz0 Þ ¼ W

dz0 eðz0 z0 Þ=l

0 2pl Z 1 0 0 dkkGðkjz 0 ÞJ 0 ðkAðz0 ÞÞ

Z

1

sz

dze

Wðk; z; wjz0 Þ

ðB7Þ

0

ðA11Þ

Applying this to Eq. (B1), we ﬁnd

0

¼ 1, this reduces to For l I0 ðq; #jq0 ; #0 ; z0 Þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Z z0 Z 1 1 0 0 0 dz0 eðz0 z0 Þ dkkGðkjz q2 þ q20 2qq0 cosð# #0 Þ ¼ 0 ÞJ 0 k 2p 0 0

ðk; s; wjz0 Þ wQ Wðk; 0; wjz0 Þ þ Q ðsw þ Q ÞW c sz0 Þ ¼ ðW 0 ðk; sjz0 Þ þ e 2

Dividing by Q(sw + Q) and integrating over w(1, 1), we ﬁnd

ðA12Þ

which is similar to Eq. (51) with a modiﬁed argument of the Bessel function and could have been derived from geometrical consider ¼ 0, we can reduce Eq. (A11) to the convenient form, ations. For l

ðB8Þ

Z

0

w0 Wðk; 0; w0 jz0 Þ þ W0 ðk; sjz0 Þ ðsw þ Q Þ 1 0 ðk; sjz0 Þ þ esz0 Þ ¼ cKðsÞðW 0

dw

ðB9Þ

782

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

where

Z

1 KðsÞ ¼ 2

1

1

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ! 2 dw 1 1 þ s2 k pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ log 2 2 Q ðsw þ Q Þ 2 s2 k 1 s2 k

I0 ðq; 0jz0 Þ ¼

0

0

dw

1

w0 Wðk; 0; w0 jz0 Þ 0 ðk; sjz0 Þ þ esz0 ÞVðsÞ ¼ esz0 þ ðW ðsw0 þ Q 0 Þ

ðB11Þ

ðB12Þ

Now set s = Q/w in Eq. (B11) and

VðQ =wÞWðk; 0; wjz0 Þ ¼

c 2Q

Z

1

dw

0

0

0

0

w Wðk; 0; w jz0 Þ c Qz0 =w þ e 2wQ w0 Q wQ 0 ðB13Þ

This is a singular integral equation for W(k, 0, w|z0) from which we may obtain the scalar surface ﬂux from

W0 ðk; 0; jz0 Þ ¼

Z

0 ðkjz0 Þ dkkJ0 ðkqÞJ 0 ðkq0 ÞG

dwWðk; 0; wjz0 Þ

ðB14Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 w 1þk n¼ Q

ðC1Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 0 ðkjz0 Þ ¼ A0 ðe m2 þk2 z0 ez0 Þ G pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Z 2 2 c 1 dtAðtÞ e 1þk t z0 =t ez0 þ 2 0

ðC2Þ

m2 ð1 m2 Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ A0 ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 m2 þ k Hð1= m2 þ k2 ; kÞðc þ m2 1Þð1 m2 þ k2 Þ

ðC3Þ

gðc; tÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ AðtÞ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 2 1 þ k t2 Hðt= 1 þ k t2 ; kÞðt 1 þ k t2 Þ

ðC4Þ

Now as q ? q0, the dominant term in Eq. (C2) is

c z0 e 2

Z

1

dtAðtÞ

ðC5Þ

0

Inserting (C5) in (C1) we get

I0sing ðq; 0jz0 Þ ¼

Eq. (B13) may be transformed into a more friendly form by introducing the new variable

c z0 e 4p

Z

1

0

dkkJ0 ðkqÞJ 0 ðkq0 Þ

gðc; tÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 2 2 1 þ k t Hðt= 1 þ k t2 ; kÞðt 1 þ k t2 Þ ðC6Þ

ðB15Þ

For large k, the H function becomes unity and Eq. (C6) can be rewritten

it then becomes

c 2

0

1

0

KðnÞUðnÞ þ

1

where

But from Eqs. (B6) and (B7)

c Wðk; 0; wjz0 Þ ¼ ½W0 ðk; Q =wjz0 Þ þ eQz0 =w 2wQ

Z

ðB10Þ

If we set V(s) = 1 cK(s), then Eq. (B9) becomes

Z

1 2p

Z

1

dn0

0

gðn0 Þn0 Uðn0 Þ c ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ez0 2 n n0 2n 1 þ k

pﬃﬃﬃﬃﬃﬃﬃﬃ

1þk2 =n

ðB16Þ

where

I0sing ðq; 0jz0 Þ

c z0 e 4p

Z

1

dtgðc; tÞ

Z

0

1

dkk 0

J 0 ðkqÞJ 0 ðkq0 Þ

ðC7Þ

2

1 t2 þ k t2

Now the integral over k can be carried out to give

UðnÞ ¼ g 2 ðnÞWðk; 0; ngðnÞjz0 Þ

ðB17Þ

1 gðnÞ ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 þ k ð1 n2 Þ

ðB18Þ

I0sing ðq; 0jz0 Þ ( pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Z K 0 ðq 1 t 2 =tÞI0 ðq0 1 t2 =tÞ; c z0 1 dt e gðc; tÞ p ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 4p 0 t K 0 ðq0 1 t 2 =tÞI0 ðq 1 t2 =tÞ;

q > q0 q < q0 ðC8Þ

and

KðnÞ ¼ 1

c 2

Z

1

dn0

1

! gðn0 Þ Q n p ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ ¼ V ¼ V 2 w n n0 1þk

ðB19Þ

Eq. (B16) is amenable to the FN and other methods of solution (Ganapol, 2008a). Solutions of this type of equation have also been described by Williams (1968a,b) by analytic methods. Finally we note that Eqs. (B4) and (B5) may be combined to produce an integral equation for the scalar ﬂux W0(k, z|z0) in the form

W0 ðk; zjz0 Þ ¼

c 2

Z

1

0.07

0.06

0

dz Kðjz z0 j; kÞ½W0 ðk; z0 jz0 Þ þ dðz0 z0 Þ

Without singularity + limit term limit term transport without singularity

ðB20Þ

0

0.05

where

KðzÞ ¼

We note that this expression is asymmetric about q = q0 as expected. Eq. (C8) applies very near the singularity q = q0. Indeed, it is easy to show that when q is very close to q0

Z

1

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

dw pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ e 2 w 1 þ k w2

1þk2 w2 z=w

without singularity+limit term

ðB21Þ

0.04

Eq. (B20) is amenable to solution by the methods described by Naz and Loyalka (2008).

0.03

0

transport

limit term

Appendix C. An estimate of the ﬂux in the neighbourhood of ring and point sources

0.02

C.1. Values of q leading to singularities

0.01 2.90

In Section 3 of the paper we have from Eq. (70) the ring source case, if we set SðkÞ ¼ J 0 ðkq0 Þ, then

without singularity

2.95

3.00

ρ

3.05

3.10

Fig. 3. Flux for a ring source in the neighbourhood of the singularirty: c = 0.99, q0 = 3.0, and z0 = 1.0.

783

M.M.R. Williams / Annals of Nuclear Energy 36 (2009) 767–783

c z0 1 e pﬃﬃﬃﬃﬃﬃﬃﬃﬃ E1 ðjq q0 jÞ 4p 2 qq0 c z0 1 1 e log 4p 2q0 jq q0 j

But in all cases, SðkÞ can be written as

I0sing

SðkÞ ¼ 1 þ a1 k2 þ a2 k4 þ ðC9Þ

It is useful to see how Eq. (C8) blends in with the complete solution (C2). What we have done is to remove the singular component from Eq. (C2) so that it reads

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Z 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 nonsing ðkjz0 Þ ¼ A0 e m2 þk2 z0 ez0 þ c dtAðtÞe 1þk t z0 =t G 2 0 ðC10Þ This is inserted into Eq. (C1) and solved numerically and, as expected, no problems are encountered at q = q0. Then we evaluate Eq. (C8) and add it to the surface ﬂux arising from Eq. (C10). Fig. 3 shows the results over the range 2.90–3.10 with q0 = 3. The non-singular component is almost a constant background and the singular term of Eq. (C8) is an asymmetric spike. When added, they show how the true solution is approached. A special case of Eq. (C8) is that of the point source as discussed in Section 2 of the text where SðkÞ ¼ 1. This case may be obtained from Eq. (C8) if q0 = 0, when we ﬁnd

I0sing ðq; 0jz0 Þ

c z0 e 4p

Z

1 0

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dt gðc; tÞK ð q 1 t2 =tÞ 0 t2

ðC11Þ

Bearing in mind that the dominant part of the integrand of Eq. (C11) is when t ? 0, we may write

I0sing ðq; 0jz0 Þ

c z0 e 4p

Z

1 0

dt c z0 K 0 ðq=tÞ ¼ e 4pq t2

Z

1

dsK 0 ðsÞ

q

ðC12Þ Hence when q ? 0, using

R1 0

dsK 0 ðsÞ ¼ p=2, we ﬁnd

c 8

qI0sing ðq; 0jz0 Þ ¼ ez0 lim q!0

ðC13Þ

This result is consistent with the numerical results obtained from Section 9 and for c = 0.99 and z0 = 1.0, the limiting value is 0.045525.... It is curious that the ring singularity changes from a logarithmic one to an inverse power one as q0 ? 0, but this is born out by the numerical results as well as the analysis. C.2. Behaviour of the ﬂux at large values of q from the source Consider the general expression for the ﬂux for an arbitrary source as deﬁned by Eq. (70) of the text, viz:

I0 ðq; 0jz0 Þ ¼

1 2p

Z

1

0

0 ðkjz0 Þ dkkJ0 ðkqÞSðkÞG

ðC14Þ

Now set kq = w, whence

I0 ðq; 0jz0 Þ ¼

Z

1 2pq

2

0

1

0 ðw=qjz0 Þ dwwJ0 ðwÞSðw=qÞG

ðC15Þ

ðC16Þ

So Eq. (C15) becomes

I0 ðq; 0jz0 Þ ¼

1 2pq

Z

2

0

1

a1 w2 a2 w4 0 ðw=qjz0 Þ dwwJ0 ðwÞ 1 þ 2 þ 4 þ G

q

q

ðC17Þ hence the leading term is independent of the source shape and so the asymptotic behaviour, as q ? 1, is also independent of source shape. As seen in Section 4 of the text, Eqs. (79) and (83), the case for c = 1 leads to an asymptotic behaviour of O(1/q3). References Cassell, J.S., Williams, M.M.R., 2007a. Radiation transport and internal reﬂection in a two region, turbid sphere. J. Quant. Spectrosc. Rad. Transfer 104, 400. Cassell, J.S., Williams, M.M.R., 2007b. An integral equation for radiative transport in an inﬁnite cylinder with Fresnel reﬂection. J. Quant. Spectrosc. Rad. Transfer 105, 12. Davison, B., 1957. Neutron Transport Theory. Oxford University Press. Davis, P.J., Rabinowitz, P., 1967. Numerical Integration. Blaisdell Publishing Co., NY. Elliott, J.P., 1952. The transport theory of a point source of neutrons on the free surface of a half-space. Atomic Energy Research Establishment, Harwell Report AERE T/R 972 (United Kingdom, Ministry of Supply). Elliott, J.P., 1955. Milne’s problem with a point source. Proc. Roy. Soc. A228, 424– 433. Ganapol, B.D., Kornreich, D.E., 2008 (Interior beam searchlight semi-analytical benchmark). In: International Conference on the Physics of Reactors ‘‘Nuclear Power: a Sustainable Source”. Casino-Kursaal Conference Centre, Interlaken, Switzerland. Ganapol, B.D., 2008a (Analytical Benchmarks for Nuclear Engineering Applications: Case Studies in Neutron Transport Theory). OECD, NEA, Paris. Ganapol, B.D., 2008b. Private communication. Gradshteyn, I.S., Ryzhik, I.M., 1994. Table of Integrals Series and Products, ﬁfth ed. Academic Press. Grosjean, C.C., 1965. On the series expansion of certain types of Fourier integrals in the neighbourhood of the origin. Bull. Soc. Math. Belg. XVII, 251. Lam, S.K., Leonard, A., 1971. Milne’s problem for two-dimensional transport in a quarter space. J. Quant. Spectrosc. Rad. Transfer 11, 893. Loyalka, S.K., 2008. Private communication. Loyalka, S.K., Williams, M.M.R., in press. Computation of some transport theory expressions for a line source in an inﬁnite half-space. Ann. Nucl. Energy. Naz, S., Loyalka, S.K., 2008. One speed criticality problems for a bare slab and sphere: some benchmark results II. Ann. Nucl. Energy 35, 2426. Noble, B., 1958. Methods Based on the Wiener–Hopf Technique. Pergamon Press, Oxford. Rybicki, G., 1971. The searchlight problem with isotropic scattering. J. Quant. Spectrosc. Rad. Transfer 11, 827. Williams, M.M.R., 1968a. Diffusion length and criticality problems in two and three dimensional one-speed neutron transport theory. I. Rectangular coordinates. J. Math. Phys. 9, 1873. Williams, M.M.R., 1968b. Diffusion length and criticality problems in two and three dimensional, one-speed neutron transport theory. II. Circular cylindrical coordinates. J. Math. Phys. 9, 1885. Williams, M.M.R., 1971. Mathematical Methods in Particle Transport Theory. Butterworths, London. Williams, M.M.R., 1982. The three-dimensional transport equation with applications to energy deposition and reﬂection. J. Math. Phys. A: Math. Gen. 15, 965. Williams, M.M.R., 2007a. The transport and diffusion theory of a line source in an inﬁnite half-space with internal reﬂection. Ann. Nucl. Energy 34, 910. Williams, M.M.R., 2007b. The searchlight problem in radiative transfer with internal reﬂection. J. Phys. A: Math. Theory 40, 6407.

Three-dimensional transport theory: An analytical solution of an internal beam searchlight problem-I

Three-dimensional transport theory: An analytical solution of an internal beam searchlight problem-I

Recommend Documents