- Email: [email protected]
A probabilistic study of the influence of parameter uncertainty on solutions of the neutron transport equation
Progress in Nuclear Energy 52 (2010) 580–588
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
A probabilistic study of the influence of parameter uncertainty on solutions of the neutron transport equation Matthew Eaton, 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 2BP, UK
a b s t r a c t Keywords: Uncertainty Neutron transport Stochastic
The influence of uncertainty in the scattering and absorption cross sections on the solution of the neutron transport equation is studied using polynomial chaos theory. The uncertainty is defined by means of uniform and log-uniform probability distributions. By expanding the neutron flux in a series of polynomial chaos functions we may reduce the stochastic transport equation to a set of coupled deterministic equations, analogous to those that arise in multi-group neutron transport theory, with the effective multi-group transfer scattering coefficients containing information about the uncertainty. This procedure enables existing transport theory computer codes to be used, with little modification, to solve the problem. Applications are made to a transmission problem and a constant source problem in a slab. In addition, we also study the rod model for which exact analytical solutions are readily available. In all cases, numerical results in the form of mean, variance and sensitivity are given which illustrate how absorption and scattering cross section uncertainty influences the solution of the transport equation. 2010 Elsevier Ltd. All rights reserved.
1. Introduction An important source of error in calculations employing the neutron transport equation arises from uncertain data values. These errors originate from lack of knowledge about the absorption and scattering cross sections. It is vital, therefore, to know how changes in these parameters affect the outcome of any calculations. For example, what are the uncertainties in the reflection and transmission coefficients of a slab in relation to uncertainties in the absorption and scattering cross sections? When the uncertainties are relatively small, these matters can be been handled by perturbation theory (Williams, 1986; Cacuci, 2003). However there is another effective method of dealing with such a problem which does not have the limitations of perturbation theory and which may be used for arbitrarily large parameter uncertainties; namely, polynomial chaos expansions. In order to proceed, therefore, we deliberately introduce uncertainty into the parameters of the transport equation, thereby converting it into a stochastic integrodifferential equation. This is done by assuming that the parameters are functions of a random variable which, in turn, allows us to prescribe a probability distribution for the values of the parameters. The polynomial chaos expansion (PCE) (Ghanem and Spanos, 1991;
* Corresponding author. E-mail address: [email protected] (M.M.R. Williams). 0149-1970/$ – see front matter 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.pnucene.2010.01.002
Emery, 2005; Wiener, 1938) is then employed to transform the stochastic equation into a set of coupled deterministic transport equations for the coefficients in the PCE. By a further transformation, these coupled equations may be cast into a form that is analogous to a multi-group neutron (or radiation) transport formulation but with probabilistic matrix elements taking the place of the group-to-group energy exchange cross sections. The only difference between the equations derived below from those of multi-group theory is the presence of coupled boundary conditions. This transformation therefore enables us to use well-established computer codes to solve the associated multi-group equations. In this paper we will give examples based upon the onedimensional transport equation for a slab. Two cases are considered: 1) a constant source in a bare slab and 2) the transmission of neutrons through a plane slab with an isotropic source incident on one face. We will solve the transport equation numerically but, because certain analytic solutions exist for the cases we have chosen, we will be able to compare the PCE method with the exact result for mean and variance and hence estimate convergence rates. We also solve the problems using the rod model, for which exact solutions are also available over the complete range of distance. Some additional calculations are carried out, which use the exact transport equation, by a modification of the computer code EVENT (DeOliveira, 1986). The general conclusion is that the PCE method shows considerable promise and can deal with a variety of probability distribution functions describing parameter uncertainty.
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
It is important to note that a number of other fields of technology have benefited from the use of PCE to study the effects of data uncertainty. For example, the uncertainty of modes of vibration in structural mechanics due to errors in elastic constants (Li and Ghanem, 1998), the flow of groundwater in the geo-sphere due to rock fractures (Ghanem, 1998) and uncertainties in fluid flow due to both geometric uncertainty and basic data (Le Maitre et al., 2001). The extension to neutron transport is an obvious extension, as we hope to convince the reader.
D
The method to be described works just as well for the energy dependent transport equation but for simplicity in this introductory paper we use the one speed equation for the angular neutron flux fðr; U; xÞ, which may be written as (Davison, 1957)
U:Vfðr; U; xÞ þ SðxÞfðr; U; xÞ ¼ Ss ðxÞ
Z
dU0 f U:U0 ; x f r; U0 ; x þ Sðr; U; xÞ
N X
fn ðr; UÞFn ðxÞ
¼ hFm ðxÞSðr; U; xÞi Z N X fn r; U0 Fm ðxÞSs ðxÞf U:U0 ; x Fn ðxÞ : þ dU0
ð1Þ
(2)
where Fn ðxÞ are the polynomial chaos functions which are to be determined by the nature of the problem and obey the orthogonality relation
hFn ðxÞFm ðxÞih
R
2 dxpðxÞFn ðxÞFm ðxÞ ¼ dnm Nm
(3)
pðxÞ being an appropriate probability weight function of the random variables x. Eq. (3) also defines the normalisation coefficient Nm. For full details of the concept of polynomial chaos functions, the reader is referred to Ghanem and Spanos (1991). However, in the simplest case of one random variable, x, the polynomials form a complete set of orthogonal functions with the weight function given by the appropriate probability distribution (Xiu and Karniadakis, 2002a,b). For more than one random variable, the sequence of polynomials can be generalised as explained in Ghanem and Spanos (1991) or Williams (2006). The expansion coefficients fn ðr; UÞ are to be determined in the manner described below. Inserting Eq. (2) into Eq. (1) we find
U:V
N X
fn ðr; UÞFn ðxÞ þ SðxÞ
n¼0
¼ S s ðxÞ
N X
fn ðr; UÞFn ðxÞ
Let us now define
dU0 f ðU:U0; xÞ
N X
fn r; U0 Fn ðxÞ þ Sðr; U; xÞ:
hFm ðxÞSðxÞFn ðxÞi Nn Nm
Fm ðxÞSs ðxÞf U:U0 ; x Fn ðxÞ Bmn U:U0 ¼ Nn Nm
(6)
(7)
Sm ðr; UÞ ¼ hFm ðxÞSðr; U; xÞi=Nm
b m ðr; UÞ ¼ Nm f ðr; UÞ: f n
(8)
We observe that the matrices A and B are symmetric. The angular brackets in Eqs. (4)–(7) are defined by
hFm SFn i ¼
Z
dxpðxÞSðxÞFm ðxÞFn ðxÞ
(9)
where pðxÞ is the weight function associated with PCE polynomials Fn . It should be noted that the evaluation of this matrix element depends on the nature of the statistics of SðxÞ and x. For example, if S obeys a prescribed probability distribution p1 ðSÞ and suppose that we have decided that the weight function associated with the PCE polynomial Fn ðxÞ is pðxÞ; then if these two distributions span different spaces the integral (9) as it stands needs interpretation. For simplicity we consider the case where there is only one random variable x, then we can define a uniformly distributed random variable u˛ð0; 1Þ, and write (Carter and Cashwell, 1975)
u ¼
Z S S0
p 1 S 0 dS 0
(10)
Thus as u varies from zero to unity, so the value of S will vary from its minimum of S0 to its maximum of S1 and be weighted according to p1 ðSÞ. But we must also refer x to the uniform variable u which we do by setting
u ¼
Z x 0 0 p x dx :
(11)
x0
From (10) we have the inverse S ¼ f(u) and from (11) the inverse x ¼ g(u). In some cases explicit expressions can be obtained for f and g, but in general this is not so. Now we transform Eq. (9) by setting du ¼ pðxÞdx, whence
hFm SFn i ¼
Z
1 0
duf ðuÞFm ðgðuÞÞFn ðgðuÞÞ:
(12)
With the matrix elements defined, we may now write Eq. (5) in the form
b m ðr; UÞ þ U:V f
n¼0
Z
ð5Þ
n¼0
n¼0
Z
fn ðr; UÞhFm ðxÞSðxÞFn ðxÞi
b m as We also define the modified flux component f
In Eq. (1) the total cross section SðxÞ ¼ Sa ðxÞ þ Ss ðxÞ, where Sa and Ss are the absorption and scattering cross sections, respectively. f ðU:U0; xÞ is the differential scattering probability per unit solid angle. Sðr; U; xÞ is an independent source term which itself may contain uncertainties. The xðx1 ; x2 ; .xN Þ are independent random variables which will describe the data uncertainty. It should be noted that we will not be considering multiplying systems in this paper. The introduction of fissile material transforms the equation into an eigenvalue problem and this raises a host of difficulties. Treatment of these difficulties is deferred to a future publication, but some discussion is given in the conclusions in Section 6. Let us assume that the random flux can be expanded in a PCE in the form (Ghanem and Spanos, 1991; Williams, 2006)
fðr; U; xÞ ¼
N X n¼0
Amn ¼
2. Theory
E
U:Vfm ðr; UÞ F2m ðxÞ þ
581
N X
b n ðr; UÞ Amn f
n¼0
(4)
n¼0
Now multiplying Eq. (4) by pðxÞFm ðxÞ and averaging over all the random variables x, leads to
¼
N Z X
b 0 dU0 Bmn U:U0 f n r; U þ Sm ðr; UÞ:
n¼0
In matrix form, Eq. (13) becomes
(13)
582
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
b þ Af b ¼ Bf b þ S: U:V f
(14)
To solve this equation, it is useful to diagonalise the matrix A as follows 1
A ¼ L LL
(15)
where L is a NN diagonal matrix with diagonal elements li and L is defined by the eigenvalue problem
AR ¼ lR
(16)
where L is an NN matrix with column vectors given by each eigenvector R of (16) corresponding to the eigenvalue li. Inserting A as given by (15) in Eq. (14) and re-arranging gives 1
1
U:Vj þ Lj ¼ L BLj þ L S
(17)
b ¼ Lj. In discrete form Eq. (17) becomes where we have set f
U:Vji ðr; UÞ þ li ji ðr; UÞ ¼
N Z X
dU0 Tij U:U0 jj r; U0 þ Qi ðr; UÞ
j¼0
(18)
It should be noted that where it is typographically more elegant, we use an over-bar to denote an average instead of C.D. 3. Boundary conditions Let us consider the conditions that must hold between two dissimilar media. We know that the angular intensity is continuous at an interface so that we may write for two adjacent media 1 and 2, at the boundary point rs,
fð1Þ ðrs ; U; xÞ ¼ fð2Þ ðrs ; U; xÞ; n:Us0
where n is a unit vector normal to the interface. Now we insert the PCE’s for each side, which in general will be different because the cross sections are different and we may wish to choose different chaos polynomials. We then find N X
0
Tij U:U and
Qj ¼
¼ L1 B U:U0 L
(19)
ij
N X
fnð2Þ ðrs ; UÞFnð2Þ ðxÞ:
(25)
n¼0
Now multiply by Fð1Þ m ðxÞ and average to get N D E X ð1Þ 2 ð1Þ ð1Þ ð2Þ fm ðrs ; UÞ ¼ fð2Þ Nm n ðrs ; UÞ Fm ðxÞFn ðxÞ
(26)
n¼0
L1 S :
Using the relation between f and j from Eq. (21), viz
(20)
j
It is important to note that if the probability distribution p1 ðSÞ, describing the distribution of S(x), is also the weight function for the orthogonality condition of the Fn(x), then the matrix T will be diagonal. The diagonal elements of T are equal to the eigenvalues of the diagonalised matrix B, analogous to Eq. (16). This means that Eq. (18) decouples and the solutions only differ because of different sources and boundary conditions. If this condition is not met, or one of the parameters has a different distribution, e.g. log-uniform, then off-diagonal matrix elements will appear and there is no decoupling. Generally, these off-diagonal elements are significantly smaller in magnitude than those of the diagonal ones. Eq. (18) is in the form of standard multi-group neutron transport equations with the stochastic matrix elements related to the groupto-group transfer cross sections. Thus it is possible to use existing multi-group computer codes to solve for the PCE coefficients. To recover the components of the flux in Eq. (2) we use Eq. (8) and the b ¼ Lj as defined above, just below Eq. (17), to get relation f
fi ðr; UÞ ¼
ð1Þ fð1Þ n ðrs ; UÞFn ðxÞ ¼
n¼0
with i ¼ 0,1,2,.N, and we have defined
(24)
fð1Þ m ðrs ; UÞ ¼
N 1 X ð1Þ
Nm
(21)
(27)
ð2Þ ð2Þ
(28)
and
fð2Þ m ðrs ; UÞ ¼
N 1 X ð2Þ
Nm
Lmj jj ðrs ; UÞ
j¼0
and inserting (27) and (28) into (26) we find ð1Þ
Nm
N X
ð1Þ ð1Þ
Lmj jj ðrs ; UÞ¼
N D X
ð2Þ Fð1Þ m Fn
n¼0
j¼0
N E 1 X
ð2Þ ð2Þ Lnj jj ðrs ; UÞ ð2Þ Nn j¼0
(29)
which on re-arrangement gives N X
ð1Þ ð1Þ
Lmj jj ðrs ; UÞ¼
j¼0
N 1X Lij jj ðr; UÞ Ni
ð1Þ ð1Þ
Lmj jj ðrs ; UÞ
j¼0
N X
jð2Þ ðrs ; UÞ j
j¼0
N X
1 ð2Þ L ð1Þ ð2Þ nj n¼0 Nm Nn
D
ð2Þ Fð1Þ m Fn
E
(30)
or in matrix form
j¼0
From Eq. (2) we may now calculate the mean and variance of the angular intensity and any higher moments of interest, viz: the mean is
hfðr; UÞihf0 ðr; UÞ ¼ f0 ðr; UÞ
(22)
and the variance
E
D
s2f ¼ ðfðr; UÞ fðr; UÞÞ2 ¼
N X
D E f2n ðr; UÞ F2n :
(23a)
For the scalar flux 2
2
sf0 ¼ ðf0 ðrÞ f0 ðrÞÞ
E
¼
N X
D E f20n ðrÞ F2n
(23b)
where
f0n ðrÞ ¼
Z
dUfn ðr; UÞ:
(23c)
(32)
for medium 2, which implies that, in this case only, each individual component of the PCE expansion satisfies the condition
jð2Þ n ðrs ; UÞ ¼ 0 : n:U < 0
n¼1
(31)
If we decide to retain the same PCE in each region as indeed may (1,2) ¼ L(2). For a bare be quite acceptable, then Fnð1Þ ¼ Fð2Þ n and Z surface the boundary condition reduces to
Lð2Þ jð2Þ ðrs ; UÞ ¼ 0 : n:U < 0
n¼1
D
Lð1Þ jð1Þ ¼ Z ð1;2Þ jð2Þ :
(33)
In Eqs. (32) and (33), n is an outward pointing unit normal to the surface. For the case where a source f(U;x), with uncertainty, is incident on the surface the condition becomes
fðrs ; U; xÞ ¼ f ðU; xÞ; n:U < 0
(34)
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
Expanding in PCE as in Eq. (2), and using the orthogonality relation, we find
1 b m ðrs ; UÞ: hFm ðxÞf ðU; xÞi ¼ f Nm
(35)
b ¼ Lj, Eq. (35) becomes Then from the relation f
583
The PCE expansion of Eqs. (37) and (38) is given by reference to Eqs. (18) and (21) where we have chosen the Legendre polynomials for the polynomial chaos expansion because the random variable x˛(0,1). Thus we write
f ðxÞ ¼
N X
f n ðxÞPn ð2x 1Þ
(44a)
n¼0
jm ðrs ; UÞ ¼
N X
L1
n¼0
1 mn Nn
hFn ðxÞf ðU; xÞi
(36)
Thus in order to use the standard multi-group codes we must change the boundary conditions as shown in Eqs. (31) and (36). For vacuum boundaries, no changes are needed.
where the weight function p(x) ¼ 1 for such polynomials (Xiu and Karniadakis, 2002a,b). After diagonalisation,
f i ðxÞ ¼
N N pffiffiffiffiffiffiffiffiffiffiffiffiffi X 1X Lij j Lij j 2i þ 1 j ðxÞ ¼ j ðxÞ: Ni j ¼ 0 j¼0
(44b)
The expansion coefficients are then given by 4. The rod model It is useful to illustrate the above theory using a transport equation for which exact solutions can be obtained. For the case of the slab, the rod model is appropriate. In this model, the neutrons are supposed to travel either in the positive or negative x-directions. Thus the transport Eq. (1), may be written in slab geometry as (setting m ¼ 1)
1 dfþ ðxÞ 1 þ SðxÞfþ ðxÞ ¼ Ss ðxÞ fþ ðxÞ þ f ðxÞ þ S dx 2 2
(37)
and
1 df ðxÞ 1 þ SðxÞf ðxÞ ¼ Ss ðxÞ fþ ðxÞ þ f ðxÞ þ S dx 2 2
(38)
1 djþ 1X þ i ðxÞ þ li jþ Tij jj ðxÞ þ j j ðxÞ þ Qi i ðxÞ ¼ dx 2 j 2 and
1 dj ðxÞ 1X þ i Tij jj ðxÞ þ j þ li j i ðxÞ ¼ j ðxÞ þ Qi dx 2 j 2
a) uniform distribution
1
PðSÞ ¼
dJ þ Sa f ¼ S dx
SðxÞ ¼ Smin þ ðSmax Smin Þx
(40)
where f ¼ fþ þ f and J ¼ fþ f . Eqs. (39) and (40) may then be combined to give
d2 f SSa f þ SS ¼ 0 dx2
(46)
(47)
b) log-uniform distribution
PðSÞ ¼
logðSmax =Smin Þ
S
; Smin S Smax
(48)
and, for sampling, from Eq. (10)
(42)
and
fðaÞ ¼ JðaÞ:
; Smin S Smax
or, for sampling, from Eq. (10)
(41)
The associated boundary conditions for a surface source at x ¼ 0 is fþ ð0Þ ¼ f and for a bare surface at x ¼ a is f ðaÞ ¼ 0. In terms of f and J these become
fð0Þ þ Jð0Þ ¼ 2f
Smax Smin
where x is a random number uniformly distributed in the range (0,1).
and
df þ SJ ¼ 0 dx
(45b)
The parameters in Eqs. (37) and (38) have been modelled as follows. For example if the parameter is uncertain over a range (Smin,Smax), then
where f refer to the angular neutron flux in the positive and negative directions, respectively. These equations may be reduced to the following physically more appealing form
(39)
(45a)
(43)
We shall see that the constant source problem and the incident beam problem can be solved exactly from Eqs. (41)–(43). Strictly speaking, Sa ðxÞ and Ss ðxÞ should be written as Sa ðx1 Þ and Ss ðx2 Þ where x1 and x2 are independent random variables. However, to simplify the analysis in this introductory paper, we will assume that x1 ¼ x2 ¼ x. Physically, this means that the S’s are perfectly correlated and that any uncertainty in one is directly reflected in the other. This is not a necessary approximation, but it avoids the need to use second order polynomial chaos functions at this stage. It does not affect the principle of the method.
SðxÞ ¼ Smin
Smax Smin
x
(49)
where again x is a random number uniformly distributed in the range (0,1). Other distributions may be used. Note that if all of the parameters follow Eq. (47) then, as discussed earlier, the Tij are diagonal, i.e. proportional to the Kronecker delta dij. It is also important to note, that if the cross section uncertainty depends on several parameters, for example in a resonance this may be the peak value s0(x1) and the width Gn ðx2 Þ, then the associated polynomial chaos function would be of the form Fn ðx1 ; x2 Þ. The use of such polynomials is explained in Williams (2006) and Ghanem and Spanos (1991), but essentially the matrix elements in Eq. (9) would become
hFm SFn i ¼
Z
1 0
dx1
Z
1 0
dx2 Fm ðx1 ; x2 ÞSðs0 ðx1 Þ; Gn ðx2 ÞÞFn ðx1 ; x2 Þ
This complicates the issue but does not affect the general approach, as we have intimated above.
584
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
For the bare slab containing a constant source, the boundary conditions are jþ i ð0Þ ¼ ji ðaÞ ¼ 0. But for the case of an incident source on the face x ¼ 0 they are
j i ðaÞ ¼ 0 1 jþ : i ð0Þ ¼ f L
(50)
i0
We also note in Eqs. (45a) and (45b) that Qi ¼ SðL1 Þi0 . The intensity can be reconstructed from Eqs. (2) and (21), and we find for the mean intensity
D
E
fþ ðxÞ þ f ðxÞ ¼ fþ 0 ðxÞ þ f0 ðxÞ
hfðxÞi ¼
2 1 þ fn ðxÞ þ f n ðxÞ 2n þ 1 n¼1 N X
N pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2n þ 1 Lnm j m ðxÞ
R ¼
Eqs. (45a) and (45b) may be evaluated numerically by a second order, ordinary differential equation solver. In this work we used the IMSL routine DBVPFD. Results will be given below for the mean and variance. 4.1. Analytic solutions The exact solutions of Eq. (41) for the two cases are 1) Constant source (Q ¼ 1)
1
Sa
ewx þ ewðaxÞ 1 1 þ a þ ð1 aÞewa
!
1 a2 1 e2wa ð1 þ aÞ2 ð1 aÞ2 e2wa
:
(59)
4aewa 2
ð1 þ aÞ ð1 aÞ2 e2wa
:
(60)
1
Sa
1
ewx
(52)
:
(53)
ð1 þ aÞ2 ð1 aÞ2 e2wa
2ewx : 1þa
f0 ð0Þ ¼ pffiffiffiffiffiffipffiffiffiffi pffiffiffiffiffiffi Sa S þ Sa
f0 ð0Þ ¼
2a ð1 þ aÞewx þ ð1 aÞewð2axÞ ð1 þ aÞ2 ð1 aÞ2 e2wa
" 2S
(62)
(63)
Ss
rffiffiffiffiffiffi# 1
Sa S
(64)
This result is also equal to the rod value of Eq. (55) when x ¼ 0 and simplifies to
pffiffiffiffi
2 S f0 ð0Þ ¼ pffiffiffiffi pffiffiffiffiffiffi S þ Sa
The current for this problem is
and for the half-space
15
Eq. (62) is the transport theory equivalent of Eq. (53) when a ¼ N. Fortuitously, Eqs. (62) and (53) are equal when a ¼ N and simplify to
(54)
(55)
S
Sa
The exact value is:
:
For the half-space
1 JðxÞ ¼ f0 ðxÞ ¼
Ss
3
2) Albedo problem
2) Transmission problem (unit incident source)
ð1 þ aÞewx ð1 aÞewð2axÞ
f0 ð0Þ ¼
2sffiffiffiffiffiffi 14 S
1
a 1 ewa : Sa 1 þ a þ ð1 aÞewa
fðxÞ ¼ 2
Clearly as a/N, R ¼ 1, T ¼ 0, i.e. all particles incident will be reemitted. Again we remind the reader that, in Eqs. (55)–(61), S ¼ S(x) and so R ¼ R(x) and T ¼ T(x), i.e. explicit functions of the random variable x. These forms for R and T are particularly useful as they can be calculated by PCE and compared with the exact stochastic averages of R and T from Eqs. (59) and (60) by quadrature. We also note the exact results, with isotropic scattering, for the surface scalar intensity and albedo for the half-space solution of the exact transport equation (Davison, 1957), viz:
The exact value is, for unit source:
In the general case, the surface intensity is
1
(61)
1) Constant source problem
!
1þa
aS 2 ; T ¼ : 2 þ aS 2 þ aS
(51)
where w2 ¼ SSa and a ¼ w/S. All of these parameters are random and depend on x. For an infinite half-space where a/N,
fðxÞ ¼
(58)
Using (54) and (56) we find
T ¼
m¼0
fð0Þ ¼
R ¼ f ð0Þ ¼ 12ðfð0Þ Jð0ÞÞ T ¼ fþ ðaÞ ¼ 12ðfðaÞ þ Jð0ÞÞ
For the half-space, T ¼ 0; R ¼ ð1 aÞ2 =ð1 þ aÞ2 . For zero absorption, (59) and (60) reduce to
f n ðxÞ ¼
fðxÞ ¼
(57)
and
where
fðxÞ ¼
2aewx 1þa
Also of interest for Eqs. (54) and (56) are the reflection R and transmission T coefficients defined by
R ¼
and for the variance
s2f ðxÞ ¼
JðxÞ ¼
(56)
(65)
It is not apparent from the literature that the surface scalar intensity for the rod model has been shown to be exact before for both constant source and albedo half-space problems.
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588 Table 1 Mean and relative standard deviation of the cross sections for uniform (U) and loguniform (LU) distributions. Type
Sa
Ss
S
sSa =Sa
sSs =Ss
sS =S
U LU
0.03 0.02485
0.7000 0.6805
0.7300 0.7054
0.3848 0.4551
0.1649 0.1692
0.1740 0.1795
The exact reflection coefficient R for the albedo problem does not reduce to a simple form but can be written as
pffiffiffiffiffiffiffiffiffiffiffiZ 1 Jð0Þ ¼ R ¼ 1 2 1 c dmmHðm; cÞ
2.0
deterministic current Deterministic flux Mean current (uniform) Mean flux (uniform)
1.5
1.0
(66)
0
0.5
where c ¼ Ss =S and H(m;c) is the Chandrasekhar H function (Chandrasekhar, 1960). Thus Eq. (66) can be calculated to any desired degree of accuracy by quadrature. Eq. (66) is the transport theory equivalent of Eq. (59) when a ¼ N. The corresponding value for the rod model is from Eq. (59)
Jð0Þ ¼
585
0
0
1
2
3
4
5
x
1a : 1þa
(67)
Fig. 1. Mean flux and current for transmission of neutrons.
This, however, is not equal to Eq. (66). The stochastic average of any of the quantities defined by Eqs. (51)–(67) can be obtained by direct quadrature in the form
n f ¼
Z
1
n
dxf ðxÞ
(68)
0
where f is the quantity whose stochastic moment is needed. We shall be concerned with the mean value Cf D and the variance s2f ¼ Cf 2 D Cf D2 , which will be compared with the PCE equivalents. 5. Numerical examples and discussion We shall illustrate the above techniques by two simple examples, viz: a) Transmission through a slab b) Constant source in a slab
5.1. Data We now give the data employed in the calculations and discuss some of their properties. For examples a) and b), the absorption and scattering cross sections are required. Here we consider the data described by either a uniform or a log-uniform distribution, as in Eqs. (46)–(49), with the limits
cm1 ;
cm1
Sa min ¼ 0:01 Samax ¼ 0:05 Ss min ¼ 0:5 cm1 ; Ssmax ¼ 0:9 cm1
uniform 1 1 S ¼ ðS0 þ S1 Þ; s2S ¼ ðS1 S0 Þ2 2 12 loguniform ðS1 S0 Þ ðS S0 ÞððS1 þ S0 ÞlogðS1 =S0 Þ2ðS1 S0 ÞÞ S¼ ; s2 ¼ 1 logðS1 =S0 Þ S log2 ðS1 =S0 Þ (70) It is seen that the relative standard deviation arising from the log-uniform distribution is always greater than that from the uniform distribution. These fluctuations will be reflected in the solution of the associated transport equation as will be demonstrated below. A useful measure of the sensitivity of the system to parameter uncertainties is the ‘sensitivity coefficient’, defined by the ratio of the relative standard deviation of the parameter of interest to the relative standard deviation of the uncertain parameter. For example, if we wish to know how sensitive a parameter R is to a given uncertainty in a cross section S, then we must calculate
0.20
(69)
0.15
The statistical moments of these coefficients are given by the following expressions and their numerical values are shown in Table 1. Note that we have set S0 ¼ Smin and S1 ¼ Smax.
0.10
Relstd current (loguni) Relstd flux (loguni) Relstd current (uni) Relstd flux (uni)
Table 2 Reflection and transmission coefficients for the rod model in a slab of thickness 5 cm.
0.05 Type
R
RelstdR
SR
T
RelstdT
ST
U LU ULU
0.5652 0.5695 0.5758
0.02407 0.02717 0.02784
0.138 0.158 –
0.3034 0.3198 0.3136
0.1930 0.1845 0.1859
1.109 1.031 –
U ¼ uniform distribution for Ss and Sa ; LU ¼ log-uniform distribution for Ss and Sa . ULU ¼ uniform distribution for Ss and log-uniform pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi for Sa . Relative standard deviation ¼ Relstd ¼ variance=mean: SR and ST are the sensitivity coefficients of R and T to cross section uncertainty.
0
0
1
2
3
4
x Fig. 2. Relative standard deviation for transmission of neutrons.
5
586
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
Table 3 Polynomial chaos expansion coefficients for R and T log-uniform distribution, a ¼ 5.0 cm.
Table 5 Surface flux, relative standard deviation and sensitivity for the albedo problem. Uniform distribution (rod model).
n
Rn exact
Rn PCE
Rn exact
Tn PCE
a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0) Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0) S (a)
1 2 3 4 5 6 7 8 9 10
0.56947 2.51912E2 1.17658E2 1.18645E3 1.88353E5 1.87656E5 2.46500E6 5.98277E8 3.30194E8 6.24389E9
0.56947 2.51916E2 1.17658E2 1.18647E3 1.88255E5 1.87647E5 2.46050E6 5.97112E8 2.79565E8 6.03360E9
0.31981 0.10217 4.06331E3 2.24357E5 9.64129E5 2.28575E5 2.17992E6 4.64831E8 4.67760E8 7.27465E9
0.31981 0.10217 4.06367E3 2.25064E5 9.64043E5 2.28587E5 2.17771E6 4.64831E8 4.39209E8 7.77338E9
5 10 20 50 100 N
1.5652 1.6432 1.6680 1.6712 1.6712 1.6712
0.2972 0.1306 2.909E2 3.429E4 2.097E7 0
0.1930 0.3594 0.7166 1.6112 2.5899 –
SR ¼ ðsR =RÞ=ðsS =SÞ
(71)
This is a useful measure as it shows the fractional change in the property of interest due to a given fractional change in a particular cross section. 5.2. Transmission through a slab We compare in Table 2 the average reflection, R, and transmission T, coefficients for the slab using three types of probability distribution. The slab is of thickness 5 cm. In Fig. 1 we show the mean flux and current for this problem for uniform and log-uniform distributions. To be more specific we define the flux and the current as given by Eqs. (39) and (40). The results are compared with the corresponding deterministic values with average parameters. The difference between the two cases is barely noticeable graphically, which is not surprising since the range of values in (69) is not great. Of more interest is a measure of the fluctuations, viz: the relative standard deviation for the U and LU cases, which we show in Fig. 2. Here we note that the fluctuations in the flux decrease initially and then increase towards the other face. The current fluctuations increase monotonically as distance increases. The general tendency for uncertainty to increase with depth is not unexpected because the effect of variations in the cross sections compound with distance from the source. The results in Table 2 and Figs. 1 and 2 are from Eqs. (45a) and (45b) using 20 terms in the PCE. In fact, only 6 terms are needed to get 4 significant figure accuracy for the variance, but it is instructive to examine the relative magnitudes of the expansion coefficients. The accuracy of the results is obtained from the exact analytical results in Eqs. (59) and (60) in the form given by Eq. (68). As a measure of the accuracy of the convergence of the PCE series, we have shown in Table 3 the PCE coefficients for R and T for the loguniform distribution which is a more severe test than the uniform case. The exact expansion coefficients are obtained by representing Eqs. (59) and (60) in terms of polynomial chaos expansions. It is clear that convergence is very rapid, but it should be noted that the values given for the PCE entries were from calculations using 20 terms in the series, only 10 of which are shown here. In Table 4 we give some exact results for the surface flux and reflection coefficients for a half-space as defined by Eqs. (65)–(67).
Table 4 Exact transport theory values for a half-space surface flux and reflection coefficient. Type
U LU
Rod
Exact
Cf0 ð0ÞD
Relstd
R
Relstd
R
Relstd
1.6712 1.6938
2.0669E2 2.1760E2
0.6356 0.6597
5.7850E2 5.9629E2
0.6712 0.6938
5.1463E2 5.3123E2
1.5689 1.6432 1.6619 1.6629 1.6629 1.6629
8.691E3 7.712E3 1.825E2 2.065E2 2.067E2 2.067E2
0.3034 0.1392 3.805E2 1.667E3 2.181E5 0
0.050 1.109 0.044 2.066 0.105 4.118 0.119 9.260 0.119 14.88 0.119 –
The value of Cf0 ð0ÞD is the same for the rod and exact model. The rod model value of the reflection coefficient R underestimates it and overestimates the relative standard deviation. A feature of some interest is the variation of the surface fluxes, relative standard deviations and sensitivities as a function of slab width. In Tables 5–8 we give values for these quantities for the albedo problem with uniform and log-uniform distributions of the cross sections, calculated by the rod model and the computer code EVENT (DeOliveira, 1986), modified to include the pseudo-multigroup transfer coefficients. The columns headed f0;det refer to the cases where the average coefficients are used directly in the deterministic solutions; this is the crudest form of approximation. There are three matters for discussion regarding Tables 5–8. Firstly, the difference between the rod model and exact equation, secondly the general behaviour of the mean and variance and finally the accuracy of the results. There is a difference between the two transport models as would be expected. Both models are in reasonable agreement at the source face at x ¼ 0, but at the exit face there are significant differences, both in mean value and in variance. This is to be expected since uncertainties accumulate with distance and the rod model is clearly not effective in diffusing the neutrons around due to the backward-forward nature of the motion. The important measure of uncertainty is the sensitivity S. We note that this is relatively small at the entry face because the source is deterministic and the neutrons have not travelled deeply into the medium. However, at the exit face, we observe a substantial increase in sensitivity as the slab thickness increases. The difference between the two models for sensitivity is however not great and both predict large values differing at most by a factor of 2 but often by much less. To some extent the accuracy of the EVENT results may be assessed from the half-space values for the surface fluxes and variances because for this case we have the exact analytical results, as given by Eq. (65). We add the caveat here that the EVENT results for a half-space are not as precise as those for a finite slab. The reason being that the half-space requires a very thick slab and the consequent additional number of mesh points makes the resolution not as sharp as for the finite slab. We elaborate further on this matter below. In practice, for the half-space, the mean and mean square intensity are only found to 3 significant figures whereas the variance, which is the difference between two nearly equal numbers, is of order 104 for the data used, and thus will be inaccurate. This does not occur for the thinner slabs, i.e. less than 50 mean free paths. For this reason the results given in Table 9 show the mean and the mean square separately and their
Table 6 Surface flux, relative standard deviation and sensitivity for the albedo problem. Uniform distribution (EVENT). a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0) Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0)
S (a)
5 10 20 50
1.6343 1.6593 1.6631 1.6637
0.1491 4.083E2 3.303E3 3.405E6
7.701E2 0.1862 0.4193 1.1345
0.443 1.070 2.410 6.520
1.6355 1.6592 1.6608 1.6620
3.517E3 9.013E3 1.061E2 1.072E2
0.6972 0.4184 0.1412 7.994E3
0.020 0.052 0.061 0.062
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588 Table 7 Surface flux, relative standard deviation and sensitivity for the albedo problem. Loguniform distribution (rod model). a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0)
Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0) S (a)
5 10 20 50 100 N
1.5695 1.6572 1.6890 1.6938 1.6938 1.6938
0.3198 0.1553 4.789E2 2.646E3 3.974E5 0
0.3136 0.1465 3.777E2 7.096E4 9.460E7 0
0.1845 0.3309 0.6178 1.2430 1.8856 –
1.5733 1.6573 1.6821 1.6839 1.6839 1.6839
9.859E3 8.336E3 1.902E2 2.174E2 2.176E2 2.1760E2
0.055 1.028 0.046 1.843 0.106 3.442 0.121 7.144 0.121 10.50 0.121 –
associated errors are compared with the exact results. The errors are less than 0.5% and 1%, respectively, and show that the method is viable. The number of terms in the chaos expansion in EVENT is 5. Additional terms make little difference to the results. The mean values from EVENT are accurate to within about 0.5%. However the relative standard deviation has an error of about 50%. There are good reasons for this discrepancy which concern the structure of EVENT. The code uses a finite element, spherical harmonics representation of the angular flux based upon a variational formulation of the second order even parity form of the transport equation (Vladimirov, 1963). Linear finite elements are used in space with the resolution of the mesh tuned in order to reach spatial convergence without making the cells too thin optically. The caveat on not making the cells too small is because the even parity equation performs poorly for cases where the optical thickness of the element is small (Ackroyd et al., 1986). If therefore the mesh resolution is made too fine the simulations may fail to converge. A spherical harmonics expansion of order N ¼ 99 (Bell and Glasstone, 1970) was used to ensure angular convergence and double precision was also employed. The problem of accuracy lies therefore not only with the limitations just noted, but also because the point of interest is at the surface where the flux varies very steeply. In addition, the variance depends on the difference between two surface quantities and this leads to error accumulation. However, the purpose of this paper is not so much to obtain accuracy but to assess feasibility. It is clear that existing multi-group neutron transport codes can be used for this problem thereby obviating the need for significant development work. The accuracy will be much greater for parameters that do not depend on surface values; for example reaction rates in local regions.
5.3. Constant source in a slab
Table 8 Surface flux, relative standard deviation and sensitivity for the albedo problem. Loguniform distribution (EVENT). a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0) Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0)
S (a)
5 10 20 50
1.6476 1.6792 1.6849 1.6858
0.1644 5.125E2 5.378E3 5.612E5
7.719E2 0.1795 0.3674 0.8061
0.430 1.000 2.047 4.491
4.726E3 1.038E2 1.206E2 1.222E2
0.7399 0.4795 0.1872 1.384E2
Table 9 Comparison of EVENT values at surface of half-space with exact analytic results. Prob distn
Mean (EVENT)
Mean (EXACT)
%
Meansquare (EVENT)
Relstd (EXACT)
%
U LU
1.6637 1.6869
1.6712 1.6938
0.45 0.41
2.7682 2.8460
2.7942 2.8704
0.93 0.85
Table 10 Mean flux, relative standard deviation and sensitivity coefficient of surface flux for constant source problem in a slab of thickness 5 cm. Type
Mean flux
Relstd
S
U LU ULU
2.2311 2.2753 2.2733
0.047849 0.046203 0.046380
0.275 0.258 –
In Table 11 we give the exact surface values for the flux arising from a constant source in an infinite half-space. In this case we note a significant difference for the mean intensity between uniform (U) and log-uniform (LU) distributions. The relative standard deviation however is not so much affected. A feature of some interest is the variation of the surface fluxes and relative standard deviations as a function of slab width. In Table 12 we give values for these quantities for the constant source problem. It is interesting to note from Table 12 that, for the source problem, the flux and variance depend strongly on the slab width; increasing as the width increases. For the 5 cm slab the Relstd is about 5%, whereas for the slab of infinite thickness it is about 35%. The sensitivity S also increases rapidly with slab thickness. 5.4. Some comments on the meaning of statistical moments Occasionally in the literature the significance of the mean and its variance are misunderstood. In this respect, and with relation to the above theory, it is important to note that the stochastic properties obtained from the above analysis do not refer to a specific system but to an ensemble of such systems. For example, one particular system will have one set of cross sections and will respond deterministically to those cross sections. Thus CfD and s2f give a measure of the uncertainty arising from a very large number (theoretically infinite) of systems whose cross sections vary
0.08
For completeness, we also give results for the case of a uniform constant source in a bare slab. Only the rod model is used but the results are instructive. For this case we show, in Table 10, the value of the flux at the surface and its relative standard deviation. Only small differences are noted. These may be seen more clearly in Fig. 3 for the relative standard deviation across the slab. Clearly, for the flux, it is symmetric about the centre and for the current anti-symmetric because the current changes sign at the mid-point.
1.6494 1.6785 1.6823 1.6823
587
0.026 0.058 0.067 0.068
0.04
0
Relstd Current (loguni) Relstd flux (loguni) Relstd current (uni) relstd flux (uni)
-0.04
-0.08
0
1
2
3
4
x Fig. 3. Relative standard deviation for constant source problem.
5
588
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
Table 11 Exact transport theory values for a half-space surface flux.
Table 12 Surface flux, relative standard deviation and sensitivity for the source problem. Uniform distribution.
Type
Cf0 ð0ÞD
Relstd
S
U LU
6.2824 7.2045
0.3450 0.3391
1.983 1.894
according to the distribution p(S). To put it otherwise, the mean and average tell us only that any prediction about the system of interest will be uncertain by the amounts indicated by the stochastic moments; it is not known which of this myriad of ensemble members actually corresponds to our system. An illustrative example arises by considering the simple model equation, viz:
dfðt; xÞ ¼ aðxÞfðt; xÞ dt
(72)
from which we have
fðt; xÞ ¼ fð0Þexp½aðxÞt
(73)
a
Cf0 ð0ÞD
f0;det ð0Þ
Relstd (0)
S
5 10 20 50 100 N
2.2311 3.7922 5.3591 6.2244 6.2815 6.2824
2.2317 3.7703 5.1502 5.6127 5.6184 5.6184
0.04785 0.1139 0.2261 0.3373 0.3497 0.3500
0.275 0.655 1.300 1.939 2.010 2.011
into a form which can be solved using existing computer codes with only small modifications. Finally, we comment on the eigenvalue problem that arises when fission is included. In that case, the eigenvalue keff will also be a random variable and it is necessary to expand it in terms of a PCE. However, the resulting equations are non-linear and require special root-finding techniques that are beyond the scope of the present paper. However, in part II of this work we shall describe a new method which avoids these non-linearities (Williams, 2010).
The mean value of f is then
hfðt; xÞi ¼ fð0Þ
Z
b
dxpðxÞexp½aðxÞt
(74)
Acknowledgements
a
For a uniform distribution, as given by Eq. (47), we can evaluate the integral in (67) to get
hfðt; xÞi ¼
i fð0Þ h a0 t ea1 t e tða1 a0 Þ
Dr. Eaton would like to acknowledge the support of The Royal Academy of Engineering and EPSRC under their Research Fellowship Scheme.
(75)
In actual practice, we can only measure one system with specific values of coefficients determined by a value x ¼ xi. Thus the solution would be as in Eq. (73) with x ¼ xi. This of course varies markedly from Eq. (75) and stresses the point made above. 6. Summary and conclusions The purpose of this work has been to develop a method for assessing the influence of data uncertainty on some important parameters in neutron transport calculations. We have concentrated on surface quantities such as reflection and transmission coefficients as well as surface fluxes for some simple, one-dimensional problems. The essence of the method is the representation of the data by random variables, the statistics of which reflect the nature of the uncertainty. For example, a piece of data may be known between upper and lower limits with a certain probability distribution. In this case we have used uniform and log-uniform distributions but we could equally well have used any other distribution. The method is particularly useful when the data uncertainty is large and the usual techniques of perturbation theory fail. It also avoids the often unphysical representation of errors by Gaussian variables which, for large uncertainty, can lead to negative values for a physical quantity. The use of polynomial chaos is adopted which is a well-proven technique in structural mechanics and fluid dynamics, but has only recently found its way into neutron transport and radiative transfer calculations (Emery, 2005; Williams, 2006, 2007). It may well be, however, that this paper is the first to use it for assessing data uncertainty in neutron transport and to cast the resulting equations, which may be quite complex,
References Ackroyd, R.T., Issa, J.G., Riyait, N.S., 1986. Treatment of voids in finite element transport methods. Prog. Nucl. Energy 18, 85. Bell, G.I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand, Reinhold Company. Cacuci, D.G., 2003. Sensitivity and Uncertainty Analysis. Chapman and Hall. Carter, L.L., Cashwell, E.D., 1975. Particle Transport with the Monte Carlo Method. USAEC. Report. Chandrasekhar, S., 1960. Radiative Transfer. Dover Publications. DeOliveira, C.R.E., 1986. An arbitrary geometry finite element method for multigroup neutron transport with anisotropic scattering. Prog. Nucl. Energy 18, 227. Davison, B., 1957. Neutron Transport Theory. Oxford University Press. Emery, A.F., 2005. Some thoughts on solving the radiative transfer equation in stochastic media using polynomial chaos and wick products as applied to radiative equilibrium. J. Quant. Spectrosc. Radiat. Transf. 93, 61. Ghanem, R.G., Spanos, P.D., 1991. Stochastic Finite Elements. Dover Publications. Ghanem, R.G., 1998. Scales of fluctuation and the propagation of uncertainty in random porous media. Water Resour. Res. 34, 2123. Li, R., Ghanem, R., 1998. Adaptive polynomial chaos expansions applied to statistics of extremes in non-linear random vibration. Prob. Eng. Mech. 13, 125. Le Maitre, O.P., Knio, O.M., Najm, H.N., Ghanem, R.G., 2001. A stochastic projection method for fluid flow. J. Comp. Phys. 173, 481. Vladimirov, V.S., 1963. Mathematical Problems in the One-Velocity Theory of Particle Transport. Atomic Energy of Canada Ltd. Report AECL 1661. Wiener, N., 1938. The homogeneous chaos. Am. J. Math. 60, 897. Williams, M.L., 1986. Perturbation theory for nuclear reactor analysis. In: Ronen, Y. (Ed.), Handbook of Nuclear Reactor Calculations, vol. III. CRC Press. Williams, M.M.R., 2010. A method for solving stochastic eigenvalue problems. Appl. Math. and Comput. 215, 3906. Williams, M.M.R., 2006. Polynomial chaos functions and stochastic differential equations. Ann. Nucl. Energy 33, 774. Williams, M.M.R., 2007. Polynomial chaos functions and neutron diffusion. Nucl. Sci. Eng. 155, 109. Xiu, D., Karniadakis, G.E., 2002a. Modelling uncertainty in steady state diffusion problems via generalized polynomial chaos. Comput. Methods Appl. Mech. Eng. 191, 4927. Xiu, D., Karniadakis, G.E., 2002b. The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM J. Sci. Comput. 24, 619.
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
A probabilistic study of the influence of parameter uncertainty on solutions of the neutron transport equation Matthew Eaton, 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 2BP, UK
a b s t r a c t Keywords: Uncertainty Neutron transport Stochastic
The influence of uncertainty in the scattering and absorption cross sections on the solution of the neutron transport equation is studied using polynomial chaos theory. The uncertainty is defined by means of uniform and log-uniform probability distributions. By expanding the neutron flux in a series of polynomial chaos functions we may reduce the stochastic transport equation to a set of coupled deterministic equations, analogous to those that arise in multi-group neutron transport theory, with the effective multi-group transfer scattering coefficients containing information about the uncertainty. This procedure enables existing transport theory computer codes to be used, with little modification, to solve the problem. Applications are made to a transmission problem and a constant source problem in a slab. In addition, we also study the rod model for which exact analytical solutions are readily available. In all cases, numerical results in the form of mean, variance and sensitivity are given which illustrate how absorption and scattering cross section uncertainty influences the solution of the transport equation. 2010 Elsevier Ltd. All rights reserved.
1. Introduction An important source of error in calculations employing the neutron transport equation arises from uncertain data values. These errors originate from lack of knowledge about the absorption and scattering cross sections. It is vital, therefore, to know how changes in these parameters affect the outcome of any calculations. For example, what are the uncertainties in the reflection and transmission coefficients of a slab in relation to uncertainties in the absorption and scattering cross sections? When the uncertainties are relatively small, these matters can be been handled by perturbation theory (Williams, 1986; Cacuci, 2003). However there is another effective method of dealing with such a problem which does not have the limitations of perturbation theory and which may be used for arbitrarily large parameter uncertainties; namely, polynomial chaos expansions. In order to proceed, therefore, we deliberately introduce uncertainty into the parameters of the transport equation, thereby converting it into a stochastic integrodifferential equation. This is done by assuming that the parameters are functions of a random variable which, in turn, allows us to prescribe a probability distribution for the values of the parameters. The polynomial chaos expansion (PCE) (Ghanem and Spanos, 1991;
* Corresponding author. E-mail address: [email protected] (M.M.R. Williams). 0149-1970/$ – see front matter 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.pnucene.2010.01.002
Emery, 2005; Wiener, 1938) is then employed to transform the stochastic equation into a set of coupled deterministic transport equations for the coefficients in the PCE. By a further transformation, these coupled equations may be cast into a form that is analogous to a multi-group neutron (or radiation) transport formulation but with probabilistic matrix elements taking the place of the group-to-group energy exchange cross sections. The only difference between the equations derived below from those of multi-group theory is the presence of coupled boundary conditions. This transformation therefore enables us to use well-established computer codes to solve the associated multi-group equations. In this paper we will give examples based upon the onedimensional transport equation for a slab. Two cases are considered: 1) a constant source in a bare slab and 2) the transmission of neutrons through a plane slab with an isotropic source incident on one face. We will solve the transport equation numerically but, because certain analytic solutions exist for the cases we have chosen, we will be able to compare the PCE method with the exact result for mean and variance and hence estimate convergence rates. We also solve the problems using the rod model, for which exact solutions are also available over the complete range of distance. Some additional calculations are carried out, which use the exact transport equation, by a modification of the computer code EVENT (DeOliveira, 1986). The general conclusion is that the PCE method shows considerable promise and can deal with a variety of probability distribution functions describing parameter uncertainty.
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
It is important to note that a number of other fields of technology have benefited from the use of PCE to study the effects of data uncertainty. For example, the uncertainty of modes of vibration in structural mechanics due to errors in elastic constants (Li and Ghanem, 1998), the flow of groundwater in the geo-sphere due to rock fractures (Ghanem, 1998) and uncertainties in fluid flow due to both geometric uncertainty and basic data (Le Maitre et al., 2001). The extension to neutron transport is an obvious extension, as we hope to convince the reader.
D
The method to be described works just as well for the energy dependent transport equation but for simplicity in this introductory paper we use the one speed equation for the angular neutron flux fðr; U; xÞ, which may be written as (Davison, 1957)
U:Vfðr; U; xÞ þ SðxÞfðr; U; xÞ ¼ Ss ðxÞ
Z
dU0 f U:U0 ; x f r; U0 ; x þ Sðr; U; xÞ
N X
fn ðr; UÞFn ðxÞ
¼ hFm ðxÞSðr; U; xÞi Z N X fn r; U0 Fm ðxÞSs ðxÞf U:U0 ; x Fn ðxÞ : þ dU0
ð1Þ
(2)
where Fn ðxÞ are the polynomial chaos functions which are to be determined by the nature of the problem and obey the orthogonality relation
hFn ðxÞFm ðxÞih
R
2 dxpðxÞFn ðxÞFm ðxÞ ¼ dnm Nm
(3)
pðxÞ being an appropriate probability weight function of the random variables x. Eq. (3) also defines the normalisation coefficient Nm. For full details of the concept of polynomial chaos functions, the reader is referred to Ghanem and Spanos (1991). However, in the simplest case of one random variable, x, the polynomials form a complete set of orthogonal functions with the weight function given by the appropriate probability distribution (Xiu and Karniadakis, 2002a,b). For more than one random variable, the sequence of polynomials can be generalised as explained in Ghanem and Spanos (1991) or Williams (2006). The expansion coefficients fn ðr; UÞ are to be determined in the manner described below. Inserting Eq. (2) into Eq. (1) we find
U:V
N X
fn ðr; UÞFn ðxÞ þ SðxÞ
n¼0
¼ S s ðxÞ
N X
fn ðr; UÞFn ðxÞ
Let us now define
dU0 f ðU:U0; xÞ
N X
fn r; U0 Fn ðxÞ þ Sðr; U; xÞ:
hFm ðxÞSðxÞFn ðxÞi Nn Nm
Fm ðxÞSs ðxÞf U:U0 ; x Fn ðxÞ Bmn U:U0 ¼ Nn Nm
(6)
(7)
Sm ðr; UÞ ¼ hFm ðxÞSðr; U; xÞi=Nm
b m ðr; UÞ ¼ Nm f ðr; UÞ: f n
(8)
We observe that the matrices A and B are symmetric. The angular brackets in Eqs. (4)–(7) are defined by
hFm SFn i ¼
Z
dxpðxÞSðxÞFm ðxÞFn ðxÞ
(9)
where pðxÞ is the weight function associated with PCE polynomials Fn . It should be noted that the evaluation of this matrix element depends on the nature of the statistics of SðxÞ and x. For example, if S obeys a prescribed probability distribution p1 ðSÞ and suppose that we have decided that the weight function associated with the PCE polynomial Fn ðxÞ is pðxÞ; then if these two distributions span different spaces the integral (9) as it stands needs interpretation. For simplicity we consider the case where there is only one random variable x, then we can define a uniformly distributed random variable u˛ð0; 1Þ, and write (Carter and Cashwell, 1975)
u ¼
Z S S0
p 1 S 0 dS 0
(10)
Thus as u varies from zero to unity, so the value of S will vary from its minimum of S0 to its maximum of S1 and be weighted according to p1 ðSÞ. But we must also refer x to the uniform variable u which we do by setting
u ¼
Z x 0 0 p x dx :
(11)
x0
From (10) we have the inverse S ¼ f(u) and from (11) the inverse x ¼ g(u). In some cases explicit expressions can be obtained for f and g, but in general this is not so. Now we transform Eq. (9) by setting du ¼ pðxÞdx, whence
hFm SFn i ¼
Z
1 0
duf ðuÞFm ðgðuÞÞFn ðgðuÞÞ:
(12)
With the matrix elements defined, we may now write Eq. (5) in the form
b m ðr; UÞ þ U:V f
n¼0
Z
ð5Þ
n¼0
n¼0
Z
fn ðr; UÞhFm ðxÞSðxÞFn ðxÞi
b m as We also define the modified flux component f
In Eq. (1) the total cross section SðxÞ ¼ Sa ðxÞ þ Ss ðxÞ, where Sa and Ss are the absorption and scattering cross sections, respectively. f ðU:U0; xÞ is the differential scattering probability per unit solid angle. Sðr; U; xÞ is an independent source term which itself may contain uncertainties. The xðx1 ; x2 ; .xN Þ are independent random variables which will describe the data uncertainty. It should be noted that we will not be considering multiplying systems in this paper. The introduction of fissile material transforms the equation into an eigenvalue problem and this raises a host of difficulties. Treatment of these difficulties is deferred to a future publication, but some discussion is given in the conclusions in Section 6. Let us assume that the random flux can be expanded in a PCE in the form (Ghanem and Spanos, 1991; Williams, 2006)
fðr; U; xÞ ¼
N X n¼0
Amn ¼
2. Theory
E
U:Vfm ðr; UÞ F2m ðxÞ þ
581
N X
b n ðr; UÞ Amn f
n¼0
(4)
n¼0
Now multiplying Eq. (4) by pðxÞFm ðxÞ and averaging over all the random variables x, leads to
¼
N Z X
b 0 dU0 Bmn U:U0 f n r; U þ Sm ðr; UÞ:
n¼0
In matrix form, Eq. (13) becomes
(13)
582
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
b þ Af b ¼ Bf b þ S: U:V f
(14)
To solve this equation, it is useful to diagonalise the matrix A as follows 1
A ¼ L LL
(15)
where L is a NN diagonal matrix with diagonal elements li and L is defined by the eigenvalue problem
AR ¼ lR
(16)
where L is an NN matrix with column vectors given by each eigenvector R of (16) corresponding to the eigenvalue li. Inserting A as given by (15) in Eq. (14) and re-arranging gives 1
1
U:Vj þ Lj ¼ L BLj þ L S
(17)
b ¼ Lj. In discrete form Eq. (17) becomes where we have set f
U:Vji ðr; UÞ þ li ji ðr; UÞ ¼
N Z X
dU0 Tij U:U0 jj r; U0 þ Qi ðr; UÞ
j¼0
(18)
It should be noted that where it is typographically more elegant, we use an over-bar to denote an average instead of C.D. 3. Boundary conditions Let us consider the conditions that must hold between two dissimilar media. We know that the angular intensity is continuous at an interface so that we may write for two adjacent media 1 and 2, at the boundary point rs,
fð1Þ ðrs ; U; xÞ ¼ fð2Þ ðrs ; U; xÞ; n:Us0
where n is a unit vector normal to the interface. Now we insert the PCE’s for each side, which in general will be different because the cross sections are different and we may wish to choose different chaos polynomials. We then find N X
0
Tij U:U and
Qj ¼
¼ L1 B U:U0 L
(19)
ij
N X
fnð2Þ ðrs ; UÞFnð2Þ ðxÞ:
(25)
n¼0
Now multiply by Fð1Þ m ðxÞ and average to get N D E X ð1Þ 2 ð1Þ ð1Þ ð2Þ fm ðrs ; UÞ ¼ fð2Þ Nm n ðrs ; UÞ Fm ðxÞFn ðxÞ
(26)
n¼0
L1 S :
Using the relation between f and j from Eq. (21), viz
(20)
j
It is important to note that if the probability distribution p1 ðSÞ, describing the distribution of S(x), is also the weight function for the orthogonality condition of the Fn(x), then the matrix T will be diagonal. The diagonal elements of T are equal to the eigenvalues of the diagonalised matrix B, analogous to Eq. (16). This means that Eq. (18) decouples and the solutions only differ because of different sources and boundary conditions. If this condition is not met, or one of the parameters has a different distribution, e.g. log-uniform, then off-diagonal matrix elements will appear and there is no decoupling. Generally, these off-diagonal elements are significantly smaller in magnitude than those of the diagonal ones. Eq. (18) is in the form of standard multi-group neutron transport equations with the stochastic matrix elements related to the groupto-group transfer cross sections. Thus it is possible to use existing multi-group computer codes to solve for the PCE coefficients. To recover the components of the flux in Eq. (2) we use Eq. (8) and the b ¼ Lj as defined above, just below Eq. (17), to get relation f
fi ðr; UÞ ¼
ð1Þ fð1Þ n ðrs ; UÞFn ðxÞ ¼
n¼0
with i ¼ 0,1,2,.N, and we have defined
(24)
fð1Þ m ðrs ; UÞ ¼
N 1 X ð1Þ
Nm
(21)
(27)
ð2Þ ð2Þ
(28)
and
fð2Þ m ðrs ; UÞ ¼
N 1 X ð2Þ
Nm
Lmj jj ðrs ; UÞ
j¼0
and inserting (27) and (28) into (26) we find ð1Þ
Nm
N X
ð1Þ ð1Þ
Lmj jj ðrs ; UÞ¼
N D X
ð2Þ Fð1Þ m Fn
n¼0
j¼0
N E 1 X
ð2Þ ð2Þ Lnj jj ðrs ; UÞ ð2Þ Nn j¼0
(29)
which on re-arrangement gives N X
ð1Þ ð1Þ
Lmj jj ðrs ; UÞ¼
j¼0
N 1X Lij jj ðr; UÞ Ni
ð1Þ ð1Þ
Lmj jj ðrs ; UÞ
j¼0
N X
jð2Þ ðrs ; UÞ j
j¼0
N X
1 ð2Þ L ð1Þ ð2Þ nj n¼0 Nm Nn
D
ð2Þ Fð1Þ m Fn
E
(30)
or in matrix form
j¼0
From Eq. (2) we may now calculate the mean and variance of the angular intensity and any higher moments of interest, viz: the mean is
hfðr; UÞihf0 ðr; UÞ ¼ f0 ðr; UÞ
(22)
and the variance
E
D
s2f ¼ ðfðr; UÞ fðr; UÞÞ2 ¼
N X
D E f2n ðr; UÞ F2n :
(23a)
For the scalar flux 2
2
sf0 ¼ ðf0 ðrÞ f0 ðrÞÞ
E
¼
N X
D E f20n ðrÞ F2n
(23b)
where
f0n ðrÞ ¼
Z
dUfn ðr; UÞ:
(23c)
(32)
for medium 2, which implies that, in this case only, each individual component of the PCE expansion satisfies the condition
jð2Þ n ðrs ; UÞ ¼ 0 : n:U < 0
n¼1
(31)
If we decide to retain the same PCE in each region as indeed may (1,2) ¼ L(2). For a bare be quite acceptable, then Fnð1Þ ¼ Fð2Þ n and Z surface the boundary condition reduces to
Lð2Þ jð2Þ ðrs ; UÞ ¼ 0 : n:U < 0
n¼1
D
Lð1Þ jð1Þ ¼ Z ð1;2Þ jð2Þ :
(33)
In Eqs. (32) and (33), n is an outward pointing unit normal to the surface. For the case where a source f(U;x), with uncertainty, is incident on the surface the condition becomes
fðrs ; U; xÞ ¼ f ðU; xÞ; n:U < 0
(34)
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
Expanding in PCE as in Eq. (2), and using the orthogonality relation, we find
1 b m ðrs ; UÞ: hFm ðxÞf ðU; xÞi ¼ f Nm
(35)
b ¼ Lj, Eq. (35) becomes Then from the relation f
583
The PCE expansion of Eqs. (37) and (38) is given by reference to Eqs. (18) and (21) where we have chosen the Legendre polynomials for the polynomial chaos expansion because the random variable x˛(0,1). Thus we write
f ðxÞ ¼
N X
f n ðxÞPn ð2x 1Þ
(44a)
n¼0
jm ðrs ; UÞ ¼
N X
L1
n¼0
1 mn Nn
hFn ðxÞf ðU; xÞi
(36)
Thus in order to use the standard multi-group codes we must change the boundary conditions as shown in Eqs. (31) and (36). For vacuum boundaries, no changes are needed.
where the weight function p(x) ¼ 1 for such polynomials (Xiu and Karniadakis, 2002a,b). After diagonalisation,
f i ðxÞ ¼
N N pffiffiffiffiffiffiffiffiffiffiffiffiffi X 1X Lij j Lij j 2i þ 1 j ðxÞ ¼ j ðxÞ: Ni j ¼ 0 j¼0
(44b)
The expansion coefficients are then given by 4. The rod model It is useful to illustrate the above theory using a transport equation for which exact solutions can be obtained. For the case of the slab, the rod model is appropriate. In this model, the neutrons are supposed to travel either in the positive or negative x-directions. Thus the transport Eq. (1), may be written in slab geometry as (setting m ¼ 1)
1 dfþ ðxÞ 1 þ SðxÞfþ ðxÞ ¼ Ss ðxÞ fþ ðxÞ þ f ðxÞ þ S dx 2 2
(37)
and
1 df ðxÞ 1 þ SðxÞf ðxÞ ¼ Ss ðxÞ fþ ðxÞ þ f ðxÞ þ S dx 2 2
(38)
1 djþ 1X þ i ðxÞ þ li jþ Tij jj ðxÞ þ j j ðxÞ þ Qi i ðxÞ ¼ dx 2 j 2 and
1 dj ðxÞ 1X þ i Tij jj ðxÞ þ j þ li j i ðxÞ ¼ j ðxÞ þ Qi dx 2 j 2
a) uniform distribution
1
PðSÞ ¼
dJ þ Sa f ¼ S dx
SðxÞ ¼ Smin þ ðSmax Smin Þx
(40)
where f ¼ fþ þ f and J ¼ fþ f . Eqs. (39) and (40) may then be combined to give
d2 f SSa f þ SS ¼ 0 dx2
(46)
(47)
b) log-uniform distribution
PðSÞ ¼
logðSmax =Smin Þ
S
; Smin S Smax
(48)
and, for sampling, from Eq. (10)
(42)
and
fðaÞ ¼ JðaÞ:
; Smin S Smax
or, for sampling, from Eq. (10)
(41)
The associated boundary conditions for a surface source at x ¼ 0 is fþ ð0Þ ¼ f and for a bare surface at x ¼ a is f ðaÞ ¼ 0. In terms of f and J these become
fð0Þ þ Jð0Þ ¼ 2f
Smax Smin
where x is a random number uniformly distributed in the range (0,1).
and
df þ SJ ¼ 0 dx
(45b)
The parameters in Eqs. (37) and (38) have been modelled as follows. For example if the parameter is uncertain over a range (Smin,Smax), then
where f refer to the angular neutron flux in the positive and negative directions, respectively. These equations may be reduced to the following physically more appealing form
(39)
(45a)
(43)
We shall see that the constant source problem and the incident beam problem can be solved exactly from Eqs. (41)–(43). Strictly speaking, Sa ðxÞ and Ss ðxÞ should be written as Sa ðx1 Þ and Ss ðx2 Þ where x1 and x2 are independent random variables. However, to simplify the analysis in this introductory paper, we will assume that x1 ¼ x2 ¼ x. Physically, this means that the S’s are perfectly correlated and that any uncertainty in one is directly reflected in the other. This is not a necessary approximation, but it avoids the need to use second order polynomial chaos functions at this stage. It does not affect the principle of the method.
SðxÞ ¼ Smin
Smax Smin
x
(49)
where again x is a random number uniformly distributed in the range (0,1). Other distributions may be used. Note that if all of the parameters follow Eq. (47) then, as discussed earlier, the Tij are diagonal, i.e. proportional to the Kronecker delta dij. It is also important to note, that if the cross section uncertainty depends on several parameters, for example in a resonance this may be the peak value s0(x1) and the width Gn ðx2 Þ, then the associated polynomial chaos function would be of the form Fn ðx1 ; x2 Þ. The use of such polynomials is explained in Williams (2006) and Ghanem and Spanos (1991), but essentially the matrix elements in Eq. (9) would become
hFm SFn i ¼
Z
1 0
dx1
Z
1 0
dx2 Fm ðx1 ; x2 ÞSðs0 ðx1 Þ; Gn ðx2 ÞÞFn ðx1 ; x2 Þ
This complicates the issue but does not affect the general approach, as we have intimated above.
584
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
For the bare slab containing a constant source, the boundary conditions are jþ i ð0Þ ¼ ji ðaÞ ¼ 0. But for the case of an incident source on the face x ¼ 0 they are
j i ðaÞ ¼ 0 1 jþ : i ð0Þ ¼ f L
(50)
i0
We also note in Eqs. (45a) and (45b) that Qi ¼ SðL1 Þi0 . The intensity can be reconstructed from Eqs. (2) and (21), and we find for the mean intensity
D
E
fþ ðxÞ þ f ðxÞ ¼ fþ 0 ðxÞ þ f0 ðxÞ
hfðxÞi ¼
2 1 þ fn ðxÞ þ f n ðxÞ 2n þ 1 n¼1 N X
N pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2n þ 1 Lnm j m ðxÞ
R ¼
Eqs. (45a) and (45b) may be evaluated numerically by a second order, ordinary differential equation solver. In this work we used the IMSL routine DBVPFD. Results will be given below for the mean and variance. 4.1. Analytic solutions The exact solutions of Eq. (41) for the two cases are 1) Constant source (Q ¼ 1)
1
Sa
ewx þ ewðaxÞ 1 1 þ a þ ð1 aÞewa
!
1 a2 1 e2wa ð1 þ aÞ2 ð1 aÞ2 e2wa
:
(59)
4aewa 2
ð1 þ aÞ ð1 aÞ2 e2wa
:
(60)
1
Sa
1
ewx
(52)
:
(53)
ð1 þ aÞ2 ð1 aÞ2 e2wa
2ewx : 1þa
f0 ð0Þ ¼ pffiffiffiffiffiffipffiffiffiffi pffiffiffiffiffiffi Sa S þ Sa
f0 ð0Þ ¼
2a ð1 þ aÞewx þ ð1 aÞewð2axÞ ð1 þ aÞ2 ð1 aÞ2 e2wa
" 2S
(62)
(63)
Ss
rffiffiffiffiffiffi# 1
Sa S
(64)
This result is also equal to the rod value of Eq. (55) when x ¼ 0 and simplifies to
pffiffiffiffi
2 S f0 ð0Þ ¼ pffiffiffiffi pffiffiffiffiffiffi S þ Sa
The current for this problem is
and for the half-space
15
Eq. (62) is the transport theory equivalent of Eq. (53) when a ¼ N. Fortuitously, Eqs. (62) and (53) are equal when a ¼ N and simplify to
(54)
(55)
S
Sa
The exact value is:
:
For the half-space
1 JðxÞ ¼ f0 ðxÞ ¼
Ss
3
2) Albedo problem
2) Transmission problem (unit incident source)
ð1 þ aÞewx ð1 aÞewð2axÞ
f0 ð0Þ ¼
2sffiffiffiffiffiffi 14 S
1
a 1 ewa : Sa 1 þ a þ ð1 aÞewa
fðxÞ ¼ 2
Clearly as a/N, R ¼ 1, T ¼ 0, i.e. all particles incident will be reemitted. Again we remind the reader that, in Eqs. (55)–(61), S ¼ S(x) and so R ¼ R(x) and T ¼ T(x), i.e. explicit functions of the random variable x. These forms for R and T are particularly useful as they can be calculated by PCE and compared with the exact stochastic averages of R and T from Eqs. (59) and (60) by quadrature. We also note the exact results, with isotropic scattering, for the surface scalar intensity and albedo for the half-space solution of the exact transport equation (Davison, 1957), viz:
The exact value is, for unit source:
In the general case, the surface intensity is
1
(61)
1) Constant source problem
!
1þa
aS 2 ; T ¼ : 2 þ aS 2 þ aS
(51)
where w2 ¼ SSa and a ¼ w/S. All of these parameters are random and depend on x. For an infinite half-space where a/N,
fðxÞ ¼
(58)
Using (54) and (56) we find
T ¼
m¼0
fð0Þ ¼
R ¼ f ð0Þ ¼ 12ðfð0Þ Jð0ÞÞ T ¼ fþ ðaÞ ¼ 12ðfðaÞ þ Jð0ÞÞ
For the half-space, T ¼ 0; R ¼ ð1 aÞ2 =ð1 þ aÞ2 . For zero absorption, (59) and (60) reduce to
f n ðxÞ ¼
fðxÞ ¼
(57)
and
where
fðxÞ ¼
2aewx 1þa
Also of interest for Eqs. (54) and (56) are the reflection R and transmission T coefficients defined by
R ¼
and for the variance
s2f ðxÞ ¼
JðxÞ ¼
(56)
(65)
It is not apparent from the literature that the surface scalar intensity for the rod model has been shown to be exact before for both constant source and albedo half-space problems.
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588 Table 1 Mean and relative standard deviation of the cross sections for uniform (U) and loguniform (LU) distributions. Type
Sa
Ss
S
sSa =Sa
sSs =Ss
sS =S
U LU
0.03 0.02485
0.7000 0.6805
0.7300 0.7054
0.3848 0.4551
0.1649 0.1692
0.1740 0.1795
The exact reflection coefficient R for the albedo problem does not reduce to a simple form but can be written as
pffiffiffiffiffiffiffiffiffiffiffiZ 1 Jð0Þ ¼ R ¼ 1 2 1 c dmmHðm; cÞ
2.0
deterministic current Deterministic flux Mean current (uniform) Mean flux (uniform)
1.5
1.0
(66)
0
0.5
where c ¼ Ss =S and H(m;c) is the Chandrasekhar H function (Chandrasekhar, 1960). Thus Eq. (66) can be calculated to any desired degree of accuracy by quadrature. Eq. (66) is the transport theory equivalent of Eq. (59) when a ¼ N. The corresponding value for the rod model is from Eq. (59)
Jð0Þ ¼
585
0
0
1
2
3
4
5
x
1a : 1þa
(67)
Fig. 1. Mean flux and current for transmission of neutrons.
This, however, is not equal to Eq. (66). The stochastic average of any of the quantities defined by Eqs. (51)–(67) can be obtained by direct quadrature in the form
n f ¼
Z
1
n
dxf ðxÞ
(68)
0
where f is the quantity whose stochastic moment is needed. We shall be concerned with the mean value Cf D and the variance s2f ¼ Cf 2 D Cf D2 , which will be compared with the PCE equivalents. 5. Numerical examples and discussion We shall illustrate the above techniques by two simple examples, viz: a) Transmission through a slab b) Constant source in a slab
5.1. Data We now give the data employed in the calculations and discuss some of their properties. For examples a) and b), the absorption and scattering cross sections are required. Here we consider the data described by either a uniform or a log-uniform distribution, as in Eqs. (46)–(49), with the limits
cm1 ;
cm1
Sa min ¼ 0:01 Samax ¼ 0:05 Ss min ¼ 0:5 cm1 ; Ssmax ¼ 0:9 cm1
uniform 1 1 S ¼ ðS0 þ S1 Þ; s2S ¼ ðS1 S0 Þ2 2 12 loguniform ðS1 S0 Þ ðS S0 ÞððS1 þ S0 ÞlogðS1 =S0 Þ2ðS1 S0 ÞÞ S¼ ; s2 ¼ 1 logðS1 =S0 Þ S log2 ðS1 =S0 Þ (70) It is seen that the relative standard deviation arising from the log-uniform distribution is always greater than that from the uniform distribution. These fluctuations will be reflected in the solution of the associated transport equation as will be demonstrated below. A useful measure of the sensitivity of the system to parameter uncertainties is the ‘sensitivity coefficient’, defined by the ratio of the relative standard deviation of the parameter of interest to the relative standard deviation of the uncertain parameter. For example, if we wish to know how sensitive a parameter R is to a given uncertainty in a cross section S, then we must calculate
0.20
(69)
0.15
The statistical moments of these coefficients are given by the following expressions and their numerical values are shown in Table 1. Note that we have set S0 ¼ Smin and S1 ¼ Smax.
0.10
Relstd current (loguni) Relstd flux (loguni) Relstd current (uni) Relstd flux (uni)
Table 2 Reflection and transmission coefficients for the rod model in a slab of thickness 5 cm.
0.05 Type
R
RelstdR
SR
T
RelstdT
ST
U LU ULU
0.5652 0.5695 0.5758
0.02407 0.02717 0.02784
0.138 0.158 –
0.3034 0.3198 0.3136
0.1930 0.1845 0.1859
1.109 1.031 –
U ¼ uniform distribution for Ss and Sa ; LU ¼ log-uniform distribution for Ss and Sa . ULU ¼ uniform distribution for Ss and log-uniform pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi for Sa . Relative standard deviation ¼ Relstd ¼ variance=mean: SR and ST are the sensitivity coefficients of R and T to cross section uncertainty.
0
0
1
2
3
4
x Fig. 2. Relative standard deviation for transmission of neutrons.
5
586
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
Table 3 Polynomial chaos expansion coefficients for R and T log-uniform distribution, a ¼ 5.0 cm.
Table 5 Surface flux, relative standard deviation and sensitivity for the albedo problem. Uniform distribution (rod model).
n
Rn exact
Rn PCE
Rn exact
Tn PCE
a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0) Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0) S (a)
1 2 3 4 5 6 7 8 9 10
0.56947 2.51912E2 1.17658E2 1.18645E3 1.88353E5 1.87656E5 2.46500E6 5.98277E8 3.30194E8 6.24389E9
0.56947 2.51916E2 1.17658E2 1.18647E3 1.88255E5 1.87647E5 2.46050E6 5.97112E8 2.79565E8 6.03360E9
0.31981 0.10217 4.06331E3 2.24357E5 9.64129E5 2.28575E5 2.17992E6 4.64831E8 4.67760E8 7.27465E9
0.31981 0.10217 4.06367E3 2.25064E5 9.64043E5 2.28587E5 2.17771E6 4.64831E8 4.39209E8 7.77338E9
5 10 20 50 100 N
1.5652 1.6432 1.6680 1.6712 1.6712 1.6712
0.2972 0.1306 2.909E2 3.429E4 2.097E7 0
0.1930 0.3594 0.7166 1.6112 2.5899 –
SR ¼ ðsR =RÞ=ðsS =SÞ
(71)
This is a useful measure as it shows the fractional change in the property of interest due to a given fractional change in a particular cross section. 5.2. Transmission through a slab We compare in Table 2 the average reflection, R, and transmission T, coefficients for the slab using three types of probability distribution. The slab is of thickness 5 cm. In Fig. 1 we show the mean flux and current for this problem for uniform and log-uniform distributions. To be more specific we define the flux and the current as given by Eqs. (39) and (40). The results are compared with the corresponding deterministic values with average parameters. The difference between the two cases is barely noticeable graphically, which is not surprising since the range of values in (69) is not great. Of more interest is a measure of the fluctuations, viz: the relative standard deviation for the U and LU cases, which we show in Fig. 2. Here we note that the fluctuations in the flux decrease initially and then increase towards the other face. The current fluctuations increase monotonically as distance increases. The general tendency for uncertainty to increase with depth is not unexpected because the effect of variations in the cross sections compound with distance from the source. The results in Table 2 and Figs. 1 and 2 are from Eqs. (45a) and (45b) using 20 terms in the PCE. In fact, only 6 terms are needed to get 4 significant figure accuracy for the variance, but it is instructive to examine the relative magnitudes of the expansion coefficients. The accuracy of the results is obtained from the exact analytical results in Eqs. (59) and (60) in the form given by Eq. (68). As a measure of the accuracy of the convergence of the PCE series, we have shown in Table 3 the PCE coefficients for R and T for the loguniform distribution which is a more severe test than the uniform case. The exact expansion coefficients are obtained by representing Eqs. (59) and (60) in terms of polynomial chaos expansions. It is clear that convergence is very rapid, but it should be noted that the values given for the PCE entries were from calculations using 20 terms in the series, only 10 of which are shown here. In Table 4 we give some exact results for the surface flux and reflection coefficients for a half-space as defined by Eqs. (65)–(67).
Table 4 Exact transport theory values for a half-space surface flux and reflection coefficient. Type
U LU
Rod
Exact
Cf0 ð0ÞD
Relstd
R
Relstd
R
Relstd
1.6712 1.6938
2.0669E2 2.1760E2
0.6356 0.6597
5.7850E2 5.9629E2
0.6712 0.6938
5.1463E2 5.3123E2
1.5689 1.6432 1.6619 1.6629 1.6629 1.6629
8.691E3 7.712E3 1.825E2 2.065E2 2.067E2 2.067E2
0.3034 0.1392 3.805E2 1.667E3 2.181E5 0
0.050 1.109 0.044 2.066 0.105 4.118 0.119 9.260 0.119 14.88 0.119 –
The value of Cf0 ð0ÞD is the same for the rod and exact model. The rod model value of the reflection coefficient R underestimates it and overestimates the relative standard deviation. A feature of some interest is the variation of the surface fluxes, relative standard deviations and sensitivities as a function of slab width. In Tables 5–8 we give values for these quantities for the albedo problem with uniform and log-uniform distributions of the cross sections, calculated by the rod model and the computer code EVENT (DeOliveira, 1986), modified to include the pseudo-multigroup transfer coefficients. The columns headed f0;det refer to the cases where the average coefficients are used directly in the deterministic solutions; this is the crudest form of approximation. There are three matters for discussion regarding Tables 5–8. Firstly, the difference between the rod model and exact equation, secondly the general behaviour of the mean and variance and finally the accuracy of the results. There is a difference between the two transport models as would be expected. Both models are in reasonable agreement at the source face at x ¼ 0, but at the exit face there are significant differences, both in mean value and in variance. This is to be expected since uncertainties accumulate with distance and the rod model is clearly not effective in diffusing the neutrons around due to the backward-forward nature of the motion. The important measure of uncertainty is the sensitivity S. We note that this is relatively small at the entry face because the source is deterministic and the neutrons have not travelled deeply into the medium. However, at the exit face, we observe a substantial increase in sensitivity as the slab thickness increases. The difference between the two models for sensitivity is however not great and both predict large values differing at most by a factor of 2 but often by much less. To some extent the accuracy of the EVENT results may be assessed from the half-space values for the surface fluxes and variances because for this case we have the exact analytical results, as given by Eq. (65). We add the caveat here that the EVENT results for a half-space are not as precise as those for a finite slab. The reason being that the half-space requires a very thick slab and the consequent additional number of mesh points makes the resolution not as sharp as for the finite slab. We elaborate further on this matter below. In practice, for the half-space, the mean and mean square intensity are only found to 3 significant figures whereas the variance, which is the difference between two nearly equal numbers, is of order 104 for the data used, and thus will be inaccurate. This does not occur for the thinner slabs, i.e. less than 50 mean free paths. For this reason the results given in Table 9 show the mean and the mean square separately and their
Table 6 Surface flux, relative standard deviation and sensitivity for the albedo problem. Uniform distribution (EVENT). a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0) Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0)
S (a)
5 10 20 50
1.6343 1.6593 1.6631 1.6637
0.1491 4.083E2 3.303E3 3.405E6
7.701E2 0.1862 0.4193 1.1345
0.443 1.070 2.410 6.520
1.6355 1.6592 1.6608 1.6620
3.517E3 9.013E3 1.061E2 1.072E2
0.6972 0.4184 0.1412 7.994E3
0.020 0.052 0.061 0.062
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588 Table 7 Surface flux, relative standard deviation and sensitivity for the albedo problem. Loguniform distribution (rod model). a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0)
Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0) S (a)
5 10 20 50 100 N
1.5695 1.6572 1.6890 1.6938 1.6938 1.6938
0.3198 0.1553 4.789E2 2.646E3 3.974E5 0
0.3136 0.1465 3.777E2 7.096E4 9.460E7 0
0.1845 0.3309 0.6178 1.2430 1.8856 –
1.5733 1.6573 1.6821 1.6839 1.6839 1.6839
9.859E3 8.336E3 1.902E2 2.174E2 2.176E2 2.1760E2
0.055 1.028 0.046 1.843 0.106 3.442 0.121 7.144 0.121 10.50 0.121 –
associated errors are compared with the exact results. The errors are less than 0.5% and 1%, respectively, and show that the method is viable. The number of terms in the chaos expansion in EVENT is 5. Additional terms make little difference to the results. The mean values from EVENT are accurate to within about 0.5%. However the relative standard deviation has an error of about 50%. There are good reasons for this discrepancy which concern the structure of EVENT. The code uses a finite element, spherical harmonics representation of the angular flux based upon a variational formulation of the second order even parity form of the transport equation (Vladimirov, 1963). Linear finite elements are used in space with the resolution of the mesh tuned in order to reach spatial convergence without making the cells too thin optically. The caveat on not making the cells too small is because the even parity equation performs poorly for cases where the optical thickness of the element is small (Ackroyd et al., 1986). If therefore the mesh resolution is made too fine the simulations may fail to converge. A spherical harmonics expansion of order N ¼ 99 (Bell and Glasstone, 1970) was used to ensure angular convergence and double precision was also employed. The problem of accuracy lies therefore not only with the limitations just noted, but also because the point of interest is at the surface where the flux varies very steeply. In addition, the variance depends on the difference between two surface quantities and this leads to error accumulation. However, the purpose of this paper is not so much to obtain accuracy but to assess feasibility. It is clear that existing multi-group neutron transport codes can be used for this problem thereby obviating the need for significant development work. The accuracy will be much greater for parameters that do not depend on surface values; for example reaction rates in local regions.
5.3. Constant source in a slab
Table 8 Surface flux, relative standard deviation and sensitivity for the albedo problem. Loguniform distribution (EVENT). a
Cf0 ð0ÞD f0;det ð0Þ Relstd (0) Cf0 ðaÞD
f0;det ðaÞ
Relstd (a) S (0)
S (a)
5 10 20 50
1.6476 1.6792 1.6849 1.6858
0.1644 5.125E2 5.378E3 5.612E5
7.719E2 0.1795 0.3674 0.8061
0.430 1.000 2.047 4.491
4.726E3 1.038E2 1.206E2 1.222E2
0.7399 0.4795 0.1872 1.384E2
Table 9 Comparison of EVENT values at surface of half-space with exact analytic results. Prob distn
Mean (EVENT)
Mean (EXACT)
%
Meansquare (EVENT)
Relstd (EXACT)
%
U LU
1.6637 1.6869
1.6712 1.6938
0.45 0.41
2.7682 2.8460
2.7942 2.8704
0.93 0.85
Table 10 Mean flux, relative standard deviation and sensitivity coefficient of surface flux for constant source problem in a slab of thickness 5 cm. Type
Mean flux
Relstd
S
U LU ULU
2.2311 2.2753 2.2733
0.047849 0.046203 0.046380
0.275 0.258 –
In Table 11 we give the exact surface values for the flux arising from a constant source in an infinite half-space. In this case we note a significant difference for the mean intensity between uniform (U) and log-uniform (LU) distributions. The relative standard deviation however is not so much affected. A feature of some interest is the variation of the surface fluxes and relative standard deviations as a function of slab width. In Table 12 we give values for these quantities for the constant source problem. It is interesting to note from Table 12 that, for the source problem, the flux and variance depend strongly on the slab width; increasing as the width increases. For the 5 cm slab the Relstd is about 5%, whereas for the slab of infinite thickness it is about 35%. The sensitivity S also increases rapidly with slab thickness. 5.4. Some comments on the meaning of statistical moments Occasionally in the literature the significance of the mean and its variance are misunderstood. In this respect, and with relation to the above theory, it is important to note that the stochastic properties obtained from the above analysis do not refer to a specific system but to an ensemble of such systems. For example, one particular system will have one set of cross sections and will respond deterministically to those cross sections. Thus CfD and s2f give a measure of the uncertainty arising from a very large number (theoretically infinite) of systems whose cross sections vary
0.08
For completeness, we also give results for the case of a uniform constant source in a bare slab. Only the rod model is used but the results are instructive. For this case we show, in Table 10, the value of the flux at the surface and its relative standard deviation. Only small differences are noted. These may be seen more clearly in Fig. 3 for the relative standard deviation across the slab. Clearly, for the flux, it is symmetric about the centre and for the current anti-symmetric because the current changes sign at the mid-point.
1.6494 1.6785 1.6823 1.6823
587
0.026 0.058 0.067 0.068
0.04
0
Relstd Current (loguni) Relstd flux (loguni) Relstd current (uni) relstd flux (uni)
-0.04
-0.08
0
1
2
3
4
x Fig. 3. Relative standard deviation for constant source problem.
5
588
M. Eaton, M.M.R. Williams / Progress in Nuclear Energy 52 (2010) 580–588
Table 11 Exact transport theory values for a half-space surface flux.
Table 12 Surface flux, relative standard deviation and sensitivity for the source problem. Uniform distribution.
Type
Cf0 ð0ÞD
Relstd
S
U LU
6.2824 7.2045
0.3450 0.3391
1.983 1.894
according to the distribution p(S). To put it otherwise, the mean and average tell us only that any prediction about the system of interest will be uncertain by the amounts indicated by the stochastic moments; it is not known which of this myriad of ensemble members actually corresponds to our system. An illustrative example arises by considering the simple model equation, viz:
dfðt; xÞ ¼ aðxÞfðt; xÞ dt
(72)
from which we have
fðt; xÞ ¼ fð0Þexp½aðxÞt
(73)
a
Cf0 ð0ÞD
f0;det ð0Þ
Relstd (0)
S
5 10 20 50 100 N
2.2311 3.7922 5.3591 6.2244 6.2815 6.2824
2.2317 3.7703 5.1502 5.6127 5.6184 5.6184
0.04785 0.1139 0.2261 0.3373 0.3497 0.3500
0.275 0.655 1.300 1.939 2.010 2.011
into a form which can be solved using existing computer codes with only small modifications. Finally, we comment on the eigenvalue problem that arises when fission is included. In that case, the eigenvalue keff will also be a random variable and it is necessary to expand it in terms of a PCE. However, the resulting equations are non-linear and require special root-finding techniques that are beyond the scope of the present paper. However, in part II of this work we shall describe a new method which avoids these non-linearities (Williams, 2010).
The mean value of f is then
hfðt; xÞi ¼ fð0Þ
Z
b
dxpðxÞexp½aðxÞt
(74)
Acknowledgements
a
For a uniform distribution, as given by Eq. (47), we can evaluate the integral in (67) to get
hfðt; xÞi ¼
i fð0Þ h a0 t ea1 t e tða1 a0 Þ
Dr. Eaton would like to acknowledge the support of The Royal Academy of Engineering and EPSRC under their Research Fellowship Scheme.
(75)
In actual practice, we can only measure one system with specific values of coefficients determined by a value x ¼ xi. Thus the solution would be as in Eq. (73) with x ¼ xi. This of course varies markedly from Eq. (75) and stresses the point made above. 6. Summary and conclusions The purpose of this work has been to develop a method for assessing the influence of data uncertainty on some important parameters in neutron transport calculations. We have concentrated on surface quantities such as reflection and transmission coefficients as well as surface fluxes for some simple, one-dimensional problems. The essence of the method is the representation of the data by random variables, the statistics of which reflect the nature of the uncertainty. For example, a piece of data may be known between upper and lower limits with a certain probability distribution. In this case we have used uniform and log-uniform distributions but we could equally well have used any other distribution. The method is particularly useful when the data uncertainty is large and the usual techniques of perturbation theory fail. It also avoids the often unphysical representation of errors by Gaussian variables which, for large uncertainty, can lead to negative values for a physical quantity. The use of polynomial chaos is adopted which is a well-proven technique in structural mechanics and fluid dynamics, but has only recently found its way into neutron transport and radiative transfer calculations (Emery, 2005; Williams, 2006, 2007). It may well be, however, that this paper is the first to use it for assessing data uncertainty in neutron transport and to cast the resulting equations, which may be quite complex,
References Ackroyd, R.T., Issa, J.G., Riyait, N.S., 1986. Treatment of voids in finite element transport methods. Prog. Nucl. Energy 18, 85. Bell, G.I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand, Reinhold Company. Cacuci, D.G., 2003. Sensitivity and Uncertainty Analysis. Chapman and Hall. Carter, L.L., Cashwell, E.D., 1975. Particle Transport with the Monte Carlo Method. USAEC. Report. Chandrasekhar, S., 1960. Radiative Transfer. Dover Publications. DeOliveira, C.R.E., 1986. An arbitrary geometry finite element method for multigroup neutron transport with anisotropic scattering. Prog. Nucl. Energy 18, 227. Davison, B., 1957. Neutron Transport Theory. Oxford University Press. Emery, A.F., 2005. Some thoughts on solving the radiative transfer equation in stochastic media using polynomial chaos and wick products as applied to radiative equilibrium. J. Quant. Spectrosc. Radiat. Transf. 93, 61. Ghanem, R.G., Spanos, P.D., 1991. Stochastic Finite Elements. Dover Publications. Ghanem, R.G., 1998. Scales of fluctuation and the propagation of uncertainty in random porous media. Water Resour. Res. 34, 2123. Li, R., Ghanem, R., 1998. Adaptive polynomial chaos expansions applied to statistics of extremes in non-linear random vibration. Prob. Eng. Mech. 13, 125. Le Maitre, O.P., Knio, O.M., Najm, H.N., Ghanem, R.G., 2001. A stochastic projection method for fluid flow. J. Comp. Phys. 173, 481. Vladimirov, V.S., 1963. Mathematical Problems in the One-Velocity Theory of Particle Transport. Atomic Energy of Canada Ltd. Report AECL 1661. Wiener, N., 1938. The homogeneous chaos. Am. J. Math. 60, 897. Williams, M.L., 1986. Perturbation theory for nuclear reactor analysis. In: Ronen, Y. (Ed.), Handbook of Nuclear Reactor Calculations, vol. III. CRC Press. Williams, M.M.R., 2010. A method for solving stochastic eigenvalue problems. Appl. Math. and Comput. 215, 3906. Williams, M.M.R., 2006. Polynomial chaos functions and stochastic differential equations. Ann. Nucl. Energy 33, 774. Williams, M.M.R., 2007. Polynomial chaos functions and neutron diffusion. Nucl. Sci. Eng. 155, 109. Xiu, D., Karniadakis, G.E., 2002a. Modelling uncertainty in steady state diffusion problems via generalized polynomial chaos. Comput. Methods Appl. Mech. Eng. 191, 4927. Xiu, D., Karniadakis, G.E., 2002b. The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM J. Sci. Comput. 24, 619.