- Email: [email protected]
Method for designing the cell of a heterogeneous reactor
METHOD FOR DESIGNING THE CELL OF A RETERUGENED~SREACTOR* KU~~~TS~V
E. S,
MO8OOW (Received
7
December
1964)
ONEof the basic problems in the theory of nuclear reactor design is that of calculating the cell of a heterogeneous reactor (see e. g. fI] 1. To solve this problem, various approximate methods have been developed, most commonly based on P,- and P,- approximations (see [21). It has recently been discavered that more accurate methods are needed, baaed on approximate solution of the kinetic equation. Such methods include jh;)+o-called Srmethod of Carlson, and V.S. Vladimirov’ s method (see
We shall assume that the statement of the problem is based on isolation of the cell in the form of an infinite __ __circular cylinder in acoordante with the Wigner - Seitz method ( [ll , [21). In this case, in the one-group approximation and on the assumption that the neutron scattering is isotropic, the kinetic equation may be written in the form
where A&r, y, rr1 is the neutron phase density, f(r) is the neutron source distribution along the cylinder radius, a(r) = E,(r) + X,frf and 1, are the macroscopic scatter and absorption cross-sections), P(r) = Es(r), r Is the radius vector of a point in the cylindrical
l
2,.
Yychisl.
Mat. mat. Fiz. 5, 3, 488 - 599, 1989, 128
(1,
The cell
of
a heterogeneous
reactor
129
coordinste system, y = cos 5, where 5 is the angle between the neutron velocity and the cylinder generator, and I.I= cos Q, Whre Q is the azimuth angle. The
neutron velocity
is taken equal to unity.
For a complete solution of the problem we also require the following definitions: 1. the neutron density n(r) = which
can
s N(r* 699) dJa*
be transformed to
or
(2) 2. the neutron flux density F(r)=
5 N(r,i2)cos(v,~)dQ,
where v is the normal vector to the surface element dQ, *positive” side of the element. We can write F(r) as
F(r) =
F+(r)
-
defining the
F-(r),
where
F+(r) = F-(r)
=
SN(r,
1 N;r. +
-4)
dl)cos(v, IWS(V,
Q)dQ, 42)
I&.
The integration is performed in both cases over the positive hemisphere of unit radius. Since the angle between v and Q is obtuse in the second case, to achieve the condition P(r) > 0 we have to take the absolute valUe of COS(V, - P).
The relationship
130
E.S.
Kuznetsoc
F(r)=!j N(r,
g, cp)coscpsin2~d@+
0 --n
or
is also easily
established
for F(r).
For a strict statement of the problem there is some lack of definiteness in the choice of the boundary condition on the outer boundary of the cell. It is assumed in diffusion theory that as many neutrons leave any given cell as enter it. This condition implies that the neutron flux density is zero on the outer boundary of the cell. In terms of kinetic theory, it can be written as F(R)
= 0
or F+(R)
-F-(R)
= 0.
(3)
We do not usually have definite enough information on the angular distribution of the neutrons incident on the cell face from outside. We shall start from the very simple assumption that the neutron emission incident on the cell boundary from outside has an isotropic angular distribution. It is interesting to observe that communications have lately appeared in the American literature (see e.g. [31) in which criticism is made of a condition on the Wigner - Seitz cell boundary involving the reflection of neutrons in accordance with Snellius’s law, which has fairly recently come into use. In the case of thin moderators this form of the boundary condition leads to a strong rise in the disadvantage factor El], which would seem to be bound up with the artificial isolation of the cell in accordance with the Wigner - Seitz method. Some authors (see [31) have proposed using an isotropic distribution of the neutrons reflected from the cell boundary, which might be achieved by adding an optically extremely thick outer lager consisting of purely scattering material. This assumption leads to the formula
F-(R)= d?,
The cell
of
a heterogeneous
reactor
131
where H is the as yet ~~0~ intensity of the neutron radiation incident from outside. We obtain from (3) and the lctst formula the boundary condition H = ; F-+(R)
(4)
Bearing in mind the application of the method to problems of optics, we shall in future use the rather more general condition
H=
p-F+(R), n:
where A is a given constant, satisfying the inequality 0 < A < 1. Obviously, A denotes the albedo. For A = 1 we obtain our boundary condition (4). We now perform the following
As a result,
substitution
in equation (1):
(1) becomes
where
(8)
132
E.S.
The solution
Kutnetsov
of (7) must obviously m(r> Y, P> =
flak
be a linear
function of H
Y, 4.4 + ERi (f, y, p) .
(9)
Substituting this in (7) and equating separately the terms in which H is present, and absent, on the right- and left-hand sides, we arrive at a system of two integro-differential equations
which only differ in the first term on the right-hand sides. We shall in future describe the first of these equations as the fundamental, and the second as the auxiliary, equation. Equation (7) has to be solved under the usual boundary condition, according to which no neutrous are incident on the surface of the cylinder fro? outside. For, with r = R, t_he second of relationships (61 yields H = MR., y. p) I- H, whence i’V(R, y, cl) = 0 ( - 15
CIS 01.
On substituting for fl from (9) in this, u so,. H8,(R. y, ~1 = 0 (-1s Since &,(R,
ye 14 > 0, k,(R,
we have h$(R,
y, p) +
y, p) > 0, H > 0 for any y and p, it
follows from this that &,,(R. y, IA) = 0, k,(R, y, IJ) = 0 ( - 12 ~1(I 0). Hence equation (10) must also be solved under the assumption that no neutrons are incident on the cylinder surface from outside. Let us turn to determining the constant H. For r = R the first
of
The
relationships
cell
of
a heterogeneous
133
reactor
(6) yields
Hence, using the definition P(R)
of F’(R), =
P+(R)
we get -I- HT,
where
‘&e quantity W(Ii)
F+(R) in (II)
= F&3+(IS)+ ~~~~(~),
needs no explanation. Substituting
itself
depends on H. As is easily seen,
where a new notation is introduced which
Consequently, F+(R) = Fe+(R) + ff [T + PI+(R)].
this in (5). we obtain the equation for H H =~[~~+(R)+HIT+~,+(R)ll,
whence
Since we have for the problems (10) Fop(R) = pi-(B) PO+(R) = PO(R)
snd F*+(R)
In the case of a cell which case
= Pi(R),
calculation,
= 0,
and we finally
then
get
we have to put A = 1 here, in
The above implies the following order for solving the problem of cell calculation:
E.S. Kuznetsov
134
1. evaluation of the function O(F), appearing in the first term of the right-hand side of the auxiliary integro-differential equation. by using (8). We have developed for this purpose a method based on the employment of special quadrature formulae, and have designed a suitable electronic computer programme; 2. approximate solution
of the integro-differential
equations (10)
for the ~~0~ functions iii;(r, 5, &.t) and R,( r, 5, u). This represents the fundamental stage of the work. It is carried out by using a special finite - difference method, based in its essential features on the familiar Carlson scheme; a computer programmewas designed permitting high accuracy of working. Apart from fiO and %, themselves, this progrsmmealso gives the values of &,(r),
and
&(F)
F*(F).
Fl(r);
3. evaluation of the constant H from (13) (in particular, for A = 1, from (14)). The relevant programmecombines evaluation of the integral T from formula (12)) of the constant H, and of the function n(r) from (15); 4. evaluation of the neutron density n(t) after substituting expressions (6) in (2) n(r)
=
from the formula obtained
R (f) + 4nHo (r)
or
n(r)
=
20(r)
+ W%(r)
+ ~TW.J(F)],
(15)
where
5. evaluation of the angular distributions by (9), become
from formulae (6).
which,
The
cell
of
a heterogeneous
135
reactor
+exp[.A__j a(s)sds ]J )/I-$ ,)/s”-rql-p2)
(-1
Q cl.<
0).
Some particular cases of these formulae may be noted. For F-= 0 we get the angular distribution on the cylinder axis. It must obviously be dependent only on y, and not on p. Both formulae here give the seme result \-
NW,Y)=mp(o,y)+H[m*(o,Y,+*=P~ For
&
[We)]
(W
get
r = R we
a(s)s ds
- @(I - P2) N(&Y,P) ==H
(-1 <
I1
to<
CL<
11,
(17)
G 0).
In an approximate calculation from (16) and (17) we obtain the values of the functions only at discrete points yj, Tie where yj denotes the Gaussian base-point number j. Since the point y = 1 (3 = 0) is not included in the Gaussian base-points, a special calculation has to be carried out for it on the basis of theoretical considerations. For y = 1 equation (1) gives
W, 1, ~1)= -
PC’) fia(r)
* sId,,s
N(r,
0
_i
Y,
II-
cL)dcL + f(r)
CL2
or
W, 1, IL)= In particular,
fWW,
f(r) a(r) *
for r = 0 we get
NOA 1,~) = B(O)nW 4na (0)
or, by (15),
I
4na (r)
1
f(O)
40) ’
a(r) ’
136
E.S.
NW, I, CL) =
*i$,
Kuznetsov
+s
{~o(O)+H[~i(O)+4~o(0)1)
The quantities appearing known from calculations.
in this
formula are either
given,
or are
The method described can be extended to the multigroup cell calculation. We base the solution of the problem in this case on the system of equations
II-r”( Ir~+qf-f$)
+%(r)Ng@,y, CL)=Q&)
(18)
(g = 1,2,..., G). Here g is the number of the group, G the total number of groups, ag(r) is the macroNg(r, y, p) is the group g neutron phase density, scopic absorption cross-section, allowing for scattering of all types from group and capture, pgg ‘(r) is the matrix of neutron transitions g’ to group g, allowing for elastic and non-elastic scatter. Since no restrictions, apart from the condition that they be positive, are imposed on the elements of the matrix pgg’, it takes account of transitions both from higher groups to lower, and from lower to higher (we shall assign a higher number to a group with higher neutron energy). This enables us to introduce into the problem processes bound up with thermolizations
(g= 42, . ..t G), where
n,(r)=&y jNB(W,P) dl.c v - Y2 -i
(g=1,2,:..,G)
0
is the group g neutron density, of neutrons of the same group. The boundary conditions
h’g(f4 Y, cl) = where the
Hg are
constants
for Hg
and fg(r)
(18)
is the given
source
intensity
are
&=I,2
,...,
G; p<(I),
depending only on the group number, and R is
The
the cylinder
of
cell
* heterogeneous
13:
reactor
radius.
By analogy
we shall
with (51, 4 HR = ;
assume that the conditions
Fs+(R)
(g =
are also satisfied in the multigroup signifying the albedo and satisfying P,+uq=
1, 2, . . . , G),
where the A, are constants
case,
0 4 A, 4 1 (g = 1, 2, . . . , 0
~)cos(~,s-~)~sz
~&vg(z-q
(1%
and
(g=1,2,.:.,G).
+ As distinct from the one-group case, in which the neutron velocity can be taken as unity, we cannot neglect in the multigroup case the different group velocities eg 1, 2, . . . , 0. We therefore have to pass to a general definition of the neutron flux, corresponding to the group number g
(g=
F,(r) = ggS445
ww~, wm
(g=
1,2, . . . . G).
+ This conclusion nected with F,(F) F,(r)
obviously by =
Fg+(r)
In the case of cell be put equal to unity.
also -
design
applies
to F,+(r)
all
, con-
1,2, . . .. G).
(g=
Fg-(r)
and F,-(F)
the A, in boundary condition
(19) must
Using the substitution
N,(r,
Y,
r
___ 1 -VI- Y2[ s1s” 1
~1)= mg(r, Y, CL)+ Hg exp - -
ag(s)sds
rJg
g=1,2,...,G),
(O
iJrg(c yI CL)= flg(r, Y, cl) + Hg exp
(-1 GpcLOo; we reduce
the fundamental
1
1
-
___ 11 -
g=
system (18)
+ +(I - 11”)
1 y2sTl/s”- rqa- p2)‘. R
ag(s) sas
1,2, . . . . G),
to the form
(20)
138
E.S.
Kurnetsov
where
Further, we put
(g=1,2
,...,
G).
(22)
g’=l obtain, for determining the unknownfunctions md fl,+, g(r, y, p) C + 1 systems of integro-differential We now
fio,g(r, y, p) equations
(g=GL..,G), +: $Br i ipo,zPAw+-2 4
_P
~l-V,(p~+~-ff$f-)+ug(r)iVg~,g(r,p,y)=
==. (g’=
For the functions n,(r) tions
I,& . . . . G; g=
1,2, . . . . G).
and the constants F,+(R), we have the equa-
The
cell
of
a heterogeneous
%tr) = b(r) + 4nHgogb), tg =
139
reactor
Fg+(R)
=
F,+(R) +
H,Tg
1, 2, . . ., 9,
where
Using relationships
ng(r)=
(22),
fio,g(r)+
we obtain
from this
~Hgriig.,g(r)+4I1HgWg(T)
(g =
1,2,.
. . , G)
g’=i
and
Fg+(fi)
= Fbtg tfi) + 5
g’=i
H,J%
g(R) + HgTg
(g=
I,?
. . . . G).
In order to find the Hg, it only remains for us to substitute the boundary conditions (19) and solve the resulting system
H,=--for
Ag [b(R)+i H&,,(R)+&T,] al g’?=i
(23)
(23)
in
(g = 1,2, . . . , G)
the H,. For example,
in the case of
two groups we have the system
It might be possible to start from the hypothesis that the neutron radi“reflected” from the cell boundary, is independent of the neutron ation, velocity, i.e. independent of the group number. In this case the constants Hg in (20) and (21) would have to be replaced by the same constant H, and instead of (22) we should write K(r, As a result,
Y, I4 =
we obtain
No,&,
Y, cl) + Hmi,s (r, Y,
two systems of
integro-differential
P). equations
E.S. Kurnetsov
140
(g=
172,. . ..G).
(g =
1, 2, . . .,
Furthermore, we hsve
n&T(r)=
50,&J
+
H[Ri,g(r)
+ 4nog(r)]
G)
and
Fg+(JO =
Fo,~+(R)
= H{&,,+(R) + TV} (s = b2,. . ., G).
We denote by A, the size of the velocity the g-th group. We now have from (24)
interval
(24)
corresponding to
Putting
i; “FB-(R)Ag = -&$ Fg+(R)Ag = H, g=i g--1 where 0 < A S 1, we obtain for H the formula
1-
$2g=1 {h*,(R)+ T,) A,
Calculations using the method described were csrried out for a number of examples in the case of the one-group problem. The solutions were obtained by computer, using special programmes. We shall not dwell here
The
cell
of
reactor
a heterogeneous
141
on the approximate method of solving the neutron transfer equation cl), on which the compilation of the programmeswas based, since this is of no importance for the purposes of the present paper. It is sufficient to say that the approximate solution was based on an improved version of Carlson’ s-SN-method (see [Zl). Experience of computer computation shows that evaluation of the function w(r) for all Fk from formula (121 takes up a great deal of time. The evaluation of this function, and of the interval T, can be avoided if the problem is not reduced to zero boundary conditions*. For this purpose we’ put in (If
N(r7Y7 PL)=
--
wr>Y, CL)-I- HNl(F, Y, cl),
where &,(r,
y, 1.4 is the solution
of (X) under the zero boundary condi-
tion,
y, cl) is the solution
of (1) for f(r)
N,(r.
E 0, but under the
boundary condition Nr(fi, y, IA = 1 for directions and g is an unknownconstant.
inward to the cylinder,
As a consequence,
F(R) = &l(R) + ml(R)
= 0
and
IT Since
pression
PO(R)
=
z--z
Fe(R)
(14) for i,
FofR) Fi (4
ax-Pi+(R)
and Pi+(R)
= Pi(R) i- T,
PO(R) we arrive at ex-
where i = H.
It must be mentioned, however, that the accuracy of this method depends on the fineness of the mesh used in the finite - difference part of the computation, whereas the direct evaluation of w(r) and T involves only the accuracy of the quadrature formulae employed and can be realized independently of the choice of mesh. A similar modification of the method is likewise.possible in the multigroup case. We shall first consider two examples concerned with the one-zone problem with a constant solution known in advance (these examples were selected by V.A. Chnyanov from considerations bound up with the diffusion approximation). *
This remark follows of
reducing
from discussion with A.V. Voronkov as to means the volume of machine computations.
142
iiuznetsou
E.S.
In this case the solution of the problem for an infinite circular cylinder can be reduced to solution of the integral equation
0
0
where K(r, p) is a kernel whose concrete form will not be quoted here (see e.g. 141). In the case of constant f and p, we have R
R ?z(r),= 4nffofr)+W
5 0
W,P)P%fB~
(25)
~(P)~(~,P}P~P* 0
is easily shown that, if the constants H, f, o, and p are in a certain ratio, this equation has a constant solution n(r) = 4df. For, on substituting this value of n(r) in the right-hand side of (25), we get
It
P~PdP+4nw~
n(r) = 4nHo(r)+4nfjKo-,
w,
P)PdP.
0
0
We now use the relationship R c&p,
s
p)pdp
=
1-
w(r).
0
which we established for a = 1.
in another paper (see [41, formula (2.21)):
Hence,
R s
fl(r,p)pdp
=
1 -o(r).
0
Now,
n(r) = 4n(f + flH) + 4n(fZ - f - pH)o(r).
be constant,
we have to satisfy
In order for n(r) to
the condition H = f t pH,
whence
H=
f/Cl - p,. For instance, with f = 1, p = 0.9, we get H = 10. Notice that this result is independent of the choice of cylinder radius R. A concrete calculation was performed for R = 1 and for R = 3, using the’ integro-differential equations (10) and further formulae relating to the cell calculation in the one-group approximation. In both cases the following values were taken for the physical parameters: a = 1.0, p = 0.9, f = 1. The finite - difference method used for solving equations (IO) was characterized by the following data. The interval 0 S r < R was divided in both cases into a hundred equal parts. Of these
The
ceZZ
of
reactor
a heterogeneous
143
hundred intervals, the third, counting from the cylinder boundary, was in turn divided into five equal parts, the second into ten equal parts, and the first into ten equal parts. Thus there were altogether K = 122 unequal intervals. The interval - 14 p 61 was divided into I = 64 equal intervals. As regards the variable y, its interval 0 < y < 1 was characterized by J = 7 Gaussian base-points. The computation was performed by successive approximations with a calculation such that the remainder term of the iterational process did not exceed E = 10e4. The results of the computations were with
R=l,
with
R =
H = 10.0227597; H = 10.0137399.
3,
We took the following group Of examples on cell calculation The initial physical data are contained in Table 1. TABLE 1. Physical
from [51.
parameters -
Jvb
R
RI
CM
PI
a1 c&t-’
CM
CM-’
k I 5
1.30 1.30 1.30 1.50 1.70
0.7221 0.7221 0.7221 0.7221 0.7221
11.90 11.28 7.93 7.93 7.93
0.3991
0.3991 0.3991 0.3991 0.3991
0.63721 0.3721 0.3725 0.3925 0.3925
0.3717882 0.3717882 0.3921882 0.3921882 0.3921882
4.88299 4.652288 3.541005 3.606925 3.672845
In all five examples it is a question of two-zone cells, consisting of an interior uranium block of radius RI and an outer moderator zone of radius R2. The “optical radius” r = alR1 + a2(R2 - RI) is given in the last column of the table. Table 2 gives some data concerning the organization difference computation of the five versions mentioned TABLE 2. NO. Jv"P
; 1 5
K
116 116 116
I
64 lf
J
:: 7
Finite
- difference
Of
NO.
I--
intervals in zone I
;: ;: 34
Of
-
meshes
t--
intervals of zone I1
E 75 70 66
of the finite above.
Remarks In all versions the last two intervals of zOne I and the first and last intervals of zone 11, were subdivided into five parts each
144
Kuznetsov
E.S.
TABLE 3. Computational
results -
n,
e
H
&I
0.86692
i EEE 199:2125 160.7575 2 133.5051 5
0.87980 0.94051
0.95196 0.96013
1.9208 1.9052 1.8096 1.940 2.0722
1.891 1.877 1.791 1.931 2.070
1.898
The computational results are given in Table 3, which includes, in addition to H, the thermal efficiency 8, calculated from the formula
and the “disadvantage
factor”
ll = &/nl
RI
where
s
s
n(r)&
0
Ri2
the thermal neutrons, RI
n(r)rdr
iii =
for
’
52 = ;i2 _ Ri2 *
The last column but one gives the values of the disadvantage factor calculated by Amouyal et al. (see [51) by a special method, different from ours, while the last column gives the values of the disadvantage factor computed by Carlson in the Se-approximation, for variant No. 1. A comparison of our calculations with the data of Amouyal reveals that our figures are somewhat higher while at the same time remaining fairly strictly parallel.
The
cefl
of
reactor
a heterogeneous
145
In addition, we calculated the neutron density distribution as a function of the radius r for all five versions. The results are shown in the figure, which gives the ratio n( t-)/n(O). The zone boundaries in the individual versions are marked by vertical lines, with the corresponding figures on the abscissa axis. The last example, set up by L.V. Maiorov, is a three-zone cell and is characterized by the following data: 0 < r < R, = 1.75. (natural uranium); f - 0, a1 = 0.62046, p1 = 0.427; R, = 1. 95 (aluminium) If- =O, a2 = 0.098474, pz= 0.08442; Rl = 1.75
Angular distributions
5
3
2
6
I,
222.78 22%:: 217:38 212.84 206.47 197.73 185.6G 168.46 142.34 107.21
202.33 201.85 200.26 197.49 193.37 187.60 179.67 168.72 153.09 129.36 98.098
103.24 :FZ! 1oO:so 98.733 95.833 91.838 86.293 78.344 66.194 49.758
74.815 74.599 73.938 72.796 71.102 68.739 65.518 61.118 54.974 46.067 35.618
As a result of the computations we obtained
55.822 55.642 55.096 54.139 52.734 50.788 48.162 44.629 39.820 33.190 26.279
H = 9.60955629,
5.8818 5.8707 5.8371 5.7735 5.6911 5.5682 5.3984 5.1625 4.8247 4.3157 3.6654 0 =
0.8139390.
The variation of the neutron density (normalized as in the previous examples) is shown in curve 6 of the figure. In addition to the cell characteristics enumerated above, we computed the angular distributions of the phase densities on the cylinder axis, which are of interest to experimental workers. As is easily seen, on the cylinder axis the phase densities depend only on the variable y = cos 5,
146
Kuznetsou
E.S.
which considerably simplifies the computational the computations are given in Table 4.
formulae.
The results
of
Analysis of these figures reveals that the shape of the curves N = N(y) give a good reflection of the cell properties. In particular, an extremely sensitive characteristic of the physical and geometrical properties of the cell is provided by the ratio of the values of the neutron phase density corresponding to a direction along the cylinder axis, and corresponding to a direction perpendicular to the axis: N(O)/%(l). This ratio sidered:
has the following
values
in the examples that we have con-
TABLE 5 No. of
variant
2
1
2.0630 --
2.0779
~(O)/~(l)
3
4
2.0749 -
2.1004
5
6
2.1242
1.6047
The author wishes to thank A. I. Vaskin, E. I. Acknowledgements: Panfilov, 0. V. Sysoev, E.G. Natrusov and G.P. Churkin for cooperation in the computer programming and in carrying out the computations, and V.A. Chuyanov and L. V. Maiorov for preparing the individual physical examples and for participating in the discussion of the results. TransEated
by D.E. Brown
REFERENCES 1.
GLASSTONE, 8. snd EDLUND, M.C. The theory, Van Nostrand.
2.
MARCHUK. G.I. Methods of calculation of nuclear reactors rascheta yadernykh reaktorov), Moscow, Gosatomizdat,
3.
4.
POMRANING,G.C. The Wigner - Seitz calculational method, Nucl. Sci.
of
cell,
nuclear
reactor
(letodg 1961.
a discussion and a simple 17. 3, 311 - 314, 1963.
Engn.
KUZNETSOV. E.S. On the temperature distribution in an infinite cylinder and in a sphere with non-monochromatic radiant equilibrium Xh.
5.
elements
vpchisl.
Hat.
mat.
Fiz.
2,
2,
217
-
240,
1962.
AMOUYAL. A., BENOIST. P. and HOROWITZ, J. Nouvelle m&hode de deftermination du facteur d'utllisatlon thermique d'une cellule, J. nucl. Energy, 6, l/2. 79 - 98. 1957.
E. S,
MO8OOW (Received
7
December
1964)
ONEof the basic problems in the theory of nuclear reactor design is that of calculating the cell of a heterogeneous reactor (see e. g. fI] 1. To solve this problem, various approximate methods have been developed, most commonly based on P,- and P,- approximations (see [21). It has recently been discavered that more accurate methods are needed, baaed on approximate solution of the kinetic equation. Such methods include jh;)+o-called Srmethod of Carlson, and V.S. Vladimirov’ s method (see
We shall assume that the statement of the problem is based on isolation of the cell in the form of an infinite __ __circular cylinder in acoordante with the Wigner - Seitz method ( [ll , [21). In this case, in the one-group approximation and on the assumption that the neutron scattering is isotropic, the kinetic equation may be written in the form
where A&r, y, rr1 is the neutron phase density, f(r) is the neutron source distribution along the cylinder radius, a(r) = E,(r) + X,frf and 1, are the macroscopic scatter and absorption cross-sections), P(r) = Es(r), r Is the radius vector of a point in the cylindrical
l
2,.
Yychisl.
Mat. mat. Fiz. 5, 3, 488 - 599, 1989, 128
(1,
The cell
of
a heterogeneous
reactor
129
coordinste system, y = cos 5, where 5 is the angle between the neutron velocity and the cylinder generator, and I.I= cos Q, Whre Q is the azimuth angle. The
neutron velocity
is taken equal to unity.
For a complete solution of the problem we also require the following definitions: 1. the neutron density n(r) = which
can
s N(r* 699) dJa*
be transformed to
or
(2) 2. the neutron flux density F(r)=
5 N(r,i2)cos(v,~)dQ,
where v is the normal vector to the surface element dQ, *positive” side of the element. We can write F(r) as
F(r) =
F+(r)
-
defining the
F-(r),
where
F+(r) = F-(r)
=
SN(r,
1 N;r. +
-4)
dl)cos(v, IWS(V,
Q)dQ, 42)
I&.
The integration is performed in both cases over the positive hemisphere of unit radius. Since the angle between v and Q is obtuse in the second case, to achieve the condition P(r) > 0 we have to take the absolute valUe of COS(V, - P).
The relationship
130
E.S.
Kuznetsoc
F(r)=!j N(r,
g, cp)coscpsin2~d@+
0 --n
or
is also easily
established
for F(r).
For a strict statement of the problem there is some lack of definiteness in the choice of the boundary condition on the outer boundary of the cell. It is assumed in diffusion theory that as many neutrons leave any given cell as enter it. This condition implies that the neutron flux density is zero on the outer boundary of the cell. In terms of kinetic theory, it can be written as F(R)
= 0
or F+(R)
-F-(R)
= 0.
(3)
We do not usually have definite enough information on the angular distribution of the neutrons incident on the cell face from outside. We shall start from the very simple assumption that the neutron emission incident on the cell boundary from outside has an isotropic angular distribution. It is interesting to observe that communications have lately appeared in the American literature (see e.g. [31) in which criticism is made of a condition on the Wigner - Seitz cell boundary involving the reflection of neutrons in accordance with Snellius’s law, which has fairly recently come into use. In the case of thin moderators this form of the boundary condition leads to a strong rise in the disadvantage factor El], which would seem to be bound up with the artificial isolation of the cell in accordance with the Wigner - Seitz method. Some authors (see [31) have proposed using an isotropic distribution of the neutrons reflected from the cell boundary, which might be achieved by adding an optically extremely thick outer lager consisting of purely scattering material. This assumption leads to the formula
F-(R)= d?,
The cell
of
a heterogeneous
reactor
131
where H is the as yet ~~0~ intensity of the neutron radiation incident from outside. We obtain from (3) and the lctst formula the boundary condition H = ; F-+(R)
(4)
Bearing in mind the application of the method to problems of optics, we shall in future use the rather more general condition
H=
p-F+(R), n:
where A is a given constant, satisfying the inequality 0 < A < 1. Obviously, A denotes the albedo. For A = 1 we obtain our boundary condition (4). We now perform the following
As a result,
substitution
in equation (1):
(1) becomes
where
(8)
132
E.S.
The solution
Kutnetsov
of (7) must obviously m(r> Y, P> =
flak
be a linear
function of H
Y, 4.4 + ERi (f, y, p) .
(9)
Substituting this in (7) and equating separately the terms in which H is present, and absent, on the right- and left-hand sides, we arrive at a system of two integro-differential equations
which only differ in the first term on the right-hand sides. We shall in future describe the first of these equations as the fundamental, and the second as the auxiliary, equation. Equation (7) has to be solved under the usual boundary condition, according to which no neutrous are incident on the surface of the cylinder fro? outside. For, with r = R, t_he second of relationships (61 yields H = MR., y. p) I- H, whence i’V(R, y, cl) = 0 ( - 15
CIS 01.
On substituting for fl from (9) in this, u so,. H8,(R. y, ~1 = 0 (-1s Since &,(R,
ye 14 > 0, k,(R,
we have h$(R,
y, p) +
y, p) > 0, H > 0 for any y and p, it
follows from this that &,,(R. y, IA) = 0, k,(R, y, IJ) = 0 ( - 12 ~1(I 0). Hence equation (10) must also be solved under the assumption that no neutrons are incident on the cylinder surface from outside. Let us turn to determining the constant H. For r = R the first
of
The
relationships
cell
of
a heterogeneous
133
reactor
(6) yields
Hence, using the definition P(R)
of F’(R), =
P+(R)
we get -I- HT,
where
‘&e quantity W(Ii)
F+(R) in (II)
= F&3+(IS)+ ~~~~(~),
needs no explanation. Substituting
itself
depends on H. As is easily seen,
where a new notation is introduced which
Consequently, F+(R) = Fe+(R) + ff [T + PI+(R)].
this in (5). we obtain the equation for H H =~[~~+(R)+HIT+~,+(R)ll,
whence
Since we have for the problems (10) Fop(R) = pi-(B) PO+(R) = PO(R)
snd F*+(R)
In the case of a cell which case
= Pi(R),
calculation,
= 0,
and we finally
then
get
we have to put A = 1 here, in
The above implies the following order for solving the problem of cell calculation:
E.S. Kuznetsov
134
1. evaluation of the function O(F), appearing in the first term of the right-hand side of the auxiliary integro-differential equation. by using (8). We have developed for this purpose a method based on the employment of special quadrature formulae, and have designed a suitable electronic computer programme; 2. approximate solution
of the integro-differential
equations (10)
for the ~~0~ functions iii;(r, 5, &.t) and R,( r, 5, u). This represents the fundamental stage of the work. It is carried out by using a special finite - difference method, based in its essential features on the familiar Carlson scheme; a computer programmewas designed permitting high accuracy of working. Apart from fiO and %, themselves, this progrsmmealso gives the values of &,(r),
and
&(F)
F*(F).
Fl(r);
3. evaluation of the constant H from (13) (in particular, for A = 1, from (14)). The relevant programmecombines evaluation of the integral T from formula (12)) of the constant H, and of the function n(r) from (15); 4. evaluation of the neutron density n(t) after substituting expressions (6) in (2) n(r)
=
from the formula obtained
R (f) + 4nHo (r)
or
n(r)
=
20(r)
+ W%(r)
+ ~TW.J(F)],
(15)
where
5. evaluation of the angular distributions by (9), become
from formulae (6).
which,
The
cell
of
a heterogeneous
135
reactor
+exp[.A__j a(s)sds ]J )/I-$ ,)/s”-rql-p2)
(-1
Q cl.<
0).
Some particular cases of these formulae may be noted. For F-= 0 we get the angular distribution on the cylinder axis. It must obviously be dependent only on y, and not on p. Both formulae here give the seme result \-
NW,Y)=mp(o,y)+H[m*(o,Y,+*=P~ For
&
[We)]
(W
get
r = R we
a(s)s ds
- @(I - P2) N(&Y,P) ==H
(-1 <
I1
to<
CL<
11,
(17)
G 0).
In an approximate calculation from (16) and (17) we obtain the values of the functions only at discrete points yj, Tie where yj denotes the Gaussian base-point number j. Since the point y = 1 (3 = 0) is not included in the Gaussian base-points, a special calculation has to be carried out for it on the basis of theoretical considerations. For y = 1 equation (1) gives
W, 1, ~1)= -
PC’) fia(r)
* sId,,s
N(r,
0
_i
Y,
II-
cL)dcL + f(r)
CL2
or
W, 1, IL)= In particular,
fWW,
f(r) a(r) *
for r = 0 we get
NOA 1,~) = B(O)nW 4na (0)
or, by (15),
I
4na (r)
1
f(O)
40) ’
a(r) ’
136
E.S.
NW, I, CL) =
*i$,
Kuznetsov
+s
{~o(O)+H[~i(O)+4~o(0)1)
The quantities appearing known from calculations.
in this
formula are either
given,
or are
The method described can be extended to the multigroup cell calculation. We base the solution of the problem in this case on the system of equations
II-r”( Ir~+qf-f$)
+%(r)Ng@,y, CL)=Q&)
(18)
(g = 1,2,..., G). Here g is the number of the group, G the total number of groups, ag(r) is the macroNg(r, y, p) is the group g neutron phase density, scopic absorption cross-section, allowing for scattering of all types from group and capture, pgg ‘(r) is the matrix of neutron transitions g’ to group g, allowing for elastic and non-elastic scatter. Since no restrictions, apart from the condition that they be positive, are imposed on the elements of the matrix pgg’, it takes account of transitions both from higher groups to lower, and from lower to higher (we shall assign a higher number to a group with higher neutron energy). This enables us to introduce into the problem processes bound up with thermolizations
(g= 42, . ..t G), where
n,(r)=&y jNB(W,P) dl.c v - Y2 -i
(g=1,2,:..,G)
0
is the group g neutron density, of neutrons of the same group. The boundary conditions
h’g(f4 Y, cl) = where the
Hg are
constants
for Hg
and fg(r)
(18)
is the given
source
intensity
are
&=I,2
,...,
G; p<(I),
depending only on the group number, and R is
The
the cylinder
of
cell
* heterogeneous
13:
reactor
radius.
By analogy
we shall
with (51, 4 HR = ;
assume that the conditions
Fs+(R)
(g =
are also satisfied in the multigroup signifying the albedo and satisfying P,+uq=
1, 2, . . . , G),
where the A, are constants
case,
0 4 A, 4 1 (g = 1, 2, . . . , 0
~)cos(~,s-~)~sz
~&vg(z-q
(1%
and
(g=1,2,.:.,G).
+ As distinct from the one-group case, in which the neutron velocity can be taken as unity, we cannot neglect in the multigroup case the different group velocities eg 1, 2, . . . , 0. We therefore have to pass to a general definition of the neutron flux, corresponding to the group number g
(g=
F,(r) = ggS445
ww~, wm
(g=
1,2, . . . . G).
+ This conclusion nected with F,(F) F,(r)
obviously by =
Fg+(r)
In the case of cell be put equal to unity.
also -
design
applies
to F,+(r)
all
, con-
1,2, . . .. G).
(g=
Fg-(r)
and F,-(F)
the A, in boundary condition
(19) must
Using the substitution
N,(r,
Y,
r
___ 1 -VI- Y2[ s1s” 1
~1)= mg(r, Y, CL)+ Hg exp - -
ag(s)sds
rJg
g=1,2,...,G),
(O
iJrg(c yI CL)= flg(r, Y, cl) + Hg exp
(-1 GpcLOo; we reduce
the fundamental
1
1
-
___ 11 -
g=
system (18)
+ +(I - 11”)
1 y2sTl/s”- rqa- p2)‘. R
ag(s) sas
1,2, . . . . G),
to the form
(20)
138
E.S.
Kurnetsov
where
Further, we put
(g=1,2
,...,
G).
(22)
g’=l obtain, for determining the unknownfunctions md fl,+, g(r, y, p) C + 1 systems of integro-differential We now
fio,g(r, y, p) equations
(g=GL..,G), +: $Br i ipo,zPAw+-2 4
_P
~l-V,(p~+~-ff$f-)+ug(r)iVg~,g(r,p,y)=
==. (g’=
For the functions n,(r) tions
I,& . . . . G; g=
1,2, . . . . G).
and the constants F,+(R), we have the equa-
The
cell
of
a heterogeneous
%tr) = b(r) + 4nHgogb), tg =
139
reactor
Fg+(R)
=
F,+(R) +
H,Tg
1, 2, . . ., 9,
where
Using relationships
ng(r)=
(22),
fio,g(r)+
we obtain
from this
~Hgriig.,g(r)+4I1HgWg(T)
(g =
1,2,.
. . , G)
g’=i
and
Fg+(fi)
= Fbtg tfi) + 5
g’=i
H,J%
g(R) + HgTg
(g=
I,?
. . . . G).
In order to find the Hg, it only remains for us to substitute the boundary conditions (19) and solve the resulting system
H,=--for
Ag [b(R)+i H&,,(R)+&T,] al g’?=i
(23)
(23)
in
(g = 1,2, . . . , G)
the H,. For example,
in the case of
two groups we have the system
It might be possible to start from the hypothesis that the neutron radi“reflected” from the cell boundary, is independent of the neutron ation, velocity, i.e. independent of the group number. In this case the constants Hg in (20) and (21) would have to be replaced by the same constant H, and instead of (22) we should write K(r, As a result,
Y, I4 =
we obtain
No,&,
Y, cl) + Hmi,s (r, Y,
two systems of
integro-differential
P). equations
E.S. Kurnetsov
140
(g=
172,. . ..G).
(g =
1, 2, . . .,
Furthermore, we hsve
n&T(r)=
50,&J
+
H[Ri,g(r)
+ 4nog(r)]
G)
and
Fg+(JO =
Fo,~+(R)
= H{&,,+(R) + TV} (s = b2,. . ., G).
We denote by A, the size of the velocity the g-th group. We now have from (24)
interval
(24)
corresponding to
Putting
i; “FB-(R)Ag = -&$ Fg+(R)Ag = H, g=i g--1 where 0 < A S 1, we obtain for H the formula
1-
$2g=1 {h*,(R)+ T,) A,
Calculations using the method described were csrried out for a number of examples in the case of the one-group problem. The solutions were obtained by computer, using special programmes. We shall not dwell here
The
cell
of
reactor
a heterogeneous
141
on the approximate method of solving the neutron transfer equation cl), on which the compilation of the programmeswas based, since this is of no importance for the purposes of the present paper. It is sufficient to say that the approximate solution was based on an improved version of Carlson’ s-SN-method (see [Zl). Experience of computer computation shows that evaluation of the function w(r) for all Fk from formula (121 takes up a great deal of time. The evaluation of this function, and of the interval T, can be avoided if the problem is not reduced to zero boundary conditions*. For this purpose we’ put in (If
N(r7Y7 PL)=
--
wr>Y, CL)-I- HNl(F, Y, cl),
where &,(r,
y, 1.4 is the solution
of (X) under the zero boundary condi-
tion,
y, cl) is the solution
of (1) for f(r)
N,(r.
E 0, but under the
boundary condition Nr(fi, y, IA = 1 for directions and g is an unknownconstant.
inward to the cylinder,
As a consequence,
F(R) = &l(R) + ml(R)
= 0
and
IT Since
pression
PO(R)
=
z--z
Fe(R)
(14) for i,
FofR) Fi (4
ax-Pi+(R)
and Pi+(R)
= Pi(R) i- T,
PO(R) we arrive at ex-
where i = H.
It must be mentioned, however, that the accuracy of this method depends on the fineness of the mesh used in the finite - difference part of the computation, whereas the direct evaluation of w(r) and T involves only the accuracy of the quadrature formulae employed and can be realized independently of the choice of mesh. A similar modification of the method is likewise.possible in the multigroup case. We shall first consider two examples concerned with the one-zone problem with a constant solution known in advance (these examples were selected by V.A. Chnyanov from considerations bound up with the diffusion approximation). *
This remark follows of
reducing
from discussion with A.V. Voronkov as to means the volume of machine computations.
142
iiuznetsou
E.S.
In this case the solution of the problem for an infinite circular cylinder can be reduced to solution of the integral equation
0
0
where K(r, p) is a kernel whose concrete form will not be quoted here (see e.g. 141). In the case of constant f and p, we have R
R ?z(r),= 4nffofr)+W
5 0
W,P)P%fB~
(25)
~(P)~(~,P}P~P* 0
is easily shown that, if the constants H, f, o, and p are in a certain ratio, this equation has a constant solution n(r) = 4df. For, on substituting this value of n(r) in the right-hand side of (25), we get
It
P~PdP+4nw~
n(r) = 4nHo(r)+4nfjKo-,
w,
P)PdP.
0
0
We now use the relationship R c&p,
s
p)pdp
=
1-
w(r).
0
which we established for a = 1.
in another paper (see [41, formula (2.21)):
Hence,
R s
fl(r,p)pdp
=
1 -o(r).
0
Now,
n(r) = 4n(f + flH) + 4n(fZ - f - pH)o(r).
be constant,
we have to satisfy
In order for n(r) to
the condition H = f t pH,
whence
H=
f/Cl - p,. For instance, with f = 1, p = 0.9, we get H = 10. Notice that this result is independent of the choice of cylinder radius R. A concrete calculation was performed for R = 1 and for R = 3, using the’ integro-differential equations (10) and further formulae relating to the cell calculation in the one-group approximation. In both cases the following values were taken for the physical parameters: a = 1.0, p = 0.9, f = 1. The finite - difference method used for solving equations (IO) was characterized by the following data. The interval 0 S r < R was divided in both cases into a hundred equal parts. Of these
The
ceZZ
of
reactor
a heterogeneous
143
hundred intervals, the third, counting from the cylinder boundary, was in turn divided into five equal parts, the second into ten equal parts, and the first into ten equal parts. Thus there were altogether K = 122 unequal intervals. The interval - 14 p 61 was divided into I = 64 equal intervals. As regards the variable y, its interval 0 < y < 1 was characterized by J = 7 Gaussian base-points. The computation was performed by successive approximations with a calculation such that the remainder term of the iterational process did not exceed E = 10e4. The results of the computations were with
R=l,
with
R =
H = 10.0227597; H = 10.0137399.
3,
We took the following group Of examples on cell calculation The initial physical data are contained in Table 1. TABLE 1. Physical
from [51.
parameters -
Jvb
R
RI
CM
PI
a1 c&t-’
CM
CM-’
k I 5
1.30 1.30 1.30 1.50 1.70
0.7221 0.7221 0.7221 0.7221 0.7221
11.90 11.28 7.93 7.93 7.93
0.3991
0.3991 0.3991 0.3991 0.3991
0.63721 0.3721 0.3725 0.3925 0.3925
0.3717882 0.3717882 0.3921882 0.3921882 0.3921882
4.88299 4.652288 3.541005 3.606925 3.672845
In all five examples it is a question of two-zone cells, consisting of an interior uranium block of radius RI and an outer moderator zone of radius R2. The “optical radius” r = alR1 + a2(R2 - RI) is given in the last column of the table. Table 2 gives some data concerning the organization difference computation of the five versions mentioned TABLE 2. NO. Jv"P
; 1 5
K
116 116 116
I
64 lf
J
:: 7
Finite
- difference
Of
NO.
I--
intervals in zone I
;: ;: 34
Of
-
meshes
t--
intervals of zone I1
E 75 70 66
of the finite above.
Remarks In all versions the last two intervals of zOne I and the first and last intervals of zone 11, were subdivided into five parts each
144
Kuznetsov
E.S.
TABLE 3. Computational
results -
n,
e
H
&I
0.86692
i EEE 199:2125 160.7575 2 133.5051 5
0.87980 0.94051
0.95196 0.96013
1.9208 1.9052 1.8096 1.940 2.0722
1.891 1.877 1.791 1.931 2.070
1.898
The computational results are given in Table 3, which includes, in addition to H, the thermal efficiency 8, calculated from the formula
and the “disadvantage
factor”
ll = &/nl
RI
where
s
s
n(r)&
0
Ri2
the thermal neutrons, RI
n(r)rdr
iii =
for
’
52 = ;i2 _ Ri2 *
The last column but one gives the values of the disadvantage factor calculated by Amouyal et al. (see [51) by a special method, different from ours, while the last column gives the values of the disadvantage factor computed by Carlson in the Se-approximation, for variant No. 1. A comparison of our calculations with the data of Amouyal reveals that our figures are somewhat higher while at the same time remaining fairly strictly parallel.
The
cefl
of
reactor
a heterogeneous
145
In addition, we calculated the neutron density distribution as a function of the radius r for all five versions. The results are shown in the figure, which gives the ratio n( t-)/n(O). The zone boundaries in the individual versions are marked by vertical lines, with the corresponding figures on the abscissa axis. The last example, set up by L.V. Maiorov, is a three-zone cell and is characterized by the following data: 0 < r < R, = 1.75. (natural uranium); f - 0, a1 = 0.62046, p1 = 0.427; R, = 1. 95 (aluminium) If- =O, a2 = 0.098474, pz= 0.08442; Rl = 1.75
Angular distributions
5
3
2
6
I,
222.78 22%:: 217:38 212.84 206.47 197.73 185.6G 168.46 142.34 107.21
202.33 201.85 200.26 197.49 193.37 187.60 179.67 168.72 153.09 129.36 98.098
103.24 :FZ! 1oO:so 98.733 95.833 91.838 86.293 78.344 66.194 49.758
74.815 74.599 73.938 72.796 71.102 68.739 65.518 61.118 54.974 46.067 35.618
As a result of the computations we obtained
55.822 55.642 55.096 54.139 52.734 50.788 48.162 44.629 39.820 33.190 26.279
H = 9.60955629,
5.8818 5.8707 5.8371 5.7735 5.6911 5.5682 5.3984 5.1625 4.8247 4.3157 3.6654 0 =
0.8139390.
The variation of the neutron density (normalized as in the previous examples) is shown in curve 6 of the figure. In addition to the cell characteristics enumerated above, we computed the angular distributions of the phase densities on the cylinder axis, which are of interest to experimental workers. As is easily seen, on the cylinder axis the phase densities depend only on the variable y = cos 5,
146
Kuznetsou
E.S.
which considerably simplifies the computational the computations are given in Table 4.
formulae.
The results
of
Analysis of these figures reveals that the shape of the curves N = N(y) give a good reflection of the cell properties. In particular, an extremely sensitive characteristic of the physical and geometrical properties of the cell is provided by the ratio of the values of the neutron phase density corresponding to a direction along the cylinder axis, and corresponding to a direction perpendicular to the axis: N(O)/%(l). This ratio sidered:
has the following
values
in the examples that we have con-
TABLE 5 No. of
variant
2
1
2.0630 --
2.0779
~(O)/~(l)
3
4
2.0749 -
2.1004
5
6
2.1242
1.6047
The author wishes to thank A. I. Vaskin, E. I. Acknowledgements: Panfilov, 0. V. Sysoev, E.G. Natrusov and G.P. Churkin for cooperation in the computer programming and in carrying out the computations, and V.A. Chuyanov and L. V. Maiorov for preparing the individual physical examples and for participating in the discussion of the results. TransEated
by D.E. Brown
REFERENCES 1.
GLASSTONE, 8. snd EDLUND, M.C. The theory, Van Nostrand.
2.
MARCHUK. G.I. Methods of calculation of nuclear reactors rascheta yadernykh reaktorov), Moscow, Gosatomizdat,
3.
4.
POMRANING,G.C. The Wigner - Seitz calculational method, Nucl. Sci.
of
cell,
nuclear
reactor
(letodg 1961.
a discussion and a simple 17. 3, 311 - 314, 1963.
Engn.
KUZNETSOV. E.S. On the temperature distribution in an infinite cylinder and in a sphere with non-monochromatic radiant equilibrium Xh.
5.
elements
vpchisl.
Hat.
mat.
Fiz.
2,
2,
217
-
240,
1962.
AMOUYAL. A., BENOIST. P. and HOROWITZ, J. Nouvelle m&hode de deftermination du facteur d'utllisatlon thermique d'une cellule, J. nucl. Energy, 6, l/2. 79 - 98. 1957.