Moment-based probability tables for angular anisotropic scattering

Annals of Nuclear Energy 38 (2011) 1125–1132

Contents lists available at ScienceDirect

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

Moment-based probability tables for angular anisotropic scattering Nicolas Martin a, Joachim Miss b, Alain Hébert a,⇑ a b

Institut de Génie Nucléaire, École Polytechnique de Montréal, P.O. Box 6079, Station ‘‘Centre-Ville’’ Montréal, Québec, Canada H3C 3A7 IRSN/DSU/SEC/LERD, 31, Avenue de la Division Leclerc, BP 17, 92262 Fontenay-aux-Roses, France

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 30 September 2010 Received in revised form 8 December 2010 Accepted 10 December 2010 Available online 22 January 2011

This paper presents a strategy for taking into account anisotropy scattering into a Monte Carlo algorithm relying on the subgroup method and developed in the DRAGON lattice code. For the sake of consistency, we limited our Monte Carlo code to the same cross-section libraries available for deterministic methods. However Legendre moments for the transfer cross-sections cannot be directly used during the Monte Carlo random walk, due to the presence of non-positive parts into the distributions. The discrete angle method is proposed to deal with this limitation, following an approach initially introduced in the MORET multigroup Monte Carlo code. We selected a moment approach, originally employed to compute probability tables for resonant cross-sections, to derive consistent sums of Dirac distributions conserving Legendre moments of the angular distributions. A detailed analysis of the applicability of the moment approach is here mandatory. When the moment technique fails due to incoherent Legendre moments, the discrete angle technique is substituted by legacy semi-analytical methods. We illustrate the proposed method using critical benchmarks coming from the ICSBEP handbook by comparison toward SN and other Monte Carlo results. The impact of the anisotropy scattering is also discussed on a PWR MOX assembly case. Ó 2010 Elsevier Ltd. All rights reserved.

Keywords: Anisotropy scattering Probability tables Subgroup method Monte Carlo

1. Introduction Taking into account scattering anisotropy into multigroup Monte Carlo codes is a difficult topic, due to the fact that scattering cross-sections are given in terms of Legendre moments:

rs ðE0

L X 2‘ þ 1 rs;‘ ðE0 2 ‘¼0

E; lÞ ¼

EÞP‘ ðlÞ;

ð1Þ

with

rs;‘ ðE0

EÞ ¼

Z

1

dlrs ðE0

E; lÞP ‘ ðlÞ:

ð2Þ

1

We denote here L the Legendre expansion order of the scattering cross-sections, and P‘(l), ‘ = 0, . . . , L the corresponding Legendre polynomials. These moments are perfectly suitable for deterministic methods such as the SN, PN techniques and the method of characteristics, or any transport method where the source term can be written as an angular expansion in real spherical harmonics. However for Monte Carlo methods, this formalism suffers from some mathematical inconsistency due to possible negative values of the scattering cross-sections appearing in

⇑ Corresponding author. Tel.: +1 514 340 4711x4519. E-mail address: [email protected] (A. Hébert). 0306-4549/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2010.12.013

limited Legendre expansions. The resulting angular distribution can be negative in some part of the range [1, 1]. Consequently, directly sampling into such probability density functions during a Monte Carlo random walk is impractical. Different techniques have so far been proposed for taking into account anisotropy scattering into Monte Carlo codes that use homogeneous, multigroup cross-sections. Most of them are mentioned in the reference paper of Brockmann (1981), notably the discrete angle method where angular scattering probability densities are replaced by Dirac distributions. This method was originally derived for the MORSE Monte Carlo code (Straker et al., 1970), and subsequently implemented in the multigroup version of the MORET code (Miss et al., 2007). This approach is routinely used in the APOLLO2-MORET way of the CRISTAL system (Gomit et al., 2003), the French criticality–safety package. Another class of sampling technique relies on equally or non-equally step function representation of the scattering law. Non-equally step function representation has been fully studied in Mao et al. (1998). This study is intended in the framework of a Monte Carlo code relying on the subgroup (or multiband) method for the whole energy range developed into the DRAGON lattice code (Marleau et al., 1992). This method has been first mentioned by Cullen (1974) in the USA, and Nikolaev (1976) in Russia. In this formalism, crosssections in each resonant group are divided into several bands, or subgroups, each of them associated with a weight. In our case, the moment approach proposed by P. Ribon (Ribon and Maillard,

1126

N. Martin et al. / Annals of Nuclear Energy 38 (2011) 1125–1132

1986) is selected for computing cross-section probability tables in resonant groups. The group structure used in this study is detailed in Hébert (2009) and comprises 295 groups. For the time being, only isotropic or transport corrected cross-sections were used for PWR or CANDU-6 lattice calculations (Martin and Hébert, 2011), hence mitigating the range of possible applications. The development of a consistent angular anisotropy model was therefore essential. In the case of a Monte Carlo code based on isotopic (both multigroup or subgroup) cross-sections, it is however feasible to use exactly the same data as in point-wise Monte Carlo codes. This approach is for instance applied in the Monte Carlo code TART (Cullen, 2005) and mentioned in Diop et al. (2010) for the TRIPOLI4 code. In both cases, cross-section probability tables are used for sampling the neutron free path, the collided isotope and the type of reaction, while the point-wise anisotropic laws are used in case of a scattering event, leading in practice to a continuous-energy neutron random walk. Another interesting option described in Ribon and Coste-Delclaux (2008) consists in computing the multigroup scattering laws by a preliminary continuous-energy Monte Carlo calculation, where the sampled cosine is stored into bins for each integral group-to-group transfer. Yet in our study, we have only considered data available through the DRAGON input libraries formatted by the GROUPR module of NJOY code (Macfarlane and Muir, 2000) and an additional in-house Dragr module (Hébert and Saygin, 1992), where angular probability density functions in the LAB frame are given in terms of Legendre series expansions. The objective here is to propose into the DRAGON lattice code a Monte Carlo algorithm relying on the same input data as in existing deterministic flux solvers based on the method of characteristics or the discrete ordinate method. Due to our experience in computing probability tables for cross-sections, we selected the discrete angle method to model angular anisotropy scattering in our Monte Carlo code. Furthermore, one of the main objective of the proposed subgroup Monte Carlo method was the reduction of the overall computational cost of the simulation, compared to a point-wise representation. The discrete angular technique thus introduces a satisfactory compromise between accuracy and computational resource expenses. 2. Angular anisotropy scattering in multigroup or subgroup approximations

The probability density function (pdf) for the scattering cosine can be easily derived from the Legendre expansions. In the case of multigroup cross-sections, we can write 0

g

ðlÞ ¼ rgs;0 pg0

g fg 0

gð

lÞ;

ð3Þ

where  pg0  fg0

rg

0

g

¼ rs;0g , probability of a transfer g0 g. s;0 0 PL 2‘þ1 rgs;‘ g g ð lÞ ¼ ‘¼0 2 f‘;g 0 g P ‘ ðlÞ, with f‘;g 0 g ¼ g 0 g , pdf for the g

rs;0

scattering cosine l corresponding to the transfer g0

In the case of the subgroup method, cross-sections for resonant isotopes are on the form of a probability table: {xk, rq,k}k2[1,K], with q = total, absorption, fission, scattering, etc. In our approach, the scattering matrix coefficients are available through a set of L Legendre moments for each transfer from a subgroup k of a group g toward a group g0 , denoted in the remainder of this paper g0 {g, k}:

rgs

0

fg;kg

ðlÞ ¼

L X 2‘ þ 1 g0 rs;‘ 2 ‘¼0

fg;kg

P‘ ðlÞ;

ð4Þ

Similarly to Eq. (3), we can write

rgs

0

fg;kg

¼ rg;k s;0 pg 0

fg;kg fg 0

fg;kg ð

lÞ;

ð5Þ

where 0

 pg0

fg;kg

¼

rgs;0

lÞ ¼

 fg0

fg;kg ð

fg;kg

, probability of a transfer g0

rg;k s;0

PL

2‘þ1 ‘¼0 2 f‘;g 0

lÞ, with f‘;g0

fg;kg P ‘ ð

{g,k}. fg;kg

¼

g rs;‘

0

fg;kg

0

fg;kg

g rs;0

, pdf

for the scattering cosine l corresponding to the transfer g0 {g,k}.

With the above relations, the explicit subgroup-to-subgroup transfer coefficients are unknowns, but can be approximated by (Cullen, 1986): 0

0

rsfg ;k g

fg;kg

¼ x0k rsg

0

fg;kg

:

ð6Þ

Note also that the multigroup scattering law can be deduced from the probability table coefficients, as

rgs

0

g

¼

K X

xk rgs

0

fg;kg

:

ð7Þ

k¼1

In particular, there is no transfer of a subgroup k of a group g toward a group k0 of a group g0 if there is no transfer from g to g0 . The scattering event can therefore be summarized as follow in the probability table case: (1) Selection of the secondary group g0 using pg0 fg;kg . (2) Sampling the l cosine with fg0 fg;kg ðlÞ. (3) Sampling of the new subgroup k0 in group g0 .

2.1. Multigroup cross-sections

rsg

2.2. Cross-section probability tables

g.

During the multigroup Monte Carlo random walk, any scattering event in group g is processed in two steps: (1) First, the secondary energy group is sampled using pg0 g . (2) fg0 g ðlÞ is then used to determine the scattering angle l in the LAB frame.

However the direct sampling in the truncated Legendre series fg0 g ðlÞ or fg0 fg;kg ðlÞ is inapplicable, due to the fact that the angular distribution is often negative in the interval [1, 1]. It is illustrated in Fig. 1 for the scattering law from group 6 to 6 in the SHEM-295 mesh f6<6(l) (self-scattering occurring in the energy range [9.04 MeV, 10 MeV]), for 1H in H2O as generated by NJOY99. The modification of the density function by taking its absolute value yields to negative weights, and is therefore of poor interest (Brockmann, 1981). It is instead more convenient to rewrite the angular distribution into an equivalent one that conserves its associated moments. It is at the basis of both the equally/non-equally probable step function representation and the discrete angle method. 2.3. The discrete angle method In the discrete angle method, Legendre moments are converted into discrete angles, each of them associated with a weight. This formalism is equivalent to the one proposed for cross-section

1127

N. Martin et al. / Annals of Nuclear Energy 38 (2011) 1125–1132

0.8 P1 approximation P3 approximation

0.7

P5 approximation

Probability density

0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Scattering cosine (LAB) Fig. 1. lLAB pdf for 6

probability tables, and is based on the replacement of the angular probability density by a set of Dirac distributions respecting a generalized Gaussian quadrature. The transfer index is omitted for clarity in the remainder of this section. We write f(l) the angular pdf in such a way that the L + 1 Legendre moments lead into K ¼ Lþ1 discrete angles {xk, lk}k2[1,K] conserving the L + 1 moments, 2 which imply two sets of equations:

f ðlÞ ’

K X

xk dðl  lk Þ;

ð8Þ

k¼1

and

Mn ¼

K X

xk lnk with 0 6 n 6 L;

Mn ¼

‘¼n;n2;...

dlf ðlÞl

n

with 0 6 n 6 L:

ð10Þ

1

Similarly to the cross-section case, probability tables can be computed by several methods. In the original approach programmed in the MORSE code, a recursive technique is used to compute a set of orthogonal polynomials to the angular pdf that conserve the Legendre moments. However such method may be numerically instable, due to the use of an iterative procedure. Initially derived for computing probability tables for cross-sections, the moment approach of P. Ribon, often referred to as the CALENDF method, can likely be applied to the computation of probability tables for the scattering cosine distribution. It has first been suggested and tested at IRSN by Le Cocq (1998) and is now daily used in the multigroup Monte Carlo code MORET. The mathematical procedure employed for computing angular probability tables in MORET is identical to the cross-section case and thoroughly described in Ribon and Maillard (1986). The CALENDF method is a straightforward procedure that lead to probability tables where the weights xk are positives and base points lk are in the physical support [1, 1], if and only if moments Mn are consistent (Ribon and Coste-Delclaux, 2008). The issue here is that the quality of the (simple) moments are directly linked to those of the Legendre moments computed from NJOY. We have indeed (Weisstein, 2010):

ð2‘ þ 1Þn! P ‘ ðlÞ;  ðn  ‘Þ !ð‘ þ n þ 1Þ!! 2

ð11Þ

ð2‘ þ 1Þn! f‘ :  ðn  ‘Þ !ð‘ þ n þ 1Þ!!

ð12Þ

1

X

Mn ¼

‘¼n;n2;...

2ðn‘Þ=2

1 2

The maximum Legendre order allowed in our study is L = 7, leading to the computation of the following moments for each transfer:

M0 ¼ 1;

ð13Þ

M1 ¼ f1 ;

ð14Þ

M3

1

2ðn‘Þ=2

which implies that

M2

We recall the definition of Mn , nth order moment of f(l):

X

ln ¼

ð9Þ

k¼1

Z

6 transfer, SHEM-295 groups.

M4 M5 M6 M7

1 ¼ ½f0 þ 2f 2 ; 3 1 ¼ ½3f 1 þ 2f 3 ; 5 1 ¼ ½7f þ 20f 2 þ 8f 4 ; 35 0 1 ½27f 1 þ 28f 3 þ 8f 5 ; ¼ 63 1 ½33f 0 þ 110f 2 þ 72f 4 þ 16f 6 ; ¼ 231 1 ½143f 1 þ 182f 3 þ 88f 5 þ 16f 7 : ¼ 429

ð15Þ ð16Þ ð17Þ ð18Þ ð19Þ ð20Þ

However Legendre moments can be badly computed by NJOY for some transfers, resulting to non-physical probability tables as it is illustrated in Ribon and Coste-Delclaux (2008). Consequently, special care should be taken when computing probability tables with Legendre moments. In order to avoid incoherent probability tables, we preliminarily test the applicability of the moment method. The calculation of angular probability tables with a moment approach is connected to the Hausdroff moment problem, fully described in Akhiezer (1965). Moments shall therefore fulfill strict mathematical conditions, leading in practice to a set of preliminary checks exposed below. We denote moments with the help of the expected value E:

Mn ¼

Z

1 1

  dlln f ðlÞ ¼ E ðf ðlÞÞn ;

ð21Þ

1128

N. Martin et al. / Annals of Nuclear Energy 38 (2011) 1125–1132

Central moments must be checked. They are defined as n

Un ¼ E½ðf ðlÞ  E½f ðlÞ 

ð22Þ

Simple relations exist between simple and central moments, given below for the first orders:

U2 ¼ M2  M21 ;

ð23Þ

U3 ¼ M3  3M1 M2 þ 2M31 ;

ð24Þ

U4 ¼ M4  4M1 M3 þ 6M21 M2  3M41 :

ð25Þ

We define also moments translated to the +1 and 1 abscissas by

Vn ¼

Z

1

dlln ðf ðlÞ  1Þ;

ð26Þ

1

Wn ¼

dlln ðf ðlÞ þ 1Þ:

ð27Þ

Moments shall meet the following criteria: (1) Moments of even order 2n must be strictly positives M2n > 0. (2) Moments of order 2n must be lower than moments of order 2n2: M2n < M2n2 . (3) Central moments of even order 2n must be strictly positives: U2n > 0. (4) Central moments of even order 2n must be lower than central moments of order 2n2: U2n < U2n2 . (5) Moments translated to the +1 abscissa of even order 2n must be strictly positives: V2n > 0. (6) Moments of order 2n + 1 odd must be lower in absolute value than moments of order 2n even: jM2nþ1 j < M2n . (7) Moments translated to the +1 abscissa of odd order must be strictly negatives: V2nþ1 < 0. (8) Moments translated to the 1 abscissa must be strictly positives: Wn > 0. (9) The following determinants must be positives:

  Mn1  Mn  M M n1

n2

 Mn1  Mn2  > 0; Mn3  Mn2 

  Mn1 þ Mn  M þM n1

ð28Þ

n2

 Mn1 þ Mn2  > 0: M þM  n3

ð29Þ

n2

In the case of one of these condition is not fulfilled, the order n is diminished by 2 and the process is repeated until tests are successful. The CALENDF technique is then applied with the L + 1 = 2K consistent moments. We first consider the Stieltjes generating function of the moments: 1

dl

1

2K1 X f ðlÞ ¼ z‘ M‘ þ Oðz2K Þ; 1  zl ‘¼0

ð30Þ

A Padé approximant of order [K1, K] can then be introduced to repP ‘ resent the term 2K1 ‘¼0 z M‘ , leading to

FðzÞ ¼

2K1 X

z‘ M‘ ¼

‘¼0

a0 þ a1 z þ    þ aK1 zK1 ; b0 þ b1 z þ    þ bK1 zK1 þ zK

ð31Þ

so that

zK þ

K1 X k¼0

BM B K1 B . B . @ . M1

MK1 MK2 .. . M0

0 1 0 .C B b C B B . . . M0 C CB 1 C B .. C C C C B : .. .. CB .. C ¼ B .C C . . A@ . A B @ .. A . . . M2K bK1 0 ...

M1

10

b0

1

ð33Þ

A LDLT factorization of the coefficient matrix is then performed, leading to the bk coefficients. Padé approximants of a Stieltjes function have simple real poles, consequently base points lk are the roots of the following polynomial:

ð34Þ

! b k zk

2K1 X ‘¼0

! z‘ M‘

¼

2K1 X k¼0

ak zk :

In the original method of P. Ribon, the weights are calculated using a root mean square procedure over all the L + 1 moments. We chose here to use the approach proposed in Hébert and Coste (2002), where only K moments in the interval 0 6 ‘ 6 L1 are selected. We 2 obtain the following system

0 B B B B @

1

1

...

l1

l2

.. .

.. .

... .. .

1

10

1

0

M0

1

.. .

B B CB . C ¼ B CB . C B A@ . A @

M1 ...

C C C: C A

x1 C B lK C CB x2 C

lK1 lK1 . . . lK1 1 2 K

xK

ð32Þ

ð35Þ

MðL1Þ=2

In this case, a Vandermonde matrix appears, for which an analytical solution exists for xk, 1 6 k 6 K:

xk ¼

ðL1Þ=2 X ð1ÞK1 ck;‘ M‘ ; K Q ðl‘  lk Þ ‘¼0

ð36Þ

‘¼1 ‘–k

where ck,‘ are the coefficients of the following polynomial:

ck;0 þ ck;1 z þ . . . þ ck;ðL1Þ=2 zðL1Þ=2 ¼

and

Z

MK

k¼1

1

1

FðzÞ ¼

0

K Y ðz  lk Þ ¼ b0 þ b1 z þ . . . þ bK1 zK1 þ zK :

and

Z

By identification of the coefficients of the power of z between K and 2K1, we obtain the following system where the bk are the unknowns:

K Y

ðz  l‘ Þ:

ð37Þ

‘¼1 ‘–k

The probability table obtained for the transfer 1 1 for 1H in H2O isotope is displayed in Fig. 2 when a P5 approximation is used. We obtain three discrete angles representing the forward peaked scattering law. In practice, the discrete angle method is used only when the available number of consistent moments is greater than or equal to 6, i.e., at least P5 scattering cross-section expansions. Using only one or two angles (P1 or P3 expansions) can indeed lead to ray effects in the neutron random walks. We thus kept the anisotropy scattering strategy implemented in the MORET code, where semi-analytical techniques are used in those cases. The basic idea of such techniques is to add an analytic part to the Dirac distributions in order to mitigate the ray effects. The semi-linear Coveyou technique detailed in Lux and Koblinger (1991) can preserve the P1 Legendre expansion, while the semi-continuous method proposed in Lux (1982) ensure the conservation of the P3 Legendre moments. All the above relations have been implemented in the Monte Carlo module MC: of the DRAGON framework. To summarize, taking into account anisotropy scattering depends on the consistency of the Legendre moments and can be summarized as follows:  Consistent Pn scattering cross-sections with n greater or equal to 5: the discrete angle method is used with nþ1 angles. 2  P3 scattering cross-sections: the Lux method is applied, leading to two discrete angles and a continuous part.

1129

Probability density

N. Martin et al. / Annals of Nuclear Energy 38 (2011) 1125–1132

0.6

P approximation 5

0.5 0.4 0.3 0.2 0.1 0 −0.1 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

0.4

0.6

0.8

1

Scattering cosine (LAB) 0.7

Probability

0.6 0.5 0.4 0.3 0.2 0.1 0 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

Scattering cosine base point Fig. 2. Probability table for f1

k (LAB)

l), 1H in H2O, SHEM-295 groups.

1(

 P1 scattering cross-sections: the Coveyou method is used, leading to one discrete angle and a continuous part.

 HEU-MET-FAST-004: a water-reflected, highly enriched uranium sphere.

This approach has been thoroughly validated in the case of multigroup, homogeneous cross-sections. It is routinely used in the MORET code for criticality–safety studies. We have keeped the same strategy in the case of isotopic cross-section probability tables.

All those benchmarks have high neutron leakages and a predominant fast flux spectrum for which anisotropy scattering is important, making them valuable tests for the proposed angular anisotropy scattering. The impact of the angular anisotropy treatment on the keff is discussed, by varying the order of the Legendre expansion n. We compare the results produced by the Monte Carlo code, denoted by MC:, toward those produced by the discrete ordinate solver SNT: of the DRAGON code and to the multigroup Monte Carlo MORET code.

3. Numerical results In this section, we discuss some results produced by the proposed Monte Carlo method on different benchmarks, selected in order to test the scattering anisotropy strategy. The first test cases are configurations coming from criticality–safety studies for which angular anisotropy scattering cannot be neglected. The second case is a PWR assembly with three different types of MOX fuels. Typically, angular anisotropy influences neutronics parameters such as the migration area and the albedo:  The modification of the migration area impacts the neutron leakage rate from the fissile zone. This effect is strengthened when the geometry itself tends to increase the leakage rate.  The existence of privileged scattering directions can modify the number of neutrons going to the reflector, and consequently, the number of neutrons coming back to the fissile zone. This effect is increased when the dimension of the fissile zone is small compared to the reflector’s size. 3.1. ICSBEP benchmarks The selected benchmarks are well-known critical experiments coming from the ICSBEP handbook (ICSBEP, 2002) and model spheres of fissile metal with or without reflector, namely:  HEU-MET-FAST-001: the Govida critical assembly, an unreflected critical uranium sphere.  PU-MET-FAST-001: the Jezebel critical assembly, an unreflected plutonium sphere.  PU-MET-FAST-005: a tungsten-reflected, plutonium sphere.

3.1.1. SNT: options Due to the spherical harmonics representation of the source term in the SNT: module, it is possible to take into account anisotropy scattering up to any arbitrary order. As a consequence, the 1D spherical SN results are computed using exactly the same scattering Legendre expansions than in the Monte Carlo case. The following options have been chosen to generate the SN solutions:  A spatial discretization with a mesh spacing of 0.1 cm.  An angular Gauss-Legendre quadrature of order 64. A lattice calculation is performed for each case with a two-step scheme. The first step is a preliminary self-shielding correction of the resonant cross-sections using a SN subgroup method, with the same parameters as for the main flux calculation. The subgroup method is detailed in Hébert (2009) and uses exactly the same cross-section probability tables as for the subgroup Monte Carlo code MC:. The impact of the self-shielding correction on the cross-sections is not very important in these cases, but has been kept for consistency with the subgroup Monte Carlo results. For both the subgroup Monte Carlo algorithm and the discrete ordinate method, the same input library is used. The chosen library is processed by the NJOY release 99.252 + upnea with the additional DRAGR module and based on the JEFF3.1 evaluation. The SHEM295 group structure is used with all scattering cross-sections available up to the P7 order.

1130

N. Martin et al. / Annals of Nuclear Energy 38 (2011) 1125–1132

Table 1 HEU-MET-FAST-001 and PU-MET-FAST-001 benchmarks. Legendre order

keff SN

jDkeffj (pcm)

keff MC:

jDkeffj (pcm)

keff MORET

jDkeffj (pcm)

HEU-MET-FAST-001: keff = 1.000 ± 100 pcm (read as 10 ) P0 1.1024 10240 P1 0.9904 960 P3 0.9952 480 0.9952 480 P5 P7 0.9952 480

1.1065 ± 66 pcm 0.9971 ± 67 pcm 0.9974 ± 63 pcm 0.9964 ± 69 pcm 0.9962 ± 69 pcm

10650 290 260 360 380

1.1023 ± 57 pcm 0.9953 ± 51 pcm 0.9957 ± 52 pcm 0.9964 ± 52 pcm 0.9958 ± 50 pcm

10230 470 430 360 420

PU-MET-FAST-001: keff = 1.000 ± 200 pcm P0 1.0989 P1 0.9911 P3 0.9977 P5 0.9977 0.9977 P7

1.0998 ± 74 pcm 0.9985 ± 72 pcm 0.9967 ± 79 pcm 0.9972 ± 74 pcm 0.9984 ± 72 pcm

9980 150 330 280 160

1.1001 ± 52 pcm 0.9986 ± 48 pcm 0.9978 ± 51 pcm 0.9976 ± 47 pcm 0.9968 ± 51 pcm

10010 140 220 240 320

5

9890 890 230 230 230

Table 2 PU-MET-FAST-005 and HEU-MET-FAST-004 benchmarks. Legendre order

jDkeffj (pcm)

keff MC:

jDkeffj (pcm)

keff MORET

jDkeffj (pcm)

PU-MET-FAST-005: keff = 1.000 ± 130 pcm P0 1.1529 P1 0.9761 P3 1.0022 P5 0.0031 P7 1.0031

keff SN

15,290 2390 220 310 310

1.1545 ± 74 pcm 0.9962 ± 73 pcm 1.0036 ± 73 pcm 1.0038 ± 69 pcm 0.9988 ± 68 pcm

15450 380 360 380 120

1.1571 ± 73 pcm 0.9976 ± 73 pcm 1.0038 ± 73 pcm 1.0026 ± 73 pcm 1.0028 ± 73 pcm

15,710 240 380 260 280

HEU-MET-FAST-004: keff = 0.9985 ± 224 pcm P0 1.2030 P1 0.9826 P3 1.0010 P5 1.0012 P7 1.0013

20,450 1590 250 270 280

1.2040 ± 59 pcm 1.0102 ± 58 pcm 0.9957 ± 58 pcm 1.0028 ± 58 pcm 1.0023 ± 58 pcm

20550 1170 280 430 380

1.2035 ± 62 pcm 0.9992 ± 60 pcm 0.9964 ± 54 pcm 0.9959 ± 61 pcm 0.9955 ± 57 pcm

20,500 70 210 260 300

3.1.2. MORET options Additionally, the multigroup Monte Carlo MORET code have been used to produce keff for each case. Here, the 295-group macroscopic cross-sections are produced by the DRAGON code, using the discrete ordinate lattice calculation with the same input JEFF3.1 library. The MORET code processes angular anisotropy in the same way than the MC: code, except that Legendre transfer cross-sections, and consequently angular scattering pdfs, are for homogenized mixtures and not for isotopes.

3.1.3. Interpretations Results are displayed in Table 1 for the unreflected experiments HEU-MET-FAST-001 and PU-MET-FAST-001, and in Table 2 for the reflected cases PU-MET-FAST-005 and HEU-MET-FAST-004. As expected, the relative error diminishes when the Legendre order is increased. In every case, isotropic (P0) scattering conducts to large errors for the computed keff compared to the experimental value. The Coveyou model used in both MC: and MORET codes show generally a very good agreement with the experimental keff, which is not the case for the SN result with P1 scattering model. It was already stated in Le Cocq (1998) and reveals the effectiveness of this semi-analytical method for taking into account scattering anisotropy in most materials. The only case where P3 approximation is mandatory for the MC: module is the HEU-MET-FAST-007, where the Hydrogen anisotropy scattering is badly represented with the Coveyou technique. We also denote the accuracy of the discrete angle method in every case for P5 and P7 approximations (respectively three and four discrete angles) compared to the SN cases and the experimental values. In summary, the difference between the experimental keff and the multigroup/subgroup with discrete angle scattering Monte Carlo keff is similar to those reported in the ICSBEP handbook for continuous-energy Monte Carlo codes.

3.2. PWR MOX assembly case In order to test the discrete angle method for lattice calculation, we selected a PWR assembly of 17  17 pin cells with three different MOX fuels, enriched to 2.9%, 4.4% and 5.6% in fissile Pu, respectively. The geometry of this assembly is illustrated in Fig. 3. Reflective boundary conditions are applied. In this benchmark, we compared the flux homogenized for each pin cell and condensed to two groups, one fast and one thermal with a limit set at 0.625 eV, between the proposed subgroup Monte Carlo method

Fig. 3. MOX assembly geometry.

1131

N. Martin et al. / Annals of Nuclear Energy 38 (2011) 1125–1132

Percent difference for the fast flux

Percent difference for the fast flux 16

16 8

14 12

6

12

10

4

10

8

2

6

−2

2

−1

8 −2

0

5

10

4

0

15

16

−4 0

5

10

15

Percent difference for the thermal flux 1.5

14 1

12

−3

2

Percent difference for the thermal flux 16

1.5

14

1

12 0.5

10 8

0

6 −0.5

4

0.5

10 8

0

6

−0.5

4

2 0

0

6 0

4

0

1

14

−1 0

5

10

15

−1

2 0

0

Fig. 4. Percent difference, isotropic scattering.

and the point-wise Monte Carlo code SERPENT (Leppänen, 2009). The computed keff are very similar, with typically keff = 1.16477 ± 0.00034 for the continuous-energy Monte Carlo simulation. Two distinct calculations have been performed using the MC: code:  P0 isotropic scattering, leading to a keff = 1.16557 ± 0.00033 and a absolute difference of 80 pcm toward the point-wise result.  P5 (discrete angle technique with three angles) scattering, leading to a keff = 1.16534 ± 0.00033 and a absolute difference of 57 pcm inferior to 1r. We then compare the homogenized flux per pin cell obtained with isotropic scattering with the one produced by the continuous-energy simulation. The percent difference between both codes is plotted in Figs. 4 and 5 for the P0 and P5 case, respectively. In this case, some large discrepancies appear for the fast flux in water holes (up to 8%). However the application of the discrete angle method with P5 scattering cross-sections eliminate this effect, as illustrated at the upper part of Fig. 5. 4. Conclusion A description of the mathematical formalism invoked for computing consistent probability tables for the angular scattering distribution has been given in the context of a subgroup Monte Carlo method. In case of incoherent Legendre moments, the discrete angle technique can be legitimately substituted by rigorous semianalytic methods. Pn scattering cross-sections can be therefore

5

10

15

Fig. 5. Percent difference, P5 approximation.

employed to take into account anisotropy scattering in both stochastic or deterministic transport solvers of the DRAGON lattice code. More importantly, the subgoup Monte Carlo method employed with the discrete angle method leads to results with the same level of accuracy as continuous-energy Monte Carlo codes. Acknowledgments This work was supported by grants from the Natural Science and Engineering Research Council of Canada and from the Fonds Québécois de la Recherche sur la Nature et les Technologies. References Brockmann, H., 1981. Treatment of anisotropic scattering in numerical neutron transport theory. Nucl. Sci. Eng. 77, 377–414. Straker, E.A., Stevens, P.N., Irving, D.C., Cain, V.R., 1970. The MORSE Code – A Multigroup Neutron and Gamma-Ray Monte Carlo Transport Code, Report ORNL-4585. Miss, J., Jacquet, O., Bernard, F., Forestier, B., 2007. Massive evolution in the MORET Monte Carlo code. A new multigroup-continuous energy version. In: ICNC 2007, Saint Petersburg, Russia, May 28–June 01. Gomit, J.M., Cousinou, P., Diop, C., De Grado, F., Gantenbein, F., Grouiller, J.P., Marc, A., Mijuin, D., Toubon, H., 2003. CRISTAL V1: criticality package for burn up credit. In: ICNC 2003 Conference, Tokaimura, Japan. Mao, L., Both, J.P., Nimal, J.C., 1998. Transfer matrix treatments in TRIMARAN-II, nonequally probable step function representation in multigroup Monte Carlo. Nucl. Sci. Eng. 130, 226. Marleau, G., Hébert, A., Roy, R., 1992. New computational methods used in the lattice code dragon. In: Int. Top. Mtg. on Advances in Reactor Physics, Charleston, USA, March 8–11.

1132

N. Martin et al. / Annals of Nuclear Energy 38 (2011) 1125–1132

Cullen, D.E., 1974. Application of the probability table method to multigroup calculations of neutron transport. Nucl. Sci. Eng. 55, 387–400. Nikolaev, M.N., 1976. Comments on the probability table method. Nucl. Sci. Eng. 61, 286. Ribon, P., Maillard, J.M., 1986. Les tables de probabilité. Application au traitement des sections efficaces pour la neutronique. Commissariat à l’Energie Atomique, CEA-N-2485. Hébert, A., 2009. Development of the subgroup projection method for resonance self-shielding calculations. Nucl. Sci. Eng. 162, 56–75. Martin, N., Hébert, A., 2011. A Monte Carlo lattice code with probability tables and optimized energy meshes. Nucl. Sci. Eng. 167, 1–19. Cullen, D.E., 2005. TART2005: A Coupled Neutron–photon 3-D, Combinatorial Geometry, Time Dependent Monte Carlo Transport Code, Lawrence Livermore National Laboratory, UCRL-SM-218009. Diop, C.M., Petit, O., Jouanne, C., Coste-Delclaux, M., 2010. Adjoint Monte Carlo neutron transport using cross-section probability table representation. Ann. Nucl. Energy 37, 1186–1196. Ribon, P., Coste-Delclaux, M., 2008. Angular Anisotropy Representation by Probability Tables. In: International Conference on the Physics of Reactors, September 14–19, Interlaken, Switzerland. Macfarlane, R.E., Muir, D.W., 2000. NJOY99.0 Code System for Producing Pointwise and Multigroup Neutron and Photon Cross Sections from ENDF/B Data, Los Alamos National Laboratory, PSR-480/NJOY99.0.

Hébert, A., Saygin, H., 1992. Development of DRAGR for the Formatting of DRAGON Cross-section Libraries, paper presented at the Seminar on NJOY-91 and THEMIS for the Processing of Evaluated Nuclear Data Files, NEA Data Bank, Saclay, France, April 7–8. Cullen, D.E., 1986. CRC Handbook of Nuclear Reactors Calculations, Vol. 1. Yigal Ronen CRC Press. Le Cocq, A., 1998. Contributions to the Development of Monte Carlo Techniques for Criticality Studies. Ph.D. Thesis, Paris-sud University (in French). Weisstein, E.W., 2010, Legendre Polynomial, From MathWorld – A Wolfram Web Resource. . Akhiezer, N.I., 1965. The Classical Moment Problem and some related questions in analysis. Oliver & Boyd Ltd., Edinburgh and London. Lux, I., Koblinger, L., 1991. Monte-Carlo Particle Transport Methods: Neutron and Photon Calculations. CRC Press. Lux, I., 1982. Semicontinuous selection of scattering angles from low-order pn scattering densities. Nucl. Sci. Eng. 82, 332–337. Hébert, A., Coste, M., 2002. Computing moment-based probability tables for selfshielding calculations in lattice codes. Nucl. Sci. Eng. 142, 245–257. ICSBEP, 2002. International Handbook of Evaluated Criticality Safety Benchmark Experiments, NEA/NSC/DOC(95)03/0I-VII, Organisation for Economic Cooperation and Development, Nuclear Energy Agency. Leppänen, J., 2009. PSG2/Serpent: a Continuous-energy Monte Carlo Reactor Physics Burnup Calculation Code. Methodology, User’s Manual, Validation Report.