A nodal collocation approximation for the multi-dimensional PL equations – 2D applications

Annals of Nuclear Energy 35 (2008) 1820–1830

Contents lists available at ScienceDirect

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

A nodal collocation approximation for the multi-dimensional P L equations – 2D applications M. Capilla a, C.F. Talavera a, D. Ginestar a, G. Verdú b,* a b

Departamento de Matemática Aplicada, Universidad Politécnica de Valencia, Camino de Vera 14, E-46022 Valencia, Spain Departamento de Ingeniería Química y Nuclear, Universidad Politécnica de Valencia, Camino de Vera 14, E-46022 Valencia, Spain

a r t i c l e

i n f o

Article history: Received 18 November 2007 Received in revised form 28 April 2008 Accepted 29 April 2008 Available online 9 July 2008

a b s t r a c t A classical approach to solve the neutron transport equation is to apply the spherical harmonics method obtaining a finite approximation known as the PL equations. In this work, the derivation of the PL equations for multi-dimensional geometries is reviewed and a nodal collocation method is developed to discretize these equations on a rectangular mesh based on the expansion of the neutronic fluxes in terms of orthogonal Legendre polynomials. The performance of the method and the dominant transport Lambda Modes are obtained for a homogeneous 2D problem, a heterogeneous 2D anisotropic scattering problem, a heterogeneous 2D problem and a benchmark problem corresponding to a MOX fuel reactor core. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction For reactor calculations, the neutron multi-group transport problem is usually treated using the classical diffusion approximation, but this approximation does not provide good results for complex fuel assemblies and other applications as fine mesh (pin by pin) geometry. Neutron transport equation is generally discretized using the SL method for the angular dependency, but this method exhibits ray-effects problems (Alcouffe et al., 1979). To eliminate ray-effects, the PL approximation to the transport equation is used. This approximation is well established (Davison, 1957; Weinberg and Wigner, 1958; Clark and Hansen, 1964) and is obtained by expanding the angular dependence of the neutron flux and the nuclear cross-sections in terms of spherical harmonics up to a finiteorder L. The exact transport solution is recovered as L goes to infinity. In three-dimensional geometries (3D) the number of PL equations is LðL þ 1Þ=2. These equations are complicated and can also be formulated as second-order equations, but the coupling involves not only the angular moments but also mixed spatial derivatives of these moments. This problem led to propose the simplified P L (SP L ) approximation in the early 1960s (Gelbard, 1968). For the SP L approximation, it is proposed to replace the second derivatives in the one-dimensional planar geometry P L equations by three-dimensional Laplacian operators avoiding the complexity of the full spherical harmonics approximation. In the 90s theoretical basis for the observed accuracy of the SP L equations in the multi-dimensional applications were provided * Corresponding author. Tel.: +34 963877630; fax: +34 963877639. E-mail addresses: [email protected] (M. Capilla), [email protected] (C.F. Talavera), [email protected] (D. Ginestar), [email protected] (G. Verdú). 0306-4549/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2008.04.008

(Brantley and Larsen, 2000), showing that these equations are high-order asymptotic solutions of the transport equation when diffusion theory is the leading-order approximation (Larsen et al., 1996). In this way, the SP L equations are a practical way to solve transport problems in physical systems that can be considered as nearly plane and where the neutrons flow along the direction orthogonal to the plane. But when it is applied to general 3D transport problems the accuracy of the solution cannot be increased with larger-order L of approximation and it is possible to obtain worse results with the SP L approximation than with the diffusion approximation (Tomasevic´ and Larsen, 1996). This fact motivated the use of the SP L approximation as a preconditioning procedure and originated the development of new approximations as the AL  SP 2L1 approximation and the mixed P1  DP 0 diffusion theory (Brantley and Larsen, 2000). The best solution is obtained from the P L approximation, but its implementation in a computer program is complicated. Fletcher (1983) obtained a solution of the P L equations eliminating the fields with odd-order and discretizing the resulting equations using the finite differences method and the finite elements method. In this work, we review the multi-dimensional PL equations and develop a nodal collocation method for the discretization of these equations, which is based on the expansion of the fields in terms of orthonormal Legendre polynomials. Particularly, we will focus on the eigenvalue problem known as the Lambda modes transport problem. The nodal collocation method approximates the initial differential eigenvalue problem by a generalized algebraic eigenvalue problem, from which the k-effective and the stationary neutron flux distribution of the system can be computed, being able also to obtain the subcritical eigenvalues and their corresponding eigenmodes. The method presented here generalizes the method for 1D geometries presented in a previous work (Capilla et al., 2005) to be able to treat

1821

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

multi-dimensional problems. To test the performance of the method, a homogeneous 2D eigenvalue problem, a heterogeneous 2D anisotropic scattering eigenvalue problem, a heterogeneous 2D eigenvalue problem with strong spatial heterogeneity and a MOX fuel benchmark problem have been studied.

2. The transport equation and the P L equations The time independent transport equation in terms of the directional flux distribution is given by Stacey (2001)

~r ~ Uð~ ~; EÞ þ Rt ð~ ~; EÞ X r; X r; EÞUð~ r; X Z

Z

~ 0 ; E0 ! X ~; EÞUð~ ~0 ; E0 Þ ~0 Rs ð~ r; X r; X dX Z vp ðEÞ Z 0 ~0 Uð~ ~0 ; E0 Þ þ Sð~ ~; EÞ; þ dE mRf ð~ r; E0 Þ dX r; X r; X 4p ¼

dE0

ð1Þ

~; EÞ is the neutronic flux at location ~ where Uð~ r; X r with energy E and ~; Rt is the total cross-section; Rs direction given by the unit vector X ~0 ; E0 Þ to ðX ~; EÞ; Rf is the fission is the scattering cross-section from ðX cross-section; m is the average number of neutrons per fission, vp is ~; EÞ is a neutron source. We will concentrate the spectrum and Sð~ r; X on the development of a nodal collocation method for the transport equation eigenvalue problem

~r ~ Uð~ ~; EÞ þ Rt ð~ ~; EÞ X r; X r; EÞUð~ r; X Z

Z

~ 0 ; E0 ! X ~; EÞUð~ ~0 ; E0 Þ ~0 Rs ð~ r; X r; X dX Z Z 1 vp ðEÞ ~0 Uð~ ~0 ; E0 Þ þ dE0 mRf ð~ r; E0 Þ dX r; X k 4p ¼

dE0

ð2Þ

known as the Lambda modes transport problem. Nevertheless, the same kind of method is also valid for other kind of problems such as stationary problems with neutron sources described by Eq. (1). In practical applications, to eliminate the dependence of energy in Eq. (1) or (2), an energy multi-group approximation is used. In order to facilitate the notation we will consider the monoenergetic version of these equations. For the extension of the nodal colloca~Þ is replaced by a coltion method to the multi-group case, Uð~ r; X umn vector of components corresponding to each energy group, Rt is a diagonal matrix and Rs , mRf are also adequate matrices. The spherical harmonics method to approximate the transport ~Þ, as equation consists of expanding the angular neutron flux, Uð~ r; X

~Þ ¼ Uð~ r; X

1 X þl X

~ m ~ /m l ðrÞY l ð Þ;

X

ð3Þ

l¼0 m¼l

~ Ym l ð Þ

where X are the (normalized) spherical harmonics (Gradshteyn and Ryzhik, 1980). It is also assumed that scattering depends only ~X ~0 , and the scattering cross-section on the relative angle, l0 ¼ X may be expanded as the following series of Legendre polynomials:

Rs ð~ r; l0 Þ ¼ ¼

1 X 2l þ 1 Rsl ð~ rÞPl ðl0 Þ 4p l¼0 1 X þl X



~ m ~0 Rsl ð~ rÞY m l ðXÞY l ðX Þ:

ð4Þ

l¼0 m¼l

Usually, if x, y and z are the Cartesian coordinates, the following ‘‘complex” coordinates are introduced:

1 xþ1 ¼  pffiffiffi ðx þ iyÞ; 2 since the identity

~¼ ~r X

rffiffiffiffiffiffiffi þ1 4p X k o Y 3 k¼1 1 oxk

1 x1 ¼ pffiffiffi ðx  iyÞ; 2

x0 ¼ z;

ð5Þ

is satisfied. Substituting Eqs. (3) and (4) in Eq. (2) and integrating R ~ Y m  ðX ~Þ results in the following equations for with respect to dX l the coefficients or moments /m l :

 1=2 mþ1 o/lþ1 1 ðl þ m þ 2Þðl þ m þ 1Þ  pffiffiffi ð2l þ 3Þð2l þ 1Þ ox1 2  1=2 m1 ! o/lþ1 ðl  m þ 2Þðl  m þ 1Þ þ ð2l þ 3Þð2l þ 1Þ oxþ1  1=2 mþ1 1 ðl  mÞðl  m  1Þ o/l1 þ pffiffiffi ð2l þ 1Þð2l  1Þ ox1 2  1=2 m1 ! ðl þ mÞðl þ m  1Þ o/l1 þ ð2l þ 1Þð2l  1Þ oxþ1  1=2 m o/lþ1 ðl þ m þ 1Þðl  m þ 1Þ þ ð2l þ 3Þð2l þ 1Þ ox0  1=2 m ðl þ mÞðl  mÞ o/l1 þ ð2l þ 1Þð2l  1Þ ox0 1 m m þ Rt /l ¼ Rsl /l þ d0l d0m mRf /00 ; l ¼ 0; 1; . . . ; m ¼ l; . . . ; þl; k ð6Þ where d0l is the Kronecker delta. In these equations it is understood that terms involving coefficients /m l with invalid indices l and m are zero. To obtain a finite approximation, the series in Eq. (3) is truncated at some value l ¼ L and the approximate Eq. (6) obtained are known as the P L equations. In the following we will only consider L to be an odd integer. Using the Eq. (6) with odd index l, the following relations between odd l coefficients /m l and derivatives of even l coefficients can be obtained:

/m l ¼ 

þ1 þ1 o/mr Dl X Dl X o/mr Aþr;lm ð1Þr lþ1   Ar;lm l1 þ oxr C l r¼1 oxr C l r¼1

ð7Þ

with l ¼ 1; 3; . . . ; L, where Dl ¼ 1=ðRt  Rsl Þ and the following constants have been introduced:

1 Aþ1;lm ¼ pffiffiffi ½ðl þ m þ 2Þðl þ m þ 1Þ1=2 ; 2 Aþ0;lm ¼ ½ðl þ m þ 1Þðl  m þ 1Þ1=2 ; 1 A1;l;m ¼ pffiffiffi ½ðl  mÞðl  m  1Þ1=2 ; 2 A0;lm ¼ ½ðl þ mÞðl  mÞ1=2 ; C þl ¼ ½ð2l þ 3Þð2l þ 1Þ1=2 ;

Aþþ1;lm ¼ Aþ1;l;m ;

Aþ1;lm ¼ A1;l;m ;

ð8Þ

C l ¼ C þl1 :

~Þ, is a real Eq. (7) are complex, but the angular neutronic flux, Uð~ r; X  m m function and the complex moments satisfy /m ¼ ð1Þ /l . Introl ducing the real moments, defined as

 m 1 m nm ¼ ð/l þ ð1Þm /m l ¼ Re /l l Þ; 2   1 gml ¼ Im /ml ¼ ð/ml  ð1Þm /m l Þ; 2i the following set of real relations is obtained from Eq. (7):

nm l

! ogmþ1 onm 1 1 onmþ1 þ lþ1 lþ1 þ Aþ0;lm lþ1 ¼ Dl þ A1;lm pffiffiffi  ox oy oz Cl 2 !! ogm1 1 onm1 lþ1 þAþþ1;lm pffiffiffi  lþ1 ox oy 2 ! mþ1 ogmþ1 1 1 onl1 onm þ  A1;lm pffiffiffi  l1 þ A0;lm l1 Cl ox oy oz 2

ð9Þ

1822

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

ogm1 1 onm1 l1 Aþ1;lm pffiffiffi þ l1 ox oy 2

!!! ;

The equivalent set of real conditions is

m ¼ 0; 1; . . . ; l;

! ogm onmþ1 1 1 ogmþ1 þ lþ1 lþ1 m gl ¼ Dl þ A1;lm pffiffiffi þ þ Aþ0;lm lþ1 ox oy oz Cl 2 !! onm1 1 ogm1 lþ1 þAþþ1;lm pffiffiffi  lþ1 ox oy 2 ! ogm 1 1 ogmþ1 onmþ1 l1 þ l1 þ A0;lm l1 þ  A1;lm pffiffiffi Cl ox oy oz 2 !!! m1 m1 1 ogl1 on ; m ¼ 1; 2; . . . ; l: Aþ1;lm pffiffiffi  l1 ox oy 2

1 m 1 m   ðB þ Bm ðB þ ð1Þm Bm Þ ¼ 0; l Þ ¼ l 2 l 2 l 1 1   ðBm  Bm ðBm  ð1Þm Bm Þ ¼ 0; l Þ ¼ l 2i l 2i l

that results in the following expressions for the real moments nm l and gm l 0

L1 X l X

n0l þ

0

nm l þ

Em l ¼ 0;

gml þ

ð11Þ

l0 ¼0 l even

B þ BM 0 n00 þ @ lm;l 0 l

L1 X l0 X

l0 X m0 ¼1 m0 6¼m

0C C 2Mþlm;l0 m0 nm l0 A ¼ 0;

ð17Þ

0

2M lm;l0 m0 gm ¼0 l0

l0 ¼0 m0 ¼1 6 m l0 even m0 ¼ þ1 X

Em l ¼ 

ð1Þrþr

0

r;r 0 ¼1

for l ¼ 1; 3; . . . ; L; m ¼ 1; 2; . . . ; l, where

Aþr;lm Aþr0 ;lþ1;mr o o mrr0 Dlþ1 / oxr oxr0 lþ2 C þl C þlþ1

0

m Mlm;l0 m0 ¼ Hm l H l0 2

Aþr;lm Ar0 ;lþ1;mr o o mrr0 Dlþ1 /l þ  ox ox C C r r0 0 l lþ1 r;r ¼1 !  þ o mrr0 r 0 Ar;lm Ar0 ;l1;mr o þ ð1Þ Dl1 / oxr oxr0 l C l C þl1 þ1 X



ð1Þr

2

Z 0

1

  sin ððm  m0 Þp=2Þ sin ððm þ m0 Þp=2Þ  m  m0 m þ m0 0

m Pm l ðxÞP l0 ðxÞdx;

m 6¼ m0 ;

ð18Þ

and

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2l þ 1 ðl  mÞ! ¼ : 4p ðl þ mÞ!

þ1 X Ar;lm Ar0 ;l1;mr o o mrr0 Dl1 /l2    ox ox C C r r0 l l1 r;r 0 ¼1

Hm l

0 þ ðRt  Rsl Þ/m l  dl0 dm0 mRf /0 :

In the above equations odd l components are related with even l components. Vacuum boundary conditions involving only even l components are then derived by substituting Eq. (10) in Eq. (17). It is convenient to define the vector X of odd l components nm l , m gl and the vector X of even l components. X has Lþ1 ðL þ 2Þ compo2 nents and X has Lþ1 L components. 2 For an external surface normal to the left side of axis X, Eq. (10) can be rewritten in matrix form as

Eq. (11) are complex, but the diffusive equations for the real mom ments nm l and gl , are given by m m Em l ¼ 0; 2; 4; . . . ; L  1; m ¼ 0; 1; . . . ; l; l þ ð1Þ El ¼ 0; m Þ ¼ 0; l ¼ 2; 4; . . . ; L  1; m ¼ 1; . . . ; l: iðEl  ð1Þm Em l

ð12Þ

As an example, the detailed PL equations for XY geometry are written in Appendix. 2.1. Boundary conditions

~Þ ¼ 0 for X ~~ Uð~ r; X n 6 0;

ð13Þ

where ~ n is the outwardly directed unitary normal vector to the external surface. This condition can be approximated by setting Marshak’s conditions (Stacey, 2001) 

~ ~~ ~ Ym l ðXÞUðr; XÞdX ¼ 0;

o X; ox

X ¼ DN 1

ð19Þ

where N 1 is a rectangular matrix and

We will consider vacuum boundary conditions. That is, the incoming neutronic flux is considered to be zero on an external surface of the region where the P L equations are considered. This is written as

~~ n60 X

1

0

L1 X 0

where

Z

0

2M þl0;l0 m0 nm l0 ¼ 0;

m0 ¼1

l ¼0 l0 even

ð10Þ

These relations are generalizations of the ‘‘Fick’s law”. Replacing Eq. (7) into Eq. (6) with even l, the following diffusive equations for even l coefficients, /m l , result

l ¼ 0; 2; 4; . . . ; L  1;

ð16Þ

l ¼ 1; 3; 5; . . . ; L ðoddÞ:

ð14Þ

D ¼ diag



 1 ; Rt  Rsl

l ¼ 1; 3; 5; . . .

Vacuum boundary conditions (17) in terms of the vectors X and X are written as the matrix expression

X þ MX X ¼ 0;

ð20Þ

X

where M is a rectangular matrix. On eliminating the odd l moments there results

MX X  DN 1

o X ¼ 0: ox

ð21Þ

To show how these conditions can be imposed, let us consider, for example, the boundary surface normal to the left side of axis X and the P L approximation (similar relations are obtained for other external surfaces). For this surface Marshak’s conditions are written as the complex conditions

In the same way, for an external surface normal to the right side of axis X a similar expression to Eq. (21) is obtained, changing the sign in the sum. For external surfaces normal to the left side of axes Y and Z relation (19) is expressed, respectively, as

Bm l ¼ 0;

X ¼ DN 2

l ¼ 1; 3; . . . ; L;

ð15Þ

where

Bm l ¼

Z p=2 p=2

d/

Z p 0

0



~ sinðhÞdhY m l ðXÞ

L l X X 0

l ¼0

m0 ¼l0

0

0

m ~ /m l0 Y l0 ðXÞ:

o X oy

and X ¼ DN3

o X; oz

ð22Þ

where N 2 and N3 are rectangular matrices. Analogous relations to Eq. (17) can be obtained for surfaces normal to Y and Z axes. The corresponding equations are shown in Appendix.

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

If the problem is symmetric with respect to a plane, reflective boundary conditions can be imposed, and odd l components /m l are zero at the symmetry plane.

Since PL equations (11) have a diffusive form, the discretization of these equations can be done using a nodal collocation method, previously used for the neutron diffusion equation (Hébert, 1987; Verdú et al., 1994), generalizing the method presented in Capilla et al. (2005) for multi-dimensional rectangular geometries. The first step to discretize the P L equations is to divide the region where these equations have to be solved into N prismatic nodes of the form

h i h i h i Ne ¼ xm1 ; xmþ1  yn1 ; ynþ1  zp1 ; zpþ1 : 2

2

2

Z 1=2 M X 1 du1 Pr1 ðu1 ÞP0k1 ðu1 Þ Dxe Dye k ;k ¼0 1=2 1 2 Z 1=2 du2 Pr2 ðu2 ÞP0k2 ðu2 Þxe;k1 k2 r3 : 

2

2

In Fig. 1, a scheme of the discretization for a generic node e and the convention used to label the six neighbouring elements that share a surface with element e is shown. For a generic node e the change of variables

  1 1 x xm1 þ xmþ1 ; 2 2 2 Dxe   1 1 yn1 þ ynþ1 ; y u2 ¼ 2 2 2 Dye   1 1 zp1 þ zpþ1 ; z u3 ¼ 2 2 2 Dz e

In the first case, integration is straightforward using the orthonormality properties of Pk ðuÞ. In case (ii), let us first consider node e to be an interior node. Adjacent nodes are then related imposing continuity conditions ~Þ, that implies continuity of all components /m for the flux Uð~ r; X l in expansion (3). The six neighbouring nodes which share surface with node e are labelled according to Fig. 1. For example, if the boundary between two adjacent nodes falls in the YZ plane, see Fig. 2, we consider the interface between node e and node e1 . Conm 1 1 tinuity conditions are then /m l;e ðu1 ¼  2Þ ¼ /l;e1 ðu1 ¼ þ 2Þ, l ¼ 0; 1; . . . ; L, that results in the expressions

    1 1 X e  ; u2 ; u3 ¼ X e1 þ ; u2 ; u3 ; 2 2     1 oX e 1 1 oX e1 1  ; u2 ; u3 ¼ De1 N1 þ ; u2 ; u3 ;  De N 1 Dxe ou1 Dxe1 ou1 2 2

u1 ¼

ð29Þ

ð23Þ

where Dxe ¼ xmþ1  xm1 , Dye ¼ ynþ1  yn1 , Dze ¼ zpþ1  zp1 , trans2 2 2 2 2 2 forms the node e into the cube of volume one







1 1 1 1 1 1 Nu;e ¼  ;   ;   ; : 2 2 2 2 2 2

ð24Þ

The nodal collocation method assumes that on each node the nuclear cross-sections are constant. Eq. (12) are transformed by means of the change of variables (23) for each node e. Furthermore, if X e ðu1 ; u2 ; u3 Þ denotes the previously defined vector of even l real m components nm l and gl that appear in Eq. (12), it is assumed that vector X e can be expanded in terms of orthonormal Legendre polynomials Pk ðuÞ (Capilla et al., 2005) up to a certain fixed order M, M X

X e ðu1 ; u2 ; u3 Þ ¼

where odd-order fields are removed using Eq. (19) (‘‘Fick’s laws”). On replacing series (25) into Eq. (29) and integrating with respect R 1=2 R 1=2 to 1=2 du2 1=2 du3 Pr2 ðu2 ÞPr3 ðu3 Þ there results equivalent conditions for X e;r2 r3 ðu1 Þ:

    1 1 ¼ X e1 ;r2 r3 þ ; X e;r2 r3  2 2     1 dX e;r2 r3 1 1 dX e1 ;r2 r3 1 ¼ De1 N1  þ : De N1 2 2 Dxe du1 Dxe1 du1

xe;k1 k2 k3 Pk1 ðu1 ÞPk2 ðu2 ÞPk2 ðu3 Þ;



1 1 uj 2  ; : 2 2

ð25Þ

The series (25) is then inserted into Eq. (12) and equations for the coefficients xe;k1 k2 k3 are derived multiplying by the weight function

r 1 ; r2 ; r 3 ¼ 0; 1; . . . ; M

and integrating over the cube (24). In performing this process three kind of terms appear:

Fig. 1. Position of the neighbouring elements of node e.

(i) ‘‘diagonal terms”, with no derivatives. (ii) Terms with second derivatives, as for example,

F e;x;r1 r2 r3 ¼ 

Dl;e Dx2e

Z

1=2

du1 Pr1 ðu1 Þ

1=2

o2 X e;r2 r3 ðu1 Þ; ou21

ð26Þ

where X e;r2 r3 ðu1 Þ is

X e;k2 k3 ðu1 Þ ¼

M X

xe;k1 k2 k3 Pk1 ðu1 Þ

ð27Þ

k1 ¼0

(Analogous terms appear for the other spatial directions). (iii) Terms with cross derivatives, as for example,

ð30Þ

Eqs. (26) and (30) describe a one-dimensional problem for variable u1 . Continuity conditions (29) are then rewritten as (Capilla et al., 2005)

k1 ;k2 ;k3 ¼0

W r1 r2 r3 ¼ Pr1 ðu1 ÞPr2 ðu2 ÞPr3 ðu3 Þ;

ð28Þ

1=2

3. The nodal collocation method

2

Ge;x;y;r1 r2 r3 ¼

1823

Fig. 2. Boundary between two adjacent nodes.

1824

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

  X W e;x M1 dX e;r2 r3 1 ¼ MðM þ 1ÞDxe D1  ck ½ð1Þk xe;kr2 r3  xe1 ;kr2 r3 ; e 2 du1 2 k¼0   M 1 X W e;x 1 ¼ ck ð1  Dxe D1 Þð1Þk xe;kr2 r3 X e;r2 r3  e 2 2 k¼0

W e;x þDxe D1 ð31Þ x e1 ;kr 2 r 3 ; e 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi kðkþ1Þ where the coefficient ck ¼ 2k þ 1ð1  MðMþ1Þ Þ and the ‘‘coupling  factor” W e;x is the square matrix solution of

W e;x N1 D1 e 2

¼ diag

De1 N1 : Dxe De1 þ Dxe1 De

ð32Þ

This procedure is repeated for node e2 adjacent to node e in forward X direction, and a ‘‘coupling factor” W þ e;x is defined as in Eq. (32) replacing e by e2 and e1 by e. Notice that Eq. (32) simplifies for isoDe De W tropic scattering to the diagonal form 2e;x ¼ diag Dxe De þD1xe De and the 1 1 components of X e;r2 r3 in Eq. (31) are decoupled. From Eq. (31) the computation of integral (26) is the vector analogue of the procedure described in Capilla et al. (2005), and results in the expression

F e;x;r1 r2 r3 ¼

M 1 h X

D ;M

D ;M

D ;M



D ;M

W e;x MðM þ 1Þ ck cr ð1Þr ; D xe 2

! W e;x W þe;x MðM þ 1Þ ck cr ð1Þrþk þ Dxe 2 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dl;e þ 2r þ 1 2k þ 1ð1 þ ð1Þrþk Þ 2 Dxe 8  rðrþ1Þ > < 1  MðMþ1Þ kðk þ 1Þ if k < r;    kðkþ1Þ > : 1  MðMþ1Þ rðr þ 1Þ if k P r;

ð34Þ

¼ D1xe De Fx and Fx is the

  1 MX þ MðM þ 1Þ De N1 Fx ¼ M X : Dxe

ð39Þ

On comparing Eq. (38) with Eq. (31), it results that Eq. (33) is valid for all nodes taking into consideration that, for boundary nodes, the coupling factors are X e;x and that there is no adjacent node in the direction of the boundary. A similar procedure is applied to axis Y and Z. Lastly, in case (iii), cross derivatives are readily computed using the relation

Z

1=2

duPr ðuÞP0k ðuÞ ¼



0; k < r þ 1; pffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r þ 1 2k þ 1ð1  ð1Þrþ1 Þ;

W e;x MðM þ 1Þ ¼ ck cr ð1Þk : D xe 2

ð35Þ

and then the following relations occur:

M X

Ge;x;y;r1 r2 r3 ¼

ð36Þ

ð37Þ

M X

Ae;x;y;r1 k1 ;r2 k2 xe;k1 k2 r3 ;

ð41Þ

1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r 1 þ 1 2k1 þ 1ð1  ð1Þr1 þk1 Þ Dxe Dy e pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2r 2 þ 1 2k2 þ 1ð1  ð1Þr2 þk2 Þ:

ð42Þ

k1 ¼r 1 þ1 k2 ¼r 2 þ1

where

Finally, once an appropriate ordering of the indices is chosen, the previous procedure approximates Eq. (12) by a generalized eigenvalue problem

AX ¼

If the node is adjacent to the vacuum boundary as, for example, in the case shown in Fig. 3, then Marshak vacuum boundary conditions (21) are used. After performing the change of variable (23) and integrating, there results

    1 1 dX e;r2 r3 1  ¼0 M X X e;r2 r3  De N1  2 2 Dxe du1

k P r þ 1:

For example, term (28) results in the expression

Ae;x;y;r1 k1 ;r2 k2 ¼

þ

D ;M Ae2l ;r1 ;e;k

2

ð38Þ

ð40Þ

where the indices e1 and e2 correspond to adjacent nodes in X direction and

Ae;rl 1 ;e;k ¼

X e;x

Ae;rl 1 ;e1 ;k xe1 ;kr2 r3 þ Ae;rl 1 ;e;k xe;kr2 r3 þ Ae2l ;r1 ;e;k xe2 ;kr2 r3 ; ð33Þ

D ;M

where the ‘‘boundary coupling factor” is solution of

1=2

i

k¼0

Ae;rl 1 ;e1 ;k ¼ 

  X Xe;x M1 dX e;r2 r3 1 ¼ MðM þ 1ÞDxe D1  ck ð1Þk xe;kr2 r3 ; e 2 du1 2 k¼0   M   1 X Xe;x 1 ¼ ck 1  Dxe D1 ð1Þk xe;kr2 r3 ; X e;r2 r3  e 2 2 k¼0

1 BX; k

ð43Þ

m where the vector X has components nm l;e;k1 k2 k3 , gl;e;k1 k2 k3 and A, B are matrices of dimension

N  G  NLeg  N even ;

ð44Þ

where N is the number of nodes, G is the number of energy groups, NLeg ¼ Md is the number of Legendre coefficients, with M the order in Legendre series (25) and d the spatial dimension and finally Neven ¼ Lþ1 ðL þ 2Þ is the number of even l component, being L the or2 der of the PL approximation. To increase numerical efficiency, a serendipity approximation can be consistently introduced in the nodal collocation method (Hébert, 1987). This approximation reduces m the number of unknowns to the set fnm l;e;k1 k2 k3 ; gl;e;k1 k2 k3 ; k1 ¼ 0; . . . ; M  1; k2 ¼ 0; . . . ; M  1  k1 ; k3 ¼ 0; . . . ; M  1  k1  k2 g and replaces Eqs. (33) and (41) by

F e;x;r1 r2 r3 ¼

M1k X2 k3

D ;Mk k3

ðAe;rl 1 ;e1 ;k2

xe1 ;kr2 r3

k¼0 D ;Mk k3

þ Ae;rl 1 ;e;k 2 Ge;x;y;r1 r2 r3 ¼

M1 X

D ;Mk k3

xe;kr2 r3 þ Ae2l ;r1 ;e;k2

M1k X1

xe2 ;kr2 r3 Þ;

Ae;x;y;r1 k1 ;r2 k2 xe;k1 k2 r3 :

ð45Þ

k1 ¼r 1 þ1 k2 ¼r 2 þ1

Fig. 3. Node with vacuum boundary.

The dimension of the eigenvalue problem (43) changes in Eq. (44) because the number of Legendre coefficients N Leg is now

1825

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

NLeg

8 for 1D; > < M; for 2D; ¼ 12 MðM þ 1Þ; > :1 MðM þ 1ÞðM þ 2Þ; for 3D: 6

6.5

ð46Þ

P1 P3 P5 S16

6

4. Numerical results The nodal collocation method developed above has been implemented into a computer code written in FORTRAN 90, which solves the eigenvalue problem (43) for an arbitrary PL approximation, with L odd. In this section, we first consider a homogeneous one-group and isotropic scattering eigenvalue problem. Second, we study a heterogeneous one-group eigenvalue problem with anisotropic scattering. Third, we consider a heterogeneous 2D eigenvalue problem with strong spatial heterogeneity. Finally, we will treat a more realistic MOX reactor problem. In all problems, we compare the results obtained with the nodal collocation method to results from the discrete-ordinates code TWODANT (Alcouffe et al., 1995), which solves the transport problem using a fine spatial mesh and an angular ðSN Þ quadrature set.

We consider a homogeneous square region of width 2 cm and height 2 cm, with vacuum boundary conditions. The nuclear cross-sections for this problem are Rt ¼ 1 cm1 , Rs ¼ 0:9 cm1 and mRf ¼ 0:25 cm1 . In Table 1, we show the fundamental eigenvalue and those corresponding to the subcritical modes. We have taken the parameters ðN x  N y ; MÞ ¼ ð8  8; 3Þ (N x and N y are the number of intervals in each direction, and M is the orthonormal Legendre polynomials order in Eq. (25)). Also we compare the first eigenvalue obtained using the successive P L approximations with a reference solution from the discrete-ordinates TWODANT code (Alcouffe et al., 1995) using the S16 quadrature set, with 150  150 spatial mesh and a convergence criterion of epsi ¼ 107 . In Fig. 4 we plot the S16 scalar flux distribution and the P L solutions with L ¼ 1; 3; 5 obtained with the nodal collocation method, along the horizontal line y ¼ 1 cm. In all the calculations, the normalization condition is N X G X

V core

e¼1 g¼1

j/00;eg j

5 4.5 4 3.5 3 2.5 2

0

0.5

1

1.5

2

x (cm) Fig. 4. P 1 , P 3 , P 5 and S16 scalar flux along y ¼ 1 cm for the homogeneous eigenvalue problem.

problem it is necessary to consider a large L in the PL approximation and either a large number of nodes or a large order M in the polynomial expansion to obtain reasonable results.

4.1. Homogeneous eigenvalue problem

1

Scalar Flux

5.5

mRf ;eg DV e ¼ 1;

where V core is the core area, N is the number of nodes (N ¼ N x  Ny ), G is the number of energy groups and DV e ¼ Dxe Dye for node e. We observe from this figure that, since the homogeneous region considered has small dimensions, the angular dependence in the solution is large and it is necessary to consider a large order L in the PL equations to obtain an accurate solution for the transport problem. To investigate the convergence of the nodal collocation method, we have calculated the fundamental eigenvalue keff for different values of ðN x  N y ; MÞ and for successive P L approximations for odd L. The obtained results are shown in Table 2. These results were computed without making use of the serendipity approximation because the problem size is small. We observe that for this

4.2. Heterogeneous eigenvalue problem with anisotropic scattering The second application is a heterogeneous one-group eigenvalue problem with anisotropic scattering cross-sections (Hébert, 2006). The two-dimensional geometry benchmark is described in Fig. 5, where units are in cm, and the cross-sections data for the three materials are given in Table 3. The domain is surrounded by a vacuum boundary condition. In Table 4 we present the four dominant eigenvalues calculated with the nodal collocation method for the P1 , P 3 , P 5 and P7 approx-

Table 2 Eigenvalue keff of P L equations for a homogeneous square region, computed for several values of ðN x  N y ; MÞ ðN x  N y ; MÞ

P1

P3

P5

P7

P9

P 11

ð1  1; 3Þ ð1  1; 4Þ ð1  1; 5Þ ð1  1; 6Þ ð1  1; 7Þ ð2  2; 3Þ ð2  2; 4Þ ð2  2; 5Þ ð4  4; 3Þ ð4  4; 4Þ ð8  8; 3Þ

0.3331 0.3328 0.3328 0.3328 0.3328 0.3328 0.3328 0.3328 0.3328 0.3328 0.3328

0.3866 0.3801 0.3766 0.3764 0.3763 0.3802 0.3766 0.3763 0.3776 0.3763 0.3766

0.3923 0.3857 0.3821 0.3817 0.3815 0.3854 0.3819 0.3815 0.3830 0.3814 0.3818

0.3935 0.3868 0.3832 0.3828 0.3826 0.3865 0.3830 0.3825 0.3841 0.3824 0.3829

0.3939 0.3873 0.3837 0.3832 0.3830 0.3870 0.3834 0.3829 0.3845 0.3828 0.3833

0.3941 0.3875 0.3839 0.3834 0.3832 0.3872 0.3836 0.3831 0.3847 0.3830 0.3835

Table 1 Four dominant eigenvalues for the homogeneous one-group eigenvalue problem Eigenvalue

P1

P3

P5

P7

P9

P 11

S16

keff 2nd 3rd 4th

0.3328 0.1248 0.1248 0.0768

0.3766 0.1683 0.1683 0.1054

0.3818 0.1788 0.1788 0.1150

0.3829 0.1812 0.1812 0.1181

0.3833 0.1819 0.1819 0.1192

0.3835 0.1822 0.1822 0.1196

0.3891

Fig. 5. Description of the two-dimensional geometry.

1826

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

Table 3 Material cross-sections for the two-dimensional benchmark Material

Rt ðcm1 Þ

Rs0 ðcm1 Þ

Rs1 ðcm1 Þ

mRf ðcm1 Þ

1 2 3

0.025 0.025 0.075

0.013 0.024 0.000

0.000 0.002 0.000

0.0155 0.000 0.000

Table 4 Dominant eigenvalues for the two-dimensional benchmark Eigenvalue

P1

P3

P5

P7

S16

keff 2nd 3rd 4th

0.9644 0.6392 0.6218 0.4374

0.9925 0.7126 0.6996 0.5379

0.9937 0.7190 0.7066 0.5519

0.9940 0.7202 0.7078 0.5545

0.9940

100

P1 P3 S16

90 80

Scalar Flux

70 60 50 40 30 20 10 0 200

250

300

350

400

x (cm) Fig. 6. P 1 , P 3 and S16 scalar flux along y ¼ 160 cm for the two-dimensional benchmark (200 < x < 400 cm).

imations. To perform the P L calculations, we have considered that the two-dimensional domain is subdivided into 10  8 squared nodes, with dimensions of 40  40 cm, and the Legendre polynomial order of M ¼ 6. Also, the serendipity approximation is used. The PL results are compared against the calculation with TWODANT for the right-superior quarter of the reactor, obtained with 800  640 mesh points, S16 quadrature set and a convergence criterion of 107. In Fig. 6 we display the scalar flux distributions obtained from the P1 , P 3 and S16 calculations along the central line y ¼ 160 cm. The P3 curve shows very good agreement with the reference S16 solution. Table 5 shows the variation of the keff eigenvalue as M (the Legendre polynomials order) increases. Again, we have considered one node per squared region with side length of 40 cm. When calculations are performed with a finer mesh, the same results are obtained with lower Legendre polynomial order. In Hébert (2006), a

divergence of the power iteration is observed at all orders of the spherical harmonics approximations. This fact forced to replace the vacuum boundary condition by a zero-flux boundary condition. In this situation our method provides a keff of 0.9936 with size mesh 40  40 cm, M order 6 and the P 5 approximation. The reference value provided by the author, keff ¼ 0:991902  76:9 pcm, was taken as the average of a mixed-primal and mixed-dual P5 approximation with analytical integration and 6  6 mesh splitting. We can conclude that our results are very satisfactory. We must to remark that with our method we can reach converged results with finer mesh using lower order for M. 4.3. Heterogeneous eigenvalue problem We consider now a one-group, isotropic-scattering problem presented in Brantley and Larsen (2000). The system is comprised of two materials as shown in Fig. 7 (units are in cm), with vacuum boundary conditions. The materials cross-sections for this problem are given in Table 6. This system is optically thin and possesses strong spatial heterogeneity. Table 7 contains the four dominant eigenvalues obtained from the nodal collocation method for the P L approximations with L ¼ 1; 3; 5; 7; 9. We have considered 20  20 nodes and M ¼ 5 polynomials in the nodal expansion and the serendipity approximation. We have also included the reference transport solution obtained using TWODANT for the right-superior quarter of the reactor in Fig. 7. To obtain the reference solution we have used the S16 quadrature set and a 150  150 spatial mesh. This problem was solved using the SP 3 approximation in Brantley and Larsen (2000). The relative error of the first eigenvalue (k eff ) was 0.932%. The P 3 result decreases this error to 0.831%. In Figs. 8 and 9 we plot PL with L ¼ 1; 3; 5 and S16 scalar flux distributions along the horizontal lines y ¼ 14:5 cm and y ¼ 18:9 cm, respectively. In both cases, the P3 and P5 fluxes are more accurate than the diffusion flux. Along the y ¼ 14:5 cm, the problem has a one-dimensional character. The system becomes increasingly two-dimensional in character as y increases. When the material interfaces become more two-dimensional (or multi-dimensional), the assumptions used in the derivation of the SP 3 approximations become less reasonable, as it was already pointed out in Brantley

Fig. 7. Heterogeneous eigenvalue problem geometry. Table 5 Eigenvalue keff computed for several values of M and 40  40 cm nodes M

P1

P3

P5

P7

3 4 5 6

0.9655 0.9645 0.9644 0.9644

0.9910 0.9923 0.9924 0.9925

0.9918 0.9936 0.9936 0.9937

0.9919 0.9939 0.9939 0.9940

Table 6 Heterogeneous problem material cross-sections Material

Rt ðcm1 Þ

Rs ðcm1 Þ

mRf ðcm1 Þ

M F

1.0 1.5

0.93 1.35

0.0 0.24

1827

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830 Table 7 Computed eigenvalues for heterogeneous problem

8

P1

P3

P5

P7

P9

S16

6

keff 2nd 3rd 4th

0.7768 0.6718 0.6665 0.5904

0.7994 0.6973 0.6941 0.6184

0.8033 0.7012 0.6983 0.6226

0.8041 0.7020 0.6992 0.6235

0.8043 0.7023 0.6995 0.6237

0.8061

4

7

2nd Eigenvector

Eigenvalue

P1 P3 P5 S16

6

P1 P3 P5

2 0 -2 -4 -6 -8

4

2

4

6

8

10

12

14

16

18

20

18

20

x (cm)

3

9

2

8

1

7

0 10

12

14

16

18

20

x (cm) Fig. 8. Scalar flux for P L approximations and S16 reference result, along y ¼ 14:5 cm for the heterogeneous one-group problem (10 6 x 6 20 cm).

3rd Eigenvector

Scalar Flux

5

P1 P3 P5

6 5 4 3 2 1

1.6

P1 P3 P5 S16

1.4

0

4

6

8

10

12

14

16

x (cm)

1.2

Scalar Flux

2

Fig. 10. Eigenvectors for the subcritical modes of the heterogeneous eigenvalue problem (y ¼ 14:5 cm, 0 6 x 6 20 cm).

1 0.8 0.6 0.4 0.2 0 10

12

14

16

18

20

x (cm) Fig. 9. Scalar flux for P L and S16 along y ¼ 18:9 cm for the heterogeneous one-group problem (10 6 x 6 20 cm).

and Larsen (2000). However, in this situation the P3 and P5 approximations give reasonable accurate results (see Fig. 9). In Fig. 10 we show the normalized eigenvectors obtained for the second and third subcritical modes, using the P 1 , P 3 and P5 approximations. We observe that to obtain a converged result for these subcritical modes it is also necessary to consider higher-order approximations than the neutron diffusion equation.

Fig. 11. MOX benchmark problem geometry.

4.4. MOX benchmark problem The fourth problem we have considered is a modification of the MOX benchmark problem presented in Brantley and Larsen (2000). We consider the application of the PL approximation to a more ‘‘realistic” problem with mixed-oxide fuels. The core configuration

is composed of 7  7 fuel assemblies of two types (MOX=UO2 ) arranged as it is shown in Fig. 11. The core is surrounded by a reflector. Each assembly has dimensions of 21:42  21:42 cm. The UO2 fuel assembly has the cross-sections of material ‘‘U” in Cavarec et al. (1994), and each MOX fuel assembly has the cross-sections of

1828

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

Table 8 Two-group assembly cross-sections (cm1 ) for MOX benchmark problem, with g ¼ 1 for the fast energy group and g ¼ 2 for the thermal energy group Assembly type

Group

Rt

mRf

Rs;1!g

Rs;2!g

MOX fuel

1 2 1 2 1 2

0.550 1.060 0.570 1.100 0.611 2.340

0.0075 0.450 0.005 0.125 0.000 0.000

0.520 0.015 0.540 0.020 0.560 0.050

– 0.760 – 1.000 – 2.300

UO2 fuel Reflector

171.4 2.21

1.38

0.553

100 –0.276

Table 9 Dominant eigenvalues for MOX benchmark problem

–1.11

Eigenvalue

P1

P3

S16

keff 2nd 3rd 4th

0.9929 0.9667 0.9667 0.9398

0.9925 0.9665 0.9665 0.9399

0.9925

–1.93

–2.76

21.42 21.42

100

171.4

Fig. 12. Spatial power map for the 2nd mode of the MOX benchmark. Table 10 Relative assembly-averaged core power distributions and percent relative errors for MOX problem (right-superior quadrant of core in Fig. 11) S16 P1 P3

0.6062 0.6041 0.6065 1.7399 1.7458 1.7380 1.7862 1.7869 1.7885 2.7083 2.7233 2.7055

(0.346) (0.049) (0.339) (0.109) (0.039) (0.129) (0.554) (0.103)

0.5545 0.5523 0.5545 1.1674 1.1654 1.1684 2.2038 2.2133 2.2014 1.7862 1.7869 1.7885

(0.398) (0.000) (0.171) (0.086) (0.431) (0.109) (0.039) (0.129)

0.3847 0.3822 0.3848 0.8476 0.8448 0.8478 1.1674 1.1654 1.1684 1.7399 1.7458 1.7380

0.1760 0.1743 0.1761 0.3847 0.3822 0.3848 0.5545 0.5523 0.5545 0.6062 0.6041 0.6065

(0.650) (0.026) (0.330) (0.023) (0.171) (0.086) (0.339) (0.109)

171.4

(0.966) (0.057)

2.15

(0.650) (0.026)

1.35

(0.398) (0.000)

0.538

100

(0.346) (0.049)

–0.269

–1.08 Table 11 Comparison of the P1 , P3 and S16 eigenvalue keff and assembly core power errors ðmax ; Þ for MOX benchmark problem Method

keff

% Relative error

max

P1 P3 S16

0.9929 0.9925 0.9925

0.040 0.00 Ref.

0.966 0.129 Ref.

(%)

  (%) 0.385 0.068 Ref.

material ‘‘P3” (PuO2 with high enrichment in Pu) in Cavarec et al. (1994). The materials two-group cross-sections considered are shown in Table 8. The results for the four dominant eigenvalues computed with the PL approximations are presented in Table 9. These results have been obtained using one node per assembly (N x ¼ N y ¼ 9), 7 Legendre polynomials and the serendipity approximation. In Table 10 we compare the P 1 and P 3 relative assembly-averaged core power distribution to S16 transport values. Each triplet in the table corresponds to a core assembly in Fig. 11 from rows 2–5 and columns 4–7 (the right-superior 4  4 quadrant of core). The reference result, S16 , was obtained with a 132  132 mesh per assembly. Also the Table shows percent relative differences for each assembly when the power distribution is computed with the P1 and P 3 approximations. We compare the fundamental eigenvalue and the relative assembly core-power distribution to the reference calculation. The comparison is summarized in Table 11. We give the maximum relative error (max ) and the mean relative error ð Þ in the average

–1.88

–2.69

21.42 21.42

100

171.4

Fig. 13. Spatial power map for the 3rd mode of the MOX benchmark.

Table 12 Eigenvalue keff and matrix dimensions of the eigenvalue problem for MOX benchmark Nx  Ny

99

18  18

M

3 4 5 6 7 3 4 5 6

P1

P3

keff

Dim

keff

Dim

0.9893 0.9912 0.9923 0.9927 0.9929 0.9917 0.9926 0.9929 0.9929

972 1620 2430 3402 4536 3888 6480 9720 13608

0.9893 0.9910 0.9919 0.9923 0.9925 0.9915 0.9923 0.9925 0.9925

3888 6480 9720 13608 18144 15552 25920 38880 54432

assembly power of the fundamental mode. Also the Table shows the relative error in keff . In Table 10 the P3 errors are, in general, lower than the errors obtained using the diffusion approximation. The same occurs with the mean relative error and the relative error in keff . The diffusion approximation yields an eigenvalue relative

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

error of 0.04%. The P 3 approximation reduces the error in keff one order of magnitude. In Figs. 12 and 13 we display the P3 spatial core power distribution for the second and third subcritical modes, which are azimuthal degenerate modes. In Table 12, we compare keff calculated for different values of N x , N y and M. Also we list the matrix dimension of the generalized eigenvalue problem given by Eq. (43). Again, from the computational point of view, it is more favourable increasing the polynomial order M than increasing the number of nodes in the discretization of the reactor core. 5. Conclusions A nodal collocation approximation for the multi-dimensional P L equations has been developed dividing the region where the equations have to be solved into rectangular nodes and assuming that the neutronic fields can be expanded in terms of orthonormal polynomials. This method has been implemented into a computer code that solves the Lambda Modes transport problem. This code has been applied to solve four problems: a small slab homogeneous eigenvalue problem, a heterogeneous anisotropic scattering eigenvalue problem, a very heterogeneous isotropic eigenvalue problem of small size and a realistic MOX reactor problem. The obtained results have been compared with the results from the discrete-ordinates code TWODANT (S16 quadrature set). In the first case, a low order PL approximation gives inaccurate results, since the homogeneous region considered has a small dimension, the angular dependence in the solution is large and it is necessary to consider a large order L to obtain results similar to the transport solution. The second problem, with anisotropic scattering, is a variant of the 2D benchmark treated in Hébert (2006), where the zero-flux boundary condition is replaced by a vacuum boundary condition. The dimensions of this benchmark are larger than the other benchmarks. The P 3 approximation shows a good agreement with the reference S16 result. The third problem has been proposed in Brantley and Larsen (2000) and it possesses a strong spatial heterogeneity. It is de-

1829

signed to obtain a problem where the assumptions to be made for the validity of the SP 3 are not fulfilled. It is shown that the P3 and P5 approximations give rather accurate results for this problem. The fourth problem considered is a modification of a benchmark problem corresponding to a MOX reactor. The results show that the P3 approximation is enough to obtain results of the spatial power distribution very close to the S16 solution obtained with TWODANT. This shows the applicability of the multi-dimensional PL equations for the analysis of realistic reactors with mixed-oxide fuels. In summary, we believe that with the availability of fast computers the implementation of multi-dimensional P L equations in the neutronic modules of reactor core simulators becomes feasible and attractive enough to be able to analyse cores with different kind of fuels, and that the use of the spherical harmonics method to solve the neutron transport equation will become more widely used in the nuclear engineering field. Acknowledgements The authors would like to express their gratitude to the reviewers for their suggestions and helpful comments. Partial support to perform this work has been received from the spanish Ministerio de Educación y Ciencia under project ENE2005-09219-C02-00, the Generalitat Valenciana under Project GV06/095, and the Vicerrectorado de Investigación Desarrollo e Innovación de la Universidad Politécnica de Valencia under Project PPI-05-05-5701. Appendix A A.1. P L equations for XY geometry If the problem does not depend on z, the flux U clearly verifies ¼ 0 and UðhÞ ¼ Uðp  hÞ, so only l þ m even moments /m l remain  2 in the series (3) and N even ¼ 12 ðL þ 1Þ . Then the first set of equations (12) yields the expressions

o/ oz

 

   o o o o m2 o o o o m2   C þ2 þ g D  D n D þ D lþ1 lþ1 lþ1 lþ1 2 ox ox oy oy lþ2 ox oy oy ox lþ2     o o o o m o o o o mþ2 n þ C þ2 n Dlþ1 þ Dlþ1 Dlþ1  Dlþ1  C þ2 0 þ2 ox ox oy oy lþ2 ox ox oy oy lþ2   

 o o o o o o o o m2 n Dlþ1 þ Dlþ1 Dlþ1  Dlþ1 gmþ2  C 0þ  2 ox oy oy ox lþ2 ox ox oy oy l

    o o o o o o o o m Dlþ1 þ Dlþ1 Dlþ1 þ Dlþ1 n gm2 þ C 0þ þ 0 ox oy oy ox l ox ox oy oy l 

   o o o o mþ2 o o o o Dlþ1  Dlþ1 nl Dlþ1 þ Dlþ1  C 0þ  gmþ2 þ2 l ox ox oy oy ox oy oy ox 

   o o o o m2 o o o o 0 þ gm2 Dl1  Dl1 nl Dl1 þ Dl1  C 2 ox ox oy oy ox oy oy ox l     o o o o m o o o o mþ2 0 n n  C þ C 0 D þ D D  D l1 l1 l1 l1 0 þ2 ox ox oy oy l ox ox oy oy l

    o o o o o o o o m2 Dl1 þ Dl1 Dl1  Dl1 n gmþ2 þ C 2  2 ox oy oy ox l ox ox oy oy l2

    o o o o o o o o m Dl1 þ Dl1 Dl1 þ Dl1 n gm2  C 2 þ 0 ox oy oy ox l2 ox ox oy oy l2 

   o o o o mþ2 o o o o Dl1  Dl1 nl2  Dl1 þ Dl1 þC 2 gmþ2 þ2 ox ox oy oy ox oy oy ox l2 0 þ ðRt  Rsl Þnm l ¼ mRf dl0 dm0 n0

ðA:1Þ

1830

M. Capilla et al. / Annals of Nuclear Energy 35 (2008) 1820–1830

for l ¼ 0; 2; 4; . . . ; L  1; L odd, l þ m even. The second set of equations (12) is the same as Eq. (A.1) interchanging the moments nlm and glm . The following constants are used:

C þ2 2;lm C þ2 0;lm

1 ½ðl þ 4  mÞðl þ 3  mÞðl þ 2  mÞðl þ 1  mÞ ¼ 2 ð2l þ 3Þ½ð2l þ 1Þð2l þ 5Þ1=2

C þ2 2;l2;mþ2 ;

¼

þ2 C þ2 þ2;lm ¼ C 2;l;m ;

C 2 0;lm

¼

;

¼

Bm l 

C 0 þ2;lm

ðA:2Þ

;

¼

rffiffiffi 21 5 Rt

where Ra ¼ Rt  Rs0 and D ¼



o2 ox2

0

D n22

o2 2 2 g oxoy 2

!

l ¼2 l0 even

l ¼ 1; 3; . . . ; L odd;

0

2M þl0;l0 m0 gm ¼ 0; l0

m0 ¼1

L1 X l X 0

0

2M þlm;l0 m0 gm ¼ 0; l0

m0 ¼1

l ¼2 l0 even m0 6¼m

8 9 > > 0 > > L1 < l = X X m6¼0   0 m0 gl  M lm;l0 0 nl0 þ 2M lm;l0 m0 nl0 ¼0 > > > m0 ¼1 l0 ¼0 > : ; 0 0 l even

m 6¼m

for l ¼ 1; 3; . . . ; L, m ¼ 1; 2; . . . ; l, where

Ml0;l0 0 n0l0 ¼ 0;

L1 X

Mlm;l0 m n0l0 ¼ 0;

gm6l ¼0 

L1 X

ðA:8Þ

M lm;l0 m gm ¼0 l0

l0 ¼0 l0 even

for l ¼ 1; 3; . . . ; L, m ¼ 1; 2; . . . ; l, where 0

0

¼0 þ nm6 l

L1 X

m Mlm;l0 m ¼ Hm l H l0 4p

0

0

l ¼ 1; 3; . . . ; L odd;

l ¼2 l0 even

which are evaluated and give the relations L1 X l X

 ~~ sin hdhY m l Uðr; XÞ ¼ 0;

0

ðA:4Þ

n0l þ

Z p=2 0

0

¼0 þ nm6 l

The condition of zero incoming flux on the left side of Y axis is approximated by the expressions

0

d/

l ¼2 l0 even

 2  oyo 2 .

 ~~ sin hdhY m l Uðr; XÞ ¼ 0;

ðA:6Þ

0

A.2. Marshak’s vacuum boundary conditions

Z p

!

which are evaluated and give the relations

n0l þ

! rffiffiffi 1 5 2 2  2 o2 2 0 0 0 pffiffiffi Dn0  Dn þ Rt n2 þ D n2  2 g ¼ 0; 21Rt 2 7 Rt 3 oxoy 2 3 5 Rt rffiffiffi 3 1 1 2  0 Dn2 þ Rt n22  pffiffiffiffiffiffi D n00 þ D n2 ¼ 0;  7Rt 2 7 3 R t 30R rffiffiffiffiffiffi t 2 rffiffiffi 2 3 2 1 o 2 2 o Dg22  Rt g22  ðA:3Þ n0 þ n0 ¼ 0; 7 Rt 15 Rt oxoy 0 7Rt 3 oxoy 2

d/

m 6¼ m0 :

0

ðA:7Þ

1 mRf n00 ; k

Z p

Z 2p

C 0 2;l;m :

1 1  Dn0 þ Ra n00 þ pffiffiffi Dn02  3Rt 0 3 5Rt

Bm l 

0

m Pm l ðxÞP l0 ðxÞdx;

2 C 2 2;lm ¼ C 2;l;m ;

If, for example, L ¼ 3 then P3 equations for XY geometry are obtained (Fletcher, 1983). Assuming constant cross-sections and isotropic scattering (Rsl ¼ 0 for l > 0), the equations are

¼

1

1  ð1Þmm 1  ð1Þmþm  m  m0 m þ m0

The condition of zero incoming flux on the left side of Z axis is approximated by the expressions

C þ2 0;l2;m ;

1 ½ðl þ 2 þ mÞðl þ 1 þ mÞðl  mÞðl  1  mÞ 2 ð2l þ 3Þð2l þ 1Þ

C 0þ 2;l;m ;

Z 0

1 ½ðl þ 2 þ mÞðl þ 1 þ mÞðl  mÞðl  1  mÞ1=2 C 0 ; 2;lm ¼ 2 ð2l þ 1Þð2l  1Þ ðl þ 2Þðl þ 1Þ þ m2 lðl  1Þ þ m2 C 0þ ; C 0 ; 0;lm ¼ 0;lm ¼ ð2l þ 3Þð2l þ 1Þ ð2l þ 1Þð2l  1Þ C 0þ þ2;lm

m0 Hm l H l0 2

2

1=2

C 0þ 2;lm ¼

¼

1=2

1 ½ðl þ 1  mÞðl þ 1 þ mÞðl þ 2  mÞðl þ 2 þ mÞ1=2 ¼ ; 2 ð2l þ 3Þ½ð2l þ 1Þð2l þ 5Þ1=2

C 2 2;lm

0

Mlm;l0 m0

ðA:5Þ

Z 0

1

m Pm l ðxÞP l0 ðxÞdx:

ðA:9Þ

For the right side of Y and Z axis the same expressions as Eqs. (A.5) and (A.8) result, changing the sign in the sums. References Alcouffe, R.E., Larsen, E.W., Miller Jr., W.F., Wienke, B.R., 1979. Computational efficiency of the numerical methods for the multigroup discrete ordinates neutron transport equations: the slab geometry case. Nucl. Sci. Eng. 71, 111. Alcouffe, R.E., Baker, R.S., Brinkley, F.W., Marr, D.R., O’Dell, R.D., Walters, W.F., 1995. DANTSYS: a diffusion accelerated neutral particle transport code system. Los Alamos National Laboratory, LA-12969-M. Brantley, P.S., Larsen, E.W., 2000. The simplified P 3 approximation. Nucl. Sci. Eng. 134, 1–21. Capilla, M., Talavera, C.F., Ginestar, D., Verdú, G., 2005. A nodal collocation method for the calculation of the lambda modes of the P L equations. Ann. Nucl. Energ. 32, 1825–1853. Cavarec, C., Perron, J.F., Verwaerde, 1994. Benchmark calculations of power distributions within assemblies. Nuclear Energy Agency Committee on Reactor Physics, HT-12/94006 A, NEA/NSC/DOC (94), 28. Clark Jr., M., Hansen, K.F., 1964. Numerical Methods of Reactor Analysis. New York Academic, New York. Davison, B., 1957. Neutron Transport Theory. Oxford University Press, London. Fletcher, J.K., 1983. A solution of the neutron transport equation using spherical harmonics. J. Phys. A: Math. Gen. 16, 2827–2835. Gelbard, E.M., 1968. Spherical harmonics methods: P L and double P L approximations. In: Kelberg, D., Okrent, C.N. (Eds.), Computing Methods in Reactor Physics. Gordon and Breach, New York. Gradshteyn, I.S., Ryzhik, I.M., 1980. Table of Integrals, Series and Products. In: Jeffrey, A. (Ed.). Academic Press, Orlando, Florida. Hébert, A., 1987. Development of the nodal collocation method for solving the neutron diffusion equation. Ann. Nucl. Energ. 14 (10), 527–541. Hébert, A., 2006. The search for superconvergence in spherical harmonics approximations. Nucl. Sci. Eng. 154, 134–173. Larsen, A.W., Morel, J.E., McGhee, J.M., 1996. Asymptotic derivation of the multigroup P 1 and simplified P N equations with anisotropic scattering. Nucl. Sci. Eng. 123, 328–342. Stacey, W.M., 2001. Nuclear Reactor Physics. Wiley, New York. Tomasevic´, D.I., Larsen, E.W., 1996. The simplified P 2 approximation. Nucl. Sci. Eng. 122, 309–325. Verdú, G., Ginestar, D., Vidal, V., Muñoz-Cobo, J.L., 1994. k modes of the neutron diffusion equation. Ann. Nucl. Energ. 21 (7), 405–421. Weinberg, A.M., Wigner, E.P., 1958. The Physical Theory of Neutron Chain Reactors. Chicago University Press, Chicago.