- Email: [email protected]
An exact transient stochastic solution for low-power neutron multiplication
Ann. nucL Energy, Vol. 12, No. 2, pp. 65 80, 1985 Printed in Great Britain. All rights reserved
0306-4549/8553.00+0.00 Copyright © 1985 Pergamon Press Ltd
AN EXACT TRANSIENT STOCHASTIC SOLUTION FOR LOW-POWER NEUTRON MULTIPLICATION G. T. PARKS* and J. D. LEWINS Engineering Laboratories, University of Cambridge, Trumpington Street, Cambridge CB2 1PZ, England (Received 15 August 1984)
Abstract--An exact solution for arbitrary physical data of the neutron-chain multiplication problem at low power is sought for the transient behaviour of the prompt-neutron stochastic chain to complementthe known fully-developedsolutions.Three means of inverting the equation for the generating function in order to find it explicitlyare investigatedand found to be unsatisfactory.A method of findingindividual neutron probabilities as implicitfunctions of time in a source-freeassemblyis demonstrated and typical results, for sub- and supercritical cases, are evaluated and comparison made with the 'quadratic approximation' model.
I. I N T R O D U C T I O N
There has been considerable work in the stochastic theory of multiplying chains since the original studies by Galton and Watson in the field of social anthropology [see reviews in Harris (1963) and Williams (1974)]. Solutions, both exact and approximate, have been found to a wide range of specific problems. In particular, it is well-known that an exact solution exists for the case of fission producing precisely two neutrons (e.g. Kendall, 1948), variously known as the 'linear birth and death', 'binary fission' or 'escalator' problem (Bartlett, 1979), and also for the form of the problem after truncation of the exact fission distribution at second order (Bell, 1963). Kendall demonstrated, though considering a generalized problem with time-varying coefficients, that 'no exact solution' exists except for the 'linear birth and death' case. This result is clearly applicable to the problems with constant coefficients studied here, but in what sense is there 'no solution'? This must mean 'no explicit closed solution', for there is no doubt that, in some senses, the general stochastic problem is well posed. In particular, as the finite moments have closed equations, these are in principle always available. In addition, Williams (1979) has given exact solutions for several other mathematical models of fission-neutron distributions through use of Lagrange expansions, concentrating predominantly on finding extinction probabilities. Other work in the field has considered the prompt-jump approximation and numerical solutions (Harris, 1964; Smith, 1970), and, most recently, Campbell et al. (1981) and Lewins (1981) found solutions for fully-developed subcritical assemblies including a source and for super-critical source-free reactors, respectively, both based on arbitrary physical data. As yet, however, there has been little consideration of the transient behaviour of systems, particularly for models based on empirical data. It is the form of these transients leading to the developed solutions for empirical data which is investigated here. Such solutions provide standards against which the accuracy of the various approximate models may be checked. We also seek to help develop understandingof the problem and consequently this paper is written with this didactic purpose in mind--we show some of the methods that did not work. 2. M O D E L S AND E Q U A T I O N S
The foward stochastic equations of the discrete model may be written (Lewins, 1978) in terms of the generating function F(x, t) in the lumped, low-power case as
r
OF O(X)-- x_ + 0 ~ - = k zf ~-
~x
+S[f~(x)--l]F,
* Present address : GEC Electrical Projects Ltd, Rugby, Warwickshire, England. 65
(1)
66
G . T . PARKS and J. D. LEWINS
where 1
is the mean rate for a neutron to cause fission,
Tf
1
S
L(x)
is the mean rate for a neutron to be captured, is the mean rate of source events, is the associated auxiliary generating function for such source neutrons
and
g(x)
is the auxiliary generating function for induced fission ;
S, zf and z¢ are here assumed to be constants independent of time. The main generating function F(x, t) is defined by
F(x,t) = ~ P,(t)x";
F(x,O) = x N,
(2)
n=O
where P,(t) is the probability that the reactor contains exactly n neutrons at time t and N is the initial number of neutrons at time t = 0. If F(x, t) can be determined analytically, two useful types of information are available : the ruth probability
O'F
= m! P.,(t)
(3)
= (n(n-- 1)... (n-- m + 1)>.
(4)
G~xm x=O
and the mth factorial moment
63xm x = 1
The auxiliary generating functions are defined by
g ( x ) = ~ q~x ~ and
f~(x)= ~, p~x ~,
v=O
(5)
v=O
where qv
is the probability of exactly v neutrons being released in an induced fission event
and Pv is the probability of exactly v neutrons being released in a source event. For simplicity it is generally assumed thatf~(x) = x, i.e. that each source event results in the release of exactly one neutron. This convention will be followed throughout this paper. If, in addition, it is assumed that each induced fission event causes precisely two neutrons to be released, i.e. 9(x) = x 2, the general stochastic equation reduces to OF
~t
-
k
X2 - X zf
OF +
~
-
7x
+ S(x-1)F.
(6)
This is the form of the problem known as the 'linear birth and death' approximation. The authors find this and the various other titles for this model ('escalator' and 'binary fission') either misleading or ambiguous. As this case considers singlet independent sources and doublet fissions, we suggest that in the context of neutron stochastics this equation be known as the 'doublet fission' model. Whatever its name, its solution is readily found. The independence of neutron behaviour in a low-power reactor means that the solutions for differing initial conditions may be found by superposition from two fundamental solutions : Finn = F ~ G ' ~ ,
where m n
G1
is the number of neutrons initially present, is the number of independent sources [assumed to have the same auxiliary function fs(X)], is the solution for rn = l, n = 0 (i.e. for no source, one initial neutron)
(7)
An exact solution for low-power neutron multiplication
67
and Fo
is the solution for m = O, n = 1 (i.e. for no initial neutrons, one source).
For the 'doublet fission' model /~(x- 1 ) T - ( 2 x - / ~ )
(8)
and
2-/~
(9)
Fo = [2 T _ I~_~( T _ I )x ] s/a, where 2=
1
,
#=--
%'f
1
T=exp(A-#)t.
and
(10)
2"c
If, alternatively, the following identifications are made J,=lg(2)(1)/Tf
1
I~=--+[1--gm(1)+½gtZ~(1)]/rf,
and
(11)
Z"c
where
d"9
9(")(1) = dx" x=l
=,
(12)
the resulting solutions from equations (8) and (9) are exact to the second moment ofn. As this change is equivalent to defining
g(x) =
1 + g m ( 1 ) ( x - 1)+½g(2)(1)(x - 1)2
(13)
this case is known as the 'quadratic approximation'. It is important to note that the single-source solution Fo may be shown (e.g. Bartlett, 1979) to be related to the single initial neutron solution G 1 by the convolution
F o = e x p { f l S [ G , ( x , t - - z ) - l ] d r }.
(14)
Thus only one of the fundamental solutions needs to be found directly from the general partial differential stochastic equation.
3. GENERAL NO-SOURCE CASE--SOLUTION BY LAGRANGE'S METHOD With no source present the general stochastic equation (1) reduces to
[
~F _ _ 0 ( x ) - x
~3F + 1--x I ~-x
The subsidiary generating function 9(x) is determined by physical data appropriate to the reactor fuel etc. (Diven et al., 1956). Typically it is found that no more than five neutrons are released in a fission reaction, so 5
5 g~V)(1)(1 _x)V
O(x) = ~" qvx~= ~ v=O
v!
'
(16)
q,x v -- + ~f 3¢ J Ox'
(17)
v=0
where 9~)(1) is defined by equation (12). Substituting for g(x) in equation (15) gives
= ~
v=o
68
G . T . PARKS and J. D. LEWINS
therefore OF Ot
q5 rr ( x - a o ) ( x - a l ) ( x - a 2 ) ( x - a 3 ) ( x - a 4 )
OF 0x'
(18)
where ao, a l , . . . , a4 are the roots of the polynomial equation 1
x 1-x q,x v - - + - = O .
5
-- ~ "L'f v = O
"t'f
(19)
27c
x = 1 is necessarily a solution of equation (19) and it is found that there are two further real roots and a complex conjugate pair. Hence OF Ot
q5 ( x - 1 ) ( x - a O ( x - a 2 ) ( x - a 3 ) ( x - a * zf
0F ) ~-x = 0.
(20)
For a single initial neutron F - G~, so applying Lagrange's method dt
- z f dx
1
dGa -
qs(x--1)(x--at)(x--a2)(x--a3)(x--a~)
0
(21)
therefore G1 = 4,
aconstant.
(22)
Also -- dt = zf d x / q s ( x - 1) (x - a 0 ( x - a2) (x - a3) (x - a*)
(23)
therefore dt=-dx
A ~ +
B -x--a
C
+ 1
X--a
+ 2
D x--a 3
+
D* ] --~- , x--a 3
(24)
where A, B, C and D are obtained by the usual partial fractions techniques. Hence t = ~k--[A In ( x - I ) + B In ( x - a O + C
In ( x - a 2 ) + D In ( x - a 3 ) + D * In ( x - a~')],
(25)
where ~ is a constant. Now ~b =f(~k), an arbitrary function of ~, therefore In ( x - 1)+ ... + D * In (x-a*)).
G1 = f ( t + A
Whent=0,
(26)
Gt = x a n d x = f ( A In (x-- 1) + ... + D* In (x--a*)) ;
(27)
equation (27) is of the form x = f(u(x)), where u(x) = A In ( x - 1 ) + B In ( x - a O + C
In ( x - a 2 ) + D In ( x - a 3 ) + D * In ( x - a * ) .
(28)
Given this it should be possible to model f ( ) by an appropriate series, or other, representation. Then (29)
G I = fm(t + U(X)),
where fro is the model off. Equation (29) would then allow Gx to be evaluated or manipulated as a function of t and x as required. We therefore seek to invert u(x) as x =f(u).
4. U N S U C C E S S F U L
METHODS
T o h e l p v i s u a l i z e t h e p r o b l e m u ( x ) m a y b e f o u n d f o r t y p i c a l p h y s i c a l d a t a ( T a b l e 1).Ifzf/z c = 1.4(corresponding to a sub-critical system), the roots of equation (19) are as follows: ao = 1,
al -~ 0.970,
a 2 -~ - 3 . 5 3 8
and
a 3 -- - 0 . 9 6 8 + 3 . 1 6 3 j ,
An
exact solution for low-power neutron multiplication
69
Table l. Sample 23 ~U fission data (after Diven et al., 1956) qo = 0.027 ql = 0.158 q2 = 0.339
q3 = 0.305 q4 = 0.133 q5 = 0.038
whence ~'f
A "" 0.5293 - - , qs
Obviously only the real part of
~'f
Tf
B -- --0.5374 - - , q5
"/~f
C -~ 0.0029 - - , qs
u(x) is of interest. Manipulation
D ~- (0.0026-0.0012j) - - . q5
of equation (28) gives
~ . (u(x)) = A In (1 - x ) + B In ( x - a O + C In ( x - a 2 ) + 2 R In ( x + 1 ) + R In (1 + z 2 ) + 2 I tan-1 z,
(30)
where
R = ~ e (D), I = J ~ (O) and
z = J~
(31)
(a3)/(x +
1).
Substitution into equation (30) gives the results in Table 2. This allows x to be plotted as a function of u, as in Fig. 1. Data for a typical super-critical case is given in Table 6. Three different forms of representation are investigated in the search for a suitable model for f: (a) Lagrange interpolation (b) Lagrange expansion (c) McLaurin series.
(a) Lagrange interpolation Lagrange's interpolation formula for a function may be written as
fro(u) = ~ lff(u,), o
where
l,
= fi u-u;. j=o Ul--Uj
(32)
This effectively gives a polynomial of order n in u which matchesf(u) at n + 1 points (u~, f(u3). Clearly, ifthe function is matched at a large number of points, the accuracy of the representation will be very good. It is also easy to find
Table 2. u(x) for "rf = 1.4
u(x)
ANE 12:2-C
x
( ×lO-3zf~qs/
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
23.19 26.10 29.57 33.85 39.32 46.68 57.33 74.58 108.77 218.92
70
G.T. PARKSand J. D. LEWINS 10 ~..-.-.0,-
09 08 O7 0.6 "~05 04 0.3 0.2 01
0
25
50
75
100
I 150
125
I
I
I
175
200
225
t./ (X 10-3"i'f /q5 )
Fig. 1. Plot ofx as a function ofu for -- = 1.4. Tc
Gl(x, t) from this model:
fi Uw-u,i
(33)
i=0 j=0
where
w = t + u(x).
(34)
Although it is relatively simple to write a computer program to generate Gl(X, t) from this expression, it must be remarked that information is contained not only in G1 but also in its derivatives with respect to x. The analogous expression for gG1/gx is
dG 1 = dx
ifi
__1
ifi
( w - u,)f(ui) ~x"
(35)
~=oj=o ui-- uj k=O t=O =~i
~:i ~=i~=k
Derivatives of considerably higher order than this are of interest and it is apparent that their expressions in Lagrange interpolation notation will be correspondingly complicated. Although a formal differentiating program might be possible, the difficulties in assimilating these expressions in a program are such that the Lagrange interpolation representation only seems viable in the absence of a suitable or simpler alternative.
(b) Lagrange expansion An initially more promising means of representing fasf= would appear to be through a Lagrange expansion, as suggested by Gurland (1978). Various forms of Lagrange expansion have been used in the literature in this connection although the correspondence is not always apparent (Williams, 1979). Gurland proposes that the expansion in the form forf(u(x)) becomes
x = 1 + ~ bi ui
around
u = 0,
i=1
acknowledging that the practicability of the method to find the bi depends on the nature of the previously defined coefficients 'A' of the partial fraction expression, without actually finding the b~. However, for x to be a proper stochastic variable we require 0 ~< x ~< 1 and Gurland's form is unattractive; there
An exact solution for low-power neutron multiplication
71
is no reason to suppose the bi are such as to secure this limit to x. If instead we consider
i.(u): xo+ k=l L
k~
\U--Uo/d .....
L
(36)
,
where D" = d"/dx" and taking x o = O, this reduces to
i
k=l
(u-,o),
k!
\ U - U o / 3x=O
k
= i k=l
k!
(37)
'
where ,, =
,
L
(38)
\U-Uo/J~=o
whence (W-Uo) ~
Gl(x,t)= ~ h k - -
k=l
(39)
k!
N o w consider the first coefficient of this series, h 1,
therefore
hl~~ A
x
(41)
ln(l_x)+Bln(~X)+Cln(~_2X)+Dln(~3X)+D,
ln(a*--x']]
\~-/jx=o
"
It is difficult to state precisely what value this gives for ha, and similar problems are encountered for all the other coefficients hk. It is thus found that, even though the conditions of applicability for this expansion are not contravened, it cannot be applied to this particular problem. This is not so surprising as "it is unusual for the general coefficient in the inverted series to be expressible in a convenient form" (Jeffreys and Jeffreys, 1956). We are aware of the successful use of the Lagrange expansion in obtaining the extinction probability Po(t) (e.g. Williams, 1979) and this remains a technique worth further investigation.
(c) McLaurin series One alternative method of expressing a function as a polynomial series is by a McLaurin series representation :
fro(U) = f(Uo) + (U-- UO) dd~ . . . . -t
(U--Uo) z d2f 2!
du2 . . . . +
(U-Uo) 3 d y . . . . 3! du 3 + ..-.
(42)
Now df du
d f dx dx du
and
df
dx
f(u) = x,
(43)
thus
dR du" Differentiating equation (28) :
db dx
A B C D D* + + + + x - 1 x - a 1 x - - a 2 x - - a 3 x--a*
(44)
72
G . T . PARKS and J. D. LEWINS
SO
(45)
du __ zf ( x - - 1 ) ( x - - a O ( x - - a 2 ) ( x - - a 3 ) ( x - - a ' ~ ) dx q5
and d x _ 1 (qo + q 1x + q2 x2 + q 3X3 q- q4 x4 -~-q 5x5 -- X) + 1 (1 -- X).
du
zf
(46)
re
Taking Uo = u(0), d f .=uo
qo
duu
1
= a - --
+ -- -
'~f
where
"£'c
(47)
Pdx=o,
1(o )1
(48)
d2f d ( d f ) dx du 2 -- dx duu duu = p2(x)pl(x)'
(49)
pl(x)
zr
q~x~--x
+ --zc ( l - - x ) ,
and so
where Pz(X)
dx
vqvx ~- 1 _ 1
~f
zc
(50)
Hence
d2f2 du
= fl - (PlP2)lx=O.
(51)
= Y ~ [Pl(P3Pl +PzZ)]lx=o,
(52)
.=uo
Similarly ddu 3 f3 u=uo
dV du4 . . . . = 6 = [PI(P4P~ + 4plPzPa + p3)] Ix= o, dSf 3 du5 . . . . = e = [ p l ( p s p l +
7 2 11 z 422 4 PxPzP4+ P l P z P a + PlP3+P2)]Ix=o
(53)
(54)
and du6 . . . . = ~ = [Pl(P6P~ + l lplaP2P5 + 32plp2p4 2 2 + 15plpap4 3 3 + 34plP2P3 2 2 + P2)] 5 1~=o, + 26paP2P3
(55)
where pa(x)
dp2 dx
p4(x)
--
1 5
~ v(v--1)q~,x v-2,
dp3 _ 1 ~ v ( v _ l ) ( v _ 2 ) q v x ~ _ 3 , dx "~f 3
ps(x) = dp4 _ i ~ v ( v _ l ) ( v _ 2 ) ( v _ 3 ) q ~ x ~ _ , , , dx
(56)
"of 2
(57)
(58)
~f 4
terminating at p6(x)
dps dx
120gs zr
(59)
An exact solution for low-power neutron multiplication
73
Noting that P7 = 0,it is clear that it is possible to find as many coefficients for the McLaurin series representation as are required. Therefore, (U - - UO) 3
y,(u) = ~(u-- Uo) + fl (u-- Uo)2 ~- 7 - 2
~-
(60)
3!
and
Gl(x,t) = ct(W-Uo)+fl ( w - u ° ) 2 + 7 (- w --u ° ) 3 + "", 2 3!
(61)
where again w = t + u(x). As u(x) is infinitely differentiable with respect to x, it can be seen that this representation o f f can be suitably manipulated, and thus any probability or moment may be found providing that the series is convergent in the region
of interest. It is well-known that a McLaurin series is not necessarily convergent but that it usually does converge for a specific range of values, so the task here is to determine over what time range this representation accurately reproduces probabilities. This is most easily done for a simplified case such as the doublet fission model, for which q2 = 1 and all other q~ = 0. Whence, taking zf/z~ = 1.4 again, 1.4 zf
p f(0)
2.4 p2(0)
--
2 Zf
p3(0)
"~f
P'> 3(0) = 0.
Thus, ~=
1.4
3.36 /~=----,
,
e
338.7238
11.984
7=
56.9856
2416.060
and
~
6--
~6
To study the convergence characteristics of the McLaurin series the following approximate solutions for the extinction probability Po = GI(0, t) are considered: t2
pl=~t,
t3 P3 = ,02-I-7 3.1'
p2=ctt+fl~,
t4
g,t 5
Pa=P3+ ~..,
PS =Pa q- 5~
Ct6 and
(62)
p6=p5+~..
These may be compared directly with the known exact solution which, substituting for rf/r c in equation (8), may be written as P0 =
1 - exp ( - 0.4t/~f)
(63)
1 1 - - - exp ( - 0 . 4 t / r f ) 1.4
The consequent numerical results are displayed in Table 3. When plotted (Fig. 2) these results demonstrate clearly that for t < 0.84~f the addition of further terms to the McLaurin series increases the accuracy of the representation, but for t > 0.84~f these extra terms only serve to increase the errors. Thus, even an infinite McLaurin series can only usefully model the doublet fission stochastic behaviour for 0 ~< t ~ 0 . 8 4 ~ ' f . The empirical fission data case is of course less clear-cut. Using the figures in Table 1 and rf/T c = 1.4 again, the following McLaurin series coefficients are obtained : 1.427
3.199334
e
47.516521
8.553538
and
45.883982 ~ - - -
23.145467
(J.
/4
1. I~ARKS a n d J. 1). LEWINS
. . . Table 3. M c L a u r i n series approximations to
Po(t). for the doublet fission model when -~f-
= 1.4
t -Tf
Pl
P2
P3
P4
P5
P6
Po
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0.140000 0.280000 0.420000 0.560000 0.700000 0.840000 0.980000 1.120000 1.260000
0.123200 0.212800 0.268800 0.291200 0.280000 0.235200 0.156800 0.044800 -0.100800
0.125917 0.228779 0.322728 0.419029 0.529667 0.666624 0.841885 1.067435 1.355256
0.124960 0.224980 0.303495 0.358245 0.381267 0.358902 0.271792 0.094880 -0.202588
0.124988 0.225883 0.310355 0.387149 0.469476 0.578395 0.746202 1.019822 1.464187
0.124985 0.225668 0.307908 0.373404 0.417044 0.421834 0.351415 0.140612 -0.319137
0.124985 0.225709 0.308551 0.377834 0.436591 0.487014 0.530727 0.568955 0.602647
These coefficients in conjunction with equations (62) yield the numerical results in Table 4. These may be compared with the exact probability curve found by the method described in Section 5. These figures are taken from Table 5. These curves are shown in Fig. 3. Although less well-defined than in the 'doublet fission' case, it is apparent that the solutions diverge from around t -~ ~f. Unfortunately this useful time period is much less than that of interest to the nuclear engineer, so one is forced to conclude that the McLaurin series model o f f is not suitable after all. By studying an alternative means of deriving an equivalent McLaurin series representation of the doublet fission model, the problem of its convergence is made clearer. The work above seeks to modelffor Gl(x, t) =f(t +u(x)). Because of the simplicity of the doublet fission case, instead of defining the arbitrary constant ~ in equation (25) by O=t+~ln
~
,
(64)
as is effectively done here, one could redefine a new constant ~ by x-- 1 ( x - - 1)T = e (x-")~ = - e (x-")t x-#/k x-#/k '
so now
7((x-1)q
G,(x,t)= \ ~ ]
and
where
T = e Ix-")t,
G,(x,O)=x.
(66)
16 14 12 10
P1
Q, 0 8
'°5
06 ,q_
04-
~P_ o
0.2 0 t
0.1
0.2
0.3
0,4
--02
0.5 t/'rf
0.6
07
0.8 '~.hoO9
--04 F i g . 2. M c L a u r i n series a p p r o x i m a t i o n s to
Po(t) for
(65)
the 'doublet fission' m o d e l w h e n z r = 1.4.
An exact s o l u t i o n lor l o w - p o w e r n e u t r o n m u l t i p l i c a t i o n
Table 4. McLaurin series approximations to
/3
. . . ~f Po(t)for the empirical fission data case when - -
= 1.4
Tc
t -Tf
P~
P2
P3
P4
P5
P6
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
0.142700 0.285400 0.428100 0.570800 0.713500 0.856200 0.998900 1.141600 1.284300 1.427000 1.569700 1.712400
0.126703 0.221413 0.284130 0.314853 0.313583 0.280320 0.215063 0.117813 -0.011430 -0.172667 -0.365897 -0.591120
0.128129 0.232818 0.322621 0.406091 0.491782 0.588247 0.704040 0.847715 1.027825 1.252923 1.531563 1.872298
0.128032 0.231275 0.314809 0.398279 0.431507 0.463262 0.472489 0.452699 0.395085 0.288528 0.119593 -0.127470
0.128036 0.231402 0.315772 0.402334 0.443881 0.494052 0.539040 0.582451 0.628902 0.684499 0.757308 0.857833
0.128037 0.231406 0.315818 0.402595 0.444877 0.497026 0.546538 0.599157 0.662770 0.748227 0.870206 1.048123
Writing x--1
#/2=a a n d t a k i n g x o -- O, i.e.
Yo = 1/a, o n e
(67)
x-#/2'
finds that d " f y=yo dy"
So using a McLaurin
y
and
nI(-a)"+l
(68)
(I -a)"
series representation,
n=l
=
~
-a
.=1
(69)
(1 --a)"
gI!
(l-ay)"
(70)
(l-a)"
18
1 6
1 4
12
1 0
P
a
~ 06
0
022
0.0.4
006 t/.rfO 8
I 0
'12
--0.2
-o~
~
>k,k°
--0.6
Fig. 3. M c L a u r i n series a p p r o x i m a t i o n s to
Po(t) for
the e m p i r i c a l fission d a t a case w h e n Ze = 1.4.
76
G . T . PARKS and J. D. LEWINS
This is a geometric progression and, assuming it converges, (1 --ay) (
fm(y)= --a (~_a)
1--ay~ -1 1-- l - - a ]
(71)
or
fro(Y) = 11--ay y '
(72)
whence
GI,, - 1 - - a y T 1 -yT
where
T = e ~-~)'.
(73)
'
Substituting for a and y and rearranging gives #(x- 1)T-(2x-/~)
Gain(x, t) = 2(x-- 1)T--(2x-/.0 "
(74)
Thus, comparing equation (74) with equation (8), it is seen that the McLaurin series representation admits an exact solution of the doublet fission problem, providing it converges. It is convergent if the geometric progression converges, i.e. 1 -ay
1-a
< 1.
(75)
The region of interest is 0 ~< x ~< 1 which is equivalent to 1/a >>-y >t O. Thus if R = (1 -- ay)/(1 -- a), IRlmax =
when
(76)
y = 0.
(77)
So for complete convergence l la
< 1.
(78)
Therefore : if a < 1,
equation (78) ~ a < 0
for convergence;
if a > 1,
equation (78) ~ a > 2
for convergence. ~Cf
Thus the McLaurin series only converges over the complete range of interest if ~f < 0 or - - > 2. This Tc
Tc
Tf
unfortunately excludes most of the range of values o f - - of interest to the nuclear engineer, so again one concludes Tc
that the McLaurin series representation off is not suitable for this problem, even though the formal solution to the doublet model is correct.
5. SUCCESSFUL METHOD
Suppose in the analysis of Section 3 one writes, instead of ~ = f(~,), ~ = 0(~b). Then
t +u(x) = O(GO.
(79)
O(x) = u(x)
(80)
When t = O, G1 = x, therefore
An exact solution for low-power neutron multiplication
77
and
u(G,) = t +u(x)
(81)
or
A l n \ x( G _ lI - I h + ] B In ( ~ -) + C
( G l - a 2 ~ + D In ( ~ ) +- D * l n ( G ~ - a '\~x h- a *
In \ x - - a 2 /
J =t.
(82)
A l t h o u g h this expression c a n n o t easily be inverted to give G1 as a n explicit function o f x a n d t, it m a y be used to evaluate t as a function o f x a n d G,, a n d hence give Ga as a n implicit function o f x a n d t. In particular, setting x = 0 gives the extinction p r o b a b i l i t y Po(t) implicitly from
u(Po) = t + u(O).
(83)
Consider n o w the effect of differentiating e q u a t i o n (81) with respect to x:
63G1
if(GO ~-x = u'(x),
(84)
where du
A x- 1 +
u'(x) - dx N o t i n g t h a t ~~G1 5=0 =
B
C
x-a~- +
D -+ x-a2 x-a
D* + 3 x --a '~
(85)
P~(t), setting x = 0 in e q u a t i o n (84) gives u'(Po)P~ = u'(0).
(86)
It is s h o w n a b o v e t h a t Po m a y be found as an implicit function of time from e q u a t i o n (83), so Pa m a y also be found as an implicit function of time from u'(0)
(87)
P1 = u,(Po)"
By further differentiation with respect to x a n d setting x = 0 successive probabilities P , m a y be found as implicit functions of time. F o r example, P2 =
P3 =
u"(O)-- u"(Po)P ~ 2u'(Po)
,
(88)
u"(O)- 6u"(Po)P1P2 - u'"(eo)e 3 6u'(Po)
etc.,
(89)
where
d2u
du"
u"(x)=~-
and
u. . . .
(90)
dx
To illustrate this method, results are t a b u l a t e d below for Po, P1 a n d
P2. This is d o n e for b o t h the sub-critical case,
zf r~ = 1.4, a n d for a super-critical c a s e , - = 0.7. "6c
17c
F o r r f = 1.4, as before, Zc
a 1 -~ 0.970, "gf
A ~- 0 . 5 2 9 3 - - , qs ANE 12:2-D
a 2 ~ --3.538
B -~ --0.5374 zf q5
and
a 3 ~- - 0 . 9 6 8 + 3 . 1 6 3 j ~Tf
c ~ 0.0029--, q5
D -~ (0.0026--0.0012j) zf q5
78
G.T.
PARKS and J. D. LEwlrqs
Table 5. Exact probabilities for the general stochastic case when
Table 6. Exact probabilities for the general stochastic case when
z~ = 1.4
zf = 0.7
"~e t ( x 10-3 zf~ \
Po
"Co
Pt
P2
qs/
0.00 2.91 6.38 10.66 16.13 23.49 34.14 51.39 85.58 195.73
( t x l O azf_)
Po
P~
P2
0.00 5.92 13.56 24.05 39.99 69.41 185.47
0.000 0.100 0.200 0.300 0.400 0.500 0.600
1.000 0.793 0.597 0.418 0.257 0.120 0.012
0.000 0.042 0.073 0.086 0.081 0.054 0.008
qs/
0.000 0.100 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900
1.000 0.845 0.697 0.557 0.426 0.308 0.204 0.117 0.052 0.011
0.000 0.024 0.044 0.059 0.068 0.069 0.062 0.047 0.027 0.007
For z~ = 0.7, Te
al -~ 0.614, 17f
A -~ 0 . 0 4 9 0 - - , q5
a2 -~ - 3 . 2 4 0 ~'f
B -~ - 0 . 0 6 0 4 - - , q5
and
a 3 -~ - 0 . 9 3 8 + 2 . 9 5 4 j qSf
C -~ 0 . 0 0 4 4 - - , q5
D -~ (0.0035-0.0015j) re. q5
Substitution into equation (83) etc. gives the results in Tables 5 and 6. It is interesting to compare these results with the probabilities predicted by the 'quadratic approximation' model, where 2 and/~ [as defined by equation (11)] are chosen so that the resulting solution is exact to the second moment of n. The solution for G t [equation (8)] yields the following expressions for the probabilities under consideration : Po(t) = Gl(O,t)
#
1-- T
= 2- # / 2 - - 7 z
(91)
P~(t) = Oxx ~G1 x=o = (1--kt/)t)ZT/(~/2 - T) 2
(92)
and /~2(t)
= _l ~1
(93)
2 c~x2 x= o = (1 - #/2)2(1 - T ) T / ( # / 2 - T) 3,
giving the results in Tables 7 and 8. Table 7. Probabilities predicted by the 'quadratic approximation' for
Tf = ~c
Table 8. Probabilities predicted by the 'quadratic approximation' for
~f Te
1.4
-- =0.7
.o
..
0.00
0.000
1.000
2.91 6.38 10.66 16.13 23.49 34.14 51.39 85.58 195.73
0.153 0.283 0.396 0.497 0.588 0.672 0.752 0.831 0.910
0.714 0.508 0.358 0.246 0.162 0.101 0.056 0.024 0.006
0.000 0.112 0.148 0.146 0.126 0.098 0.070 0.043 0.021 0.005
t(× I0-3~) 0.00 5.92 13.56 24.05 39.99 69.4l 185.47
15o
Pl
P2
0.000 0.196 0.341 0.453 0.544 0.619 0.677
1.000
0.573 0.330 0.183 0.092 0.035 0.002
0.000 0.165 0.165 0.122 0.073 0.032 0.002
An exact solution for low-power neutron multiplication
79
These exact and approximate probabilities are plotted in Figs 4-6. These graphs show that the exact probabilities are indeed of the form one would anticipate from the physical considerations, but show a substantial inaccuracy in the probabilities predicted by the 'quadratic approximation' during the transient development. This is consistent with the discrepancies found in the fully-developed solutions by Campbell et al. (1981) for the subcritical case and Lewins (1981) for the super-critical case. [] Exact P1 rf /Tc = t 4 10 I
~
1011
• Quadratic approximation rf / ~ : 1 4
09
09
:: Exact P1 rf
o k
o
./,/
o,1
/r c =07
• Quadratic approximation rf /Tc = 0 7
o7
•
05 04
03 F/m~ / 02
'"/~/
O11~
03
. Quadratic approximation rf / ~ ----14
02
O Exact PO0 ~ / ~ = 0 7
y 0
D Exact PO ~/Tc ----14
• Quadrat c approximation ~ / ~ : 0 7 [ 25
I 50
I 75
l 100
I 125
for
z ~ = 1.4
"['¢
I 175
J 200
I
0
25
50
75
/qs)
t (X lO-3"rf
Fig. 4. Po(t)
I 150
0 1
t (xlO
and
014~1
for
[3 Exact P2 rf /Tc = 1 4
L
016V~
Fig. 5. Pl(t)
T~, =0.7.
1~c
O 20 I 018
\
• Quadratic approximation Tf / r c ~ 14 o Exact P2 rf / r c = 0 7 • Quadratic approximation r f / ' r c : 0 7
/
oo i ..,x " oo ![ O0 2
0
100
25
50
75
100
125
t (xlO-3"q
Fig. 6. P2(t)
for
~rr= 1.4 "Ce
150
175
/qs) and
Zf=o.7. "£c
200
__
125
150
rf
/qs)
--3
~f = 1.4 Tc
and
"L'f
175
--=0.7. Tc
200
80
G . T . PARKS and J. D. LEW1NS 6. CONCLUSIONS
The success in finding the exact neutron probabilities for the problem based on empirical fission data for a source-free system is undoubtedly due to avoiding the need to find the generating function explicitly; the unsuccessful methods having failed to do this. An explicit expression for G1 is required, however, to evaluate the other basic generating function F o through the convolution integral. Unless F o is available, the system with completely general initial conditions remains insoluble exactly. Thus further work must be directed towards finding a successful means of inverting the relationship for G 1 found here. However, the implicit method we have given could be the basis of a numerical convolution. The implicit method we have given here is straightforward and lends itself to numerical computation with a digital computer. The resulting numerical values demonstrate the weakness of the quadratic approximation for this model of prompt-neutron stochastics in the transient period before the establishment of the well-developed or asymptotic solution, a conclusion that reinforces the results previously obtained for asymptotic solutions. REFERENCES
Abramowitz M. and Stegun I. A. (1964) Handbook of Mathematical Functions. NBS, Washington, D.C. Bartlett M. S. (1979) An Introduction to Stochastic Processes, 3rd edn. Cambridge Univ. Press, U.K. Bell G. 1. (1963) Ann. Phys. (N.Y.) 21,243. Campbell G. A., Chan V. and Lewins J. D. (1981) Ann. nucl. Energy 8, 1. Divert B. C., Martin H. C., Taschek R. F. and Terrel J. (1956) Phys. Rev. 101, 1012. Gurland J. (1978) A Seminar on Statistics. CNEN, Rome. Harris D. R. (1964) Naval Reactors Physics Handbook, Vol. 1 (Edited by Radkowsky A.). USAEC, Washington, D.C. Harris T. E. (1963) The Theory of Branching Processes. Springer-Verlag, Berlin. Jeffreys H. and Jeffreys B. S. (1956) Methods of Mathematical Physics, 3rd edn. Cambridge Univ. Press, U.K. Kendall D. G. (1948) Ann. math. Statist. 19, 1. Lewins J. D. (1978) Ann. nucl. Energy 5, 141. Lewins J. D. (1981) Ann. nucl. Energy 8, 15. Pacilio N., Colombino A., Mosiello R., Norelli F. and Jorio V. M. (1978) Adv. nucl. Sci. Technol. II, 67. Smith P. R. (1970) J. nucl. Energy 24, 199. Williams M. M. R. (1974) Random Processes in Nuclear Reactors. Pergamon Press, Oxford. Williams M. M. R. (1979) Ann. nucl. Energy 6, 43.
0306-4549/8553.00+0.00 Copyright © 1985 Pergamon Press Ltd
AN EXACT TRANSIENT STOCHASTIC SOLUTION FOR LOW-POWER NEUTRON MULTIPLICATION G. T. PARKS* and J. D. LEWINS Engineering Laboratories, University of Cambridge, Trumpington Street, Cambridge CB2 1PZ, England (Received 15 August 1984)
Abstract--An exact solution for arbitrary physical data of the neutron-chain multiplication problem at low power is sought for the transient behaviour of the prompt-neutron stochastic chain to complementthe known fully-developedsolutions.Three means of inverting the equation for the generating function in order to find it explicitlyare investigatedand found to be unsatisfactory.A method of findingindividual neutron probabilities as implicitfunctions of time in a source-freeassemblyis demonstrated and typical results, for sub- and supercritical cases, are evaluated and comparison made with the 'quadratic approximation' model.
I. I N T R O D U C T I O N
There has been considerable work in the stochastic theory of multiplying chains since the original studies by Galton and Watson in the field of social anthropology [see reviews in Harris (1963) and Williams (1974)]. Solutions, both exact and approximate, have been found to a wide range of specific problems. In particular, it is well-known that an exact solution exists for the case of fission producing precisely two neutrons (e.g. Kendall, 1948), variously known as the 'linear birth and death', 'binary fission' or 'escalator' problem (Bartlett, 1979), and also for the form of the problem after truncation of the exact fission distribution at second order (Bell, 1963). Kendall demonstrated, though considering a generalized problem with time-varying coefficients, that 'no exact solution' exists except for the 'linear birth and death' case. This result is clearly applicable to the problems with constant coefficients studied here, but in what sense is there 'no solution'? This must mean 'no explicit closed solution', for there is no doubt that, in some senses, the general stochastic problem is well posed. In particular, as the finite moments have closed equations, these are in principle always available. In addition, Williams (1979) has given exact solutions for several other mathematical models of fission-neutron distributions through use of Lagrange expansions, concentrating predominantly on finding extinction probabilities. Other work in the field has considered the prompt-jump approximation and numerical solutions (Harris, 1964; Smith, 1970), and, most recently, Campbell et al. (1981) and Lewins (1981) found solutions for fully-developed subcritical assemblies including a source and for super-critical source-free reactors, respectively, both based on arbitrary physical data. As yet, however, there has been little consideration of the transient behaviour of systems, particularly for models based on empirical data. It is the form of these transients leading to the developed solutions for empirical data which is investigated here. Such solutions provide standards against which the accuracy of the various approximate models may be checked. We also seek to help develop understandingof the problem and consequently this paper is written with this didactic purpose in mind--we show some of the methods that did not work. 2. M O D E L S AND E Q U A T I O N S
The foward stochastic equations of the discrete model may be written (Lewins, 1978) in terms of the generating function F(x, t) in the lumped, low-power case as
r
OF O(X)-- x_ + 0 ~ - = k zf ~-
~x
+S[f~(x)--l]F,
* Present address : GEC Electrical Projects Ltd, Rugby, Warwickshire, England. 65
(1)
66
G . T . PARKS and J. D. LEWINS
where 1
is the mean rate for a neutron to cause fission,
Tf
1
S
L(x)
is the mean rate for a neutron to be captured, is the mean rate of source events, is the associated auxiliary generating function for such source neutrons
and
g(x)
is the auxiliary generating function for induced fission ;
S, zf and z¢ are here assumed to be constants independent of time. The main generating function F(x, t) is defined by
F(x,t) = ~ P,(t)x";
F(x,O) = x N,
(2)
n=O
where P,(t) is the probability that the reactor contains exactly n neutrons at time t and N is the initial number of neutrons at time t = 0. If F(x, t) can be determined analytically, two useful types of information are available : the ruth probability
O'F
= m! P.,(t)
(3)
= (n(n-- 1)... (n-- m + 1)>.
(4)
G~xm x=O
and the mth factorial moment
63xm x = 1
The auxiliary generating functions are defined by
g ( x ) = ~ q~x ~ and
f~(x)= ~, p~x ~,
v=O
(5)
v=O
where qv
is the probability of exactly v neutrons being released in an induced fission event
and Pv is the probability of exactly v neutrons being released in a source event. For simplicity it is generally assumed thatf~(x) = x, i.e. that each source event results in the release of exactly one neutron. This convention will be followed throughout this paper. If, in addition, it is assumed that each induced fission event causes precisely two neutrons to be released, i.e. 9(x) = x 2, the general stochastic equation reduces to OF
~t
-
k
X2 - X zf
OF +
~
-
7x
+ S(x-1)F.
(6)
This is the form of the problem known as the 'linear birth and death' approximation. The authors find this and the various other titles for this model ('escalator' and 'binary fission') either misleading or ambiguous. As this case considers singlet independent sources and doublet fissions, we suggest that in the context of neutron stochastics this equation be known as the 'doublet fission' model. Whatever its name, its solution is readily found. The independence of neutron behaviour in a low-power reactor means that the solutions for differing initial conditions may be found by superposition from two fundamental solutions : Finn = F ~ G ' ~ ,
where m n
G1
is the number of neutrons initially present, is the number of independent sources [assumed to have the same auxiliary function fs(X)], is the solution for rn = l, n = 0 (i.e. for no source, one initial neutron)
(7)
An exact solution for low-power neutron multiplication
67
and Fo
is the solution for m = O, n = 1 (i.e. for no initial neutrons, one source).
For the 'doublet fission' model /~(x- 1 ) T - ( 2 x - / ~ )
(8)
and
2-/~
(9)
Fo = [2 T _ I~_~( T _ I )x ] s/a, where 2=
1
,
#=--
%'f
1
T=exp(A-#)t.
and
(10)
2"c
If, alternatively, the following identifications are made J,=lg(2)(1)/Tf
1
I~=--+[1--gm(1)+½gtZ~(1)]/rf,
and
(11)
Z"c
where
d"9
9(")(1) = dx" x=l
=
(12)
the resulting solutions from equations (8) and (9) are exact to the second moment ofn. As this change is equivalent to defining
g(x) =
1 + g m ( 1 ) ( x - 1)+½g(2)(1)(x - 1)2
(13)
this case is known as the 'quadratic approximation'. It is important to note that the single-source solution Fo may be shown (e.g. Bartlett, 1979) to be related to the single initial neutron solution G 1 by the convolution
F o = e x p { f l S [ G , ( x , t - - z ) - l ] d r }.
(14)
Thus only one of the fundamental solutions needs to be found directly from the general partial differential stochastic equation.
3. GENERAL NO-SOURCE CASE--SOLUTION BY LAGRANGE'S METHOD With no source present the general stochastic equation (1) reduces to
[
~F _ _ 0 ( x ) - x
~3F + 1--x I ~-x
The subsidiary generating function 9(x) is determined by physical data appropriate to the reactor fuel etc. (Diven et al., 1956). Typically it is found that no more than five neutrons are released in a fission reaction, so 5
5 g~V)(1)(1 _x)V
O(x) = ~" qvx~= ~ v=O
v!
'
(16)
q,x v -- + ~f 3¢ J Ox'
(17)
v=0
where 9~)(1) is defined by equation (12). Substituting for g(x) in equation (15) gives
= ~
v=o
68
G . T . PARKS and J. D. LEWINS
therefore OF Ot
q5 rr ( x - a o ) ( x - a l ) ( x - a 2 ) ( x - a 3 ) ( x - a 4 )
OF 0x'
(18)
where ao, a l , . . . , a4 are the roots of the polynomial equation 1
x 1-x q,x v - - + - = O .
5
-- ~ "L'f v = O
"t'f
(19)
27c
x = 1 is necessarily a solution of equation (19) and it is found that there are two further real roots and a complex conjugate pair. Hence OF Ot
q5 ( x - 1 ) ( x - a O ( x - a 2 ) ( x - a 3 ) ( x - a * zf
0F ) ~-x = 0.
(20)
For a single initial neutron F - G~, so applying Lagrange's method dt
- z f dx
1
dGa -
qs(x--1)(x--at)(x--a2)(x--a3)(x--a~)
0
(21)
therefore G1 = 4,
aconstant.
(22)
Also -- dt = zf d x / q s ( x - 1) (x - a 0 ( x - a2) (x - a3) (x - a*)
(23)
therefore dt=-dx
A ~ +
B -x--a
C
+ 1
X--a
+ 2
D x--a 3
+
D* ] --~- , x--a 3
(24)
where A, B, C and D are obtained by the usual partial fractions techniques. Hence t = ~k--[A In ( x - I ) + B In ( x - a O + C
In ( x - a 2 ) + D In ( x - a 3 ) + D * In ( x - a~')],
(25)
where ~ is a constant. Now ~b =f(~k), an arbitrary function of ~, therefore In ( x - 1)+ ... + D * In (x-a*)).
G1 = f ( t + A
Whent=0,
(26)
Gt = x a n d x = f ( A In (x-- 1) + ... + D* In (x--a*)) ;
(27)
equation (27) is of the form x = f(u(x)), where u(x) = A In ( x - 1 ) + B In ( x - a O + C
In ( x - a 2 ) + D In ( x - a 3 ) + D * In ( x - a * ) .
(28)
Given this it should be possible to model f ( ) by an appropriate series, or other, representation. Then (29)
G I = fm(t + U(X)),
where fro is the model off. Equation (29) would then allow Gx to be evaluated or manipulated as a function of t and x as required. We therefore seek to invert u(x) as x =f(u).
4. U N S U C C E S S F U L
METHODS
T o h e l p v i s u a l i z e t h e p r o b l e m u ( x ) m a y b e f o u n d f o r t y p i c a l p h y s i c a l d a t a ( T a b l e 1).Ifzf/z c = 1.4(corresponding to a sub-critical system), the roots of equation (19) are as follows: ao = 1,
al -~ 0.970,
a 2 -~ - 3 . 5 3 8
and
a 3 -- - 0 . 9 6 8 + 3 . 1 6 3 j ,
An
exact solution for low-power neutron multiplication
69
Table l. Sample 23 ~U fission data (after Diven et al., 1956) qo = 0.027 ql = 0.158 q2 = 0.339
q3 = 0.305 q4 = 0.133 q5 = 0.038
whence ~'f
A "" 0.5293 - - , qs
Obviously only the real part of
~'f
Tf
B -- --0.5374 - - , q5
"/~f
C -~ 0.0029 - - , qs
u(x) is of interest. Manipulation
D ~- (0.0026-0.0012j) - - . q5
of equation (28) gives
~ . (u(x)) = A In (1 - x ) + B In ( x - a O + C In ( x - a 2 ) + 2 R In ( x + 1 ) + R In (1 + z 2 ) + 2 I tan-1 z,
(30)
where
R = ~ e (D), I = J ~ (O) and
z = J~
(31)
(a3)/(x +
1).
Substitution into equation (30) gives the results in Table 2. This allows x to be plotted as a function of u, as in Fig. 1. Data for a typical super-critical case is given in Table 6. Three different forms of representation are investigated in the search for a suitable model for f: (a) Lagrange interpolation (b) Lagrange expansion (c) McLaurin series.
(a) Lagrange interpolation Lagrange's interpolation formula for a function may be written as
fro(u) = ~ lff(u,), o
where
l,
= fi u-u;. j=o Ul--Uj
(32)
This effectively gives a polynomial of order n in u which matchesf(u) at n + 1 points (u~, f(u3). Clearly, ifthe function is matched at a large number of points, the accuracy of the representation will be very good. It is also easy to find
Table 2. u(x) for "rf = 1.4
u(x)
ANE 12:2-C
x
( ×lO-3zf~qs/
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
23.19 26.10 29.57 33.85 39.32 46.68 57.33 74.58 108.77 218.92
70
G.T. PARKSand J. D. LEWINS 10 ~..-.-.0,-
09 08 O7 0.6 "~05 04 0.3 0.2 01
0
25
50
75
100
I 150
125
I
I
I
175
200
225
t./ (X 10-3"i'f /q5 )
Fig. 1. Plot ofx as a function ofu for -- = 1.4. Tc
Gl(x, t) from this model:
fi Uw-u,i
(33)
i=0 j=0
where
w = t + u(x).
(34)
Although it is relatively simple to write a computer program to generate Gl(X, t) from this expression, it must be remarked that information is contained not only in G1 but also in its derivatives with respect to x. The analogous expression for gG1/gx is
dG 1 = dx
ifi
__1
ifi
( w - u,)f(ui) ~x"
(35)
~=oj=o ui-- uj k=O t=O =~i
~:i ~=i~=k
Derivatives of considerably higher order than this are of interest and it is apparent that their expressions in Lagrange interpolation notation will be correspondingly complicated. Although a formal differentiating program might be possible, the difficulties in assimilating these expressions in a program are such that the Lagrange interpolation representation only seems viable in the absence of a suitable or simpler alternative.
(b) Lagrange expansion An initially more promising means of representing fasf= would appear to be through a Lagrange expansion, as suggested by Gurland (1978). Various forms of Lagrange expansion have been used in the literature in this connection although the correspondence is not always apparent (Williams, 1979). Gurland proposes that the expansion in the form forf(u(x)) becomes
x = 1 + ~ bi ui
around
u = 0,
i=1
acknowledging that the practicability of the method to find the bi depends on the nature of the previously defined coefficients 'A' of the partial fraction expression, without actually finding the b~. However, for x to be a proper stochastic variable we require 0 ~< x ~< 1 and Gurland's form is unattractive; there
An exact solution for low-power neutron multiplication
71
is no reason to suppose the bi are such as to secure this limit to x. If instead we consider
i.(u): xo+ k=l L
k~
\U--Uo/d .....
L
(36)
,
where D" = d"/dx" and taking x o = O, this reduces to
i
k=l
(u-,o),
k!
\ U - U o / 3x=O
k
= i k=l
k!
(37)
'
where ,, =
,
L
(38)
\U-Uo/J~=o
whence (W-Uo) ~
Gl(x,t)= ~ h k - -
k=l
(39)
k!
N o w consider the first coefficient of this series, h 1,
therefore
hl~~ A
x
(41)
ln(l_x)+Bln(~X)+Cln(~_2X)+Dln(~3X)+D,
ln(a*--x']]
\~-/jx=o
"
It is difficult to state precisely what value this gives for ha, and similar problems are encountered for all the other coefficients hk. It is thus found that, even though the conditions of applicability for this expansion are not contravened, it cannot be applied to this particular problem. This is not so surprising as "it is unusual for the general coefficient in the inverted series to be expressible in a convenient form" (Jeffreys and Jeffreys, 1956). We are aware of the successful use of the Lagrange expansion in obtaining the extinction probability Po(t) (e.g. Williams, 1979) and this remains a technique worth further investigation.
(c) McLaurin series One alternative method of expressing a function as a polynomial series is by a McLaurin series representation :
fro(U) = f(Uo) + (U-- UO) dd~ . . . . -t
(U--Uo) z d2f 2!
du2 . . . . +
(U-Uo) 3 d y . . . . 3! du 3 + ..-.
(42)
Now df du
d f dx dx du
and
df
dx
f(u) = x,
(43)
thus
dR du" Differentiating equation (28) :
db dx
A B C D D* + + + + x - 1 x - a 1 x - - a 2 x - - a 3 x--a*
(44)
72
G . T . PARKS and J. D. LEWINS
SO
(45)
du __ zf ( x - - 1 ) ( x - - a O ( x - - a 2 ) ( x - - a 3 ) ( x - - a ' ~ ) dx q5
and d x _ 1 (qo + q 1x + q2 x2 + q 3X3 q- q4 x4 -~-q 5x5 -- X) + 1 (1 -- X).
du
zf
(46)
re
Taking Uo = u(0), d f .=uo
qo
duu
1
= a - --
+ -- -
'~f
where
"£'c
(47)
Pdx=o,
1(o )1
(48)
d2f d ( d f ) dx du 2 -- dx duu duu = p2(x)pl(x)'
(49)
pl(x)
zr
q~x~--x
+ --zc ( l - - x ) ,
and so
where Pz(X)
dx
vqvx ~- 1 _ 1
~f
zc
(50)
Hence
d2f2 du
= fl - (PlP2)lx=O.
(51)
= Y ~ [Pl(P3Pl +PzZ)]lx=o,
(52)
.=uo
Similarly ddu 3 f3 u=uo
dV du4 . . . . = 6 = [PI(P4P~ + 4plPzPa + p3)] Ix= o, dSf 3 du5 . . . . = e = [ p l ( p s p l +
7 2 11 z 422 4 PxPzP4+ P l P z P a + PlP3+P2)]Ix=o
(53)
(54)
and du6 . . . . = ~ = [Pl(P6P~ + l lplaP2P5 + 32plp2p4 2 2 + 15plpap4 3 3 + 34plP2P3 2 2 + P2)] 5 1~=o, + 26paP2P3
(55)
where pa(x)
dp2 dx
p4(x)
--
1 5
~ v(v--1)q~,x v-2,
dp3 _ 1 ~ v ( v _ l ) ( v _ 2 ) q v x ~ _ 3 , dx "~f 3
ps(x) = dp4 _ i ~ v ( v _ l ) ( v _ 2 ) ( v _ 3 ) q ~ x ~ _ , , , dx
(56)
"of 2
(57)
(58)
~f 4
terminating at p6(x)
dps dx
120gs zr
(59)
An exact solution for low-power neutron multiplication
73
Noting that P7 = 0,it is clear that it is possible to find as many coefficients for the McLaurin series representation as are required. Therefore, (U - - UO) 3
y,(u) = ~(u-- Uo) + fl (u-- Uo)2 ~- 7 - 2
~-
(60)
3!
and
Gl(x,t) = ct(W-Uo)+fl ( w - u ° ) 2 + 7 (- w --u ° ) 3 + "", 2 3!
(61)
where again w = t + u(x). As u(x) is infinitely differentiable with respect to x, it can be seen that this representation o f f can be suitably manipulated, and thus any probability or moment may be found providing that the series is convergent in the region
of interest. It is well-known that a McLaurin series is not necessarily convergent but that it usually does converge for a specific range of values, so the task here is to determine over what time range this representation accurately reproduces probabilities. This is most easily done for a simplified case such as the doublet fission model, for which q2 = 1 and all other q~ = 0. Whence, taking zf/z~ = 1.4 again, 1.4 zf
p f(0)
2.4 p2(0)
--
2 Zf
p3(0)
"~f
P'> 3(0) = 0.
Thus, ~=
1.4
3.36 /~=----,
,
e
338.7238
11.984
7=
56.9856
2416.060
and
~
6--
~6
To study the convergence characteristics of the McLaurin series the following approximate solutions for the extinction probability Po = GI(0, t) are considered: t2
pl=~t,
t3 P3 = ,02-I-7 3.1'
p2=ctt+fl~,
t4
g,t 5
Pa=P3+ ~..,
PS =Pa q- 5~
Ct6 and
(62)
p6=p5+~..
These may be compared directly with the known exact solution which, substituting for rf/r c in equation (8), may be written as P0 =
1 - exp ( - 0.4t/~f)
(63)
1 1 - - - exp ( - 0 . 4 t / r f ) 1.4
The consequent numerical results are displayed in Table 3. When plotted (Fig. 2) these results demonstrate clearly that for t < 0.84~f the addition of further terms to the McLaurin series increases the accuracy of the representation, but for t > 0.84~f these extra terms only serve to increase the errors. Thus, even an infinite McLaurin series can only usefully model the doublet fission stochastic behaviour for 0 ~< t ~ 0 . 8 4 ~ ' f . The empirical fission data case is of course less clear-cut. Using the figures in Table 1 and rf/T c = 1.4 again, the following McLaurin series coefficients are obtained : 1.427
3.199334
e
47.516521
8.553538
and
45.883982 ~ - - -
23.145467
(J.
/4
1. I~ARKS a n d J. 1). LEWINS
. . . Table 3. M c L a u r i n series approximations to
Po(t). for the doublet fission model when -~f-
= 1.4
t -Tf
Pl
P2
P3
P4
P5
P6
Po
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0.140000 0.280000 0.420000 0.560000 0.700000 0.840000 0.980000 1.120000 1.260000
0.123200 0.212800 0.268800 0.291200 0.280000 0.235200 0.156800 0.044800 -0.100800
0.125917 0.228779 0.322728 0.419029 0.529667 0.666624 0.841885 1.067435 1.355256
0.124960 0.224980 0.303495 0.358245 0.381267 0.358902 0.271792 0.094880 -0.202588
0.124988 0.225883 0.310355 0.387149 0.469476 0.578395 0.746202 1.019822 1.464187
0.124985 0.225668 0.307908 0.373404 0.417044 0.421834 0.351415 0.140612 -0.319137
0.124985 0.225709 0.308551 0.377834 0.436591 0.487014 0.530727 0.568955 0.602647
These coefficients in conjunction with equations (62) yield the numerical results in Table 4. These may be compared with the exact probability curve found by the method described in Section 5. These figures are taken from Table 5. These curves are shown in Fig. 3. Although less well-defined than in the 'doublet fission' case, it is apparent that the solutions diverge from around t -~ ~f. Unfortunately this useful time period is much less than that of interest to the nuclear engineer, so one is forced to conclude that the McLaurin series model o f f is not suitable after all. By studying an alternative means of deriving an equivalent McLaurin series representation of the doublet fission model, the problem of its convergence is made clearer. The work above seeks to modelffor Gl(x, t) =f(t +u(x)). Because of the simplicity of the doublet fission case, instead of defining the arbitrary constant ~ in equation (25) by O=t+~ln
~
,
(64)
as is effectively done here, one could redefine a new constant ~ by x-- 1 ( x - - 1)T = e (x-")~ = - e (x-")t x-#/k x-#/k '
so now
7((x-1)q
G,(x,t)= \ ~ ]
and
where
T = e Ix-")t,
G,(x,O)=x.
(66)
16 14 12 10
P1
Q, 0 8
'°5
06 ,q_
04-
~P_ o
0.2 0 t
0.1
0.2
0.3
0,4
--02
0.5 t/'rf
0.6
07
0.8 '~.hoO9
--04 F i g . 2. M c L a u r i n series a p p r o x i m a t i o n s to
Po(t) for
(65)
the 'doublet fission' m o d e l w h e n z r = 1.4.
An exact s o l u t i o n lor l o w - p o w e r n e u t r o n m u l t i p l i c a t i o n
Table 4. McLaurin series approximations to
/3
. . . ~f Po(t)for the empirical fission data case when - -
= 1.4
Tc
t -Tf
P~
P2
P3
P4
P5
P6
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
0.142700 0.285400 0.428100 0.570800 0.713500 0.856200 0.998900 1.141600 1.284300 1.427000 1.569700 1.712400
0.126703 0.221413 0.284130 0.314853 0.313583 0.280320 0.215063 0.117813 -0.011430 -0.172667 -0.365897 -0.591120
0.128129 0.232818 0.322621 0.406091 0.491782 0.588247 0.704040 0.847715 1.027825 1.252923 1.531563 1.872298
0.128032 0.231275 0.314809 0.398279 0.431507 0.463262 0.472489 0.452699 0.395085 0.288528 0.119593 -0.127470
0.128036 0.231402 0.315772 0.402334 0.443881 0.494052 0.539040 0.582451 0.628902 0.684499 0.757308 0.857833
0.128037 0.231406 0.315818 0.402595 0.444877 0.497026 0.546538 0.599157 0.662770 0.748227 0.870206 1.048123
Writing x--1
#/2=a a n d t a k i n g x o -- O, i.e.
Yo = 1/a, o n e
(67)
x-#/2'
finds that d " f y=yo dy"
So using a McLaurin
y
and
nI(-a)"+l
(68)
(I -a)"
series representation,
n=l
=
~
-a
.=1
(69)
(1 --a)"
gI!
(l-ay)"
(70)
(l-a)"
18
1 6
1 4
12
1 0
P
a
~ 06
0
022
0.0.4
006 t/.rfO 8
I 0
'12
--0.2
-o~
~
>k,k°
--0.6
Fig. 3. M c L a u r i n series a p p r o x i m a t i o n s to
Po(t) for
the e m p i r i c a l fission d a t a case w h e n Ze = 1.4.
76
G . T . PARKS and J. D. LEWINS
This is a geometric progression and, assuming it converges, (1 --ay) (
fm(y)= --a (~_a)
1--ay~ -1 1-- l - - a ]
(71)
or
fro(Y) = 11--ay y '
(72)
whence
GI,, - 1 - - a y T 1 -yT
where
T = e ~-~)'.
(73)
'
Substituting for a and y and rearranging gives #(x- 1)T-(2x-/~)
Gain(x, t) = 2(x-- 1)T--(2x-/.0 "
(74)
Thus, comparing equation (74) with equation (8), it is seen that the McLaurin series representation admits an exact solution of the doublet fission problem, providing it converges. It is convergent if the geometric progression converges, i.e. 1 -ay
1-a
< 1.
(75)
The region of interest is 0 ~< x ~< 1 which is equivalent to 1/a >>-y >t O. Thus if R = (1 -- ay)/(1 -- a), IRlmax =
when
(76)
y = 0.
(77)
So for complete convergence l la
< 1.
(78)
Therefore : if a < 1,
equation (78) ~ a < 0
for convergence;
if a > 1,
equation (78) ~ a > 2
for convergence. ~Cf
Thus the McLaurin series only converges over the complete range of interest if ~f < 0 or - - > 2. This Tc
Tc
Tf
unfortunately excludes most of the range of values o f - - of interest to the nuclear engineer, so again one concludes Tc
that the McLaurin series representation off is not suitable for this problem, even though the formal solution to the doublet model is correct.
5. SUCCESSFUL METHOD
Suppose in the analysis of Section 3 one writes, instead of ~ = f(~,), ~ = 0(~b). Then
t +u(x) = O(GO.
(79)
O(x) = u(x)
(80)
When t = O, G1 = x, therefore
An exact solution for low-power neutron multiplication
77
and
u(G,) = t +u(x)
(81)
or
A l n \ x( G _ lI - I h + ] B In ( ~ -) + C
( G l - a 2 ~ + D In ( ~ ) +- D * l n ( G ~ - a '\~x h- a *
In \ x - - a 2 /
J =t.
(82)
A l t h o u g h this expression c a n n o t easily be inverted to give G1 as a n explicit function o f x a n d t, it m a y be used to evaluate t as a function o f x a n d G,, a n d hence give Ga as a n implicit function o f x a n d t. In particular, setting x = 0 gives the extinction p r o b a b i l i t y Po(t) implicitly from
u(Po) = t + u(O).
(83)
Consider n o w the effect of differentiating e q u a t i o n (81) with respect to x:
63G1
if(GO ~-x = u'(x),
(84)
where du
A x- 1 +
u'(x) - dx N o t i n g t h a t ~~G1 5=0 =
B
C
x-a~- +
D -+ x-a2 x-a
D* + 3 x --a '~
(85)
P~(t), setting x = 0 in e q u a t i o n (84) gives u'(Po)P~ = u'(0).
(86)
It is s h o w n a b o v e t h a t Po m a y be found as an implicit function of time from e q u a t i o n (83), so Pa m a y also be found as an implicit function of time from u'(0)
(87)
P1 = u,(Po)"
By further differentiation with respect to x a n d setting x = 0 successive probabilities P , m a y be found as implicit functions of time. F o r example, P2 =
P3 =
u"(O)-- u"(Po)P ~ 2u'(Po)
,
(88)
u"(O)- 6u"(Po)P1P2 - u'"(eo)e 3 6u'(Po)
etc.,
(89)
where
d2u
du"
u"(x)=~-
and
u. . . .
(90)
dx
To illustrate this method, results are t a b u l a t e d below for Po, P1 a n d
P2. This is d o n e for b o t h the sub-critical case,
zf r~ = 1.4, a n d for a super-critical c a s e , - = 0.7. "6c
17c
F o r r f = 1.4, as before, Zc
a 1 -~ 0.970, "gf
A ~- 0 . 5 2 9 3 - - , qs ANE 12:2-D
a 2 ~ --3.538
B -~ --0.5374 zf q5
and
a 3 ~- - 0 . 9 6 8 + 3 . 1 6 3 j ~Tf
c ~ 0.0029--, q5
D -~ (0.0026--0.0012j) zf q5
78
G.T.
PARKS and J. D. LEwlrqs
Table 5. Exact probabilities for the general stochastic case when
Table 6. Exact probabilities for the general stochastic case when
z~ = 1.4
zf = 0.7
"~e t ( x 10-3 zf~ \
Po
"Co
Pt
P2
qs/
0.00 2.91 6.38 10.66 16.13 23.49 34.14 51.39 85.58 195.73
( t x l O azf_)
Po
P~
P2
0.00 5.92 13.56 24.05 39.99 69.41 185.47
0.000 0.100 0.200 0.300 0.400 0.500 0.600
1.000 0.793 0.597 0.418 0.257 0.120 0.012
0.000 0.042 0.073 0.086 0.081 0.054 0.008
qs/
0.000 0.100 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900
1.000 0.845 0.697 0.557 0.426 0.308 0.204 0.117 0.052 0.011
0.000 0.024 0.044 0.059 0.068 0.069 0.062 0.047 0.027 0.007
For z~ = 0.7, Te
al -~ 0.614, 17f
A -~ 0 . 0 4 9 0 - - , q5
a2 -~ - 3 . 2 4 0 ~'f
B -~ - 0 . 0 6 0 4 - - , q5
and
a 3 -~ - 0 . 9 3 8 + 2 . 9 5 4 j qSf
C -~ 0 . 0 0 4 4 - - , q5
D -~ (0.0035-0.0015j) re. q5
Substitution into equation (83) etc. gives the results in Tables 5 and 6. It is interesting to compare these results with the probabilities predicted by the 'quadratic approximation' model, where 2 and/~ [as defined by equation (11)] are chosen so that the resulting solution is exact to the second moment of n. The solution for G t [equation (8)] yields the following expressions for the probabilities under consideration : Po(t) = Gl(O,t)
#
1-- T
= 2- # / 2 - - 7 z
(91)
P~(t) = Oxx ~G1 x=o = (1--kt/)t)ZT/(~/2 - T) 2
(92)
and /~2(t)
= _l ~1
(93)
2 c~x2 x= o = (1 - #/2)2(1 - T ) T / ( # / 2 - T) 3,
giving the results in Tables 7 and 8. Table 7. Probabilities predicted by the 'quadratic approximation' for
Tf = ~c
Table 8. Probabilities predicted by the 'quadratic approximation' for
~f Te
1.4
-- =0.7
.o
..
0.00
0.000
1.000
2.91 6.38 10.66 16.13 23.49 34.14 51.39 85.58 195.73
0.153 0.283 0.396 0.497 0.588 0.672 0.752 0.831 0.910
0.714 0.508 0.358 0.246 0.162 0.101 0.056 0.024 0.006
0.000 0.112 0.148 0.146 0.126 0.098 0.070 0.043 0.021 0.005
t(× I0-3~) 0.00 5.92 13.56 24.05 39.99 69.4l 185.47
15o
Pl
P2
0.000 0.196 0.341 0.453 0.544 0.619 0.677
1.000
0.573 0.330 0.183 0.092 0.035 0.002
0.000 0.165 0.165 0.122 0.073 0.032 0.002
An exact solution for low-power neutron multiplication
79
These exact and approximate probabilities are plotted in Figs 4-6. These graphs show that the exact probabilities are indeed of the form one would anticipate from the physical considerations, but show a substantial inaccuracy in the probabilities predicted by the 'quadratic approximation' during the transient development. This is consistent with the discrepancies found in the fully-developed solutions by Campbell et al. (1981) for the subcritical case and Lewins (1981) for the super-critical case. [] Exact P1 rf /Tc = t 4 10 I
~
1011
• Quadratic approximation rf / ~ : 1 4
09
09
:: Exact P1 rf
o k
o
./,/
o,1
/r c =07
• Quadratic approximation rf /Tc = 0 7
o7
•
05 04
03 F/m~ / 02
'"/~/
O11~
03
. Quadratic approximation rf / ~ ----14
02
O Exact PO0 ~ / ~ = 0 7
y 0
D Exact PO ~/Tc ----14
• Quadrat c approximation ~ / ~ : 0 7 [ 25
I 50
I 75
l 100
I 125
for
z ~ = 1.4
"['¢
I 175
J 200
I
0
25
50
75
/qs)
t (X lO-3"rf
Fig. 4. Po(t)
I 150
0 1
t (xlO
and
014~1
for
[3 Exact P2 rf /Tc = 1 4
L
016V~
Fig. 5. Pl(t)
T~, =0.7.
1~c
O 20 I 018
\
• Quadratic approximation Tf / r c ~ 14 o Exact P2 rf / r c = 0 7 • Quadratic approximation r f / ' r c : 0 7
/
oo i ..,x " oo ![ O0 2
0
100
25
50
75
100
125
t (xlO-3"q
Fig. 6. P2(t)
for
~rr= 1.4 "Ce
150
175
/qs) and
Zf=o.7. "£c
200
__
125
150
rf
/qs)
--3
~f = 1.4 Tc
and
"L'f
175
--=0.7. Tc
200
80
G . T . PARKS and J. D. LEW1NS 6. CONCLUSIONS
The success in finding the exact neutron probabilities for the problem based on empirical fission data for a source-free system is undoubtedly due to avoiding the need to find the generating function explicitly; the unsuccessful methods having failed to do this. An explicit expression for G1 is required, however, to evaluate the other basic generating function F o through the convolution integral. Unless F o is available, the system with completely general initial conditions remains insoluble exactly. Thus further work must be directed towards finding a successful means of inverting the relationship for G 1 found here. However, the implicit method we have given could be the basis of a numerical convolution. The implicit method we have given here is straightforward and lends itself to numerical computation with a digital computer. The resulting numerical values demonstrate the weakness of the quadratic approximation for this model of prompt-neutron stochastics in the transient period before the establishment of the well-developed or asymptotic solution, a conclusion that reinforces the results previously obtained for asymptotic solutions. REFERENCES
Abramowitz M. and Stegun I. A. (1964) Handbook of Mathematical Functions. NBS, Washington, D.C. Bartlett M. S. (1979) An Introduction to Stochastic Processes, 3rd edn. Cambridge Univ. Press, U.K. Bell G. 1. (1963) Ann. Phys. (N.Y.) 21,243. Campbell G. A., Chan V. and Lewins J. D. (1981) Ann. nucl. Energy 8, 1. Divert B. C., Martin H. C., Taschek R. F. and Terrel J. (1956) Phys. Rev. 101, 1012. Gurland J. (1978) A Seminar on Statistics. CNEN, Rome. Harris D. R. (1964) Naval Reactors Physics Handbook, Vol. 1 (Edited by Radkowsky A.). USAEC, Washington, D.C. Harris T. E. (1963) The Theory of Branching Processes. Springer-Verlag, Berlin. Jeffreys H. and Jeffreys B. S. (1956) Methods of Mathematical Physics, 3rd edn. Cambridge Univ. Press, U.K. Kendall D. G. (1948) Ann. math. Statist. 19, 1. Lewins J. D. (1978) Ann. nucl. Energy 5, 141. Lewins J. D. (1981) Ann. nucl. Energy 8, 15. Pacilio N., Colombino A., Mosiello R., Norelli F. and Jorio V. M. (1978) Adv. nucl. Sci. Technol. II, 67. Smith P. R. (1970) J. nucl. Energy 24, 199. Williams M. M. R. (1974) Random Processes in Nuclear Reactors. Pergamon Press, Oxford. Williams M. M. R. (1979) Ann. nucl. Energy 6, 43.