Numerical solution to the problem of criticality by Monte Carlo method

}'rll;Ic,l

Ill ( ; ~

tI B r l l , l l l ~

\ll rl~rl;- r~crtmd

:,';

{'zr_',Ht~,,n {~rc~, i~lc

N U M E R I C A L S O L U T I O N TO THE P R O B L E M OF CRITICALITY BY MONTE C A R L O M E T H O D JAN

K'~ N C L

Nuclear Research Institute, 250 68 t~.e2, Czechoslovakia

Abstract. In this contribution, a new method of numerical solution to the c r i t i c a l i t y problem for neutron transport equation is suggested. This method is based on results of the theory of p o s i t i v e operators developed by Krein and Rutman. The problem is solved by the Monte Carlo method and the random process is c o n s t r u c t e d in such a way that d i f f e r e n c e s between results obtained and the exact ones would be arbitrarily small. The method can a d v a n t a g e o u s l y be applied in case of n o n a n a l o g random processes. Numerical results of the method are also presented.

Introduction. Let us consider a spatial region 0 and study the c r i t i c a l i t y problem for the transport equation. Ue wish to find a real number ~ > 0 and n o n n e g a t i v e nontrivial function ~ ( m ( O ) such that the equation -4

Ay holds. on t h e

x,E,~

~2 H (~) 1 n

Zi

+ i

Here, m(O) is a linear s e t H m D x E 3 and

o ~

where

= B~

the

aq

following

~=4

notations

C~

(t)

space

of complex

functions

are

used=

... the number of different nuclei, ... m i c r o s c o p i c total cross-section, ... m i c r o s c o p i c fission c r o s s - s e c t i o n , ... the number of neutrons created by fission, fission

spectrum,

are defined

v

... coordinates of location,energy and d i r e c t i o n , ... surface of the unit sphere, ... the number of nuclei of material i per unit

...

which

respectively, volume,

244

~i

J. K', ~ ( L

... m i c r o s c o p i c Problem

differential

41) can be rewritten

= A"B V

scattering

into the form

'4"4c f

*

(2)

T h i s i s an i n t e g r a l e q u a t i o n and i t s r e g i o n o f For the sake o f b r e v i t y , we d e n o t e t h e i n t e g r a l operators A-iB and A-1C by

Ks(#'E'~ ,

cross-section.

,EI,~-~)

and

integration i s t h e s e t N. k e r n e l s c o r r e s p o n d i n g to t h e

Kf(~,E,~,

,Et,~ ,

),

respectively. Equation (2) is f r e q u e n t l y used as a s t a r t i n g p o i n t to numerical s o l u t i o n of the c r i t i c a l i t y problem by Monte Carlo method. Usually, i t is proceeded by the method of successive generations ( some d e t a i l s can be found in the paper by Frank-Kamenietzky (1978)). U n f o r t u n a t e l y , such a computation leads to a bias. I t has been shown that the e r r o r depends on the number of p a r t i c l e s present i n the generation and i t does not decrease with the number of generations used ( see e.g. book of Zolotukhin and Mayorov (1984)). Here, we suggest another way. Using i t we obtain numerical r e s u l t s the bias of which can be made a r b i t r a r i l y sma11. The method can also be applied in nonanalog processes.

Random process and the s o l u t i o n of the problem. If the spatial r e g i o n O i s bounded and c o n v e x , the l i n e a r space m(O) can be chosen as a subspace o f t h e Lebesgue space L 1 and i n such a way t h a t t h e following c o n d i t i o n s are s a t i s f i e d : i)

i£) iii)

~(~,E,~

) = O ,

x~

70,

n~

~

0

f o r any ~m(O) ( ~ means e x t e r n a l n o r m a l t o t h e b o u n d a r y r e g i o n O). The o p e r a t o r s A-1B, ( I - A-1B) -1 and A-1C are bounded. The o p e r a t o r the spectral

(Shikhov t h e form

(1973)). ~=

30 o f

the

T = ( I - A-1B) - 1 A - 1 C i s p o s i t i v e T2 b e i n g compact and radius r(T) of the operator T is positive. Under t h e s e c o n d i t i o n s ,

T~

,

~

problem

42) can be r e w r i t t e n

re(D)

into (3)

and, by t h e t h e o r y of K r e i n and Rutman ( 1 9 4 8 ) , the f o l l o w i n g assertion is true: There i s one and o n l y one e i g e n f u n c t i o n ~ , f l ~ E = l t o t h e p r o b l e m 43) ( # . # means t h e norm i n t h e space m ( O ) ) . T h i s s o l u t i o n is nonnegative, the e i g e n v a l u e ~o c o r r e s p o n d i n g to i t b e i n g p o s i t i v e and~o=r(T). N e x t , f o r any function ~ m(O) and any p o s i t i v e i n t e g e r n, t h e r e l a t i o n h o l d s . Here, ~o i s a p o s i t i v e p r o b l e m ( ] ) and T 1 i s a l i n e a r r(T I )

<

eigensolution t o the p r o b l e m o p e r a t o r such t h a t v~,.

adjoint

to

the

(5)

The criticalit', problem

245

Consider a particle in the following process ot random collissiuns in the phase space M: The particle starts its hLstury with pcsitiun and velocity randomly chosen. It colleads with elements of a mediun and it is either scattered ( i.e. its velocity is changed) or absorbed and then possibly reproduced. Velocity of the particle between any two succeeding eollissions is constant. Suc~ a process

,

will be described

"',#m ) ,

/Si

=

Here, m~ 0 i s a g i v e n i n t e g e r means

~i

i are the c o l l i s s i o n

the first collission

~of

the particle,

( xil , . . . , X n . )i ,

'

j=l,2,...,n

Dy t~e track

is

i=l,2,...,m. I a s t a g e of the t r a c k

points

point

of

the p a r t i c l e .

6a) and x i ~ M , 3

Specifically,

l xI

in the stage

( if the particle with energy E leaves direction ~ , the point x ~ ( ~ , E , ~ )

i, x i the point of absorption nc the region D at the place ~e ~0 in the is also interpreted as the point of

absorption stage

) and x~j' j # n i denotes the site of a scattering collission. Any is terminated by absorption. The particle may pass into the next stage

only by means of ~eproduction. first collission

of particle

The point

x~+l- corresponds

to the site of the

reproduced.

The random walk is described

by the following

of the first collission

quantities:

Pl(X)

... density

probability,

p(x,y)

... density of the scattering probability (i.e. p(x,y)dy is the probability that, after collission at the point x, the particle will be scattered and its new collission will happen in the neighbourhood dy of the point y),

(6b)

Pa(X)

... absorption

q (x) r q(x,y)

... probability of reproduction of the particle after its absorption, . . . d e n s i t y of t r a n s i t i o n probability ( i . e . q ( x , y ) d y i s the probability that, after a reproduction at the point x, the first collission of the particle will happen in the neighbour~ hood dy of the point y).

probability,

From now on, we will suppose that any stage /~i ,i=l,...,m will be terminated after a finite number of collissions with probability i, i.e.

lim for

fdx I

...

Jdx n p(x,x I

)

~T p ( x j , x j , l ) =

any x ~ M. Here and i n t h e f o l l o w i n g , I I f(xi)--i

it

is

0

understood

(7) that

for j~,k. f

Simultaneously with process ~-~ ~'-4

where ~'(×) = { S ( x l l P l ( X ) 0

if Pl(X)

= O,

(6), we will consider

the random ~-4

variable

~m' ~

246

J. K~ >,~ t

-Ks(Y,X)/p(x,Y) 1

u(x,y)

=lO

f

v(x)

if

p(x,y) = O,

( Pa(X)qr(x))-i

=[0

if Pa(X)qr(x) = O,

Kf(y,x)/q(x,y) w(x,y) =~0 i f q(x,y) = O, ~(x) : ( g ( x ) / p a ( x ) ~0 if Pa(X) : O, S is an element of the space m(M) and g a real function such that I:dxg(x) ~ ( x ) [ ~

e.~

for any

:~

m(O).

The quantities pl,p,pa,qr and q are supposed to be chosen in such a way that the following implications hold:

S(x) /

0

=> pl(x) / o,

Ks(Y,X) i 0 --)p(x,y) f O, Kf(y,x) / 0 ~ Pa(X)qr(x)q(x,y) # 0 g(x) 4 0 :=~ Pa(X) # O.

and

Theorem i. For any integer m ~i, the expectation M f of the random variable ~m is equal to the integral m

Jdx gfx) ~ (x)

M where

~ = Tm-I ( I - A-IB)s.

(9)

Proof. The following relation holds in general,

where the summation is taken over all tracks corresponding to a given m, P( ~m ) is the probability of the occurence of the track ~m and ~(°~'m) =~m'o Especially, using definitions (6) and (8) ~-4 ~-~

Z

M,f ,,,,,

•,v.=~

i,

5

,~,=4

,;

S

x, . .. f d d ... /,~,~. . .f d x ~

,W

~

,;

..;

"

t"1

/'1

"

A, (,,j) / / ~.=4

77~'=1

Thc criticatit', problem

: Z . . Z . . ],z~,'.. . . . .f,,,,, " .r,~Cr..

24.7

/,.2 sc,:)n 77 "-'~c,~

~.= "I

~/)

"I

#'.'l

Q.E.O.

For a given m, let us introduce

the random variable

~m' (

? : s.:)/ ,,,~-,t/

~ ,~(~.,~ X~:.).~ ;..," .,C';Z 7-;,,c,T,~7.,) .)

~

.

Proceeding similarly to the proof of Theorem find that the following theorem is true:

Theorem 2.

For any i n t e g e r

where

~

i and using condition

c,o~

(7), we

m >_l,

i s the f u n c t i o n

g i v e n by e x p r e s s i o n

(9).

Note t h a t Theorems 1 and 2 are analogous to the w e l l known theorems s t a t e d e . g . i n the book by S p a n i e r and O e l b a r d (1972) f o r the case when the t r a c k c o n s i s t s of a s i n g l e s t a g e . Now, we w i l l s i m u l a t e trials. D u r i n g any t r i a l , Theorem

3. For any ~ > 0

>"

lira

N'."~~

t h e random p r o c e s s ( 6 ) by rl m u t u a l l y i n d e p e n d e n t we w i l l r e c o r d the v a l u e s ~ ¢ o r ~ , i = l . . . . . N.

~"- J,,,,q(,,~ (,,)t=.~,) = I,~ ed-~ Z

-

,,;}cx)4(x)l>~

=0

(ll)

I,=I where P is the probability of the corresponding is the function given by expression (9).

event and

Proof. The distribution function of the variable j~,~ ( ~ ) is the same for all i=l,2,...,N. According to Theorems i and 2, the following relations hold

IN jim I

: I l l ~im]

:

I ;dx g ( x ) ¢ ( x ) I ~

~

so t h a t a s s e r t i o n ( l l ) i s a consequence of i n e q u a l i t y K c h i n c h i n ' s theorem ( see Gnedenko ( 1 9 6 9 ) ) .

(12) (12)

and of

Besides the case of bounded s p a t i a l r e g i o n O, the r e s u l t s above can a l s o be extended to i n f i n i t e b u t p e r i o d i c a l r e g i o n : L e t t h e m a t e r i a l medium be spread o v e r an i n f i n i t e s p a t i a l r e g i o n ~ ( E~ ( s u r r o u n d e d p o s s i b l y b ~ vacuum) and assume t h a t t h e r e e x i s t a bounded and £ ~ l a t i v e l y compact s e t D C ~ and t h r e e bounded v e c t o r s ~ i ' ~2 and ~3 such t h a t

248

J. K'~ NCL

Ni(~) = N i ( ~ k ~ j ) j=1,2,3;

O~,~ # O Z k b l ¢ l

;

i=l,2,...,n~

2:m~3~

,

x~O,

(i~)

k,l,m=O,l,2,...

In this case, such functions ~ ~ m(O) will play substantial role which have the property (13) with respect to their spatial variables. Then, operators A-18 and A-Ic are expressed as infinite sums ( convergent in norm) of integral operators having properties ii) and iii) and are such that the domain of integration corresponding to each o~ them is contained in M.From this it follows that relations (&) and (5) will remain true. In the next considerations, the set 0 will be either a convex and bounded region at boundary o~ which condition i) is satisfied or a subregion o~ an unbounded region ~ where the properties (13) are supposed. Theorems 1,2 and 3 enable us to solve the criticality problem (~) numerically. Indeed, using relations (~) and (9), we get ~-~"

~'~

'

(i~)

M

~:

( I - A-1B)S

where ~ and ~oare the e i g e n f u n c t i o n

and eigenvalue to the problem (3).

So

the quantities ~ a n d ?~ , i:l,2,...,rl serve as unbiased estimates of the integral expression (14). Clearly, from Eq.(14) and inequality (5) we get ~. =,~lim

(M #m

)""=

(15)

~--lim ( M # m ) - - - '

suppose that thevarianoe D(#m ) oC the random variable

#m satisfies

the inequality

0 < O(~m ) < (Random v a r i a b l e

(8) i s o b v i o u s l y bounded i f ,

for some f i n i t e

constants CI,

C 2 and C3, the inequalities S(x) ~ClP1(X),

g(x) ~C2Pa(X) ,

and ~'~ ( Y ) K s ( Y , X ) / 4 ( x )

~

~y)Kf(y,x)/~(x)

~ C3Pa(X)qr(X)q(x,y)

p(x,y)

hold on the set M. Here, we have put

Obviously,

in such a case

O(~ m ) ~ ) .

Then, by Lyapunov's

theorem

( Gnedenko (1959), §42).So, the standard e r r o r ~ ( i . e . x = l ) in computation of e x p e c t a t i o n M ~m by Monte Carlo method obeys the law

where C ~ i s a c o n s t a n t . Then, c l e a r l y q u a n t i t y ( M ~ . . ) ~ i s accuracy Em estimate of which i s given by expression .f

determined w i t h

The criticality problem

24-<)

Of course, case m ~ is not realizable in practice and the use of expression (15) depends on the value of d o m i n a n c e ratio r ( T ~ ) / ~ o . It is necessary to ensure that the number m of the realized stages be sufficiently great which is a matter of suitable choice of the functions pl,p~,p,q and q.

Applications. On t h e b a s i s o$ T h e o r e m s 1 - 3, t h e c o m p u t i n g c o d e MOCA2 has been developed. It is aimed to the numerical solution to the criticality probiem in multigroup and s p a t i a l l y three - dimensional formulation f o r VVER r e a c t o r fuel assemblies. In this code,the value of the multiplication factor k is computed according to the formula eff 4 • 4 6=4

~'=4

where N is the number of that t ~ m o ~ Z m.

trials

and m i s o

•

a suitably

chosen

The relative mean square error ~ k e f f of the quantity from the expression kef

=

•

Knowing keff,

such

kef f is calculated

-

where

are computed

integer

"

the neutron

from the vaiues

Z / flux ~

~

4

and its relative mean square error ~

according

to the formulas

P and

U n d e s i r a b l e leakage of particles in the random process is suppressed in the code. For example, any collission of a particle in a n o n f i s s i o n a b l e medium i is c o n s i d e r e d as a scattering and the weight of the particle is m u l t i p l i e d by the factor

Accuracy of the results can be i m p r o v e d u s i n g t w o o r m o r e s t e p s o f c o m p u tation: Resulting particle distribution obtained i n one s t e p s e r v e s as a starting point in the next step of calculation. Two t e s t p r o b l e m s w e r e c o m p u t e d by t h e c o d e NOCA2 and t h e r e s u l t s are shown i n t h e f o l l o w i n g table. Problem 1 (a) - d)) corresponds to the case of f u e l a s s e m b l y VVER-IOOO p l a c e d i n i n f i n i t e regular hexagonal lattice of the same a s s e m b l i e s . The h e i g h t o f t h e a s s e m b l y i s c o n s i d e r e d as i n f i n i t e and the presence of absorbing rods is supposed. Formulation of the second problem i s t h e same b u t t h e a b s o r b i n g r o d s a r e r e p l a c e d by w a t e r . The c a l c u l a t i o n s were performed in two steps with four energy groups. In both steps,we have

!50

J [<.~,",4,_'L

put m :

I 0 0 and m = 5.

o

For comparison, the same problems were computed by the method of successive generations ( code MOCA). In this casepwe have considered particles per one generation, 5 initial generations were omitted.

i00

Table Results of two test problems ( variants c) and d) correspond to the second step of computation).

MOCA2

Vat.

N.part. 600

la) tO)

i

ld) 2a)

Err.(~)

0.9585

0.19

m

0.ii

400

0.9650

0.II

900

0.9645

0.i0

1.1547

0.16

ii00

lc)

kef f

MOCA

2b)

800

1.1533

0.15

2c)

400

1.1528

0.15

2d)

II00

1.1522

0.12

N.part.

kef f

Err.(~)

124500

0.9517

0.51

177400

1.1454

0.24

References. Frank - Kamenietzky A.O. (1978) Modelling calculation by Monte Carlo method. Series Moscow (Russian). Gnedenko

8.V.

(1969)

of neutron tracks in reactor " Nuclear reactor physics" 8,

The theory of probability.

Moscow

(Russian).

Krein M.G. and Rutman M.A. (1948) Linear operators leaving in a Banach space. Usp.mat. Nauk III, N.1, 3 -95 (Russian). Spanier J. and Gelbard E.M. (1972) Monte Carlo principles transport problems. Moscow (Russian). Shikhov S.V. (Russian). Zolotukhin parameters

(1973) Problems

of mathematical

invariant

a cone

and neutron

theory of reactors.

Moscow

V.O. and Mayorov L.V. (1984) Estimate of reactor criticality by Monte Carlo method. Moscow (Russian).