Uncertainty quantification in neutron transport with generalized polynomial chaos using the method of characteristics

Annals of Nuclear Energy 45 (2012) 14–28

Contents lists available at SciVerse ScienceDirect

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

Uncertainty quantification in neutron transport with generalized polynomial chaos using the method of characteristics D.A.F. Ayres ⇑, M.D. Eaton, A.W. Hagues, M.M.R. Williams Applied Modelling and Computational Group, Department of Earth Science and Engineering, Imperial College of Science, Technology and Medicine, Prince Consort Road, London SW7 2AZ, UK

a r t i c l e

i n f o

Article history: Received 23 November 2011 Received in revised form 7 February 2012 Accepted 10 February 2012 Available online 22 March 2012 Keywords: Method of characteristics Polynomial chaos Stochastic collocation Stochastic eigenvalue

a b s t r a c t The polynomial chaos expansion has been used to solve the mono-energetic stochastic neutron transport equation with the spatial and angular components discretised using the step characteristics method. Uncertainties were assumed to arise purely from the material cross sections and a novel method for treating uncertainties in discrete, uncorrelated, material regions has been proposed. The method is illustrated by numerical and Monte Carlo simulation of the mean, variance and probability density of the scalar flux for the fixed source Reed cell problem and a critical benchmark in one dimension. For the case of the critical benchmark we compare the results from the Newton–Krylov root finding method to that of the stochastic collocation method. We find that there is no benefit in the extra computation of using the Newton–Krylov method. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Deterministic transport methods compute the solution of the neutron transport equation using input parameters which are completely specified, i.e., material cross sections and any extraneous sources are modelled with zero uncertainty. In reality, the parameters such as material cross sections, initial conditions and boundary conditions are not completely known. This uncertainty can arise from random variations in material properties, such as densities or isotopic abundances, or from heterogeneities in the materials themselves. Uncertainties in material parameters, initial conditions and boundary conditions, which then translate into uncertainties in the solution itself, can be treated using stochastic characterisations of problem inputs and unknowns. Unlike the deterministic transport equation, solutions to the stochastic transport equation are random fields and it is important to be able to study the statistical moments and probability density function (PDF) of the output quantities. The stochastic transport equation is a stochastic partial integro-differential equation (SPDE). Moments and PDFs of an SPDE may be found using a variety of methods including sampling based, perturbation and spectral expansion methods. All of these methods are used to simulate or discretise a stochastic space and are independent of the discretisation of other system variables. This work will focus on using the spectral expansion method or polynomial chaos expansion (PCE) to represent the stochastic ⇑ Corresponding author. E-mail address: [email protected] (D.A.F. Ayres). 0306-4549/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2012.02.008

component of the random neutron transport problem. The PCE method has been used extensively in areas such as structural mechanics (Schuller, 2001), heat transfer (Hien and Kleiber, 1997), fluid flow (Knio and Le Matre, 2006) and also, more recently, neutronics. Applying the PCE method to fixed source neutronics problems (Williams and Eaton, 2010a; Williams, 2007, 2006) reduces the stochastic differential to a set of coupled linear differential equations. However, for multiplying systems, the PCE formulation results in a nonlinear system of equations which is NP times larger than the original system and has NP constraints where NP is the number of terms in the PCE. There are several methods available to solve this particular problem (Zhu and Wu, 1991; Ghanem and Ghosh, 2007; Rahman, 2006; Verhoosel et al., 2006) with the PCE method and also a pseudo time dependant approach which puts the problem in a linear form (Williams, 2010a,b). The neutron transport equation is purely advective and the characteristics of the associated partial differential equation are straight lines. Along these straight lines or characteristics the partial differential reduces to a total differential. Using a constant approximation for the source variation within a characteristic line segment (Le Tellier and Hebert, 2006), the angular flux can be integrated along the characteristic to an arbitrary position within the problem domain (Askew, 1972; Yousry, 1992). The characteristic lines are chosen along directions defined by a discrete ordinate quadrature set (Carlson and Bell, 1958) where angular integrals are now represented by a sum over directions with appropriate weights (Lewis, 1993). In this paper we will characterise the uncertainty in the neutron flux and criticality constant for the characteristic form of the

15

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

neutron transport equation using the PCE method. The method of characteristics or MoC has been chosen due to its accuracy and geometrical flexibility (Askew, 1972) as a well as the fact that it has never been used for uncertainty quantification in conjunction with polynomial chaos. Based on the method reported by Williams and Eaton (Williams and Eaton, 2010a; Williams and Eaton, 2010b), we give the general theory of the mono energetic stochastic method of characteristics. We then consider two types of problem: (1) Fixed source problems in slab geometry for multiple regions with differing material and statistical properties. (2) Multiplying slab problems. For the fixed source problems we present an analytic solution for a homogeneous, purely absorbing slab and study the rate of convergence of the PCE method for different relative variances. We also show the numerical results of the PCE and Monte Carlo simulation of the Reed cell problem with uncorrelated material regions. For the multiplying system we calculate the Keffective eigenvalue and flux eigenvector in the entire parameter range for a critical bench mark problem (Sood et al., 2003) and compare with previously published results for a stochastic discrete ordinates formulation (Fichtl, 2009). 2. General theory of the Stochastic Method of Characteristics (SMoC)

where the isotropic source term is

qðr þ s0 Xm ; nÞ ¼ kðnÞmRf ðr þ s0 Xm ; nÞ

d/ðr þ s0 X; X; nÞ þ Rt ðr þ s0 X; nÞ/ðr þ s0 X; X; nÞ 0 ds ¼ Q ðr þ s0 X; X; nÞ

ð1Þ

The energy dependence has been ignored since that will simplify some of the mathematical expressions that will be derived. The n(n1, n2    nN) are independent random variables which will describe the data uncertainty. For isotropic transport, fission, scattering and external sources Eq. (1) becomes:

d/ðr þ s0 X; X; nÞ þ Rt ðr þ s0 X; nÞ/ðr þ s0 X; X; nÞ 0 ds qðr þ s0 X; nÞ ¼ 4p

/ðr þ s0 Xm ; Xm ; nÞwm

m¼1

þ RS0 ðr þ s0 Xm ; nÞ

Nl X

/ðr þ s0 Xm ; Xm ; nÞwm

m¼1

þ qext ðr þ s0 Xm ; nÞ

NP X

/ðr þ s0 X; X; nÞ ¼

/i ðr þ s0 X; XÞUi ðnÞ

kðnÞ ¼

NP X

kk Uk ðnÞ

ð7Þ

k¼0

The gPC basis functions Ui(n) obey the following orthogonality relation

hUj ðnÞUi ðnÞi ¼

Z

dnpðnÞUj ðnÞUi ðnÞ ¼ dji N2i

ð8Þ

where p(n) is an appropriate probability weight function of the random variables n and N2i is the normalisation constant. The number of basis functions in the expansion, NP, is given by the following (Ghanem and Spanos, 1991):

NP þ 1 ¼

ðp þ mÞ! p! m!

ð9Þ

In the above equation p is the maximum order of the PC polynomials and m is the number of stochastic dimensions. Fig. 1 shows the computational region of interest which is defined in cartesian geometry according to Dx, Dy and Dz. Over the computational region of interest we assume the material parameters Rt ðr þ s0 Xm ; nÞ; RS0 ðr þ s0 Xm ; nÞ and Rf(r + s0 Xm, n) are spatially uniform. Inserting Eqs. (6) and (7) into (4) yields

ð2Þ

Z qðr þ s0 X; nÞ ¼ kðnÞmRf ðr þ s0 X; nÞ /ðr þ s0 X; X; nÞdX Z 4p /ðr þ s0 X; X; nÞdX þ RS0 ðr þ s0 X; nÞ 4p

ð3Þ

The term k(n) = 1/Keffective(n) where Keffective is the multiplication factor. qext(r + s0 X, n) is an independent source term which itself may contain uncertainties. Discretising the angular variable using a discrete ordinate quadrature set yields the following expression

d/ðr þ s0 Xm ; Xm ; nÞ þ Rt ðr þ s0 Xm ; nÞ/ðr þ s0 Xm ; Xm ; nÞ 0 ds qðr þ s0 Xm ; nÞ ¼ 4p

ð6Þ

i¼0

The total neutron source which is the sum of the fission, scattering and external sources can be written explicitly as:

þ qext ðr þ s0 X; nÞ

ð5Þ

In the above expression the subscript m represents each discrete direction and Nl is the number of discrete directions. Let us assume that the flux /(r + s0 Xm, Xm, n) and the uncertain eigenvalue k(n) can be expanded in a PCE of the form (Ghanem and Ghosh, 2007)

C

This section derives the multi-dimensional form of the stochastic method of characteristics and discretises the stochastic domain using the generalised polynomial chaos (gPC) method. The characteristic form of the transport equation is derived from the integrodifferential form of the neutron transport equation (Askew, 1972; Halsall, 1980). The stochastic forward characteristic form of the neutron transport equation which is written without energy dependence is:

Nl X

ð4Þ Fig. 1. Cartesian space-angle coordinate system in three dimensions.

16

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

ARi ¼ ki Ri

NP d X /i ðr þ s0 Xm ; Xm ÞUi ðnÞ 0 ds i¼0

þ Rt ðnÞ

NP X

/i ðr þ s0 Xm ; Xm ÞUi ðnÞ ¼

i¼0

0

qðr þ s Xm ; nÞ 4p

ð10Þ

Now, projecting onto the basis Uj(n) and taking the mathematical expectation yields the following equations:

^ j ðr þ s Xm ; Xm Þ d/ þ 0 ds 0

NP X

^ i ð r þ s 0 Xm ; X m Þ Aji /

i¼0

^j ðr þ s Xm Þ q 4p

dwðr þ s0 Xm ; Xm Þ Cðr þ s0 Xm Þ þ Kwðr þ s0 Xm ; Xm Þ ¼ 0 4p ds

8 j ¼ 0; 1; . . . ; NP



¼m

kk

i¼0

k¼0



Nl X

Nl

C jki

X

^ i ðr þ s0 Xm ; Xm Þwm þ /

m¼1

^ i ðr þ s0 Xm ; Xm Þwm þ /

m¼1

NP X

Bji

i¼0 NP X hUj ðnÞqext i ðnÞi

ð12Þ

i¼0

^ j ðr þ sXm ; where we have defined the modified flux component / Xm Þ ¼ N j /j ðr þ sXm ; Xm Þ and the components

hUj ðnÞRt ðnÞUi ðnÞi Aji ¼ Ni Nj hUj ðnÞRS0 ðnÞUi ðnÞi Bji ¼ Ni Nj hUj ðnÞRf ðnÞUk ðnÞUi ðnÞi C jki ¼ Ni Nj Z

ð14Þ ð15Þ

dnpðnÞRðnÞUj ðnÞUi ðnÞ

ð16Þ

R

p1 ðR0 ÞdR0

ð17Þ

Ci ¼

1 V

As u varies from zero to unity, the value of R varies from its minimum value to its maximum value and is weighted according to p(R0 ). But we must also relate n to the uniform random variable u by setting n

pðn0 Þdn0

ð18Þ

nmin

From Eqs. (17) and (18) we have the inverses R = f(u) and n = g(u) and by setting du = p(n)dn we may transform Eq. (16) to

hUj RUi i ¼

Z

1

duf ðuÞUj ðgðuÞÞUi ðgðuÞÞ

Z

ð1  expðKi ðSm  Sin ÞÞÞ

ð24Þ

Cðr þ s0 Xm ÞdV 0

wi ðSout ; Xm Þ ¼ wi ðSin ; Xm Þ expðKi DSm Þ þ

ð20Þ

To solve (20) it is useful to diagonalise the matrix A as follows

Ci 4pKi

ð1  expðKi DSm ÞÞÞ

ð26Þ

where DSm = Sout  Sin. The average value of the modified source C is given by the following expression

C¼m

NP X

kk L1 C k L

Nl X

 m Þwm þ L1 BL wðX

m¼1

k¼0

Nl X

 m Þwm þ Cext wðX

ð27Þ

m¼1

 m is the average angular flux over the In the above expression w characteristic line segment of length DSm. The average can be calcu becomes lated from (24). The final expression for w

 ¼ wi ðSin ; Xm Þ  wi ðSout ; Xm Þ þ Ci w K i DSm 4pKi

ð28Þ

The angular flux can be recovered from the modified flux component using the expression

/i ðr; XÞ ¼

We may now express Eq. (11) in matrix form

ð25Þ

V

ð19Þ

0

0 ^ ð r þ s 0 Xm ; Xm Þ ^ d/ ^ ðr þ s0 Xm ; Xm Þ ¼ qðr þ s Xm Þ þ A/ 0 4 p ds

Ci 4 p Ki

In the above expression V is the volume of the characteristic line segment. By setting Sm = Sout and letting the origin of the characteristic r0 = (0, 0, 0), the angular flux exiting the characteristic line segment can be calculated in terms of the angular flux entering the segment,

Rmin

Z

ð23Þ

where the average source term Ci is defined as

In the above equation we note that R and n may belong to different probability spaces. In this instance we are required to map both R and n to the same probability space. It is most convenient to map them to a uniform probability space u 2 [0, 1], and write

u¼

wðr þ s0 Xm ; Xm Þwm þ Cext ðnÞ

wi ðSm ; Xm Þ ¼ wi ðSin ; Xm Þ expðKi ðSm  Sin ÞÞ

ð13Þ

C

Z

wðr þ s0 Xm ; Xm Þwm þ L1 BL

m¼1

To proceed with the integration of Fig. 22 the modified isotropic neutron source C(r + s0 Xm) is approximated by its average value C over a characteristic within the computational region of interest. This is known as the flat source approximation or step characteristics method. Other methods of source approximation are discussed in (Le Tellier and Hebert, 2006). Integrating along a characteristic from r = r0 to r = r0 + sX gives the following expression

þ

The angular brackets in Eqs. (13)–(15) are defined by

hUj ðnÞRðnÞUi ðnÞi ¼

Nl X

Nl X

m¼1

hU ðnÞqj ðr þ s Xm ; nÞi ^j ¼ j q Nj NP X

kk L1 C k L

k¼0

0

u¼

NP X

ð11Þ

with

NP X

ð22Þ

^ ¼ Lw and C ¼ L1 q ^. The modified source term where we have set / is expressed as

Cðr þ s0 Xm Þ ¼ m

0

¼

where L is an (NP + 1)  (NP + 1) matrix with column vectors given by each eigenvector Ri. Inserting A as given in (21) into (20) and re-arranging gives

NP 1 X Lij wj ðr; XÞ ¼ ðDwðr; XÞÞi Ni j¼0

ð29Þ

and the angular flux can be mapped to the modified space using

wi ðr; XÞ ¼

NP X

 1  L1 ij N j /j ðr; XÞ ¼ D /ðr; XÞ i

ð30Þ

j¼0

A ¼ LKL1

ð21Þ

where K is a (NP + 1)  (NP + 1) diagonal matrix with elements ki. L is defined by the eigenvalue problem

where D is the modified-moment operator. The mean and the variance of the angular flux in Eq. (6) can be calculated using the following relations:

17

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

h/ðr; XÞi ¼ /0 2 /

r ¼

NP X

/2i ðr; XÞh

2 ii

U

ð31Þ

P X

ð32Þ

i¼0

i¼1

In the case of a multiplying source the retrieval of the statistical moments of the multiplication constant Keffective from Eq. (7) requires a little more care. The expansion coefficients for Keffective in the form

K effective ðnÞ ¼

NP X

Gi;0;0 ki 6 6 i¼0 6 6 .. 6 . 6 6N 4 PP Gi;0;NP ki

P X

ðNj Þ2 /j ðr; XÞ ¼ ð1Þ

ð1Þ

ð2Þ

ð1Þ

ð2Þ

/i ðr; XÞhUj ðnÞUi ðnÞi

ð38Þ

i¼0

P 1 X

ð1Þ

ð34Þ

ð39Þ

ð2Þ

Ni

ð2Þ

ð40Þ

Lik wk ðr; XÞ

k¼0

P X

ð1Þ

Ljk wk ðr; XÞ ¼

P X

ð2Þ

wk ðr; XÞ

k¼0

P X

1

ð1Þ ð2Þ i¼0 N j N i

ð1Þ

ð2Þ

Lik hUj ðnÞUi ðnÞi

ð41Þ

or, in matrix form

ð35Þ

The statistical moments of Keffective can be found by inserting the expansion coefficients from Eq. (33) into Eqs. (31) and (32). 3. Boundary and interface conditions The treatment of boundary conditions and interface conditions for the stochastic neutron transport equation are outlined in the paper by (Williams and Eaton, 2010a,) and are reproduced for convenience below. We know that the angular flux for each realization, n, is continuous at the interface between two adjacent dissimilar media, so we may write at the boundary point r

/ð1Þ ðr; X; nÞ ¼ /ð2Þ ðr; X; nÞ

P 1 X

ð2Þ

/i ðr; XÞ ¼

ð1Þ

Lik wk ðr; XÞ

k¼0

k¼0

i¼0

hUi Uj Uk i hUk Uk i

ð1Þ

Ni

Inserting (39) and (40) into (38) and rearranging we find

In the above expression Gijk is the Galerkin multiplication tensor and is given by

Gijk ¼

ð37Þ

Projecting onto the gPC basis, Uj(n), and taking the mathematical expectation we find

/i ðr; XÞ ¼

3 NP P  Gi;NP ;0 ki 72 3 2 3 1 i¼0 7 K0 76 7 6.7 . 7 .. .. . 7¼6.7 76 . . 74 . 5 4 . 5 7 NP 5 K NP 0 P  Gi;NP ;NP ki

i¼0

ð2Þ

/i ðr; XÞUi ðnÞ

i¼0

ð33Þ

can be recovered by solving the following system of linear equations (Maitre and Knio, 2010): NP P

P X

Using the relation between / and w from Eq. (29)

K l Ul ðnÞ

l¼0

2

ð1Þ

/i ðr; XÞUi ðnÞ ¼

Lð1Þ wð1Þ ¼ Z ð1;2Þ wð2Þ

ð2Þ

If we decide to retain the same PCE for each region then U Ui and Z(1,2) = L(2). For a bare surface the boundary condition reduces to

Lð2Þ wð2Þ ¼ 0 :

nX<0

ð43Þ

4. Slab with non-multiplying source In this section we will only be considering non-multiplying systems, this removes the host of difficulties associated with evaluating the PC moments in the expansion of the stochastic eigenvalue. Considering only a single dimension, x, ignoring the fission term and integrating over the azimuthal direction Eqs. (26)–(28) reduce to the following

wi ðSout ; lm Þ ¼ wi ðSin ; lm Þ expðKi DSm Þ

ð36Þ

Fig. 2 illustrates the interface condition between regions with dissimilar material properties. The angular flux exiting a region penetrates into the adjacent region along the path of motion. Inserting the PCEs for each side

ð42Þ ð1Þ i =

þ

Ci ¼

Ci ð1  expðKi DSm ÞÞ 8i ¼ 0;    ; NP Ki

Nl NP   X 1X  j ðl Þwm þ Cext w L1 BL i m ij 2 j¼0 m¼1

Fig. 2. Angular flux across multiple regions reproduced from (Cacuci, 2010).

ð44Þ

ð45Þ

18

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

 i ¼ wi ðSin ; lm Þ  wi ðSout ; lm Þ þ Ci w K i DSm Ki

ð46Þ

where matrices A and B are given by Eqs. (13) and (14) and l = cos(h) where h is the angle between the characteristic line and the coordinate axis. The solution procedure for Eqs. (44)–(46) is given in the following Algorithm: Algorithm 1. Algorithm for the solution of the fixed source stochastic MoC

hUj ðnÞRðnr ÞUi ðnÞi ¼

Input: Initial angular flux and source term, ðw ; CÞ Output: Converged angular flux, ~ wn while ðk~ wðaÞ  ~ wða1Þ k > tolÞ do  ~ Call Average_Angular_Fluxðw; wðaÞ Þ Call Source_TermðC; ~ wðaÞ Þ

5. Slab with multiplying source

wi ðSout ; lm Þ ¼ wi ðSin ; lm Þ expðKi DSm Þ þ

Ci ¼ m An initial guess for the PCE angular flux w and source term, C  from must be given in order to calculate the average angular flux, w Eq. (46). Next, the average source term, C is calculated using Eq. (45) and then the angular flux is computed with a marching scheme using Eq. (44). If there is scattering in the problem then each direction will be coupled and the loop will continue until the angular flux has converged to the specified tolerance. Strictly speaking Rt(n) and RS0 ðnÞ from matrices A and B should be written Rt(n1) and RS0 ðn2 Þ where n1 and n2 are independent random variables. However, to simplify the analysis we assume n1 = n2 = n. Physically, this means that the uncertain parameters are perfectly correlated. Given that the material cross sections are uncertain in the range [Rmin, Rmax] we have modelled them as follows.

1

Rmax  Rmin RðnÞ ¼ Rmin þ ðRmax  Rmin Þn

ð47Þ n 2 ½0; 1

ð48Þ

 Log-uniform Distribution

logðRmax  Rmin Þ

ð49Þ



RðnÞ ¼ Rmin

NP X NP X

Ci ð1  expðKi DSm ÞÞ Ki

Nl NP   X X  l Þwm þ wð kk L1 C k L ðL1 BLÞij m ij

j¼0 k¼0



Nl X

ð52Þ

m¼1

j¼0

 j ðl Þwm þ Cext w i m

ð53Þ

m¼1

 i ¼ wi ðSin ; lm Þ  wi ðSout ; lm Þ þ Ci w K i DSm Ki

ð54Þ

For multiplying systems where we have NP terms in our PCE for the flux and eigenvalue, Nquad discrete directions and N spatial nodes there will be (NquadN + 1)(NP + 1) unknowns. In Eq. (52) there are (NquadN)(NP + 1) equations. The other (NP + 1) equations are derived from the normalisation condition on the eigenvalue:

/T / ¼ 1

ð55Þ

The eigenfunction / in Eq. (55) can be replaced by its PC expansion to yield

 Uniform Distribution

R  Rmax n Rmin

ð51Þ

where n = n1, n2, . . . , nr, . . . , nR and R is the number of uncorrelated regions.

Call Calc_Angular_Fluxð~ wðaÞ Þ end

pðRÞ ¼

dn1 Z Z  dn2    dnR Uj ðnÞRðnr ÞUi ðnÞpðnR Þ

This section deals with the stochastic eigenvalue problem and outlines the algorithm used to solve the resulting set of non-linear equations. Again, we will only be considering problems in slab geometry. In one dimension, Eqs. (26)–(28) reduce to the following:

~0

pðRÞ ¼

Z

n 2 ½0; 1

ð50Þ

NP X NP X ðDwi ÞT ðDwj ÞUi Uj ¼ 1 i¼0

Here D is the modified-moment operator defined by Eq. (29). Taking Eq. (56), projecting onto the basis Ul and taking the mathematical expectation yields the final (NP + 1) equations NP X NP X ðDwÞTi ðDwÞj < Ui Uj Ul >¼ dl0 i¼0

For the statistical distributions chosen to model the data uncertainties the appropriate basis polynomials are the Legendre polynomials (Xiu and Karniadakis, 2002). 4.1. Uncorrelated regions In this problem we assume that all stochastic parameters, i.e. the material cross sections, are perfectly correlated within each region and their uncertainties can be represented with the same random variable nr. Different material regions are then assumed to be uncorrelated with respect to one another and a different independent random variable is used to model the uncertainties in each. The basis functions in Eq. (6) now form a partial tensor product of the appropriate one dimensional polynomials (Ghanem and Spanos, 1991). The number of basis functions NP is given by Eq. (9). The calculation of the stochastic inner product matrices, Eqs. (13)–(15), now take the form

ð56Þ

j¼0

8l ¼ 0; . . . ; N P

ð57Þ

j¼0

To solve the (NquadN)(NP + 1) system of non-linear equations given by Eqs. (52) and (57) we use the Newton–Raphson method (Fichtl, 2009; Ghanem and Ghosh, 2007). The system of equations is first written in the form

Fð~ xÞ ¼ 0

ð58Þ

where Fl are given by ðe;aþ1Þ

F l ðSout ; lm Þ ¼ Wl 

ðe;aþ1Þ

ðSout ; lm Þ  Wl

  ðeÞ ðSin ; lm Þ exp Kl DSm

  Clðe;aÞ  ðeÞ 1  exp Kl DSm ðeÞ Kl

8e 2 f1; . . . ; Ne g ð59Þ

F NP þ1;l ¼

NP X NP X ðDwi ÞT ðDwj Þ < Ui Uj Ul > dl0 i¼0

j¼0

ð60Þ

19

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

In Eq. (59) we have defined Ne as the number of characteristic line segments, e as the line segment index and a as the Newton iteration index. The solution vector ~ x is given in the following format:

~ x ¼ ½~ /0 ; ~ /i ; . . . ; ~ /NP ; k0 ; . . . ; kNP T ~ /i ¼ ½/ ; . . . ; / ;/ ;...;/ 11

1;N quad

21

2;N quad ; . . . ; /N1 ; . . . ; /N;N quad 

ð61Þ T

ð62Þ

expanding Fð~ xÞ in a Taylor series around ~ x gives

Fð~ x þ d~ xÞ ¼ Fð~ xÞ þ J  d~ x þ Oðd~ x2 Þ

ð63Þ

where J is the Jacobian matrix

@F i J ij ¼ @xj

ð64Þ

Since the zeros of Fð~ xÞ are desired, the condition Fð~ x þ d~ xÞ ¼ 0 is imposed and, neglecting higher-order terms, the equation for the update, dx, becomes

J  d~ xðaÞ ¼ Fð~ xÞ

ð65Þ

~ xðaÞ þ d~ xðaÞ xðaþ1Þ ¼ ~

ð66Þ

As can be seen from Eq. (61) the Newton–Krylov solution vector ~ x contains the unmodified flux component /. This is because the moments of the eigenvalue inverse, k, are defined in the unmodified space. As a result, the residual of the modified angular flux given by Eq. (59) must be mapped back to real space using Eq. (29) before the update vector d~ x can be computed from Eq. (65). The solution process for the flux and eigen-moments is defined in the following algorithm Algorithm 2. Newton–Krylov algorithm for the stochastic MoC eigenvalue Input: Initial flux, eigen-moments and source term, ð~ x 0 ; CÞ n ~ Output: Converged flux and eigen-moments, x while ðk~ xðaÞ  ~ xða1Þ k > tolÞ  ~ Call Average_Fluxðw; xðaÞ Þ) Call Source_TermðC; ~ xðaÞ Þ Call Calc_Residualðreal F; C; ~ xðaÞ Þ Call Calc_Jacobian (J, ~ xðaÞ , Calc_Residual) Call GMRES_Solver (J, real_F, d~ xðaÞ ) ð a þ1Þ ð a Þ ð a Þ ~ ¼~ x þ d~ x x end

/j ðr þ sX; XÞ ¼ kj ¼

hUj ðnÞ/ðr þ sX; X; nÞi

hUj ðnÞkðnÞi

ð69Þ

N2j

To calculate a good initial guess for the residual F in Eq. (52) we must also calculate initial values for the modified average angular  This can be calculated in the following way flux w.

  j ¼ hUj ðnÞ/ðnÞi / N2j

ð70Þ

The average angular flux is then mapped to the modified space. The integrals over stochastic space, hi, are evaluated using Gaussian quadrature. For additional information regarding the application of the SCM to radiation transport problems see (Fichtl and Prinja, 2011). 6. Results We will illustrate the above techniques with three examples: 1. An analytical, fixed source, homogeneous, purely absorbing slab to demonstrate the underlying theory and study the convergence of the PCE method. 2. The Reed cell problem to demonstrate the interface condition between multiple regions with different statistical properties and uncorrelated uncertainties. 3. Keffective eigenvalue problem for a homogeneous slab. The first two of the above are fixed source problems where the external source is assumed to have zero uncertainty. Uncertainties in all of the above problems have been modelled as either uniform or log-uniform between upper and lower limits. The appropriate polynomial basis (Xiu and Karniadakis, 2002) for these distributions are the Legendre polynomials. 6.1. Problem 1: One region The one region problem is a 3 cm, purely absorbing, homogeneous slab with vacuum boundary conditions. In the absence of the scattering and fission terms Eq. (2) can be solved exactly to yield:

   qext Ra ðnÞx 1  exp  2Ra ðnÞ l    qext Ra ðnÞðx  LÞ wðx; l; nÞ ¼ 1  exp  2Ra ðnÞ l

wðx; l; nÞ ¼  is the In the above algorithm C is the modified source term, w modified average angular flux, real_F is the residual in real space, J is the Jacobian in real space and a is the Newton iteration index. To calculate the modified average flux using Eq. (54), the angular flux components of the solution vector ~ x must first be mapped to the modified space using Eq. (30). The Jacobian is recalculated for each Newton iteration but its sparsity structure remains unchanged. For a single homogeneous region, assuming the modified average flux is constant for each Newton iteration, it is possible to find an analytical solution for the Jacobian. In general, there will not exist a straight forward analytical solution and, as such, we have chosen to approximate the Jacobian using a finite difference relation as follows:

J ij 

F i ðxj þ hÞ  F i ðxj Þ h

ð68Þ

N2j

ð67Þ

To speed up the convergence of the Newton–Krylov algorithm initial values of the PC expansion coefficients for the flux and eigenvalue moments can be calculated using the Stochastic Collocation Method (SCM) (Maitre and Knio, 2010). Eqs. (6) and (7) are premultiplied by the polynomial basis functions Uj(n) and with the orthogonality relation given in Eq. (8) this yields

l>0 l<0

ð71Þ ð72Þ

with L being the width of the slab. The scalar flux can be recovered from the angular flux by integrating over angle

/ðx; nÞ ¼

Z

0

wðx; l; nÞdl þ

q

Ra ðnÞ

1

wðx; l; nÞdl

0

1 ext

¼

Z

½1  0:5fE2 ðRa ðnÞðL  xÞÞ þ E2 ðRa ðnÞxÞg

ð73Þ

where En is the exponential integral. The mean and the variance in the scalar flux as a function of position are given by the following

 /ðxÞ ¼

Z

1

pðnÞ/ðx; nÞdn

ð74Þ

n¼0

r2/ ðxÞ ¼

Z

1

 2 ðxÞ pðnÞ/2 ðx; nÞdn  /

ð75Þ

n¼0

The integrals over the stochastic dimension in the above expressions are evaluated using the adaptive Gauss–Konrod quadrature routine, qag, from the library QUADPACK. Ra(n) is assumed to vary

20

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

Table 1 Cross section data for the one region slab. Problem

Domain

Distribution

Ra

r2 =Ra

qext

A B

[0:3.0] [0:3.0]

Uniform Uniform

1.5 1.5

1/8 1/800

1.0 1.0

linearly across the stochastic space and the functional form is given by Eq. (48). The upper and lower bounds, Rmax/min can be calculated using the following relation:

pffiffiffi

Rmax=min ¼ R  3r

ð76Þ

The mean and relative standard deviation in the scalar flux have been evaluated numerically for problem A in Table 1 at 5 spatial locations with an angular discretisation of S500. The results for the mean and relative standard deviation are given in Tables 2 and 3 respectively along with the analytic values calculated from Eqs. (74) and (75). In reactor analysis a quantity of more practical interest is the reaction rate. To calculate the differential reaction rate we substitute the expression /(x) = R(x)/Ra into Eq. (2) to give

dRðx; nÞ qext Ra ðnÞ þ Ra ðnÞRðx; nÞ ¼ ds 2

ð77Þ

Inserting the PC expansions for the reaction rate and following the same procedure as before, Eq. (77) becomes the following matrix equation:

^ q dRðxÞ þ ARðxÞ ¼ 2 ds

ð78Þ

^ is now given where A is given by Eq. (13) and the modified source q by

^j ¼ q

hUðnÞRa ðnÞqext ðnÞi Nj

ð79Þ

^j can be Using the form of Ra(n) given in Eq. (48) the components q calculated exactly and are given by

qext ðRmax þ Rmin Þ 2 ext q ^1 ¼ pffiffiffi ðRmax  Rmin Þ q 6 3 ^j>1 ¼ 0 q

^0 ¼ q

ð80Þ ð81Þ ð82Þ

It is important to note that when calculating the source moments above that the external source is assumed to be determinisP tic, i.e. qext ðnÞ ¼ j¼0 Uj ðnÞqext j dj0 . Numerical calculations for the mean and relative standard deviation in the reaction rate have been performed on the same spatial mesh as the scalar flux problem above. Results for the mean and relative standard deviation are given in Tables 4 and 5 along with the analytic solution. The analytic solution to the mean and relative standard deviation in the scalar flux and reaction rate have also been evaluated at 100 discrete spatial points and compared with the deterministic value computed with the mean cross section. The results are plotted in Figs. 3 and 4. Once the coefficients in the PC expansion for the scalar flux and reaction rate have been computed, the PDF at a particular position can be found using Eq. (6). The probability space of the random variable/vector, n, is sampled and a realisation of the stochastic flux or reaction rate is computed from the PC expansion. 106 realisations are computed in this way and the output space is subdivided into 100 bins to give the PDF. To explain the phenomena exhibited in Figs. 3–5 we must first consider how the scalar flux varies as a function of absorption. It is evident that as Ra ? 1, / ? 0 and as Ra ? 0 the flux is limited by

some maximum value determined by the leakage from the system. Since this problem is purely absorbing, the variation of flux as a function of absorption is exponential i.e. we expect a higher probability of lower flux and vice versa. This can be seen in Fig. 5 where there is a strong anti-correlation between absorption rate and scalar flux. Since the PDF of the scalar flux favours lower scalar flux values, the deterministic calculation with mean cross sections in Fig. 3 is lower than the mean scalar flux from the stochastic case. Conversely, the deterministic reaction rate is larger since higher values are favoured. To explain the behaviour of the RSDs in Fig. 4 we consider multiple realisations of the scalar flux and reaction rate for different absorption cross sections, Ra. Considering how the scalar flux changes with Ra as above, we see that a spread in Ra values causes a spread in scalar flux values. But, when the scalar flux is multiplied by Ra, we see a clustering in magnitude of the reaction rate realisations. Hence, we observe a larger RSD in the scalar flux compared to the reaction rate and also an anticorrelation between the two RSDs. Fig. 6 shows the convergence rates for the mean and variance in the scalar flux and reaction rate for the different relative uncertainties in the absorption cross section as given Table 1. The L2 error in the difference between the analytic and PC solution, as given in Tables 2–5, has been plotted against various PC orders. For both the scalar flux and reaction rate larger relative uncertainties have a slower rate of convergence compared to small uncertainties. This is due to the fact that for small perturbations from the mean the output parameter space can be well approximated using a linear variation, as is the justification for perturbation theory, so we find that for low order PCE the error is small and converges quickly. The L2 error levels off at around 1011 in both cases which is due to the angular approximation. 6.2. Problem 2: The Reed cell This problem is designed to demonstrate the accuracy of the stochastic method of characteristics over more complicated geometries. Fig. 7 presents the diagram of the problem domain. The Reed cell problem is a very demanding test case even for deterministic transport codes. It consists of five regions lying within a reactor lattice cell located at the edge of a bare core. In this problem we have assumed that all macroscopic cross sections within a region are perfectly correlated since they are a property of the same material. As a result, the uncertainties of the cross sections within a region can be modelled with a single random variable. Dissimilar regions will have different variations in number density or be composed of different materials producing uncertainties which have been assumed as uncorrelated with all other regions. This situation is of particular interest in reactor physics since lattice cell codes use/ produce nuclear cross section data which is spatially constant over discrete regions. To model uncorrelated uncertainties in each material region requires an independent random variable in each, producing a total of five stochastic dimensions. Region 1 is a fuel region and is situated next to a perfect reflector. It is 2 cm in length and contains an optically thick, purely absorbing media with a uniformly distributed uncertainty in the transport cross section and a deterministic source. The next region is a 1 cm can. This is another purely absorbing region but with its data uncertainty described with a log-uniform distribution. Adjacent to the fuel can is a 2 cm void and a 3 cm moderator. The first 1 cm of the moderator has a deterministic source with scattering and transport cross sections described by a log-uniform distribution. The final 2 cm of moderator has a bare vacuum boundary with its transport and scattering cross sections described by a uniform distribution. The upper and lower bounds for the cross sections in each region with their associated statistical distributions are summarised in Table 6.

21

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28 Table 2 Mean scalar flux. x

PCE 1

PCE 2

PCE 3

PCE 4

PCE 5

Analytic

0.00 0.75 1.50 2.25 3.00

0.36163364 0.66456832 0.69427892 0.66456832 0.36163364

0.36356365 0.66679969 0.69662689 0.66679969 0.36356365

0.36367683 0.66691781 0.69674640 0.66691781 0.36367683

0.36368329 0.66692433 0.69675293 0.66692433 0.36368329

0.3636836737 0.6669247117 0.6967533127 0.6669247117 0.3636836737

0.3636836979 0.6669247359 0.6967533369 0.6669247359 0.3636836979

Table 3 Relative standard deviation in the scalar flux. x

PCE 1

PCE 2

PCE 3

PCE 4

PCE 5

Analytic

0.00 0.75 1.50 2.25 3.00

28.50216989 25.28679613 26.38431394 25.28679613 28.50216989

31.27630514 27.20292667 28.38599653 27.20292667 31.27630514

31.56155449 27.36557888 28.54920543 27.36557888 31.56155449

31.58546796 27.37753934 28.56079894 27.37753934 31.58546796

31.58726504 27.37836749 28.56158487 27.37836749 31.58726504

31.58740278 27.37842771 28.56164128 27.37842771 31.58740278

Table 4 Mean reaction rate. x

PCE 1

PCE 2

PCE 3

PCE 4

PCE 5

Analytic

0.00 0.75 1.50 2.25 3.00

0.4978183644 0.9240855461 0.9620988075 0.9240855461 0.4978183644

0.4975770107 0.9234943845 0.9613640886 0.9234943845 0.4975770107

0.4975574415 0.9234666753 0.9613332022 0.9234666753 0.4975574415

0.4975563386 0.9234654469 0.9613319414 0.9234654469 0.4975563386

0.4975562835 0.9234653905 0.9613318849 0.9234653905 0.4975562835

0.4975562806 0.9234653876 0.9613318819 0.9234653876 0.4975562806

Table 5 Relative standard deviation in the reaction rate. x

PCE 1

PCE 2

PCE 3

PCE 4

PCE 5

Analytic

0.00 0.75 1.50 2.25 3.00

0.3980985905 3.8626802012 2.6879252139 3.8626802012 0.3980985905

0.5961924326 4.2707136229 3.1783961684 4.2707136229 0.5961924326

0.6422446069 4.3058506025 3.2291990537 4.3058506025 0.6422446069

0.6489119478 4.3082387868 3.2327553946 4.3082387868 0.6489119478

0.6495994602 4.3083815833 3.2329618041 4.3083815833 0.6495994602

0.6496598538 4.3083901382 3.2329737030 4.3083901382 0.6496598538

Fig. 3. Mean scalar flux and reaction rate for the one region problem.

The following calculations were performed on a spatial mesh of 160 elements with an S4 angular discretisation. Fig. 8 shows the mean value of the scalar flux over the domain for PC expansions 1–5, Fig. 9 shows the relative standard deviation in the scalar flux for PC expansions 1–5 and Fig. 10 shows the PDFs for the scalar flux at x = 2.0 cm, x = 3.0 cm, x = 5.0 cm and x = 6.0 cm for PC 5. The data

Fig. 4. Relative standard deviation in the mean scalar flux and reaction rate for the one region problem.

for the mean and relative standard deviation are tabulated in Tables 7 and 8. All results agree to at least two decimal places with the Monte Carlo solution which was computed from 106 realisations of the deterministic system. For the number of realisations computed, the Monte Carlo solution has not fully converged. For

22

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

Fig. 6. Convergence of the mean and the variance in the scalar flux for the one region problem with different relative uncertainties.

and can regions we can see the deterministic flux is lower than the average flux. The reason for this is the same as in Problem 1; these regions are dominated by absorption and the mean of the input space corresponds to a flux in the output space which is lower than the average. The effect of absorption on the output stochastic space can be seen in Fig. 10a and b. In the other regions where scattering is prevalent the effect of absorption is lost and the output space more closely resembles the PDF of the input space Fig. 10c and e. This explains why there is little difference between the deterministic flux and mean flux in the void and moderator regions in Fig. 8. The moderator in region 4 of Fig. 7 is a slight exception since it contains a deterministic source, this explains the maximum in the PDF of Fig. 10d. A similar maximum is not observed in the PDF for region 1 since it is dominated by absorption. There is also a large variance in the void region, which itself has no uncertainty. It is the leakage current from the neighbouring region which causes this uncertainty. Fig. 5. Probability density functions of the scalar flux and absorption reaction rate at various spatial positions for the one region problem.

6.3. Problem 3: multiplying systems 106 realisations, the Monte Carlo takes approximately 3 h and 4 GB of memory to compute and increasing this by and order of magnitude becomes infeasible. In Fig. 8 the deterministic calculation of the problem using the mean value of the material cross sections is also shown. In the fuel

In this problem we demonstrate the use of the Newton–Krylov algorithm to compute the PC moments for the random angular flux and uncertain eigenvalue Keffective. The problem is the same as in the work by Fichtl (Fichtl, 2009) which will serve as a benchmark.

23

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

Fig. 7. Diagram of the one dimensional Reed cell problem.

Table 6 Cross section data for the Reed cell problem. Region

Domain

Distribution

Rt

1. 2. 3. 4. 5.

[0.0:2.0) [2.0:3.0) [3.0:5.0) [5.0:6.0) [6.0:8.0]

Uniform Log-uniform N/A Log-uniform Uniform

25.0 2.5 0.0 0.5 0.5

Fuel Can Void Moderator Moderator

2

Rt

min

75.0 7.5 0.0 1.5 1.5

0.0 0.0 0.0 0.45 0.45

qext

max

0.0 0.0 0.0 1.35 1.35

50.0 0.0 0.0 1.0 0.0

PC-1

PC-1 PC-2

Relative Standard Deviation (%)

PC-4

1.4

PC-2

70

PC-3

1.6

PC-5 MC 106

Mean (x)

RS

min

80

Deterministic

1.8

1.2 1 0.8 0.6 0.4 0.2 0

RS

max

PC-3 PC-4

60

PC-5 6

MC 10

50 40 30 20 10

0

1

2

3

4

5

6

7

8

x (cm) Fig. 8. Mean scalar flux for the uncorrelated Reed cell problem.

The problem is a critical bare homogeneous slab reactor consisting of uranium dioxide with uncertain macroscopic cross sections described by a single random variable. The problem is taken from an analytical benchmark test set originally intended for deterministic code verification (Sood petffiffiffiffial., 2003). The uncertainties in the cross 2 sections are defined by hRri ¼ 15 with the mean values given in Table 9. In this work only uniformly distributed uncertainties in the macroscopic cross sections have been considered. For results of other statistical distributions and uncertainties for this problem see the work by Fichtl (Fichtl, 2009). The following numerical calculations were performed on a spatial mesh of 100 elements with an angular discretisation of S8. Before solving the system, the convergence rates of the SCM method were investigated as a function of the convergence tolerance of the power iteration and source iteration (PI and SI) of the deterministic code. This showed how the error in the SCM residual converges as a function of the number of collocation points M. Together, these measures provide insight into choosing the optimal parameters to initialise the Newton method. Fig. 11 shows the L2 error of the

0 0

1

2

3

4

5

6

7

8

x (cm) Fig. 9. Relative standard deviation in the scalar flux for the uncorrelated Reed cell problem.

SCM residual, F, given by Eqs. (59) and (60) as a function of the tolerance of the PI and SI tolerance. We can see that the L2 error is strongly dependant upon the tolerance of the SI and PI. Since the Newton method requires a good initial guess, it is important to choose a tolerance where the SCM has completely converged. It was decided that the number of collocation points should be M P (NP + 1). From Fig. 11, we can see the optimal PI and SI tolerance is 6109. We also observe that the converged value of the L2 error decreases as NP increases. This indicates that the SCM and Galerkin projection methods will give equivalent answers for higher order PCE moments. In Figs. 12 and 13 the number of Newton iterations and time taken to converge the L2 error in the residual to a tolerance of 108 is shown for PI and SI convergence tolerances of 106 and 1012 respectively. We note that as the number of collocation points increases, the solutions converge with higher order PC expansions needing fewer iterations. From Fig. 13a and b we can see that decreasing the PI and SI tolerance from 106 to 1012 reduces the

24

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

Fig. 10. Probability density functions at various spatial positions for the uncorrelated Reed cell problem.

number of Newton iterations for PC 7 and PC 9 by an order of magnitude. This is due to the fact that the SCM initialisation for these PC expansions has converged to a value less than or equal to the convergence tolerance of the Newton solver. We can see

from Figs. 13 and 12a that M = (NP + 1) is sufficient to minimise the number of Newton iterations for a given PI and SI tolerance but to minimise the number of Newton iterations the SCM solution must be fully converged.

25

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28 Table 7 Mean scalar flux data for the Reed cell problem. x

PCE 1

PCE 2

PCE 3

PCE 4

PCE 5

Monte Carlo

1.00 2.50 4.00 5.50 7.00

1.09090909 0.05603684 1.13706491 1.82204943 0.70802565

1.09803922 0.05657704 1.13707813 1.82190388 0.70797973

1.09857035 0.05658477 1.13707249 1.82189147 0.70797211

1.09860924 0.05658547 1.13707248 1.82189126 0.70797203

1.09861207 0.05658552 1.13707247 1.82189124 0.70797201

1.09915922 0.05656677 1.13707651 1.82189607 0.70792341

Table 8 Relative standard deviation data for the Reed cell problem. x

PCE 1

PCE 2

PCE 3

PCE 4

PCE 5

Monte Carlo

1.00 2.50 4.00 5.50 7.00

28.86751346 67.80646627 9.74619467 1.49810921 12.56742385

31.94382825 75.48073410 9.73828541 1.69922967 12.59296184

32.31686570 75.85020608 9.74067999 1.70948267 12.59346571

32.35553159 75.86475749 9.74081041 1.71008246 12.59336690

32.35917741 75.86597083 9.74082033 1.71011696 12.59336682

32.34966760 75.89763410 9.74412017 1.70973508 12.59516588

Table 9 Cross section data the bare reactor problem. 0.54628 0.464338 0.0928676

Rt Rs mRf

-3

10

P=1 P=3 P=5 P=7 P=9

10-4

Residual Error

-5

10

-6

10

10-7

10-8

-9

10

10-13 10-12 10-11 10-10 10-9

10-8

10-7

10-6

10-5

10-4

10-3

10-2

10-1

Tolerance Fig. 11. L2 error in the residual of the SCM solution as a function of convergence tolerance of the PI and SI. M = (NP + 1).

From the results above, we chose the PI and SI tolerance to be 109 and set the number of collocation points M = (NP + 1). The mean in the scalar flux is given in Fig. 14, the variance in the scalar flux is given in Fig. 15 and the PDF for the stochastic Keffective in Fig. 16. Data at three spatial locations are given for the mean and RSD in Tables 10 and 11 for PCE 1, 3, 5, 7, 9 and a Monte Carlo simulation of 106 realisations. The mean tabular data agree to 5 decimal places and variance agrees to at least two with the MC simulation. The convergence of the L2 error for the entire flux vector is given in Fig. 17a. The PDF of the Keffective eigenvalue for various PCE and the MC is given in Fig. 16, the convergence of its mean and variance with respect to the MC are given in Fig. 17b. The deterministic calculation using the mean cross section values is also shown in Fig. 14 with the deterministic calculation slightly higher than the stochastic average. The decrease in the average flux compared with the deterministic calculation can be explained by considering the effect of the uncertain fission source on the output stochastic space. In this problem, the reproduction factor g = Rf/

Fig. 12. (a) The number of Newton iterations required to converge the residual to a tolerance of 1  108 as a function of M, the number of collocation points. (b) Execution time for the examples in Fig. 12a. PI and SI were converged to a tolerance of 106.

Ra = 1.1333 suggests that the output statistics will be dominated by fission. In other words, the PDF of the scalar flux will favour larger values given a uniformly distributed input space. From these

26

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28 18

PC-1 PC-3 PC-5 PC-7 PC-9 MC 106

16 14

RSDφ (%)

12 10 8 6 4 2 0 5

0

10

15

20

x (cm) Fig. 15. Variance in the scalar flux for the stochastic eigenvalue problem.

14

12

PC-1 PC-3 PC-7 PC-15 6 MC 10

P (k)

10

8

6

4

2 0.88

0.9

0.92

0.94

0.96

0.98

1

1.02

1.04

1.06

1.08

k Fig. 13. (a) The number of Newton iterations required to converge the residual to a tolerance of 1  108 as a function of M, the number of collocation points. (b) Execution time for the examples in Fig. 13a. PI and SI were converged to a tolerance of 1012.

0.14

0.1

Mean

centre due to the decreased leakage from the system. The PDFs of Keffective calculated for various PC orders are given in Fig. 16 and converge to the Monte Carlo solution as the PC order increases. For a PC expansion of order one (PC-1) the basis functions, which for Legendre polynomials are constant plus linear, are unable to capture the higher order variation in the uncertainty of Keffective. This explains the linear variation of PC-1 in Fig. 16. Higher order polynomial expansions are able to more accurately represent the non-linear variation in the response of Keffective to input parameter uncertainty. To look more closely at the convergence of the mean and variance in the scalar flux and Keffective the L2 error has been computed using the following:

Deterministic PC-1 PC-3 PC-5 PC-7 PC-9 MC 106

0.12

0.08

0.06

0.04

/ ¼ k/p ðx; nÞ  /e ðx; nÞk2 K ¼ kK p ðx; nÞ  K e ðx; nÞk2

0.02

0

Fig. 16. Probability density functions of the stochastic eigenvalue for various PC expansions.

0

5

10

15

20

x (cm) Fig. 14. Mean scalar flux for the K eigenvalue problem.

arguments it is apparent that the deterministic calculation using the average cross sections will produce a value greater than the stochastic average, as in Fig. 14. This effect is pronounced in the

In the above expressions p stands for the order of PC approximation and e stands for the exact or reference solution. The L2 error of the flux and Keffective for the PC and SCM (with a Gauss quadrature order equal to NP + 1) have been plotted as a function of NP in Fig. 17. The reference solution is a Monte Carlo simulation using 106 realisations. The mean and variance in the scalar flux and the stochastic Keffective in Fig. 17a and b converge to within a fixed value of the reference solution for both the SCM and gPC methods. Results for the

27

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28 Table 10 Mean scalar flux data for the Keffective eigenvalue problem. x

PCE 1

PCE 3

PCE 5

PCE 7

PCE 9

Monte Carlo

0.000000 5.185533 10.371065

0.019213 0.102334 0.133710

0.019223 0.102332 0.133709

0.019223 0.102332 0.133709

0.019223 0.102332 0.133709

0.019223 0.102332 0.133709

0.019225 0.102332 0.133708

Table 11 Relative standard deviation data for the Keffective eigenvalue problem. x

PCE 1

PCE 3

PCE 5

PCE 7

PCE 9

Monte Carlo

0.000000 5.185533 10.371065

15.923718 0.278977 0.928074

16.462755 0.286228 0.964104

16.463268 0.286218 0.964125

16.463269 0.286218 0.964125

16.463269 0.286218 0.964125

16.461028 0.286187 0.964059

Table 12 Probability of super critical reactor (Keffective > 1) for various PC expansions with cross pffiffiffiffi 2 section uncertainty hRri ¼ 15. SN ðSCMÞ100 results taken from (Fichtl, 2009). The 32 superscripts and subscripts stand for the spatial and angular discretisations respectively. PC order

SMOC ðNKÞ100 8

MOC ðSCMÞ100 32

MOC ðSCMÞ400 32

SN ðSCMÞ100 32

1 3 7 15 MC

0.429767 0.493562 0.496034 0.496271 0.496035

0.429767 0.493561 0.495938 0.496163

0.439722 0.505039 0.495859 0.496108

0.435864 0.501563 0.499945 0.499946 0.500370

Table 13 Mean and variance in the Keffective eigenvalue. PC order

Mean

Variance

1 2 2 4 5 6 7 8 9 MC

0.989998670207306 0.989757907180035 0.989755353069144 0.989755422422864 0.989755426319174 0.989755428569063 0.989755436041727 0.989755436038524 0.989755436037616 0.989743963004378

0.001867873785303 0.002073489493207 0.002082750786664 0.002083017496260 0.002083019452422 0.002083024143270 0.002083024543236 0.002083024535575 0.002083024534086 0.002083003766899

In the work by Fichtl the problem is solved using a discrete ordinates angular discretisation with a linear discontinuous finite element (LDFE) spatial discretisation. Results for the probability of the reactor being super critical P(Keffective > 1), are tabulated for various probabilities and uncertainties. Data for an S32 angular discretisapffiffiffiffi 2 tion for uniform uncertainties defined by hRri ¼ 15 are tabulated in

Fig. 17. Convergence of the L2 error for the mean flux, variance in the flux and Keffective for both PC and SCM. The reference value was taken to be a Monte Carlo simulation with 106 realizations.

mean and variance in Keffective are tabulate in Table 13. The large error that remains suggests that the MC with 106 realizations has not fully converged. Computing 106 realizations of the deterministic code takes approximately 30 h and to do more would be too time consuming. In both cases the convergence rate with respect to increasing polynomial order is exponential.

Table 12. Data for the SMoC solved using the Newton–Krylov method are also tabulated for an angular discretisation of S8. For comparison between the MoC and LDFE SN methods data for the MoC using the SCM with spatial discretisations of 100 and 400 elements are also given in Table 12. The Monte Carlo results for the 32 ordinate problems have not been calculate due to the time need to compute the required number of realisations. The results for SMoC ðNKÞ100 and MoC ðSCMÞ100 8 32 are very similar but show a slight difference when compared with SN ðSCMÞ100 32 . Here we have used a subscript to represent the angular discretisation and a superscript to indicate the number of spatial elements. For low 100 order PCE the MoC ðSCMÞ400 32 results are closer to the SN ðSCMÞ32 data but as the PC order increases the results converge to the MoC ðSCMÞ100 32 data. This suggests that the difference between the MoC and SN implementations is due more to the model approximation rather than the spatial discretisation.

28

D.A.F. Ayres et al. / Annals of Nuclear Energy 45 (2012) 14–28

7. Summary and conclusions In this paper we applied the pseudo-multigroup approach to uncertainty quantification outlined in (Williams and Eaton, 2010a; Williams and Eaton, 2010b) to perform uncertainty modelling for mono energetic fixed source and eigenvalue problems using the MoC method. A method for modelling problems consisting of multiple homogeneous regions of different materials with uncertainties that are independent of one another has also been proposed. The method uses an independent random variable in each material region to represent the uncertainty. The random angular flux is then expanded in a multivariate PCE with the number of stochastic dimensions equal to the number of different materials. Finally, we investigate the stochastic eigenvalue problem which arises when a fission source is included. The eigenvalue can be treated as a random variable and hence, be expanded using a PCE. The resulting set of non-linear equations were solved using a Newton–Krylov root finding method and the results of this were compared with the non-intrusive stochastic collocation method (SCM). We have used the work by Fichtl (Fichtl, 2009) as a benchmark and also to investigate the accuracy and convergence properties of applying the PCE method to stochastic neutron transport problems solved by the method of characteristics. We also found that the accuracy of the SCM for the MoC is strongly dependant on the convergence tolerance of the source and power iterations of the deterministic system. For the initialisation of the Newton–Krylov method, the optimal quadrature order for the SCM to minimise the number of non-linear iterations is M = (NP + 1) assuming the PI and SI tolerance is sufficient to completely converge the SCM solution. We also found that the SCM solution best approximates the gPC solution for higher order polynomial expansions meaning lower order approximations require more non-linear iterations to converge. This will not present a significant problem since low order polynomial expansions are usually not sufficient to appropriately model the uncertainties in the output quantities.

Acknowledgements The authors thank Professor Anil Prinja from the University of New Mexico. They would also like to thank Dr. Brian Turland and Mr. Richard Hiles from Serco for many useful discussions on this work. The authors would also like to show appreciation to the European Commission for their support through the EU FP7 Project NURISP (Nuclear Reactor Integrated Simulation Project). Dr. M.D. Eaton would like to acknowledge the support of The Royal Academy of Engineering and EPSRC under their Research Fellowship Scheme. Mr. Dan Ayres would like to acknowledge the support of EPSRC under their industrial doctorate programme as well as industrial support from NNL (National Nuclear Laboratory). Mr. Andrew

Hagues would like to thank the support of the EPSRC through their KNOO programme (EPSRC Grant EP/C549465/1). References Askew, J.R., 1972. A Characteristics Formulation of the Neutron Transport Equation in Complicated Geometries. Tech. Rep. AEEW-M 1108, United Kingdom Atomic Energy Establishment. Cacuci, D.G. (Ed.), 2010. Handbook of Nuclear Engineering, vol. 2. Springer. Carlson, B.G., Bell, G.I., 1958. Solution to the neutron transport equation by the sn method. In: 2nd U.N. Conf. on Peaceful uses of Atomic Energy, Geneva. p. 2386. Fichtl, E.D., 2009. Stochastic Methods for Uncertainty Quantification in Radiation Transport. Ph.D. Thesis, University of New Mexico. Fichtl, E.D., Prinja, A.K., 2011. The stochastic collocation method for radiation transport in random media. Journal of Quantitative Spectroscopy and Radiative Transfer 112 (4), 646–659. Ghanem, R., Ghosh, D., 2007. Efficient characterization of the random eigenvalue problem in a polynomial chaos decomposition. International Journal for Numerical Methods in Engineering 72 (4), 486–504. Ghanem, R., Spanos, P., 1991. Stochastic Finite Elements: A Spectral Approach. Springer, New York. Halsall, M.J., 1980. Cactus, a Characteristics Solution to the Neutron Transport Equations in Complicated Geometries. Tech. Rep., United Kingdom Atomic Energy Establishment. Hien, T.D., Kleiber, M., 1997. Stochastic finite element modelling in linear transient heat transfer. Computer Methods in Applied Mechanics and Engineering 144 (1–2), 111–124. Knio, O., Le Matre, O., 2006. Uncertainty propagation in CFD using polynomial chaos decomposition. Fluid Dynamics Research 38 (9), 616–640. Le Tellier, R., Hebert, A., 2006. On the integration scheme along a trajectory for the characteristics method. Annals of Nuclear Energy 33 (14–15), 1260–1269. Lewis, E.E., 1993. Computational Methods of Neutron Transport. American Nuclear Society. Maitre, O., Knio, O.M., 2010. Spectral Methods for Uncertainty Quantification. Springer. Rahman, S., 2006. A solution of the random eigenvalue problem by a dimensional decomposition method. International Journal for Numerical Methods in Engineering 67 (9), 1318–1340. Schuller, G.I., 2001. Computational stochastic mechanics – recent advances. Computers & Structures 79 (22–25), 2225–2234. Sood, A., Forster, R.A., Kent Parsons, D., 2003. Analytical benchmark test set for criticality code verification. Progress in Nuclear Energy 42 (1), 55–106. Verhoosel, C.V., Gutirrez, M.A., Hulshoff, S.J., 2006. Iterative solution of the random eigenvalue problem with application to spectral stochastic finite element systems. International Journal for Numerical Methods in Engineering 68 (4), 401–424. Williams, M.M.R., 2006. Polynomial chaos functions and stochastic differential equations. Annals of Nuclear Energy 33 (9), 774–785. Williams, M.M.R., 2007. Polynomial chaos functions and neutron diffusion. Nuclear Science and Engineering 155, 109–118. Williams, M.M.R., 2010a. A method for solving a stochastic eigenvalue problem applied to criticality. Annals of Nuclear Energy 37 (6), 894–897. Williams, M.M.R., 2010b. A method for solving stochastic eigenvalue problems. Applied Mathematics and Computation 215 (11), 3906–3928. Williams, M.M.R., Eaton, M., 2010a. A probabilistic study of the influence of parameter uncertainty on solutions of the neutron transport equation. Progress in Nuclear Energy 52 (6), 580–588. Williams, M.M.R., Eaton, M., 2010b. A probabilistic study of the influence of parameter uncertainty on solutions of the radiative transfer equation. Journal of Quantitative Spectroscopy and Radiative Transfer 111 (5), 696–707. Xiu, D., Karniadakis, G.E., 2002. The wiener–askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing 24 (2), 619–644. Yousry, Y.A., 1992. Arbitrarily high order characteristic methods for solving the neutron transport equation. Annals of Nuclear Energy 19 (10–12), 593–606. Zhu, W.Q., Wu, W.Q., 1991. A stochastic finite element method for real eigenvalue problems. Probabilistic Engineering Mechanics 6 (3–4), 228–232.