Solution of the neutron transport problem with anisotropic scattering in cylindrical geometry by the decomposition method

Annals of Nuclear Energy 36 (2009) 98–102

Contents lists available at ScienceDirect

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

Solution of the neutron transport problem with anisotropic scattering in cylindrical geometry by the decomposition method G.A. Gonçalves a, S.Q. Bogado Leite b,*, M.T. de Vilhena c a

UFRGS, Departamento de Engenharia Nuclear, Av. Osvaldo Aranha 99, 4° andar, Porto Alegre, RS 90046-900, Brazil Comissão Nacional de Energia Nuclear, Coordenação Geral de Reatores e Ciclo do Combustível, Rua General Severiano, 90, Rio de Janeiro, RJ 22294-900, Brazil c UFRGS, Departamento de Matemática Aplicada, Av. Bento Gonçalves, 9500, Porto Alegre, RS 91509-900, Brazil b

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 10 June 2008 Accepted 1 October 2008 Available online 21 November 2008

An analytical solution has been obtained for the one-speed stationary neutron transport problem, in an infinitely long cylinder with anisotropic scattering by the decomposition method. Series expansions of the angular flux distribution are proposed in terms of suitably constructed functions, recursively obtainable from the isotropic solution, to take into account anisotropy. As for the isotropic problem, an accurate closed-form solution was chosen for the problem with internal source and constant incident radiation, obtained from an integral transformation technique and the FN method. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction

where

Consider the one-speed stationary neutron transport equation with anisotropic scattering, in a bare uniform cylinder of radius R and infinite length

fl ¼ 2p

Lwðr; l; uÞ ¼ c

Z 2p Z 0

Z

1

f ðl0 ÞPl ðl0 Þdl0

ð2bÞ

with the normalization condition

1

1 f ðX0 ! XÞwðr; l0 ; u0 Þdl0 du0 þ ð1  cÞQ ðrÞ 4 p 1 ð1aÞ

1

f0 ¼ 2p

Z

1 1

f ðl0 Þdl0 ¼ 1

ð2cÞ

and using the addition theorem

subject to

1 wðR; l; uÞ ¼ G for  1 6 l 6 1 and 4p

p=2 6 u 6 3p=2;

ð1bÞ

Pl ðl0 Þ ¼ Pl ðlÞPl ðl0 Þ þ 2

l X ðl  mÞ! m 0 0 P ðlÞPm l ðl Þ cos mðu  u Þ; ðl þ mÞ! l m¼1

ð3Þ

where the invertible differential operator L is defined as



   o 1 o þ 1 wðr; l; uÞ: Lwðr; l; uÞ ¼ ð1  l2 Þ1=2 cos u  sin u or r ou ð1cÞ In the above, w(r; l, u) is the angular flux, c = (rs + mrf)/rt is the mean number of neutrons emerging from a collision, Q(r) denotes the internal source term, r is the optical distance, u and h = cos1 l are angles the neutron direction of movement makes with the radial and axial axes, respectively, as shown in Fig. 1, f(X0 ? X), assumed a function of X0 X = l0 only, is the neutron transfer function, and G is a constant. Upon expanding f in terms of Legendre polynomials

f ðX0 ! XÞ ¼ f ðl0 Þ ¼

1 4p

1 X

ð2l þ 1Þfl Pl ðl0 Þ;

ð2aÞ

an equivalent expression for Eq. (1a) is obtained as

L1 Lwðr; l; uÞ ¼ wðr; l; uÞ   Z 2p Z 1 c 1 wðr; l0 ; u0 Þdl0 du0 þ ð1  cÞQ ðrÞ ¼ L1 4p 0 4p 1 Z 2p Z 1 X 1 c F l ðl0 ; u0 ! l; uÞ þ L1 4p 0 1 l¼1  wðr; l0 ; u0 Þdl0 du0 with

F l ðl0 ; u0 ! l; uÞ ¼ ð2l þ 1Þfl "  Pl ðlÞPl ðl0 Þ þ 2

l¼0

# 0

* Corresponding author. Tel.: +55 21 21732344. E-mail address: [email protected] (S.Q. Bogado Leite). 0306-4549/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2008.10.002

ð4aÞ

 cos mðu  u Þ :

l X ðl  mÞ! m 0 P ðlÞPm l ðl Þ ðl þ mÞ! l m¼1

ð4bÞ

99

G.A. Gonçalves et al. / Annals of Nuclear Energy 36 (2009) 98–102

and consequently we can associate

z

Z 2p Z

c 4p

w1 ðr; l; uÞ ¼ L1

0

c 4p

wi ðr; l; uÞr ¼ L1

1 X

1

F l ðl0 ; u0 ! l; uÞw0 ðr; l0 ; u0 Þdl0 du0 ;

ð8aÞ

1 l¼1

Z 2p Z 0

1

1 X

F l ðl0 ; u0 ! l; uÞwi1 ðr; l0 ; u0 Þdl0 du0 ;

1 l¼0

i ¼ 2; 3; . . .

ð8bÞ

subject to

wi ðR; l; uÞ ¼ 0 for  1 6 l 6 1 and

y

0

p=2 6 u 6 3p=2:

ð8cÞ

If the series converges, the n-term partial sum of wi will be an approximate solution of w. Eqs. (8) state that wi, for i = 1, 2, . . ., can be recursively obtained from wi1, with w0(r; l, u) and w0(r) the angular flux and the neutron density, respectively, of the isotropic problem with internal sources and prescribed incident radiation. An analytical solution for the later has been proposed by Siewert and Thomas (1984) using the FN method and an integral transformation technique developed by Mitsis (1963), as summarized below.

r P x Fig. 1. Notation used – infinite cylinder.

2. A solution for the isotropic problem Following the decomposition method (Adomian, 1988), we now let

wðr; l; uÞ ¼

1 X

wn ðr; l; uÞ;

ð5Þ

In obtaining a solution for w0(r) from Eqs. (7), Siewert and Thomas first define I(R; l, u) as

n¼0

1 X

wn ðr; l; uÞ

n¼0

 Z 2p Z 1 c 1 ¼ L1 w0 ðr; l0 ; u0 Þdl0 du0 þ ð1  cÞQ ðrÞ 4p 0 4p 1 Z 2p Z 1 X 1 1 c 0 0 0 þL wn ðr; l ; u Þdl du0 4p 0 1 n¼1 Z 2p Z 1 X 1 c þ L1 F l ðl0 ; u0 ! l; uÞ 4p 0 1 l¼1 1 X

wn ðr; l0 ; u0 Þdl0 du0 :

w0 ðrÞ ¼

0

ð7aÞ

1 0

0

0

w0 ðr; l ; u Þdl du

p=2 6 u 6 3p=2;

ð7cÞ

þL

þ L1 

1 X n¼0

c 4p c 4p

0

1

Iðr; l; uÞdl du

0

1

w0 ðr; l; uÞdl du ¼ IðrÞ þ G:

ð11Þ

1

IðrÞ ¼ cK½IðrÞ

ð12aÞ

Z

Z

1

0 1

ð10bÞ

1

Z 2p Z

ðKf ÞðrÞ ¼

Z 2p Z

1

with K an integral operator defined by

wn ðr; l; uÞ ¼ w0 ðr; l; uÞ 1

0

In case of Q = G = 0, Mitsis has shown that the integral form of the transport equation can be written as

1 G for  1 6 l 6 1 and 4p

n¼0

1

ð7bÞ

1

then 1 X

Z 2p Z

and we see that

w0 ðrÞ ¼

subject to

w0 ðR; l; uÞ ¼

Z 2p Z 0

where 0

Iðr; l0 ; u0 Þdl0 du0

1

where

IðrÞ ¼

Z 2p Z

1

G 1 G dl0 du0 þ ð1  cÞQ ðrÞ  L 4p 4p 4p c ð1  cÞðQ  GÞ ¼ IðrÞ þ ; ð10aÞ 4p 4p

If we impose

w0 ðr; l; uÞ ¼ L1

0

c 4p

þ

ð6Þ

 c 1 w0 ðrÞ þ ð1  cÞQ ðrÞ ; 4p 4p

Z 2p Z

c 4p

LIðr; l; uÞ ¼

n¼0



ð9Þ

so that I(r; l, u) = 0. The transport equation, assuming isotropic scattering, can thus be written as





1 G; 4p

w0 ðr; l; uÞ ¼ Iðr; l; uÞ þ

so that Eq. (4a) becomes

1 X

Z 2p

1 n¼1 Z 1 X 1

0

1 l¼1

0

0

0

þ

0

wn ðr; l ; u Þdl du

tf ðtÞK 0 ðr=lÞI0 ðt=lÞdt

0 R

tf ðtÞK 0 ðt=lÞI0 ðr=lÞdt

r



dl

l2

;

ð12bÞ

where I0 and K0 are modified Bessel functions of first and second kind, respectively. For a problem with internal sources and incident radiation, one can readily see that I(r) satisfies

F l ðl0 ; u0 ! l; uÞ

wn ðr; l0 ; u0 Þdl0 du0

Z

r

ð7dÞ

IðrÞ ¼ K fcIðrÞ þ ð1  cÞ½Q ðrÞ  Gg:

ð13Þ

100

G.A. Gonçalves et al. / Annals of Nuclear Energy 36 (2009) 98–102

Defining U(r, l) as

Uðr; lÞ ¼

Z 0

þ

E0 ðnÞ ¼ 1  c;

tfcIðtÞ þ ð1  cÞ½Q ðtÞ  GgK 0 ðr=lÞI0 ðt=lÞdt Z

R

t fcIðtÞ þ ð1  cÞ½Q ðtÞ  GgK 0 ðt=lÞI0 ðr=lÞdt

ð14Þ

r

we see that

IðrÞ ¼

Z

1

Uðr; lÞ 0

dl

l

2

ð15Þ

;

so that I(r) and U(r, l) form a transform pair, with U(r, l) being interpreted as a pseudo-neutron distribution taking the place of w0(r; l, u). Upon differentiating Eq. (14) with respect to r and using the definition of U(r, l), it can be shown that U(r, l) satisfies

! o2 1 o 1 þ Uðr; lÞ þ cIðrÞ þ ð1  cÞ½Q ðrÞ  G ¼ 0  or 2 r or l2

ð16aÞ

E1 ðtÞ ¼ t  c=2  ct lnð1 þ 1=tÞ; 0 6 t 6 1;    c 1 E1 ðt0 Þ ¼  1 þ t20 ln 1  2 ; t0 > 1; 2 t0 Ea ðnÞ ¼ n2 Ea2 ðnÞ  c=ða þ 1Þ; Z 1 lZðR=lÞ D0 ðt0 Þ ¼ ct0 dl; t20  l2 0 Z 1 2 l ZðR=lÞ dl; D1 ðt0 Þ ¼ ct0 t20  l2 0 Z 1 2 t ZðR=tÞ  lZðR=lÞ D1 ðtÞ ¼ ZðR=tÞ  ct dl; t2  l2 0 0 6 t 6 1; D1 ðtÞ ¼ tZðR=tÞ  ct

o lK 0 ðR=lÞ Uðr; lÞjr¼R þ K 1 ðR=lÞUðR; lÞ ¼ 0: or

ð16bÞ

Uðr; lÞ ¼ l2 fAðt0 Þ½/ðt0 ; lÞ þ /ðt0 ; lÞI0 ðr=t0 Þ Z

1

 AðtÞ½/ðt; lÞ þ /ðt; lÞI0 ðr=tÞdt þ UP ðr; lÞ; ð17Þ

0

where

c 2

/ðt0 ; lÞ ¼



1

t0

t0 t0  l



t



tl

1

þ ð1  cttanh

tÞdðl  tÞ

" # N X n aa Ea ðnÞ  ð1  cÞðQ  GÞ ; NðnÞI0 ðR=nÞ a¼0

ð18cÞ

 1  2 ;

c 3 c t 2 0 t20  1 t0 h i 1 NðtÞ ¼ t ð1  cttanh tÞ2 þ ðpct=2Þ2 ;

t ZðR=tÞ  l ZðR=lÞ dl; t2  l2

Z

1

2

2

Da ðnÞ ¼ n2 Da2 ðnÞ  cn

ð19nÞ

la1 ZðR=lÞdl

ð19oÞ

with

ZðxÞ ¼ K 1 ðxÞ=K 0 ðxÞ:

ð19pÞ

In addition, for constant Q, a particular solution is

UP ðr; lÞ ¼ l2 ðQ  GÞ:

ð20Þ

Z

1

AðtÞI0 ðr=tÞdt þ Q:

ð21Þ

0

  I0 ðR=nb Þ aa Ea ðnb Þ þ Da ðnb Þ ¼ ð1  cÞðQ  GÞ I1 ðR=nb Þ a¼0

N X

Once the isotropic solution has been found, the higher order terms in the expansion of Eq. (5) are determined from Eqs. (8) and from the fact that, in an infinite cylinder, the directional deriv~ w may be written as b r ative X

~ w ¼ sin h b r X

b r ~ w þ w ¼ Zðr; l; uÞ becomes so that X

ð19bÞ ð19cÞ

ð19dÞ

pffiffiffiffiffiffiffi 1l2 r cos u

w¼e

ð19eÞ b ¼ 1; 2; . . . ; N

ð19fÞ

and Ea(n) and Da(n) are functions discussed elsewhere (Thomas and Southers, 1983). For c < 1, these functions are

Z

1 pffiffiffiffiffiffiffi 2 e 1l Zð1; r sin u; lÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d1 ¼ L1 Z; 1  l2

ð22bÞ

ð22cÞ

which defines the inverse operator L1. We now substitute I0 by a polynomial approximation in the form (Beyer, 1978)

I0 ðxÞ ¼ 1 þ   1 1 2b  1 p ; þ cos 2 2 2N

ð22aÞ

with w  wð1; r sin u; lÞ; Z  Zð1; r sin u; lÞ and 1 ¼ r cos u. The solution for the above equation, unless for an integration constant necessary for satisfying the boundary condition, is

with

n0 ¼ t0 ;

ow ; oðr cos uÞ

ow 1 Z þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi w ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi oðr cos uÞ 1  l2 1  l2

the coefficients aa are found from

nb ¼

ð19lÞ

ð19aÞ

where



ð19kÞ

3. Approximate solution for the linearly anisotropic problem

ð18bÞ

with P denoting the Cauchy principal value of the integral involving the term in parenthesis, UP is a particular solution of Eqs. (16a) and the expansion coefficients are to be determined from the boundary conditions. This has been done by Siewert and Thomas for constant Q, using the FN method, yielding

Nðt0 Þ ¼

0

w0 ðrÞ ¼ Aðt0 ÞI0 ðr=t0 Þ þ

ðt0 real for c < 1 and purely imaginary for c > 1Þ;

AðnÞ ¼

1

ð18aÞ

¼1

c /ðt; lÞ ¼ P 2

ð19jÞ

We now insert Eq. (20) into Eq. (17) and this into Eq. (15) and use Eq. (11) to obtain the desired solution for w0(r) as



with t0 the positive root of 1

ð19iÞ

0

The solution to Eqs. (16) which remains bounded at the origin, by separation of variables, is

þ

ð19hÞ

ð19mÞ Z

0 6 t 6 1;

subject to

ct0 tanh

ð19gÞ 2

r

1 X

x2m

m¼1

2 ðm!Þ2 2m

;

ð23Þ

so that the isotropic neutron density Eq. (21) becomes

w0 ðrÞ ffi c0 þ

M X m¼1

cm r 2m þ Q

ð24Þ

101

G.A. Gonçalves et al. / Annals of Nuclear Energy 36 (2009) 98–102

with

cm ¼

Aðt0 Þ ð2t0 Þ2m ðm!Þ2

þ

Z 0

1

AðtÞ ð2tÞ2m ðm!Þ2

dt;

m ¼ 0; 1; 2; . . .

ð25Þ

and where we have truncated the series at some value of m = M. Then, using Eqs. (22c) and (24) in Eq. (7a), we obtain 2 3 1 pffiffiffiffiffiffiffi r cos u M pffiffiffiffiffiffiffi Z 1l2 X c 6 e 7 2 2 w0 ðr; l; uÞ ffi cm e 1l ð12 þ r2 sin uÞm pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d15 4c0 þ 4p 1  l2 m¼1

þ

Q ; 4p ð26aÞ

Table 1 Normalized neutron density for c = 0.9, M = 5 and 200 recursive terms in the expansion of Eq. (5). r (mfp)

cf1 = 0.0

cf1 = 0.1

cf1 = 0.3

cf1 = 0.5

cf1 = 0.7

0 1 2 3 4 5 6 7 8 9 10

0.116370 0.116201 0.115660 0.114634 0.112892 0.110049 0.105488 0.098175 0.086135 0.065159 0.025981

0.117764 0.117559 0.116904 0.115673 0.113625 0.110355 0.105198 0.097067 0.084101 0.062999 0.027741

0.120195 0.119902 0.118978 0.117287 0.114576 0.110433 0.104221 0.094965 0.081169 0.060518 0.029380

0.121645 0.121257 0.120053 0.117907 0.114589 0.109746 0.102851 0.093148 0.079555 0.060528 0.033860

0.121856 0.121375 0.119900 0.117340 0.113531 0.108233 0.101106 0.091687 0.079364 0.063321 0.042490

whose solution is

" M m X X c m!ðr sin uÞ2ðmiÞ c0 þ cm 4p i!ðm  iÞ! m¼1 i¼0 # 2ij qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2i X ð2iÞ!ðr cos uÞ Q  ð1Þj : ð 1  l2 Þ j þ 4 ð2i  jÞ! p j¼0

w0 ðr; l; uÞ ffi

Table 2 Neutron densities for increasing values of M (cf1 = 0.7 and 200 recursive terms).

ð26bÞ

Keeping only terms in l = 0 and l = 1 in the summations in Eqs. (8), we obtain

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3c w1 ðr; l; uÞ ¼ L1 f1 1  l2 cos u 4p Z 2p Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  l02 cos u0 w0 ðr; l0 ; u0 Þdl0 du0 ;  0

1

0

1

ð27aÞ

Z 2p Z 1 c wi1 ðr; l0 ; u0 Þdl0 du0 wi ðr; l; uÞ ¼ L1 4p 0 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3c þ L1 f1 1  l2 cos u 4p Z 2p Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  l02 cos u0 wi1 ðr; l0 ; u0 Þdl0 du0 ;  i ¼ 2; 3; . . .

ð27bÞ

and where we have used the fact that wi(r; l, u) are even functions of l and u. Applying w0(r;l0 , u0 ) in Eq. (27a) and integrating over u0 all terms with the exception of those in odd powers of r vanish. The resulting function Zð1; r sin u; lÞ over which the L1 operator is ap2 plied is, thus, in the form Zð1; r sin u; lÞ / 1ð12 þ r 2 sin uÞn . Similar arguments are used in Eq. (27b) to obtain wi(r; l, u). The solutions to Eqs. (27) are obtained using the software Maple 9.5. 4. Numerical results for a test case and conclusion We now consider the problem of an infinite cylinder of 10 mfp radius, with constant source and a free surface boundary condition. The value of Q in this example (Q = 0.11877) was chosen in such a way as to normalize the isotropic neutron density. Using M = 5 in the polynomial expansion in Eq. (24) and 200 recursive terms in the summation in Eq. (5), we obtain the linearly anisotropic angular flux w(r; l, u) from Eqs. (26) and (27). Table 1 shows values of the corresponding neutron density, w(r), for c = 0.9, at several positions inside the cylinder and for several values of cf1. Neutron densities for cf1 = 0.0 are found in good agreement with the isotropic solution of Siewert and Thomas, indicating convergence of the polynomial expansion in Eq. (24) for M = 5. The convergence of the solution with increasing values of M is illustrated in Table 2 for cf1 = 0.7. Table 3 illustrates the convergence of the neutron density with respect to the number of recursive terms in the summation in Eq. (5), for cf1 = 0.7 and M = 7. A solution has been proposed for the one-speed stationary neutron transport problem with anisotropic scattering in an infinitely long cylindrical system by the decomposition method. Series

r (mfp)

M=5

M=7

M=9

0 1 2 3 4 5 6 7 8 9 10

0.121856 0.121375 0.119900 0.117340 0.113531 0.108233 0.101106 0.091687 0.079364 0.063321 0.042490

0.121862 0.121380 0.119905 0.117344 0.113535 0.108236 0.101106 0.091685 0.079358 0.063311 0.042476

0.121859 0.121377 0.119903 0.117343 0.113534 0.108236 0.101107 0.091687 0.079360 0.063313 0.042477

Table 3 Neutron density as function of the number of terms in the partial sum of Eq. (5) (cf1 = 0.7 and M = 7). r (mfp)

100 terms

200 terms

400 terms

0 1 2 3 4 5 6 7 8 9 10

0.121588 0.121119 0.119684 0.117186 0.113460 0.108255 0.101221 0.091882 0.079600 0.063531 0.042562

0.121862 0.121380 0.119905 0.117344 0.113535 0.108236 0.101106 0.091685 0.079358 0.063311 0.042476

0.121862 0.121380 0.119906 0.117344 0.113535 0.108236 0.101106 0.091685 0.079358 0.063311 0.042476

expansions of the angular flux have been obtained recursively from the isotropic solution using an integral transformation technique and the FN method, for cylinders with constant source and incident radiation. The results for a linearly anisotropic test case appear consistent and fast convergent with the number of recursive terms taken in the expansions. Although applied to a linearly anisotropic problem, the developed methodology is suitable to problems to any degree of anisotropy as well. It is equally applicable to cylinders with general source and boundary conditions whose isotropic solution has been obtained by other methods. This work proceeds focused in these more general cases. Acknowledgments The second and third authors are gratefully indebted to Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) of Brazil, for partial financial support to this work.

102

G.A. Gonçalves et al. / Annals of Nuclear Energy 36 (2009) 98–102

References Adomian, G., 1988. A review of the decomposition method in applied mathematics. J. Math. Anal. Appl. 135, 501–544. Beyer, W.H. (Ed.), 1978. CRC Standard Mathematical Tables, 25th ed. CRC Press, Inc., Boca Raton, Florida.

Mitsis, G.J., 1963. Transport Solutions to the Monoenergetic Critical Problems. ANL6787, Argonne National Laboratory. Siewert, C.E., Thomas, J.R., 1984. Neutron transport calculations in cylindrical geometry. Nucl. Sci. Eng. 87, 107–112. Thomas, J.R., Southers, J.D., 1983. The critical problem for an infinite cylinder. Nucl. Sci. Eng. 84, 79–82.