A boundary element–response matrix method for the multigroup criticality problems in the SP3 approximation

Annals of Nuclear Energy 42 (2012) 119–130

Contents lists available at SciVerse ScienceDirect

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

A boundary element–response matrix method for the multigroup criticality problems in the SP3 approximation V. Giusti ⇑, B. Montagnini Department of Mechanical, Nuclear and Production Engineering, Pisa University, Largo Lucio Lazzarino 2, I-56126 Pisa, Italy

a r t i c l e

i n f o

Article history: Received 22 September 2011 Received in revised form 29 November 2011 Accepted 30 November 2011 Available online 2 January 2012 Keywords: Neutron transport Simplified spherical harmonics Response matrix Boundary element method

a b s t r a c t The 3D xyz multigroup BERM-SP3 transport method is here presented. This method, based on the simplified spherical harmonics approximation (SPN), with N = 3, and the assumption of linearly anisotropic scattering, makes use of a Response Matrix (RM) solution procedure, coupled with the Boundary Element Method (BEM), by which the SP3 partial differential equations are reduced to a system of boundary integral equations in terms of partial currents. Numerical problems, all endowed with a complete set of data and the reference results obtained by well assessed transport codes such as the discrete ordinate code TORT and the Monte Carlo code MCNP, illustrate the accuracy and efficiency of the method. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Because of the recently published issue on the Simplified Spherical Harmonics (SPN) method (Holloway, 2010) and, in particular, of the state-of-art article by R.G. McClarren and the comprehensive paper by E.W. Larsen (other than several previous papers, which also give a summary of the foundations of the method), we refrain from recalling here the general SPN theory. The SPN differential equations are now widely accepted as a useful approximation of the spherical harmonics (PN) equations and have been also endowed with a set of asymptotic conditions of validity (Larsen et al. (1996); Morel et al. (1996); Brantley and Larsen (2000)) that has given a much firmer basis to the early intuitive approach by Gelbard (1960, 1961, 1962). Thus, the challenge is now mainly about the performances of the several methods that can be adopted to solve the system of 2G diffusion-like equations (where G is the number of the energy groups) to which the SPN approximation of the linear transport equation can be reduced. The case N = 3 is usually preferred for its simplicity, other than for the intrinsic limits of the approximation itself, which discourage considering, in general, the higher order approximations. As it concerns the numerical solution techniques, they are most often based on the classical finite elements, the mixed-dual elements and the nodal methods (e.g. Lemanska, 1981; Smith, 1986; Gamino, 1989, 1991; Lautard et al., 1999; Baudron and Lautard, 2007; Beckert and Grundmann, 2008; Hébert, 2010, to quote only a few of them). Another group of works relies on the Response ⇑ Corresponding author. Tel.: +39 050 2218026. E-mail address: [email protected] (V. Giusti). 0306-4549/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2011.11.021

Matrix (RM) method and even in this group quite different techniques are applied, to evaluate the RM of each region in which the physical system is subdivided. The specific RM solution procedure in the BERM-SP3 method here presented is based on the Boundary Element Method (BEM), by virtue of which the SP3 partial differential equations are transformed into a system of boundary integral equations, which is then arranged in a partial current form and solved by a moment (or weak) technique. Another RM method (Tada et al., 2010) does not make use of BEM and the region-wise SP3 differential equations are solved by means of staggered finite differences. A further method (Chiba, 2011) makes use of BEM (in its classical form; a mixed collocation-polynomial procedure is adopted in order to approximate the integral equations), but not of the RM approach, since it relies on a Newton iterative solution of the criticality problem for the whole multiplying system, taking advantage of a hierarchical domain decomposition in order to improve the computing efficiency. A third group of works is worth mentioning, which is based on a variant of the SPN method itself, the so-called AN method (Coppa and Ravetto, 1982; Ciolini et al., 2002; Giusti et al., 2010). This method, however, deserves a special treatment and will not be dealt with here. In the present paper, Section 2 introduces the SP3 multigroup differential equations and their interface and boundary conditions; a particular emphasys is given to the well known trick of considering the SP3 void boundary conditions as interface conditions with a purely absorbing outer medium. Section 3 deals with the transformation of the SP3 differential equations into both the classical and the partial current boundary integral form. The theory being essentially an extension of that developed in Cossa et al. (2010) for the

120

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

multigroup diffusion equations, the reader is referred to this paper for any detail concerning the numerical solution procedure. In the present paper, we have preferred to deal mainly with the applications and give (Section 4) a reasonably extended set of numerical results that can be obtained by the SP3 solution technique here developed. Conclusions are drawn in Section 5.

PG

g 0 ¼1

the computations. An ingenious trick (Bell and Glasstone, 1970) is that of assuming that the slowing down of the neutrons from the groups g0 < g (i.e. from the groups the energy of which is higher than that of group g) is nearly the same, if the absorption is small, as the slowing down from group g to all the groups g0 > g (i.e. to the groups with a lower energy):

2. The SP3 equations and the interface and boundary conditions

G X

The SPN equations are intrinsically approximate (in the sense that the approximation cannot be improved beyond some limit by letting N go to infinity) and, for fixed N, the accuracy of the solution is strongly problem dependent. This cannot be without consequences as regards the problems the method should or should not be applied to. For instance, applying the SPN method to a fuel assembly the cells of which are homogenized appears to be a quite reasonable practice. On the contrary, considering the same assembly as made of cells in which the fuel pins retain their own individuality should probably cause some concern and one would seek, at least now and then, for a confirmation by one of the current full transport methods, such as PN or SN. The basic limitations of the SPN (in particular, SP3) method have been here taken into consideration from the outset as regards the choice of the numerical examples, which are all referred to homogenized pin cells. We have also considered only the case of a linearly anisotropic scattering, since neglecting the anisotropy of the second order has appeared to us as a quite tolerable (even if not strictly necessary) restriction, if account is taken of the rather strong simplifications that we have admitted in the treatment of the linear anisotropic scattering itself (see below). With the notations of the paper by Beckert and Grundmann (2008), where the reader can also find a short introduction to the SP3 method, and with the above assumptions concerning the scattering, the SP3 equations for a homogeneous region Vi of the system are as follows:

g 0
h

i

~ i ðrÞ  ri / ~ i ðrÞ  2/ ~ i ðrÞ þ Di0;g D/ 0;g t;g 0;g 2;g h

i

G X

ri0;gg0 þ

g 0 ¼1

v

i g

k

!

mrif ;g0

~ i 0 ðrÞ ¼ 0; ~ i 0 ðrÞ  2/  / 0;g 2;g

ð1Þ

  G X 2 ~ i ðrÞ  9 / ~ i ðrÞ  2 þ rit;g / 0;g 5 2 2;g 5 g0 ¼1 h i ~ i 0 ðrÞ ¼ 0 ðg ¼ 1; . . . ; GÞ: ~ i 0 ðrÞ  2/  / 0;g 2;g

~ i ðrÞ Di2;g D/ 2;g

i 0;gg 0

r

þ

vig k

! i f ;g 0

mr

ð2Þ

~ i ðrÞ and / ~ i ðrÞ replace the SP3 moments Here, the unknowns / 0;g 2;g /i0;g ðrÞ and /i2;g ðrÞ as they can be derived from the plane geometry P3 equations, after the elimination of the odd-order moments and the application of the Gelbard assumptions. Namely,

~ i ðrÞ ¼ /i ðrÞ þ 2/i ðrÞ; / 0;g 2;g 0;g ~ i ðrÞ ¼ /i ðrÞ: / 2;g

ð3Þ

ri1;gg0 /i1;g0 and this would result in a remarkable increase of

ri1;gg0 /i1;g0 

As regards the physical constants, r is the total (macroscopic) cross section in group g, ri0;gg0 the zero-th moment of the differential scattering cross section from group g0 to group g, mrif ;g the ggroup fission cross section times the average number of the secondary neutrons, vg the g-group fission spectrum fraction, while k is the multiplication factor. No external (independent) sources are considered. As it concerns the diffusion coefficients Di0;g and Di2;g , a simplifying assumption introduced in the paper by Beckert and Grundmann (2008) is adopted here also. Namely, taking a full account of the group-wise first moments of the scattering cross section ri1;gg0 would imply considering an additional term of the type i t;g

ri1;g0 g /i1;g :

ð5Þ

g 0 >1

Note that this formula can be extended to the case in which upscattering processes are present, so that one has, more generally, G X

ri1;gg0 /i1;g0 

g 0 ¼1

G X

ri1;g0 g /i1;g :

ð6Þ

g 0 ¼1

Using this approximate formula, Beckert and Grundmann were able to avoid the explicit evaluation of the sum on the l.h.s. and obtain the following simple formula for the transport cross section

ritr;g ¼ rit;g 

G X

ri1;gg0 ; 2 6 g 6 G:

ð7Þ

g 0 ¼1

If g = 1, the simpler (and not subject to (6)) expression i tr;1

r

¼ rit;1  ri1;11

ð8Þ

will be here adopted. Then, one may set

Di0;g ¼

1 ; 3ritr;g

1 6 g 6 G:

ð9Þ

The diffusion coefficient Di2;g is simply given, in the present calculations, by the following formula

Di2;g ¼

9 ; 35rit;g

1 6 g 6 G;

ð10Þ

owing to our choice of neglecting the higher order scattering, which has been considered, as anticipated before, as a simplifying assumption of a minor relevance, if compared with the approximation implied by Eq. (6). Now we turn to the conditions at the interface, S say, between two adjacent regions Vi and Vj. The following conditions of continuity for the SP3 moments, let us say (somewhat inappropriately) the SP3 ‘‘fluxes’’ and ‘‘net normal currents’’:

~ i ðrS Þ ¼ / ~ j ðrS Þ; / 0;g 0;g

ð11Þ

~ j ðrS Þ; ~ i ðrS Þ ¼ / / 2;g 2;g

ð12Þ

and

 Di0;g

ð4Þ

2;g

G X

 Di2;g

~i @/ 0;g @niS ~i @/ 2;g @niS

ðrS Þ ¼ Dj0;g ðrS Þ ¼ Dj2;g

~j @/ 0;g @njS ~j @/

2;g

@njS

ðrS Þ;

ð13Þ

ðrS Þ;

ð14Þ

where niS ðnjS Þ is the outward normal of the region Vi (Vj), with niS ¼ njS , are widely accepted and we shall agree with this rule, which is corroborated by the asymptotical and variational analysis in the references of Section 1 and is common practice for the SP3 finite element or nodal methods in the case of an interface. ‘‘Partial currents’’ are introduced only if considering Marshak-like boundary conditions with an outer vacuum. In the case of a RM method, however, the rôle of representing the coupling of a region with its neighbors is always played by the ingoing and outgoing partial currents (the inverted commas will now be suppressed) and the interface

121

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

conditions must be written accordingly. The problem then arises about what definition is to be used for the partial currents. Will the following simple (and quite profitable in our response-matrix framework) Marshak-like definitions of the outward (+) and inward () normal partial currents, as suggested by the P1 theory, extended also to the second moment,

~i @/ 0;g eJ i ðrS Þ ¼ 1 / ~ i ðrS Þ  1 Di ðrS Þ; 0;g n;0;g 0;g 4 2 @niS eJ i ðrS Þ n;2;g

~i 1 ~i 1 i @/ 2;g ¼ / ðrS Þ; 2;g ðrS Þ  D2;g 4 2 @niS

ð15Þ

ð16Þ

be adequate? At least from a formal point of view they are, since they preserve, as linear combinations of the SP3 fluxes and net currents, the basic continuity conditions, Eqs. (11)–(14). Namely, if Eqs. (11) and (12) are multiplied by 1/4 and Eqs. (13) and (14) by 1/2 and the equations so obtained are summed or subtracted, the result can be written, if account is taken of the relation niS ¼ njS , as follows:

eJ i ðrS Þ ¼ eJ j ðrS Þ; n;0;g n;0;g

ð17Þ

eJ i ðrS Þ ¼ eJ j ðrS Þ; n;2;g n;2;g

ð18Þ

which are the usual conditions for the partial currents. Conversely, by summing and subtracting the ± expressions in Eqs. (15) and (16) one gets

h i ~ i ðrÞ ¼ 2 eJ iþ ðrS Þ þ eJ i ðrS Þ ; / 0;g n;0;g n;0;g h i ~ i ðrÞ ¼ 2 eJ iþ ðrS Þ þ eJ i ðrS Þ ; / 2;g n;2;g n;2;g  Di0;g 

Di2;g

~i @/ 0;g @niS ~i @/

2;g @niS

ei ðrS Þ ¼ eJ iþ n;0;g ðrS Þ  J n;0;g ðrS Þ; ðrS Þ ¼

eJ iþ ðrS Þ n;2;g j

 eJ i n;2;g ðrS Þ

Remark 1. It can be observed that more sophisticated expressions for the normal partial currents, e.g. those in Beckert and Grundmann (2008), namely

~i @/ 3 ~i 0;g eJ i ðrS Þ ¼ 1 / ~ i ðrS Þ  1 Di / ðrS Þ; ðrS Þ  n;0;g 0;g 0;g 4 2 16 2;g @niS ~i @/ 21 ~ i 2;g eJ i ðrS Þ ¼  3 / ~ i ðrS Þ  1 Di / ðrS Þ; ðrS Þ þ 2;g n;2;g 0;g 80 2 80 2;g @niS

ð23Þ ð24Þ

also allow to obtain interface conditions, as in Eqs. (17) and (18), that are equivalent to the classical conditions of Eqs. (11)–(14). ej For example, if Eq. (23) for eJ i n;0;g and J n;0;g are used in Eq. (17) and the resulting equations are summed up, the terms containing the ~i  normal derivatives drop out and we remain with ð4=3Þ/ 0;g

~ i ¼ ð4=3Þ/ ~j  / ~ j . Proceeding in the same way with Eq. (18) / 2;g 0;g 2;g ~ i ¼ ð1=7Þ/ ~j  / ~ j . Thus, after subtracting ~i  / we obtain ð1=7Þ/ 0;g 2;g 0;g 2;g i ~ j , which yields ~ the two equations, we get ð25=21Þ/ ¼ ð25=21Þ/ 0;g

0;g

Eq. (11) and, by the aid of the previous relation, Eq. (12). The procedure for obtaining Eqs. (13) and (14) is even simpler. On the other hand, Eqs. (23) and (24) should obviously be more suitable than Eqs. (15) and (16) if they are requested to provide physically accurate values of the individual partial currents and should certainly be preferred, in particular, when dealing with a boundary condition with an outer vacuum, the case that we have just excluded (except if the peripherical flux is almost negligible) by introducing the absorbing shell.

ð19Þ ð20Þ ð21Þ ð22Þ

and similarly for the V region. Then, it is immediately shown that Eqs. (17) and (18) imply Eqs. (11)–(14). In other words, at an inner interface, the SP3 solution must satisfy Eqs. (17) and (18) as well as Eqs. (11)–(14). This does not mean that definitions (15) and (16) provide a physically very realistic value of each partial current, if individually considered. Rather, even if the above partial currents may be individually poorly approximate, their sum and their difference must give, nevertheless, considerably more accurate values for ~i ; / ~ i and Di @ / ~ i =@ni , Di @ / ~ i =@ni . the quantities / 0;g 2;g 0;g S 2;g S 0;g 2;g On account of this argument we have decided: (a) when dealing with an inner interface, to use the elementary formulas (15) and (16) for the partial currents; (b) to replace the void surrounding the multiplying system, assumed to be convex, by a perfectly absorbing medium. By this way, when performing the calculation procedure, only interface conditions are to be considered. If introducing an absorbing medium of an infinite extent seems to be unpractical (which is not always true, if the boundary element technique is fully exploited), an absorbing shell with a thickness of few mean free paths will suffice to reduce the flux at its outer boundary to a very small value. The vanishing of the incoming partial currents given by the elementary formulas Eqs. (15) and (16) can then be there applied without compromising the accuracy. The introduction of the absorbing region, a well known trick also used in the case of the PN calculations (see for instance Davison and Sykes, 1957; Fletcher, 1983), may be not so costly in terms of CPU time; moreover, in many realistic reactor problems the flux at the outer reflector boundary is so small that Eqs. (15) and (16) can be directly used there to formulate the classical Marshak void conditions, without adding any absorbing shell.

3. The partial current boundary integral form of the SP3 equations The partial differential SP3 equations for multiplying systems without external sources (Eq. (1)) can be set in a group-wise block matrix form, namely

~ i ðrÞ þ D/ g

G X

~ i 0 ðrÞ ¼ 0 ðg ¼ 1; . . . ; GÞ; Q igg0 / g

ð25Þ

g 0 ¼1

h i ~ i ðrÞ; / ~ i ðrÞ ; ~ i ðrÞ ¼ / where / g 0;g 2;g " Q igg0

¼

Q i00;gg0

Q i02;gg0

Q i20;gg0

Q i22;gg0

Q i00;gg0 ¼ 

Q i20;gg0 ¼

Aigg0

# ðg; g 0 ¼ 1; . . . ; GÞ;

Aigg0

;

Q i02;gg0 ¼ 2

i 2 Agg0 ; 5 Di2;g

Q i22;gg0 ¼ 

Di0;g

g ¼ 1; . . . ; G and

Di0;g

;

i 1 Bgg0 ; 5 Di2;g

ð26Þ

ð27Þ

ð28Þ

with

Aigg0 ¼ rit;g dgg0  ri0;gg0 

vig k

mrif ;g0 ;

Bigg0 ¼ 9rit;g dgg0  4ri0;gg0  4

v

i g

k

mrif ;g0

ð29Þ ð30Þ

(the criticality parameter k will be now understood). One can avoid considering the group small blocks separately and rewrite Eq. (25) as follows:

~ i ðrÞ þ Q i / ~ i ðrÞ ¼ 0; D/ where after a renumbering of the components,

ð31Þ

122

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

 i  ~ ðrÞ; . . . ; / ~ i ðrÞ T ~ i ðrÞ ¼ / / 1 G h iT ~ i ðrÞ; / ~ i ðrÞ; . . . ; / ~ i ðrÞ; / ~ i ðrÞ ¼ / 0;1 2;1 0;G 2;G h iT ~ i ðrÞ; / ~ i ðrÞ; . . . ; / ~i ~i ¼ / 1 2 2G1 ðrÞ; /2G ðrÞ

ð32Þ

h iT ~ i ðrÞ ~ i ðrÞ ¼ / i.e. in short, / , while the 2Gx2G matrix Q i ¼ c c¼1;2G h iT Q icc0 0 is written accordingly. c;c ¼1;2G

Eq. (31) is formally identical to a multigroup diffusion equation system (where the number of the energy groups is multiplied by two); in particular, it is identical to Eq. (3) in the paper Cossa et al. (2010). The solution procedure shown in that paper has been strictly followed in the present SP3 framework and only needs to be shortly summarized. With the intent of diagonalizing the equation system one must consider the eigenvalue eqaution (the index i is momentarily suspended)

Q n ¼ kn;

ð33Þ

which is assumed to have 2G distinct real or complex eigenvalues kh. The corresponding eigenvectors nh then constitute a basis of the space C2G . The matrix N having the (normalized) nh as columns is non singular and diagonalizes Q:

N1 Q N ¼ K;

ð34Þ

where K = diag[k1, . . . , k2G]. Multiplying Eq. (31) by N

1

and setting

~ WðrÞ ¼ N1 /ðrÞ;

ð35Þ

we get

XS(rS) being the angle of aperture of the tangent cone at rS. If rS is a smooth point, then obviously XS(rS) = 2p and c(rS) = 1/2. Let now the elements of N be denoted by nch (these nch’s are the components of the eigenvectors, nh, of Q) and let nH ch be the elements of N1. Setting 2G   X   ^ cc0 r; r0 ¼ ^ h r; r0 nH 0 ; / nch w hc S S

ð42Þ

h¼1

and using Eq. (35) and the reciprocal equation, so that, in terms of the components,

wh ðrÞ ¼

2G X

nH hc0 /c0 ðrÞ;

/c ðrÞ ¼

c0 ¼1

2G X

nch wh ðrÞ;

ð43Þ

h¼1

Eq. (40) can be given the following form

cðrÞ/c ðrÞ þ

" # 2G Z ^ cc0  X   0     @/ 0 ^ cc0 r; r0 @/c0 r0 dS0 ¼ 0 r; r r /  / 0 c S S S @n0S @n0S S c0 ¼1 S

ðc ¼ 1; . . . ; 2GÞ:

ð44Þ

If the boundary values of the SP3 moments /c and their normal derivatives @/c/@nS are known, then the values of the /c’s at any point r, in particular at any interior point, can be determined by performing a simple quadrature. Let now r be taken on the boundary itself, r = rS say. Then Eq. (44) become a system of boundary integral equations, namely

cðrS Þ/c ðrS Þ þ

2G Z X

c0 ¼1 S

" # ^ cc0    0  @/c0  0  @/ 0 ^ cc0 rS ;r0 dS0 ¼ 0 r ;r r r /  / 0 S c S S S @n0S @n0S S

ðc ¼ 1; . .. ;2GÞ:

ð45Þ

with an appropriate condition at infinity, the specific form of which depends on the reality and sign of kh. While referring the reader to the above mentioned paper Cossa et al. (2010) for a more detailed analysis, it will suffice here to consider the simplest case, ~ h of Eq. kh ¼ a2h , with ah real. If we assume that the solution w (38) must be vanishing at infinity we easily obtain

This system represents the SP3 equations in their classical boundary integral form. If, for instance, the normal derivatives @/c/@nS are assigned (which corresponds to a Neumann boundary condition) then we can solve the system with respect to the /c’s. Once the boundary values of both the quantities /c and @/c/@nS are known, the moments /c at any interior point of the region can be obtained by Eq. (44), as shown above. Otherwise one can assign the boundary values of the /c’s (which corresponds to a Dirichlet boundary condition), solve Eq. (45) with respect to the normal derivatives and again proceed to the determination of the interior values of the moments and, consequently, of the physical flux. The following modified form of Eq. (45) is particularly suitable in the case of the response matrix approach. We restore the index i, which specifies the single homogeneous region, Vi, to which the boundary integral procedure is applied. Then, let eJ i; ðrS Þ be the out-

~ h ðr; r0 Þ ¼ expðah jr  r jÞ ; w 4pjr  r0 j

ward, respectively inward partial current as given by Eqs. (15) and (16) (here again, the vector components eJ i; ; eJ i; are now de-

DWðRÞ þ KWðRÞ ¼ 0;

ð36Þ

actually a set of 2G uncoupled equations in terms of the components, wh(r) say, of W(r):

Dwh ðrÞ þ kh wh ðrÞ ¼ 0;

ðh ¼ 1; . . . ; 2GÞ:

ð37Þ

Let us also consider the following inhomogeneous equations

^ h ðr; r0 Þ þ kh w ^ h ðr; r0 Þ þ dðr  r0 Þ ¼ 0; Dr w

ðh ¼ 1; . . . ; 2GÞ;

ð38Þ

0

ð39Þ

n;0;c

which has clearly the form of the Green function of the diffusion equation in an infinite, absorbing and scattering medium. The hth Eq. (37) is then written with r replaced by r0 and multiplied by ^ h ðr; r0 Þ. An integration on r0 followed by the application of the w Green’s second identity allows to shift the Laplace operator in such ^ h r20 w is replaced by the integral of a way that the integral w r h 2^ wh rr0 wh plus a surface term. Use of Eq. (38) (with r and r0 interchanged) then leads, after a few manipulations, to the following integral relationship:

# Z " ^ @ wh  0   0  ^  0  @wh  0  cðrS Þwh ðrS Þ þ r; rS wh rS  wh r; rS r dS0 @n0S @n0S S S ¼ 0 ðh ¼ 1; . . . ; 2GÞ; Z

dðrS  r0 Þ dV 0 ¼

ð40Þ

V

XS ðrS Þ 4p

;

n;2;c

~ i ’s in Eq. (32)): noted by means of a unique index c, just as for the / c

~i @/ c eJ i ðrS Þ ¼ 1 / ~ i ðrS Þ  1 Di ðrS Þ: n;c 4 c 2 c @niS

ð46Þ

We have, reciprocally,

h i ~ i ðrS Þ ¼ 2 eJ iþ ðrS Þ þ eJ i ðrS Þ ; / n;c n;c c ~i @/ c @niS

ðrS Þ ¼ 

i 1 heiþ J n;c ðrS Þ  eJ i n;c ðrS Þ : i Dc

ð47Þ ð48Þ

Substituting these expression for the moments and their normal derivatives into Eq. (45) and introducing the new kernels

where

cðrS Þ ¼

n;c

ð41Þ

^i 0  @/      cc bJ i 0 rS ; r0 ¼ 1 / ^ i 0 rS ; r0  1 rS ; r0S ; cc n;cc S S i 0i 2 @nS 4Dc

ð49Þ

the following partial current form of the SP3 boundary integral equation system is obtained

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130 2G Z X     1 bJ iþ 0 rS ; r0 eJ iþ 0 r0 dS0 ¼  1 cðrS ÞeJ i ðrS Þ cðrS ÞeJ iþ n;c ðrS Þ þ n;cc S n;c S n;c 2 2 c0 ¼1 S 2G Z X     bJ i 0 rS ; r0 eJ i 0 r0 dS0 ðc ¼ 1; . . . ; 2GÞ: þ ð50Þ n;cc S n;c S S

c0 ¼1

The convenience of this equivalent form of the boundary intei gral equations is evident: if the partial currents eJ i n;c entering V are known, the r.h.s. of Eq. (50) is also known and we can solve the equations with respect to the outgoing partial currents, eJ iþ n;c . In other words, the response of the Vi region to the injection of neutrons from the neighboring regions has been determined. The actual numerical solution of the above equations is, however, not trivial, even in the case of a Vi region with a simple, regular shape. The weak (or moment) method used in this paper, as well as in the previous papers concerning the multigroup diffusion equations (Maiani and Montagnini, 1999, 2004; Cossa et al., 2010), is based on the following expansion:

eJ i ðrS Þ ¼ n;c

M X

eJ i wm ðrS Þ; n;c;m

ð51Þ

m¼0

where the wm(rS) constitute a complete, orthonormal set of functions (here truncated after the M-th term) on the boundary surface, Si, of Vi. If this expansion is introduced into Eq. (50) and the scalar product with the weight functions wm0 is taken, then the integral equations give rise to the following (approximate) linear system 2G X M X     1 eiþ bJ iþ 0 0 rS ; r0 eJ iþ 0 0 r0 ¼  1 eJ i ðrS Þ J n;c;m ðrS Þ þ n;cc ;mm S n;c ;m S 4 4 n;c;m c0 ¼1 m0 ¼0

þ

2G X M X

c0 ¼1

    bJ i 0 0 rS ; r0 eJ i 0 0 r0 n;cc ;mm S n;c ;m S

m0 ¼0

ðc ¼ 1; . . . ; 2G;

m ¼ 0; . . . ; MÞ;

ð52Þ

where



 0

bJ i 0 0 rS ; r ¼ n;cc ;mm S

Z Z Si

Si

bJ i 0 n;cc



 0

  rS ; rS wm ðrS Þwm0 r0S dS dS0 :

ð53Þ

It is to be observed that

Z Si

1 ei cðrS ÞeJ i ; J n;c ðrS Þwm ðrS Þ dS ¼ 2 n;c;m

ð54Þ

since the set of the non-smooth boundary points is a set of zero measure. Eq. (52) can be given a more compact matrix form

Miþ~Jiþ ¼ Mi~Ji ;

ð55Þ

i±

where the M matrices have the following elements:

1 bi M i cc0 ;mm0 ¼  dcc0 dmm0 þ J n;cc0 ;mm0 : 4

ð56Þ

Upon inversion of the Mi+ matrix we obtain

~Jiþ ¼ Ri~Ji ; i

ð57Þ i+ 1

where R = [M ]

i

M

i

is the response matrix of the V region.

Remark 2. If the region Vi is a polyhedron (for instance a parallelepiped) it is convenient that the weight functions wm are splitted into products of orthonormal polynomials (e.g. Legendre polynomials) defined on each face of it. For such details see Cossa et al. (2010), where the procedures for the numerical calculation of the resulting integrals are also shown. Once the response matrices Ri have been evaluated for all the regions of the system (including the purely absorbing regions, if existing, that replace the void) the standard iterative RM method

123

is applied: if the SP3 partial currents entering the Vi regions are known, the application of the response matrices gives the SP3 partial currents exiting from them and the next step is ready to start. The k eigenvalue is updated at any step of such an outer iteration. The details of the iterative process, including the acceleration techniques, can be also found in Cossa et al. (2010).

4. Numerical results All the numerical examples shown in this section have been performed with a PC based on an Intel i7–860 CPU and equipped with 8 Gb of RAM. Moreover, in order to better exploit the multithread characteristics of the CPU, many routines of the source Fortran files take advantage of the easy parallelization offered by the OpenMP™ Fortran Compiler directives. Finally the executable BERM-SP3 file has been obtained compiling the source Fortran files Ò with the Intel Fortran Compiler, version 11.1. The examples of this section have been solved, for comparison, other than with the present BERM-SP3 code, which implements the method described in this paper, by the worldwide used Monte Carlo code MCNP, version 5 (X-5 Monte Carlo Team, 2003) and the well known 3D discrete ordinates code TORT (Rhoades and Simpson, 1997) in the S16 approximation. However, while for MCNP it was possible to use an executable file with parallel capabilities (supporting the MPI communication protocol), for TORT only a sequential version of the executable file was available. Thus, the computing times, reported hereafter, will have to be compared with caution.

4.1. Example 1: A 2D, one energy group, linearly anisotropic scattering benchmark problem This problem, proposed by Hébert, is not representative of a real situation but is designed, as stated by the author himself, in order to magnify the transport and anisotropic effects (Hébert, 2010). A homogeneous rectangular core is surrounded by a reflector which is surrounded, in its turn, by a perfectly absorbing medium. Finally, on the external surface of the absorbing medium a vacuum boundary condition is applied (see Fig. 1). MCNP, TORT and BERM-SP3 are all 3D computational codes, thus the 2D geometry of this benchmark problem had to be modified: for BERM-SP3 and TORT a finite height along the z-axis was considered by applying reflective boundary conditions on the surfaces orthogonal to this axis; for MCNP the geometry was, instead, assumed as infinite along the z-axis direction. Table 1 shows the cross sections of the three different materials of the benchmark problem. The MCNP calculation has been performed using the same cross sections reported in Table 1. The linearly anisotropic Legendre rep lÞ=2, where l is the resentation of the scattering, f ðlÞ ¼ ð1 þ 3l  its average value, has been cosine of the scattering angle and l converted into a histogram probability density function (PDF) with 1000 steps of equal area. However, it has been shown that, with this approach, MCNP can match exactly (within statistics) analyti 6 1=3. If l  > 1=3, a value of l exists below cal solutions only if l which the PDF becomes negative, a case that cannot be handled by a stochastic code as MCNP (Brown and Barnett, 2006). In such case the Legendre representation of the scattering cannot be converted into a histogram PDF without making some approximation, in order to avoid the negative values of f(l). Finally, to estimate the multiplication constant k a total number of 2500 active cycles, each one of 250,000 particles, has been used for the MCNP calculation. The MCNP value of the multiplication constant k has then been taken as reference.

124

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

4.2. Example 2: The LWR mixed-oxide benchmark of the American Nuclear Society (ANS)

Fig. 1. The geometry of the 2D, one energy group benchmark problem proposed by Hebért (dimensions are in cm).

Table 1 Cross section for the 2D, one energy group benchmark problem proposed by Hebért. Mixture

R (cm1)

Rs0 (cm1)

Rs1(cm1)

mRf (cm1)

1 2 3

0.025 0.025 0.075

0.013 0.024 0.0

0.0 0.006 0.0

0.0155 0.0 0.0

As regards TORT, we have considered it as a supplementary reference code and even, in Example 2 of Section 4.2, as the unique reference code, owing to the highly anisotropic scattering which ruled out MCNP. We have therefore always used a rather high angular approximation order (S16) and adopted a fine spatial mesh (2.0  2.0  1.0 cm3 in this example). The multiplication constants k, calculated by BERM-SP3 by expanding the partial currents into Legendre polynomials up to the Lth order, with L = 3 and 4, are compared in Table 2 with the values obtained by MCNP and TORT. The BERM-SP3 calculation performed with the L = 3 expansion provides the same value of the multiplication constant as that with L = 4, a value which differs from that obtained by MCNP by 119 pcm. As expected, the TORT value is more accurate and differs from that of MCNP by only 8 pcm. The paper by Hébert (2010) reports a TRIVAC calculation in the SP3 approximation, which gives a value of k identical, at pcm level, with that obtained by BERM-SP3. It appears, therefore, that, at least in this case, the specific SP3 numerical solution method (finite element method for TRIVAC and boundary element method for BERM-SP3) does not play a substantial rôle as regards the accuracy.

This 2D benchmark is designed to provide a simple calculational problem for a mixed-oxide lattice (Gehin and Primm, 1997). The model chosen, see Fig. 2, was that of a pressurized water reactor (PWR) with a 17  17 assembly lattice of pins, all of them with the same composition. The assembly has 24 control rod guide tubes. In the configuration here studied, these tubes are simply filled by water and thus do not contain any burnable poison rods. A central instrumentation channel is also present. The assembly is assumed to be in an infinite array of similar assemblies, so that reflective boundary conditions can be used. Moreover, no axial buckling has been assumed, i.e. the assembly has been considered infinite in the z-axis direction. Since BERM-SP3 and TORT can only handle 3D problems, a supplementary reflective boundary condition on the surfaces orthogonal to the z-axis has been applied, as it has been done in Example 1 of Section 4.1. Starting from the composition of the materials as specified in Gehin and Primm (1997), by means of suitable cell calculations performed with the SCALE6 code (ORNL, 2009), a set of 8 energy group cross sections, with linearly anisotropic scattering, as reported in Appendix A, has been derived for the homogenized fuel cells, control rod cells and instrumentation channel cell. This set of cross sections has been used for all the calculation perfomed with the BERM-SP3 and TORT.  larger than 1/3 came Unfortunately, in this problem values of l out. Thus, it was not possible to run a MCNP calculation using exactly the same set of cross sections (see Example 1, Section 4.1). The MCNP calculation was done all the same, by implementing the exact geometry (no cell homogenization) and using the most up-to-date continuous-energy cross section data libraries. The MCNP value of the multiplication constant k has been estimated using 1500 active cycles of 100,000 particles each. Of course, this very accurate calculation can be only used to show the full transport results to which the other methods should converge and is here reported only for completeness. As anticipated in Section 4.1 we have decided to make use of TORT as reference code in place of MCNP.

Table 2 The multiplication constant k for the 2D, one energy group with linearly anisotropic scattering benchmark problem. Code MCNP (Ref.) TORT S16b BERM-SP3c a b c

L

k

Dk (pcm)

CPU time (s)

3 4

0.99233a 0.99225 0.99114 0.99114

– 8 119 119

4210.0 467.0 2.2 3.3

With an estimated standard deviation of 0.00003. With meshes of 2.0  2.0  1.0 cm3. With a computational cell size of 10  10  10 cm3.

Fig. 2. The 2D geometry, after cell homogenization, of the 17  17 LWR mixedoxide assembly of the benchmark problem proposed by ANS.

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

125

Table 3 shows the values of the multiplication constant k obtained with the different numerical codes. In this case, the values obtained by BERM-SP3 with both L = 3 and L = 4 agree very well with that provided by TORT. The cell homogenization and the cross section collapsing procedures prevented the two deterministic codes from achieving a value of k with an error lower than 200 pcm with respect to that estimated by MCNP, an amount that is not to be considered, however, as exceedingly high, if account is taken of the involved approximations. Fig. 3 shows the values of the fuel cell fission power as obtained by TORT and BERM-SP3, other than the indicative results of MCNP. The maximum relative error between BERM-SP3 and TORT is only 0.43%, a fact that should be taken as a good support of the reliability of the SP3 method. The maximum discrepancy with respect to the MCNP values is also surprisingly small, i.e. 2.05% and 1.61% for TORT and BERM-SP3 respectively. 4.3. Example 3: A 3D reactor core loaded with mixed-oxide fuel assemblies This benchmark problem (Smith et al., 2005) represents a 3D reactor core made of 16 fuel assemblies (quarter core symmetry), half of which contain mixed-oxide fuel rods, and completely surrounded by a water reflector, see Fig. 4. Each fuel assembly is made of 17  17 square pin cells the side lenght of which is 1.26 cm, (for a complete description of the benchmark original geometry see Smith et al. (2005)). A vacuum boundary condition is applied on the external surface of the reflector. Moreover, in order to reduce the computational burden, an axial symmetry is assumed. Fig. 5 shows the vertical section of 1/8 of the core and the boundary conditions that have to be applied. This benchmark is a modification of a preceding OECD benchmark (Smith et al., 2003) with respect to which the height of the core has been reduced and, to replace the control rod guide composition in certain parts of the reactor, a new material describing the control rods has been defined. Three different configurations of the reactor core are proposed: Unrodded, Rodded A and Rodded B. In the Unrodded configuration the control rods are fully inserted into the axial reflector above the fuel assemblies. The second configuration, Rodded A, takes the Unrodded configuration as its starting point. In this configuration the control rods are inserted 1/3 of the way into the inner UO2 fuel assembly, as indicated by the single hatched region in Fig. 5. Similarly, the Rodded B configuration also takes the Unrodded configuration as its starting setting. In this configuration, however, the control rods are inserted 2/3 of the way into the inner UO2 assembly and 1/3 of the way into both MOX assemblies, as indicated by the single and double hatched regions in Fig. 5. The original benchmark was devised in order to test deterministic transport methods and codes on realistic problems without any spatial homogenization. However, according to what has been stated in Section 2, in this paper the different pin cells have been spatially homogenized, weighting the cross section over the volumes, with a preliminary estimate of the neutron flux within each

Fig. 3. Map of the fuel cell fission power for the 17  17 LWR mixed-oxide assembly of the benchmark problem proposed by ANS. For symmetry reasons only 1/8 of the assembly is shown. The (only indicative) MCNP values are reported in italics.

Table 3 The multiplication constant k for LWR mixed-oxide benchmark proposed by ANS. Code TORT S16a (Ref.) BERM-SP3b MCNP a b c

L

k

Dk (pcm)

CPU time (s)

3 4

1.30897 1.30890 1.30887 1.30658c

– 7 10 –

243 69 180 96,580

With meshes of 0.42  0.42  0.75 cm3. With a computational cell size coincident with that of the fuel pin cell. With an estimated standard deviation of 0.00006.

Fig. 4. Horizontal cross section of the 3D reactor core loaded with mixed-oxide fuel assemblies.

cell. The sets of cross section so obtained are reported in Appendix B. The assumption of isotropic scattering made it possible to run also MCNP calculations with the same sets of cross sections used

126

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

was evaluated. This was responsible, in the case of the MCNP calculation, of a longer computing time (about five times) with respect to that of the Unrodded and Rodded A configurations. The computing time is greater than before, since MCNP had to evaluate and store the contribution of each particle to the neutron flux for all the computational cells where the fission power had to be estimated. The MCNP values of the axially-integrated fission power, assumed as reference values, are compared in Table 5 with those obtained by BERM-SP3 and TORT. The BERM-SP3 values agree very well with those provided by MCNP, the maximum error (in correspondence of the farthest fuel cell from the centre of the reactor core along the diagonal crossing the two UO2 fuel assemblies) being of 1.3%. However, despite of the better agreement concerning the value of the multiplication constant k, the TORT values of the axially-integrated fission power show larger errors, with respect to the MCNP values, than those shown by BERM-SP3. In the case of TORT, to get more accurate values an even finer mesh than that used in the present case is probably required, obviously accepting an increasing of the computing time. Fig. 5. Vertical cross section of the 3D reactor core loaded with mixed-oxide fuel assemblies (dimensions are in cm). In the Unrodded configuration the control rods are inserted only into the upper reflector; the Rodded A configuration is characterized by the insertion of the control rods 1/3 of the way into the inner UO2 fuel assembly (single hatched region); in the Rodded B configuration the control rods are inserted 2/3 of the way into the inner UO2 fuel assembly and 1/3 of the way into the two MOX fuel assemblies (single and double hatched regions).

with the deterministic codes. For this reason, for the three configurations described above, the MCNP estimated values of the multiplication constant k, obtained using a total number of 2000 active cycles of 500,000 particles each, have been again considered as the reference values. The values of the multiplication constant k for the three configurations are summarized in Table 4. It is interesting to note that, for BERM-SP3, the L = 3 expansion of the partial currents gives a value of the multiplication constant k that differs of only few pcm from the value obtained with a L = 4 expansion. On the other hand, the computing time is reduced by a factor 1.5–1.8, depending on the configuration. For the Rodded B configuration, other than the multiplication constant k, also the fission power in the fuel cells along the two main diagonals (one of which crosses only the UO2 fuel assemblies, while the other just the two MOX fuel assemblies, see again Fig. 4),

Table 4 The multiplication constant k for a 3D reactor core loaded with mixed-oxide fuel assemblies. Code

a b c

L

k

Dk (pcm)

CPU time (m)

Unrodded configuration MCNP (Ref.) TORT S16b BERM-SP3c 3 4

1.14212a 1.14205 1.14159 1.14155

– 7 53 57

800 205 27 49

Rodded A configuration MCNP (Ref.) TORT S16b BERM-SP3c 3 4

1.12582a 1.12576 1.12496 1.12496

– 6 86 86

796 212 21 35

Rodded B configuration MCNP (Ref.) TORT S16b BERM-SP3c 3 4

1.07111a 1.07046 1.06899 1.06903

– 65 212 208

2773 206 32 48

With an estimated standard deviation of 0.00004. With meshes of 0.63  0.63  1.19 cm3. With a computational cell size coincident with that of the fuel pin cell.

5. Conclusions A multigroup SP3 boundary element–response matrix code, BERM-SP3, directly stemming from the diffusion code BERM (Cossa et al., 2010), previously developed, has been presented. As known, the boundary integral approach allows to transform the original 3D differential problem into a 2D problem of integral type, the solution of which is looked for on the boundary of the domain that is studied, clearly an important advantage with respect to, e.g. the finite element method. This achievement is counterbalanced, however, by the fact that the boundary integral method, particularly in its weak form, entails the evaluation of many multiple integrals involving the Green’s functions of the problem. Efforts have therefore been made to alleviate such a drawback, both by performing as many analytical integrations as possible and by using ad hoc recurrence formulas (see again Cossa et al. (2010)). The general architecture of the method is simple: the multiplying system (a fuel assembly, or even a reactor) is represented as an array of homogeneous regions and, for each region, the set of SP3 boundary integral equations is solved, in order to obtain the outgoing SP3 flux moments as a response to the ingoing moments coming from the adjacent regions. To this purpose, the classical boundary equations are transformed into a partial current form, which is very suitable for the application of the response matrix technique (Maiani and Montagnini, 1999, 2004). The solution procedure, as implemented in the present BERM-SP3 code, is not essentially different from that of the above mentioned BERM diffusion code, if exception is made of the void conditions (if any) at the system boundary. In this case, the vacuum is replaced by a perfect absorber (which is therefore included in the calculation) and the void conditions replaced by plain interface conditions in terms of a couple of P1-like partial currents. Although the definition of the latter quantities is entirely formal, the resulting interface continuity conditions are algebraically equivalent to the classical continuity conditions for the SP3 fluxes and net currents. The whole problem setting (partial current boundary integral equations and interface conditions) is then equivalent to the classical one. Once the response matrix of each region of the system is evaluated, the general response matrix technique can be applied, in order to determine the multiplication constant of the system and the distribution of the physical flux. The numerical examples show that the BERM-SP3 code is accurate, efficient and robust. In some cases, the results turn out to be surprisingly good. It must be recognized, however, that the numerical benchmarks reported in this paper are belonging to a class of problems for which the SP3 approximation is well suited.

127

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

Table 5 3D reactor core loaded with mixed-oxide fuel assemblies: axially-integrated fission power in the homogenized fuel cells along the two main diagonals of the four fuel assemblies. Reference values obtained by MCNP are compared with those obtained by BERM-SP3 and TORT. For the latter codes only the relative error with respect to the MCNP value is actually shown. UO2 diagonal

a b

MOX diagonal

MCNPa

TORT S16 (%)

BERM-SP3 (%)

MCNPb

TORT S16 (%)

BERM-SP3 (%)

0.06737 0.06677 0.06574 0.0631 0.0612 0.06233 0.05914 0.05486 0.05006 0.04502 0.04104 0.03548 0.03403 0.03689 0.03977 0.04073 0.03474 0.03187 0.02602 0.02295 0.01805 0.0135 0.01304 0.01628

0.6 0.3 1.2 3.7 2.6 1.6 1.3 2.2 2.9 0.7 3.6 5.7 6.0 4.6 1.8 1.5 0.2 1.2 1.8 0.5 1.1 0.5 1.1 0.2

0.3 0.2 0.3 0.5 0.4 0.4 0.3 0.4 0.4 0.1 0.2 0.1 0.1 0.4 0.6 0.6 0.7 0.7 0.9 0.9 0.9 0.7 0.5 1.3

0.03452 0.02906 0.02607 0.03351 0.03565 0.03854 0.04316 0.04526 0.05426 0.05087 0.05636 0.05276 0.05276 0.05636 0.05087 0.05426 0.04526 0.04316 0.03854 0.03565 0.03351 0.02607 0.02906 0.03452

3.1 0.7 0.1 3.3 0.7 1.4 1.1 0.9 3.3 1.3 3.4 1.5 1.5 3.4 1.3 3.3 0.9 1.1 1.4 0.7 3.3 0.1 0.7 3.1

0.4 0.1 0.1 0.1 0.2 0.1 0.3 0.1 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.1 0.3 0.1 0.2 0.1 0.1 0.1 0.4

The maximum uncertainty is 0.07% (1r). The maximum uncertainty is 0.05% (1r).

Appendix A The following tables show the eight energy group sets of cross sections, with linearly anisotropic scattering, used for the three materials of the example reported in Section 4.2. These sets have been employed for all the simulations with the BERM-SP3 and TORT codes but not  6 1=3 (see Section 4.1). with the MCNP code, owing to the violation of the condition l

Group 1

2

Mixture 1: Homogenized fuel cell Rt 2.26710E01 4.37767E01 mRf 1.32385E02 1.77845E03 v 5.98360E01 3.50904E01

3

4

5

6

7

8

8.49347E01 1.83466E03 5.07305E02

1.06286E+00 1.79311E02 5.06416E06

1. 04976E+00 3. 38421E02 5.89829E11

1.14035E+00 3.58922E02 5.35510E13

1.55597E+00 4.63059E01 9.48054E14

2.23436E+00 4.28673E01 6.78291E15

RS0 1 2 3 4 5 6 7 8

1.37891E01 0.0 0.0 0.0 0.0 0.0 0.0 0.0

6.87608E02 3.21988E01 0.0 0.0 0.0 0.0 0.0 0.0

1.48734E02 1.13782E01 7.39229E01 0.0 0.0 0.0 0.0 0.0

4.69897E05 4.30321E04 1.04130E01 8.62981E01 0.0 0.0 0.0 0.0

5.89707E07 5.46672E06 1.29351E03 1.51216E01 6.99757E01 9.20879E04 0.0 0.0

3.88884E08 3.66086E07 8.66221E05 9.89375E03 1.91378E01 5.69769E01 3.55272E03 7.34779E21

1.04247E08 1.78737E07 4.26381E05 4.87002E03 8.94810E02 4.01804E01 8.98591E01 2.98050E02

8.18126E13 2.38659E08 1.33641E05 1.53773E03 2.39499E02 8.55021E02 3.70568E01 1.94260E+00

7.17290E02 0.0 0.0 0.0 0.0 0.0 0.0 0.0

2.60752E02 1.28550E01 0.0 0.0 0.0 0.0 0.0 0.0

2.76721E03 4.89640E02 3.52157E01 0.0 0.0 0.0 0.0 0.0

2.31520E06 2.12485E05 4.52107E02 4.45630E01 0.0 0.0 0.0 0.0

1.37735E07 6.34473E07 8.60573E05 6.91123E02 3.85933E01 4.54630E04 0.0 0.0

2.28260E08 1.05540E07 4.03423E06 1.73518E03 1.00105E01 3.02211E01 1.73559E03 4.23830E21

8.61157E09 8.86090E08 2.44805E06 5.17653E04 2.52677E02 1.66938E01 4.44180E01 7.65897E03

0.0 1.90539E08 1.72578E06 8.97310E05 1.98820E03 5.05147E03 4.27903E02 4.12693E01

Mixture 2: Homogenized control rod cell 2.18664E01 4.94066E01

1.13054E+00

1.38038E+00

1.39800E+00

1.46332E+00

1.75989E+00

2.94956E+00

1.87519E02

7.65941E05

9.72551E07

6.40495E08

1.70885E08

8.49733E13

RS1 1 2 3 4 5 6 7 8

Rt

RS0 1

1.17023E01

8.22991E02

(continued on next page)

128

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

Appendix A (continued) Group 1

2

3

4

5

6

7

8

0.0 0.0 0.0 0.0 0.0 0.0 0.0

3.10563E01 0.0 0.0 0.0 0.0 0.0 0.0

1.82703E01 9.30530E01 0.0 0.0 0.0 0.0 0.0

7.36525E04 1.96973E01 1.07758E+00 0.0 0.0 0.0 0.0

9.35920E06 2.47161E03 2.71658E01 9.02115E01 8.02797E04 0.0 0.0

6.26757E07 1.65514E04 1.79445E02 3.06060E01 6.88672E01 2.06824E03 0.0

3.06105E07 8.14719E05 8.83287E03 1.45799E01 6.26596E01 1.06566E+00 2.69021E02

4.15634E08 2.55580E05 2.78903E03 3.89547E02 1.36410E01 6.71040E01 2.86496E+00

6.69920E02 0.0 0.0 0.0 0.0 0.0 0.0 0.0

4.39680E02 1.61885E01 0.0 0.0 0.0 0.0 0.0 0.0

4.47493E03 8.71617E02 5.69173E01 0.0 0.0 0.0 0.0 0.0

3.69697E06 3.64473E05 8.88003E02 6.94260E01 0.0 0.0 0.0 0.0

2.27941E07 1.07947E06 1.66234E04 1.30045E01 6.08833E01 6.20153E04 0.0 0.0

3.76710E08 1.79490E07 7.59540E06 3.24058E03 1.65452E01 4.64410E01 1.36955E03 0.0

1.41237E08 1.50777E07 4.54100E06 9.64357E04 4.12750E02 2.71343E01 6.01700E01 6.04540E03

0.0 3.31300E08 3.19411E06 1.64333E04 3.22663E03 8.08650E03 9.75053E02 5.87183E01

Mixture 3: Homogenized instrumentation channel cell 7.52986E02 1.61295E01 3.46806E01

4.12540E01

4.17651E01

4.36579E01

5.21671E01

8.62488E01

2 3 4 5 6 7 8

RS1 1 2 3 4 5 6 7 8

Rt

RS0 1 2 3 4 5 6 7 8

4.28928E02 0.0 0.0 0.0 0.0 0.0 0.0 0.0

2.63531E02 1.08749E01 0.0 0.0 0.0 0.0 0.0 0.0

5.87680E03 5.22869E02 2.89903E01 0.0 0.0 0.0 0.0 0.0

2.23649E05 2.09600E04 5.59026E02 3.26007E01 0.0 0.0 0.0 0.0

2.83488E07 2.66343E06 7.00605E04 7.73596E02 2.76193E01 2.85263E04 0.0 0.0

1.86967E08 1.78361E07 4.69171E05 5.10279E03 8.73648E02 2.15107E01 6.71924E04 0.0

5.08077E09 8.71285E08 2.30942E05 2.51177E03 4.15139E02 1.79138E01 3.22860E01 7.84192E03

7.88324E13 1.18727E08 7.24466E06 7.93109E04 1.10918E02 3.88627E02 1.91934E01 8.37666E01

2.28024E02 0.0 0.0 0.0 0.0 0.0 0.0 0.0

1.28928E02 5.21333E02 0.0 0.0 0.0 0.0 0.0 0.0

1.31752E03 2.47856E02 1.62533E01 0.0 0.0 0.0 0.0 0.0

1.07367E06 1.03755E05 2.51533E02 1.97604E01 0.0 0.0 0.0 0.0

6.59027E08 3.06708E07 4.71277E05 3.69470E02 1.73556E01 1.76790E04 0.0 0.0

1.09206E08 5.09933E08 2.15318E06 9.21563E04 4.70453E02 1.32623E01 3.85813E04 0.0

4.19050E09 4.28537E08 1.28726E06 2.74246E04 1.17523E02 7.71420E02 1.72220E01 1.72488E03

0.0 9.46307E09 9.05393E07 4.67323E05 9.18763E04 2.30379E03 2.74303E02 1.67969E01

RS1 1 2 3 4 5 6 7 8

Appendix B The following tables show the seven energy group sets of cross sections used for the materials of the example reported in Section 4.3. These sets have been employed for all the simulations with the BERM-SP3, TORT and MCNP codes. Group 1 Mixture 1: UO2 fuel cell Rt 1.70136E01 mRf 1.15899E02 v 5.87910E01

2

3

4

5

6

7

3.61474E01 1.14368E03 4.11760E01

5.25057E01 8.99074E03 3.39060E04

5.62665E01 2.56528E02 1.17610E07

4.83029E01 2.48901E02 0.0

7.55205E01 1.13474E01 0.0

1.43551E+00 2.73041E01 0.0

7.24379E02 3.02966E01 0.0 0.0 0.0 0.0 0.0

3.11386E04 5.61153E02 4.04391E01 0.0 0.0 0.0 0.0

1.58891E06 2.64805E04 9.66645E02 2.95566E01 1.02171E04 0.0 0.0

2.24899E08 2.03899E05 7.20100E03 1.79448E01 2.14763E01 1.67014E03 0.0

0.0 3.16396E06 1.12016E03 2.70394E02 2.22572E01 4.46809E01 6.15231E02

0.0 4.44099E07 2.13264E04 5.15018E03 2.59218E02 2.37875E01 1.21107E+00

3.54278E01 1.35980E03 4.11760E01

5.25488E01 9.25379E03 3.39060E04

5.69165E01 3.70984E02 1.17610E07

5.44776E01 1.72822E02 0.0

8.89776E01 3.42779E01 0.0

1.48444E+00 3.48608E01 0.0

7.18064E02 2.95451E01 0.0

3.10928E04 5.65137E02 4.04183E01

1.58998E06 2.66912E04 9.67086E02

2.25093E08 2.05522E05 7.20787E03

0.0 3.18915E06 1.12123E03

0.0 4.47635E07 2.13467E04

RS0 1 2 3 4 5 6 7

9.24944E02 0.0 0.0 0.0 0.0 0.0 0.0

Mixture 2: 4.3% MOX fuel cell 1.70149E01 1.25076E02 v 5.87910E01

Rt mRf

RS0 1 2 3

9.29259E02 0.0 0.0

129

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

Appendix B (continued) Group 1

2

3

4

5

6

7

0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0

2.97663E01 1.20606E04 0.0 0.0

1.79440E01 2.14816E01 1.98592E03 0.0

2.70415E02 2.22672E01 4.31720E01 6.16650E02

5.15058E03 2.60116E02 2.39212E01 1.20673E+00

3.23389E01 1.60984E03 4.11760E01

5.26446E01 1.34469E02 3.39060E04

5.81178E01 5.31310E02 1.17610E07

5.69713E01 2.54795E02 0.0

9.52859E01 4.52771E01 0.0

1.54939E+00 4.63444E01 0.0

7.15447E02 2.62740E01 0.0 0.0 0.0 0.0 0.0

3.12827E04 5.85465E02 4.02384E01 0.0 0.0 0.0 0.0

1.59976E06 2.77594E04 9.70209E02 2.99763E01 1.28502E04 0.0 0.0

2.26486E08 2.13747E05 7.23723E03 1.79645E01 2.16407E01 2.07503E03 0.0

0.0 3.31678E06 1.12580E03 2.70815E02 2.22799E01 4.26593E01 6.24118E02

0.0 4.65549E07 2.14337E04 5.15819E03 2.60721E02 2.40363E01 1.21055E+00

5.36609E01 5.58226E03 4.11760E01

5.65015E01 1.88797E02 3.39060E04

5.91046E01 6.35094E02 1.17610E07

5.84657E01 3.06555E02 0.0

9.85357E01 5.09760E01 0.0

1.58338E+00 5.26924E01 0.0

7.56682E02 4.84889E01 0.0 0.0 0.0 0.0 0.0

2.98907E04 4.60003E02 4.38739E01 0.0 0.0 0.0 0.0

1.52105E06 2.11391E04 9.45829E02 3.03445E01 1.33529E04 0.0 0.0

2.15205E08 1.62770E05 7.03954E03 1.79476E01 2.17701E01 2.10308E03 0.0

0.0 2.52574E06 1.09504E03 2.70569E02 2.22764E01 4.23511E01 6.26395E02

0.0 3.54518E07 2.08482E04 5.15352E03 2.60924E02 2.40829E01 1.21075E+00

3.42046E01

4.13241E01

4.10227E01

4.99498E01

8.69863E01

1.81436E+00

8.18419E02 2.56620E01 0.0 0.0 0.0 0.0 0.0

4.68575E04 8.49334E02 2.51569E01 0.0 0.0 0.0 0.0

2.42505E06 4.06529E04 1.48095E01 8.44139E02 5.20843E05 0.0 0.0

3.43390E08 3.13009E05 1.11818E02 2.74364E01 1.17324E01 1.48697E03 0.0

0.0 4.85826E06 1.73939E03 4.19586E02 3.37482E01 4.86048E01 8.54379E02

0.0 6.83095E07 3.31159E04 7.99187E03 4.02201E02 3.70457E01 1.69940E+00

3.42587E01

4.13746E01

4.10552E01

4.98110E01

8.63057E01

1.84119E+00

5.67573E02 0.0 0.0 0.0 0.0 0.0 0.0

8.18658E02 2.57074E01 0.0 0.0 0.0 0.0 0.0

4.68681E04 8.50192E02 2.51832E01 0.0 0.0 0.0 0.0

2.42560E06 4.06937E04 1.48317E01 8.43943E02 5.19224E05 0.0 0.0

3.43467E08 3.13323E05 1.11984E02 2.74660E01 1.16867E01 1.47608E03 0.0

0.0 4.86314E06 1.74200E03 4.20025E02 3.36663E01 4.82228E01 8.65241E02

0.0 6.83784E07 3.31651E04 8.00029E03 4.01234E02 3.67601E01 1.72457E+00

Mixture 7: Reflector 1.59206E01

4.12970E01

5.90310E01

5.84350E01

7.18000E01

1.25445E+00

2.65038E+00

4.44777E02 0.0 0.0 0.0 0.0 0.0 0.0

1.13400E01 2.82334E01 0.0 0.0 0.0 0.0 0.0

7.23470E04 1.29940E01 3.45256E01 0.0 0.0 0.0 0.0

3.74990E06 6.23400E04 2.24570E01 9.10284E02 7.14370E05 0.0 0.0

5.31840E08 4.80020E05 1.69990E02 4.15510E01 1.39138E01 2.21570E03 0.0

0.0 7.44860E06 2.64430E03 6.37320E02 5.11820E01 6.99913E01 1.32440E01

0.0 1.04550E06 5.03440E04 1.21390E02 6.12290E02 5.37320E01 2.48070E+00

Mixture 8: Control rod 1.92421E01

4.51705E01

7.61144E01

8.06881E01

8.29067E01

1.18711E+00

2.18305E+00

7.35858E02 3.91228E01

3.62769E04 5.53566E02

1.65984E06 2.63682E04

2.24954E08 2.03035E05

0.0 3.15056E06

0.0 4.42218E07

4 5 6 7

Mixture 3: 7.0% MOX fuel cell 1.68686E01 1.32500E02 v 5.87910E01

Rt mRf

RS0 1 2 3 4 5 6 7

9.15269E02 0.0 0.0 0.0 0.0 0.0 0.0

Mixture 4: 8.7% MOX fuel cell 1.94075E01 1.78396E02 v 5.87910E01

Rt mRf

RS0 1 2 3 4 5 6 7

1.11144E01 0.0 0.0 0.0 0.0 0.0 0.0

Mixture 5: Fission chamber channel 1.39587E01

Rt

RS0 1 2 3 4 5 6 7

5.67269E02 0.0 0.0 0.0 0.0 0.0 0.0

Mixture 6: Control rod channel 1.39642E01

Rt

RS0 1 2 3 4 5 6 7

Rt

RS0 1 2 3 4 5 6 7

Rt

RS0 1 2

1.17232E01 0.0

130

V. Giusti, B. Montagnini / Annals of Nuclear Energy 42 (2012) 119–130

Appendix B (continued) Group

3 4 5 6 7

1

2

3

4

5

6

7

0.0 0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0

6.08728E01 0.0 0.0 0.0 0.0

9.54025E02 3.67842E01 6.80436E05 0.0 0.0

7.19011E03 1.76592E01 1.78780E01 1.52821E03 0.0

1.11847E03 2.69569E02 2.18687E01 4.12872E01 5.80557E02

2.12941E04 5.13447E03 2.58982E02 2.30014E01 1.42930E+00

References Baudron, A.M., Lautard, J.J., 2007. MINOS: A simplified PN solver for core calculation. Nucl. Sci. Eng. 155, 250–263. Beckert, C., Grundmann, U., 2008. Development and verification of a nodal approach for solving the multigroup SP3 equations. Ann. Nucl. Energy 35, 75–86. Bell, G.I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand Reinhold Co. Brantley, P.S., Larsen, E.W., 2000. The simplified P3 approximation. Nucl. Sci. Eng. 134, 1–21. Brown, F., Barnett, N., 2006. A Tutorial on Using MCNP for 1-Group Transport Calculations. LA-UR-07-4594. Los Alamos National Laboratory. Chiba, G., 2011. Application of the hierarchical domain decomposition boundary element method to the simplified p3 equation. Ann. Nucl. Energy 38, 1033– 1038. Ciolini, R., Coppa, G.G.M., Montagnini, B., Ravetto, P., 2002. Simplified PN and AN methods in neutron transport. Prog. Nucl. Energy 40, 237–264, Corrigenda, 2003. Ibidem 42, 121. Coppa, G., Ravetto, P., 1982. An approximate method to study the one-velocity neutron integral transport equation. Ann. Nucl. Energy 9, 169–174. Cossa, G., Giusti, V., Montagnini, B., 2010. A boundary element–response matrix method for criticality diffusion problems in xyz geometry. Ann. Nucl. Energy 37, 953–973. Davison, B., Sykes, J.B., 1957. Neutron Transport Theory. Oxford/Clarendon Press. Fletcher, J., 1983. Solution of the multigroup neutron transport equation using spherical harmonics. Nucl. Sci. Eng. 84, 33–46. Gamino, R.G., 1989. Simplified PL nodal transport applied to two-dimensional deep penetration problems. Trans. Am. Nucl. Soc. 59, 149. Gamino, R.G., 1991. Three dimensional nodal transport using the simplified PL method. In: Proceedings of ANS Topical Meetings. Advances in Mathematics, Computations, and Reactor Physics, April 29–May 2, 1991, Pittsburgh, PA, USA. American Nuclear Society, La Grange Park, IL, USA, p. 3-1. Gehin, J.C., Primm III, R.T., 1997. American Nuclear Society. Light Water Reactor Mixed Oxide Benchmark I. Computational Physics and Engineering Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA. Gelbard, E.M., 1960. Applications of Spherical Harmonics Method to Reactor Problems. Technical Report WAPD-BT-20. Bettis Atomic Power Laboratory. Gelbard, E.M., 1961. Simplified Spherical Harmonics Equations and their Use in Shielding Problems. Technical Report WAPD-T-1182. Bettis Atomic Power Laboratory. Gelbard, E.M., 1962. Applications of Simplified Spherical Harmonics Equations in Spherical Geometry. Technical Report WAPD-TM-294. Bettis Atomic Power Laboratory. Giusti, V., Montagnini, B., Coppa, G., Dulla, S., 2010. Solution of the one-velocity 2D and 3D source and criticality problems by the boundary element–response matrix (BERM) method in the A2–SP3. Nuovo Cim. della Soc. Ital. Fisica C 33, 95– 101.

Hébert, A., 2010. Mixed-dual implementations of the simplified PN method. Ann. Nucl. Energy 37, 498–511. Holloway, J.P. (Ed.), 2010. Transport Theory Statistical Physics, vol. 39. Taylor & Francis. Larsen, E.W., Morel, J.E., McGhee, J.M., 1996. Asymptotic derivation of the multigroup P1 and simplified PN equations with anisotropic scattering. Nucl. Sci. Eng. 123, 328–342. Lautard, J.J., Schneider, D., Baudron, A.M., 1999. Mixed dual methods for neutronic reactor core calculations in the CRONOS system. In: Proceedings of International Conference on Mathematics and Computation, Reactor Physics and Environmental Analysis in Nuclear Applications, September 1999, Madrid, Spain. Senda Editorial, S.A., Isla de Saipán, 47, 28035 Madrid, Spain. pp. 814– 826. Lemanska, M., 1981. On the simplified PN method in the 2D diffusion code EXTERMINATOR. Atomkernenergie 37, 173–175. Maiani, M., Montagnini, B., 1999. A boundary element–response matrix method for the multigroup neutron diffusion equations. Ann. Nucl. Energy 26, 1341–1369. Maiani, M., Montagnini, B., 2004. A Galerkin approach to the boundary element– response matrix method for the multigroup neutron diffusion equations. Ann. Nucl. Energy 31, 1447–1475. Morel, J.E., McGhee, J.M., Larsen, E.W., 1996. A three-dimensional time-dependent unstructured tetrahedral-mesh SPN method. Nucl. Sci. Eng. 123, 319–327. ORNL, 2009. SCALE: A Modular Code System for Performing Standardized Computer Analyses for Licensing Evaluations. ORNL/TM-2005/39, version 6, vols. I-III, January 2009. Oak Ridge National Laboratory. Oak Ridge, TN, USA (Available from Radiation Safety Information Computational Center at Oak Ridge National Laboratory as CCC-750). Rhoades, W.A., Simpson, D.B., 1997. The TORT Three-Dimensional Discrete Ordinates Neutron/Photon Transport Code. ORNL/TM-13221. Oak Ridge National Laboratory. Oak Ridge, TN, USA. Smith, K.S., 1986. Multidimensional nodal transport using the simplified PL method. Trans. Am. Nucl. Soc. 52, 427–430. Smith, M.A., Lewis, E.E., Na, B.C., 2003. Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation – A 2-D/3-D Mox Fuel Assembly Benchmark (C5G7 MOX Benchmark). NEA/NSC/DOC(2003)16. OECD Nuclear Energy Agency. Smith, M.A., Lewis, E.E., Na, B.C., 2005. Benchmark on Deterministic Transport Calculations Without Spatial Homogenisation. NEA/NSC/DOC(2005)16. OECD Nuclear Energy Agency. Tada, K., Yamamoto, A., Yamane, Y., Kosaka, S., Hirano, G., 2010. Validation of neutron current formulations for the response matrix method based on the SP3 theory. Ann. Nucl. Energy 37, 22–27. X-5 Monte Carlo Team, 2003. MCNP-A General Monte Carlo N-Particle Transport Code, version 5. LA-UR-03-1987. Los Alamos National Laboratory.