- Email: [email protected]
The weighting of diffusion coefficients in cell calculations
Reactor Science and Technology fJourna! of Nuclear Energy Parts A/B) 1962, Vol. 16, pp. I to
I I.
Pergamon Press Ltd. Printed in Northern Ireland
THE WEIGHTING OF DIFFUSION COEFFICIENTS IN CELL CALCULATIONS D. c. Atomic Energy Establishment,
LESLIE
Winfrith, Dorchester,
Dorset
(Received 12 April 1961)
Abstract-The new theory of leakage in heterogeneous systems introduced in this paper is based on BENOIST’S lemma that the leakage from a region of such a system is equal to Bat*, B* being the geometric buckling and LBone-sixth of the mean square distance travelled by a neutron between birth and death, provided B1 is small. Lzis calculated exactly and directly and is shown to depend on a new type of flux [arising from a tilted source: this flux can be calculated in one lattice cell, so that it is no longer necessary to consider migration over the lattice as a whole to compute leakage. The exact calculation of t is troublesome, but a diffusion theory calculation can be carried through quite easily. From this, it appears that the radial and axial migration -areas are given by L,* = &/z, and Lz2 = 0,/C., C, being the normal cell-averaged absorption cross section. 0, is obtained by weighting the local diffusion coefficient with the usual flux, and & by weighting with a new flux derived from 6. There is thus a migration area asymmetry even in the absence of holes, and an example suggests that in some lattices’ this asymmetry may be quite substantial. The new theory may also be applied to more conventional streaming problems, and this application is discussed in the Appendix.
1, INTRODUCTION
BOTHthe absorption cross section C, and the diffusion
coefficient D vary across the lattice cell of a typical reactor. Feinberg-Galanin technique is much too complex for survey calculations, and the first step in such calculations is therefore to homogenize the lattice cell by inferring ‘paste parameters’ c, andii. za is computed by weighting the Za in each region with the local flux, but there is no general agreement as to how b should be calculated. Some authors have advocated weighting D with the local flux, others have advocated weighting E:,, and others still have advocated some very fancy prescriptions. The problem has been solved in general terms by BENOIST (1959) but his prescription contains collision probabilities which are very difficult to calculate. In this paper, a new and as yet unpublished theory of streaming due to the author is applied to the problem. On one-group diffusion theory, the leakage from a homogeneous region is given by B2L2, B2 being the geometric buckling, and L2 being one-sixth of the mean square distance travelled by a neutron between birth and death in an infinite lump of the same homogeneous medium. BENOIST (1959) has shown on the basis of one-group transport theory, that the result may be generalized to a uniform heterogeneous lattice, provided the buckling is small; it gives the leakage as F (&%3
where Lk2 is one-half the mean square of the k component of the vector joining the points of birth and death (k = x, y, z), Bk2 is the component of ihe geometric buckling in the k-direction. If B2 is large; the heterogeneity produces buckling-dependent corrections to L2. The author’s new theory gives a closed expression for L2 in a lattice, based on exact one-group transport theory, in the limit of small B2. b is then defined as L2za, since the essence of any definition of D is that it should give the leakage correctly. The paper therefore begins with an account of this new theory. 2. THE
CALCULATION
OF L”
From now on, we restrict ourselves to the usual cylindrical lattice and to a single neutron group. It is infinite in all directions, being uniform in the z-direction and periodic in the x- and y-directions. Let P(r,r’) be the flux at a point r due to a unit source of neutrons at r’ [r is the vector (x, y, z), the origin being arbitrary]; it is therefore the Green’s function for the lattice. The rate of absorption at the point r EC,(r) being the macroscopic is just &(r)P(r,r’), absorption cross section. The mean square distance travelled by neutrons from this source is therefore sss m
&(r)P(r,r’)(r
- r’)2 d%
the integral being taken over the whole of the lattice,
D. C. LESLIE
2
that is, over the whole of space. To obtain R2, the overall mean square distance for the lattice as a whole, we have to weight this according to the distribution of sources q(f) in the lattice. Thus
becomes
lengthH d3r’C,(r)P(r,r’)q(r’)(x - x’)~ NW116
Q) +
unit
B. (2.6)
B is a contribution from the boundary of the region (N,H)introduced by the change in the limits of
where the second integral is taken over one cell in the x - y plane and over unit length in the z-direction. The source function q(d) is normalized to a unit rate of neutron production in this element of volume. For example, in a metal-fuelled lattice, q will be constant in the moderator and zero in the channel and in the fuel: it will be equal to l/V,, V, being the ‘volume’ (actually an area) of moderator in one cell. In a homogeneous medium, we would put 2 = 6L2, L2being the thermal diffusion area. In a lattice, there may be asymmetries in the diffusion area, and we must put i? = 2L,2 + 2L,2 -+-2L, where
(2.2)
integration. It can be made relatively as small as we wish by taking the region to be large enough. Jgnoring this contribution, equation (2.6) may be written length
LZ21Lr2=-
N
1
2NH sss
d%$Jr){$)(r)
- 2x@(r)
N cells
+
x24dr)I
(2.7)
where (2.8a)
J$) (r) =
x’q(r’)P(r,r’) d%' sss m
and
7-g)(r) =
(2.8b)
xr2q(r’)P(r,r’) d31.'.(2.8~) US
03
unit
Similar expressions may be derived for L,"and Lz2. The significance of the suffix 0 will be made clear in (2.3) the next section. Because of the definition of P(r,r’), #0(r) is simply the definition of Lv2and Lz2being similar. In lattices the flux due to a source distribution q(r’). Since this with square, hexagonal or trigonal symmetry is the actual source distribution in the lattice, &,0(r) must be the actual flux. &j(r) is a new type of flux, L,2= L,2= L,2 (2.4) due to a ‘tilted’ source distribution x’q(r’). Figure 1 shows a comparison of a tilted source with the normal the ‘radial’ diffusion area, while Ls2 is the axial source in a one-dimensional lattice. This type of diffusion area: in general, Lr2is not equal to Lz2. flux would appear to be new, although SELENGUT Equation (2.3) represents the fate of neutrons starting from one unit volume (unit length of one cell) Tilted Conventional source source ; and migrating over the whole of space. We could, if we wished, start the neutrons from a large region ccl-Nentlonor of length H in the z-direction and containing N cells. source We would then have
7L
-Moderator --)
the neutrons still being tracked over the whole of space. We now come to the vital step in this analysis. Instead of being start&d in the region (N,H) and tracked over the whole of space, the neutrons are started everywhere, a death only being counted if it occurs inside the region (NJ-Z). Equation (2.5) then
Fuel ------+
--Moderator-+
FIG. l.-Illustration of conventional and tilted sources in a circular lattice cell. (Seclion along the x-axis.)
The weighting of diffusion coefficients in cell calculations
I -
J Fuel
+ /
Moderotor i
-
Mrxierator-
problem was almost intractable if heterogeneity was taken properly into account, since the lattice had to be treated in its entirety. Now the calculation of leakage has been put on the same footing as the calculation of thermal utilization, in that both computations can be done in a single cell, the only restriction being that the geometric buckling B2 should be small: this restriction applies equally to the calculation off. 3.
FIG. 2.--Illustration of the fluxes % and tin a circular lattice cell. (Section along the x-axis.) b = 2a; D, = 40,; No absorption.
(1960) has noticed the relevance of flux tilt to leakage problems. It is shown in the next section that x$)(r) = .$&r) - Et)(r)
(2.9)
where 4,(r) is the ‘normal’ flux defined by (2.8a) and [f)(r) is another new flux which, like q&(r), is periodic in the lattice. It is antisymmetric in the s-direction, symmetric in they-direction and uniform in thez-direction. Its behaviour in a one-dimensional lattice is illustrated in Fig. 2. It is also shown in the next section that $‘)(r), which is the flux due to a ‘curved’ source .u’2q(r’),is given by 7:‘(r) = s*&(r) L 2s$Yr)
- $‘(r)
(2.10)
where v:‘(r) is yet another new flux, which is also periodic in the lattice cell. It is symmetric in both the x’ and y’ directions (though its dependence on the two co-ordinates is different) and uniform in the z-direction. Substituting from (2.9) and (2.10) into (2.7), we find
3
SIMPLIFICATION OF THE FORMULA
BASIC
We shall now show that the flux 1’can be eliminated from the basic formula (2.12), which can be expressed in terms of ,t only. This reduces considerably the labour of using (2.12), since (2.10) shows that it is necessary to know 5 before 1’can be calculated: one step in the work can now be eliminated. For this further step we need to consider the equation satisfied by 1~. We shall take this in its exact or transport form, so that the work of this section does not involve any further approximations. BENOIST (1959) prefers to use the transport equation in its integral form, but we find the integro-differential or Boltzmann form more convenient. We write this equation in a tensor notation, in which
is the co-ordinate vector, normally denoted by r or (s, ,v, r): thus x1 is to be identified with X, xq with _r and s3 with :. Similarly {~4~} is a unit vector giving the direction cosines of the neutron path. It is often denoted by Q, and since it is a unit vector, it must satisfy
From this it will be seen that we are using the summation convention : a repeated index is to be summed over 1, 2, 3. With this notation, the Boltzmann equation for a vector flux + may be written
(2.11) ui x
(-yj,uI;)
7
c,(xj)b(sj,uk)
- 1
Since t’?)(r) is periodic, this integral can be taken over one cell only, giving as the basic formula of the new theory unit rrr Lx2 = 3
JJJ
C,(r) v?‘(r) d3r.
(2.12)
rell
Similar formulae can be derived for L,2 and Lz2. The advantage of this new procedure is that it reduces the calculation of leakage to the computation of a function in a single lattice cell. Previously the
=
C,(.Uj;U,U,')~(Xj,Ui')
d*Um’+
q(XjvU&
(3.1)
ss
Here Ct(xj) is the total cross section at the point (xj‘; &,(.x~;u,u,‘) is the scattering cross section for deflections through an angle cos-‘(u,u,‘) at the point {_uj} q(xj,uk) is the source of the vector flux I$(x~,u,). This equation brings out explicitly the dependence of the vector flux on direction as well as position. The first term on the right-hand side of (3.1) represents
D. C. LESLIE
4
the contribution to a given directed flux at a particular point from neutrons moving in other directions which suffer a suitable collision at that point. SJ d2ui’ represents an integration over all possible directions of motion of the incident neutron, and is normalized to 1. In computing reaction rates, a full knowledge of the vector flux (x,, uJ is not required. These are determined by the ‘ordinary’ flux
which is the average of the vector flux over all directions of motion, and is a function of position only. The suffix 0 in the formulae of the previous section indicate that it is this average flux which is involved. In this paper we shall assume that the scattering cross section has a simple cosine anisotropy only, so that X,(X,; I(k1lk))= E,(Xj){l
+ 3E(Xj)Uj$~)a
(3.3)
Here & is the usual isotropic scattering cross section, and p is the average of the cosine of the scattering angle, which would be zero if the scattering were purely isotropic. In a more accurate treatment there would be further terms on the right-hand side of (3.3) involving P,(u,u,‘) and perhaps higher-order spherical harmonics. However (3.3) will suffice for illustration, and the conclusions of this section are not affected by the assumption. Substituting from (3.3) into (3.1) and suppressing implicit variables, we finb
equation (3.4). We then find that 5@)satisfies uj y
+
c,5 @)- C,(@’ + P&r,+)} = -ui+
X,
ax1
This follows since axi = 6ii, the Kronecker
(3.7) delta
symbol, which is 1 if i = 1 and 0 otherwise. (3.7) is a normal Boltzmann equation, showing that 5(x) is derived from the anisotropic source -ui$, Since 4 and the cross sections are periodic in the lattice, it follows at once from (3.7) that 5(x) is also periodic. In exactly the same way, we find that Y(X)satisfies
adz)
uiaxi +
Cp
- lqvp
+ puiYi(+)} =
-2up
(3.8) v@) and the BY) also being defined by equations similar to (3.2) and (3.5). Equation (3.8) shows that Y@)must also be periodic in the lattice. We now integrate (3.8) over all directions of neutron motion, and use the definitions of VP’, Y!*) I and El”): we find I
3
a+)
awYj
+
I;&)
=
2
r”’
(3.9)
3 l
where &, is of course equal to Zc, - Z,: the second term inside-the curly brackets in (3.8) gives no contribution since the average value of ui is zero. Thus from equation (3.9), equation (2.12) may be rewritten unit
unit
L,” =
-2
Z(r) sss ! sss p(r) d% -
CC11
6
cell
t
dyr
(3.10)
where (bOis defined by equation (3.2) and the vector +i, which is also a function of position only is defined by +ifxj)
=
3.f j”iNxj4
@%2
Using Green’s theorem, the second term is equal to
(3.5)
+i is proportional to the neutron currentj,, the exact relation being jd = &$, u’being the average neutron velocity in the single group. We shall take (3.4) as the standard form of the Boltzmann equation for a flux 4 due to a source q: a more extensive discussion of this topic is given by DAVISON(1957). Thus the gux x(~), which is derived from a source qx, must satisfy
# and the xi(*) being defined by equations similar to (3.2) and (3.5); in (3.6), we have represented x as x,. In this equation we write x(*) as x1+ + P, and use
This is proportional to the current out of the unit cell of neutrons having the distribution v@). Since v,,tr) is symmetric in the x- and y-directions and uniform in the z-direction, this current is zero. Therefore the second term of (3.10) vanishes, and we have unit’ L,2 = -
f
SIJP(r)
dar.
(3.11)
cell
This is the simplification which was promised at the beginning of this section. Similar equations may be deduced for L,"and LZ2: provided that B2 is small, these equations are all exact.
The weighting of diffusion coefficients in cell calculations 4. THE
CALCULATION OF THEORY
special ‘flux’
5 BY DIFFUSION
The calculation of I# by transport theory is a formidable problem even in the simplest geometries, and the calculation of 6 is even more difficult. It therefore seems profitable to start by calculating L2 from equation (3.11) using diffusion theory, and we shall now do this. It must be remembered that diffusion theory is not very accurate in real lattices, and that a better procedure is likely to be needed. However if 4 has been calculated by diffusion theory, it seems perfectly reasonable to compute E in the same way. Strictly, we are going to use the PI procedure rather than diffusion theory, since the former can always be carried through in a consistent manner. In this procedure, we put 4 = 40 +
(4.1)
%di
and E(I) = ,$’ + n.EW. z,
5
(4.2)
Here &, and +i are defined by (3.2) and (3.5) respectively, the definition of $1 and @) being similar. Substituting from (4.1) and (4.2) into (3.7) we find
(4.7) In fact, (4.6) may be written L,” = b,. z,
(4.8)
where unit d,
=
unit D(dUr)
us
d3r
L(r) II(Sss
Cdl
d+)
(4.9)
Cdl
and unit z,
=
S,(r)
d3r.
(4.10)
sss
Cdl
Here D* is a cell-averaged diffusion coefficient, the weighting function being t,(r). The form (4.8) suggests that we should investigate the quantity Z,. As an illustration, suppose that the lattice cell is square, the half-width being a. Since the ‘unit cell’ is of unit height, we have
The integral of the second term in the inner integral is just @‘(Xl = a) - @)(x, = -a).
- &{@’ + Ff4&@) = -z&
- urujl&. (4.3)
We now multiply this equation by u, and integrate over all u. Using the obvious results
But .$’ is antisymmetric in the x,-direction and periodic in the cell, so that it must be zero on the faces x1 = a and x1 = --a of the unit cell. It follows that the second term in (4.11) is zero, and that unit
I&
d2U +u,u,dpu
= 0,
(6, being the Kronecker viously) . we find g’
=
(.&
d% = g aij
(4.12)
delta symbol defined pre-
_‘,&(A+$!g.
But the Pr approximation is just
D=
]/Wj
(4.4)
to the diffusion coefficient 1
3(% - $3
*
(4.5)
Substituting from (4.4) and (4.5) into (3.11) we find
a very suggestive formula. This formula is not exact. since it depends on the use of the PI approximation. Equation (4.6) shows that Lz2 is obtained by weighting the local diffusion coefficient D(r) with the
From this equation, it follows that Z, = Z, = Z,. Now in the Introduction, we remarked that the average absorption cross section in the lattice is obtained by weighting with the normal flux &r), that is unit uuit c, = V%Mr) d3r )/(///M) @rJ. (4.13) (IIS cell cell The numerator of this expression is the mean rate of absorption in the cell. This must be equal to the mean rate of production, which is 1, from the definition of q(r). It follows at once that z, = l/E,
(4.14)
and therefore L,” = DJZ,.
(4.15)
This is, of course, just the sort of relation which we
D. C. LESLIE
6
had hoped to find. It shows that the averaging process defined by equation (4.9) is compatible with simple definitions of migration area and diffusion coefficient. In exactly the same way, we may show that L,2 = I@,
(4.16)
I), being obtained by weighting D(r) with
holes. If diffusion theory @ valid, the axial diffusion coefficient b, is obtained by weighting the local diffusion coefficient D(r) with the normal flux, and the radial diffusion coefficient b, by weighting with the special flux c*(r) defined by equation (4.7). The corresponding migration areas are obtained by dividing the appropriate diffusion coefficients by the cell-averaged absorption cross section. 5.
In a square cell we shall obviously find (4.17)
D, = D, = b,
the radial diffusion coefficient. We would expect that the weighting function for the axial diffusion coefficient would be
this following from the equation
which is analogous to (4.4). Integrating the equation for P) analogous to (4.3) over all U, we find
f.as:“~C,J(;‘ -;$3=Dz = ay 3
a.Yi
=D
2)
(4.19)
The right-hand side of (4.19) is zero, since the reactor is uniform in the z-direction, and #+, is therefore independent of z. Combining (4.18) and (4.19), we find that $1 satisfies
CALCULATION
OF THE ASYMMETRY
Thus the axial diffusion coefficient is obtained by weighting the local diffusion coefficient with the normal flux y$(r) unit
unit
B, =
(4.20) Cdl
f3 = 0 being one of the axes of symmetry of the cell (other terms in the Fourier series being excluded by symmetry) and then assuming that g,r is the dominant term. From the remarks in Section 2 following the definition of EF’, the boundary conditions on this flux are* ” on
ap
ay@) =O On
P _x=---
-2 +;.
From this it follows that the Fourier expansion of lt’ is (2n + 1>e. More terms are now allowed, since the diagonals are no longer axes of symmetry. If we take only the first term of this series, namely j,(r) cos 6
Cdl
the axial migration area being then given by L,2 = D,/&
AREA
Having discovered that d, is not identical with b,, we must see whether the discrepancy is large enough to be important in practice. We do this by calculating both quantities for a typical lattice cell. Since we are using diffusion theory, this cell must be free from voids. We take it to consist of a fuel rod of radius a, absorption cross section C,,,, and diffusion coefficient D,. embedded in a block of moderator in the form of a square of side P. The moderator diffusion coefficient is D,, and moderator absorption is negligibly small. It is conventional in normal flux calculations to replace the square by a circle of radius b = P/15. This is equivalent to expanding the flux in a series of the form d,(r) = &g,(r) ~0s 4~~9
Eg)(r) = 0
since D is also independent of z. From this, it follows that Et) is also independent of z, and therefore that
MIGRATION
(4.21)
Thus we may conclude: There is a migration area asymmetry in reactors with cylindrical symmetry, even in the absence of
the calculation of Eg’ can be cylindricalized in the same way as that of &,. The author (unpublished) has calculated 6 in a cell with square boundaries and has * These conditions follow at once from 5 being continuous, periodic in the lattice, anti-symmetric in the x-direction and symmetric in the y-direction.
The weighting of diffusion coefficients in cell calculations
thus been able to show that as long as the fuel diameter is not more than half the lattice pitch, this cylindricalization introduces very little error. It is used here to simplify the work. With this approximation, the boundary condition on j,(r) is j,(b) = 0
so that, using the boundary condition (5.3) L,2 = nD1 ob2rg(r) dr + (Do s
x {nab(a) - na2g(u)} + 27 l&r)
~&j(r) = x+,(r) + $?(r) = r cos 0$,(r) + @j(r)
(5.2)
b, = D,+ @h(a) (Do
-
4)
- m2g(a)} -I- Zn[rg(r)
x
rb
J0
(5.7a) Similarly, from (4.20), the axial diffusion coefficient is given by 2n ‘rg(r) dr s
D, = D, $ (Do - D3
Our task is therefore to calculate the functions g(r) and h(r) (we can drop the suffices since we are only keeping one term in each series) and to average them suitably. Now in the diffusion approximation unit B,
=
’Wr)Ur)
1, SSJ
.
r”b
h,(r) cos 6.
(5.3)
dr
27r -rg(r) dr
2x
h,(b) = bg&) [and h,(O) finite].
(5.6)
whence
rather than with @‘. Following the line of argument used above, we approximate to $)(r) by
From (5.1) and (5.2) the function h,(r) must satisfy the boundary conditions
dr]
[
(5.1)
b being, as before, the outer radius of the equivalent cylindricalized cell. The general theory has been formulated in terms of the flux 5, since this is periodic in the lattice while x is not. Now that we have reduced the problem to one of calculation in a single cell, this distinction is not important, and it is convenient to work with
Dl)
J0
(5.8a)
rg(r) dr
The migration area asymmetry is therefore due to the term in curly brackets in (5.7a). These formulae can be simplified by introducing two new symbols E: the ratio of the mean flux in the fuel to the flux on the interface r = a. S: the ratio of the mean flux in the moderator to the flux on the interface. (In the usual notation, E = l/G). Using these new symbols we find
d3r
b, = D, + (Do - D3 ;‘E”+“v”E;
cell
0
(5.7b)
1
from (4.6) and (4.15) where
axf) ag) L.(r)= +o(r)f (r> = (r>ax1 ax1
b,
i$O (4 1
x1 -
from (5.2). Now
sothat 5,(r)=
-=a axI
cose---a
sin 0
ar
r
$ (30s~e + ” sin2 13-
I
D(r) = Do =D,
44
X = --
a&)
‘OE VoE 4 v,S
1; Vo=rra2;
(5.9)
ae COG
if
0 < r
if
a
(5.8b)
Vl=7r(b2-a2)
a
0 2
(5.4)
Since D(r) in fact depends only on r, and C,(r) is independent of z, we have
Now
where
D1 C (Do - 01)
=
The calculation of E and S is quite standard, but the computation of X is in general rather tedious. However if the absorption in the lattice cell is negligible, this latter computation is quite easy: provided Do + D,, this case is not trivial. Since x@)is the flux produced by a source qx,, it follows that in the diffusion approximation the isotropic component &) satisfies -gi(D(r)z) Since we are putting
i-&&)=4X1. xe) = h(r) cos 8, we find that
D. c.
8
h(r) satisfies -- l&)(r):)
(5.10)
+ (X.+$+.
LESLIE
This will be zero if D,, = D,. Since the mean free path in the fuel is presumably always less than in the moderator the asymmetry will be greatest for D,, = 0, when it is given by
If the absorption is zero, the source must also vanish and in a region where D is constant, h(r) = Ar + ;“. We shall take this form for h in the moderator region a < r < b: in the fuel region 0 < r < a the second term must be excluded, so that h(r) = Cr, 0 < r < a.
The boundary conditions are that h(r) and D dh/dr should be continuous at the interface r = a, together with equation (5.3). Since there is no absorption, the ‘normal’ flux g(r) must be constant, and we can write this equation h(b) = b&a).
Eliminating the arbitrary find that
h(a) f&a)
-=x+l=(l_v)D,+(l+v)Dl
constants A, B and C we
24
=
’
(1 + VP, + (1 - W,
Since v may be as large as 1 (equal volumes of fuel and moderator) it would seem that asymmetries as large as 50 per cent are possible. It should be noted that the asymmetry is considerably reduced if Do is not actually zero. For example, (5.13) gives, for v = a, a = 20 per cent if Do = 0. If D, = ) D,, a is only 9.9 per cent: this reduction is larger than might be expected. Asymmetries of this size would have considerable of practical consequences for the interpretation exponential experiments and in the optimization of small power reactors. There is therefore a need for a more precise calculation in which allowance is made for absorption in the fuel. For simplicity, we shall suppose that the source is zero in the fuel and flat in the moderator: we shall also ignore moderator absorption. Using the diffusion approximation in the fuel as well as in the moderator, the flux rise parameters are given by the standard formulae E
D
(1 - v)D, + (1 + v)D, "
21,&d
=
(5.14)
Kd,( KU)
.y = 1 + ; (W7)”D”E D 1
3b2 - a2 b4 b (b2 - a2)2 In ii - 4(b2 - a2)I
(5.12)
This formula has already been deduced by SELENGUT(1960) using an entirely different method. Since the physical models used by SELENGUTand by the author are both reasonable (although quite different) it is highly satisfactory that they give the same answer. The author’s method would seem to be the more general, since it covers the case of finite absorption, it can be used with more exact methods than diffusion theory, and it has revealed an asymmetry in the migration area. For this case of zero absorption, equation (5.8b) gives b, = vD, + (1 - v)D, so that the migration area asymmetry is L2 a=-4_1=f_1 L,”
(5.13)
1+ 0
(5*11)
where v = u2/b2 is the volume fraction of the cell occupied by the fuel. We can now substitute back into (5.7b), remembering that E = S = 1 since g(r) is flat. We find D
L2
orZ~_llL. L,2
b D,
(1 - v)D, + (1 + v)D, (1 - v)D, + uD, _ I = (1 + u)D, + (1 - v)D, ’ D,
(5.15) where, as usual $=3X
a0Z tr.0.
We shall also need to know the total flux rise in the moderator, which is specified by T=
(5.16)
1
We can now calculate the ‘tilted flux’ h(r), which is the solution of (5.10). In the fuel, where there is absorption but no source, the solution is h(r) = CIl(Kr)
O
the second solution K,(Kr) being excluded because it is singular at r = 0. In the moderator there is a
source but no absorption and the appropriate solution is h(r) = Ar + 5 - &qr3
acr-cb.
1
The constants A, B and C are determined
by the
The weighting of diffusion coefficients in cell calculations
continuity and boundary conditions
Jf-
cr,(Ka) = Au +
&$
(1
Ql
Kc&‘(KU)
=
A
-
;
(5.17) 1 (5.18)
-
B bg(6) = bTg(a) = Ab + h - $
and
qb3
(5.19)
1 The total source in the cell must be equal to the total absorption, so that n(b2 -
a2)q =
?ra2.
= n-(~a)~D, . E&a)
&&(a)
9
this new theory is to the calculation of cell-averaged diffusion coefficients, this averaging being done with a new kind of flux. This flux should, strictly, be calculated using transport theory and such a calculation will clearly be lengthy and difficult. To obtain some preliminary results, the flux has been calculated by diffusion theory. Provided diffusion theory is valid, the cell-averaged axial diffusion coefficient is obtained by weighting the local diffusion coefficient with the normal flux, while the averaged radial diffusion coefficient is derived from a new flux due to a ‘tilted’ source. There is thus a migration area asymmetry, even in the absence of voids.
(5.20) Eliminating C between (5.17) and (5.18), we find A--$-
where
&qa2=Z
A+;;,C
1
B
1
1 5
1
4a2 \
Z=--
‘I
[from the recurrence relations satisfied by the function I,(z)1 or (5.21)
I )
j
Finally, eliminating (5.20), we find
A
B
and
and using the relation
0
2T-t =x+1=
2
E(Kcz)~(~ -
(l+t?;+z(l-U)
u)
(5.22)
where, as before, u = a2/b2. If K = 0 and E = 1, the above equation becomes identical with (5.1 l), as it should. The formulae given above enable d, and b, to be computed from equations (5.7b) and (5.8b). Mr. H. M. Brown has carried out such a computation for the special case v = 1, Do = tDl, and has deduced the variation of the migration area asymmetry ccwith KU: his results are shown in Fig. 3. It will be seen that realistic amounts of absorption do not have much effect on the asymmetry. When the absorption is large, the statistical weight of the fuel will become very small. When this occurs the asymmetry must vanish, but the figure shows that this only happens if the absorption is extremely large. A priori, one might expect that the asymmetry would decrease monotonically to zero : the figure shows that this is not so. 6. CONCLUSIONS
I
4
6
8
IO
12
14
16
18
20
K0
FIG.
h(a) a&)
2
AND
FURTHER
WORK
A new theory of neutron currents in a heterogeneous assembly has been formulated. One application of
3.-Variation
of migration area asymmetry with fuel 1 DO 1 absorption for a lattice cell with a = b
2’%=4’
The conclusions of the previous paragraph are only valid if the individual regions of the lattice cell are large enough for diffusion theory to be applicable: they certainly do not apply if the cell consists of large numbers of very small regions, or if there are voids. It is most desirable to remove this restriction by calculating the required fluxes more accurately. This programme of work has already been started. In Appendix 1, we show how Lr2 can be calculated for a cell consisting of a circular hole surrounded by moderator, using boundary conditions worked out by CARTER and JARVIS(1960). The analysis leads to a formula already deduced by CARTER(1960), who has concluded that it is in good agreement with available evidence. At present, the theory is based on BENOIST’lemma S which relates leakage out of a large region to the mean square distance travelled by neutrons inside the region. It is hoped to reformulate the work of Section 2 so that leakage can be calculated directly. With this improvement, it may be possible to study flux matching between different zones of a reactor,
D. C. LESLIE
10
and the dependence spectrum effects.
of diffusion
coefficients
on
Acknowledgments-The author is indebted to M. P. BENOI~Tof C.E.A., Saclay, for making his work on streaming available in advance of publication and for helpful comment, to DR. C. CARTER of A.E.R.E. for many interesting discussions, and to MR. H. M. BROWNof A.E.E. for help with computation.
arguments can be used about 7 which, unlike v, is a ‘real’ flux derived from a real source qx,“). It follows that $‘(a)
C.E.A. 1354. CARTERC. (1960) Streaming Due to Holes in a Reactor. AERE R-3367. CARPER C. and JARVIS R. J. (1960) Zntuence of a Cylindrical
Press. SCHAEFERG. W. and PARKYND. M. (1958) Proceedings of the Second International Conference on the Peaceful Uses of Atomic Energy, Geneva, A/CONF/l5/P.310. United Nations, N. Y. SELENGUTD. S. (1960) Trans. amer. Nucl. Sot. 3, No. 2, 398.
q = &#.
(A.11
Since Z,, is zero in the hole, the basic result (2.12) for the radial migration area may be written &Jr)+‘(r)
and this removes the flux v, which would be difficult to calculate. Now X’“’= X’$ + 5:“’ and since the +i are zero, ,iz’ must be equal to$‘. may be written
mod
the integral being taken over the moderator only. The diffusiontheory methods of Section 4 do not apply to cells containing holes, but CARTER(1960) has shown how the transport calculation of DAVISON (1957) may be combined with diffusion theory to confirm the Behrens value of L,” for this case. It would no doubt be possible to translate CARTER’Scalculation into the language of Section 3 of this paper. The transformation of Section 3 reduces the previous equation to sss
mod
x1$‘(a)
L,a = )
X:‘(r) d%.
dS, - +
i
Thus (A.3)
sss
+
@‘(a)
i
(A.4)
mod
h0le
We shall now follow the procedure of the main text, and will calculate x1”’by diffusion theory in the moderator: equation (A.4) has been put into a form in which it is not necessary to know x(Z) in the hole. We begin by cylindricalizing the problem, replacing the square outer boundary of the cell by a circle of radius b = Zy&. of the form
The isotropic component
of the flux is then (A.5)
so that, in the diffusion approximation $
2:“’ = -30
~09 8 + 2 sine e
r
,
(A.6)
The first term of equation (A.4) is simply - 1
2VxF’(a) . aa cos 0 d0 s0
where Xj”‘(r) = -30
dh G cos 0
64.7)
(in the diffusion approximation) is the current of xtz’ in the radial direction. Substituting into equation (A.4) from (A.6) and (A.7), we find
d3r
sss
-
(A.3)
mod
xb”’ = h(r) cos f3
We now want to consider a lattice cell of a type slightly different from that treated in Section 5. It consists of a square block of moderator of side P, containing a concentric circular hole of radius a. The moderator absorption cross section and diffusion coefficient are &, and D respectively, and themoderator contains a constant source of intensity q per unit area. The conventional flux in the cell is constant and is given by
l:“‘(r) &
t:.)(r) dJr
It
sss
1
Streaming in a cell containing a circular hole
L> = -)
dS* - )
hole
Channel on a Neutron Flux which Varies Linearly in a Direction Perpendicular to the Axis of the Channel. AERE R-3366. DAVISONB. (1957) Neutron Transport Theory, Oxford University
= 4
x&“‘(a)
I
Coeflcient de Diffusion darts un Reseau Comportant de Cavites.
Ip
hole
hole
L;L = f
ASHFORDJ. R. Unpublished material. BENOISTP. (1959) Formulation GPnt?rale et Calcul Pratique du
+ p:“‘(a)} dSi = 0.
f
But $ is constant and the #Q are therefore zero everywhere. follows that
REFERENCES
APPENDlX
{x?+<(a) + 2x&“(a)
d& =
I
dSi
(A.2)
hole
second integral being taken round the channel wall: it represents the current of the flux Y out of the moderator into the hole (the contribution from the cell boundary vanishes). Now the net current of the flux 7 into the hole is zero, since there is neither absorption nor production within the hole (such
L:=D~e~rdr~~Ms’B+~sin’Bj 2=al dh z (a)
+D
COS*e
de
s0 = SOD bh(b) - ah(a) + a’ z dh (a)
1.
(A.81
Now h(r) must satisfy equation (5.10), with constant C,, D and q. It is therefore given by h(r) = AI,
+ BK,(ur)
the
4 + - r c.
(‘4.9)
where K’= 3B,&, the A and B being determined by the boundary conditions. On the outer boundary we have simply h(b) = b+
(A.10)
The weighting of diffusion coefficients in cell calculations from equation (5.3). CARTER(1960) has remarked that the most general form of condition on the inner boundary is h(a) =
5;
af
(a)
(All)
and that all the information about streaming across the hole must be contained in the value of ,!J. CARTERand JARVIS(1960) have discussed the calculation of B and have concluded that the formula 1 + 2aC,, tA.12) P= 2(1 + a&,) should be quite accurate for all sizes of hole. Eliminating A and B from equations (A.9), (A.10) and (A.11) and substituting into (A@, we find L 1 L nDb*# (1 + B + 2Bw’)M + (1 - B + 2W)N r (1 + B)M -t (1 - P)N
(A.13)
where w’ = a*]bP is the fraction of the cell volume occupied by the hole, and M = Z,(~b)K&ra) - Zl(m)Kt(Kb) N = /u?Z,‘(txa)~,(!rb)- Z&b) .
KUK,‘(KU).
Now in Section 2 it was remarked that the source 4 must be normalized to unit rate of production in the cell, so that r(b* - a*)q = 1. Combining
this with equation (A.]), we have rrDbL$ =
D (1 -
w’)C, = (1 -
1 W’)2
=-
L” 1 -
w’
where L,B = l/G = D/C, is the migration area in pure moderator. The radial streaming ratio S, is defined by S, = Lr2/L,,*, so that our final result is S,_tl
+/9+2LIw?M+(l -p+2/%‘)N (1 - w’)((l + p)6)M+ (1 - B)N}
’
(A.14)
11
This result is new. However if the moderator very weak K tends to zero, and we find
absorption
M=
1 Y (1 - W) i O(KW log Kfl) 24/w’
N=
-$
(1 + w’) +
0(K2d
log
is
Ka).
Equation (A.14) then takes the very simple form (1 + BW’) Sr = (1 - w’)(l - Pw’)’
(A.15)
This result has previously been deduced by CARTER (1960). Indeed, combining (A.15) and (A.12) we derive an equation identical to CARTER’Sequation (3.15). This provides valuable confirmation of CARTER’Sanalysis. CARTERhas had difficulty in finding material to compare with equation (A.15). There are very few relevant experiments, and some of these are not of the highest accuracy. Also before comparing this formula with experiment it is necessary to verify that the model is adequate and that effects such as anisotropic scattering and finite pile size are not important. There are also a few Monte Carlo calculations, notably those of SCHAEFERand PARKYN (1958) and ASHFORD(unpublished). These have been done in lattices of low absorption (since the effect of absorption on S, was previously unknown). The lifetime of each individual neutron in these calculations is therefore long, the number of individuals which can be tracked rather small, and the statistical uncertainty large. The number of collisions which are recorded is large but since (as BENOIST(1959) and ASHFORD(unpublished) have emphasized) there is substantial correlation between paths, this does not help: it is necessary to track the neutron from birth to death. Equation (A.14) offers an escape from this difficulty. The absorption can now be increased by a large factor to shorten the neutron life, and the results compared with this equation. If the comparison is satisfactory, the equation can be used to extrapolate to real situations of low absorption. Some work on these lines has been started by ASHFORD.
I I.
Pergamon Press Ltd. Printed in Northern Ireland
THE WEIGHTING OF DIFFUSION COEFFICIENTS IN CELL CALCULATIONS D. c. Atomic Energy Establishment,
LESLIE
Winfrith, Dorchester,
Dorset
(Received 12 April 1961)
Abstract-The new theory of leakage in heterogeneous systems introduced in this paper is based on BENOIST’S lemma that the leakage from a region of such a system is equal to Bat*, B* being the geometric buckling and LBone-sixth of the mean square distance travelled by a neutron between birth and death, provided B1 is small. Lzis calculated exactly and directly and is shown to depend on a new type of flux [arising from a tilted source: this flux can be calculated in one lattice cell, so that it is no longer necessary to consider migration over the lattice as a whole to compute leakage. The exact calculation of t is troublesome, but a diffusion theory calculation can be carried through quite easily. From this, it appears that the radial and axial migration -areas are given by L,* = &/z, and Lz2 = 0,/C., C, being the normal cell-averaged absorption cross section. 0, is obtained by weighting the local diffusion coefficient with the usual flux, and & by weighting with a new flux derived from 6. There is thus a migration area asymmetry even in the absence of holes, and an example suggests that in some lattices’ this asymmetry may be quite substantial. The new theory may also be applied to more conventional streaming problems, and this application is discussed in the Appendix.
1, INTRODUCTION
BOTHthe absorption cross section C, and the diffusion
coefficient D vary across the lattice cell of a typical reactor. Feinberg-Galanin technique is much too complex for survey calculations, and the first step in such calculations is therefore to homogenize the lattice cell by inferring ‘paste parameters’ c, andii. za is computed by weighting the Za in each region with the local flux, but there is no general agreement as to how b should be calculated. Some authors have advocated weighting D with the local flux, others have advocated weighting E:,, and others still have advocated some very fancy prescriptions. The problem has been solved in general terms by BENOIST (1959) but his prescription contains collision probabilities which are very difficult to calculate. In this paper, a new and as yet unpublished theory of streaming due to the author is applied to the problem. On one-group diffusion theory, the leakage from a homogeneous region is given by B2L2, B2 being the geometric buckling, and L2 being one-sixth of the mean square distance travelled by a neutron between birth and death in an infinite lump of the same homogeneous medium. BENOIST (1959) has shown on the basis of one-group transport theory, that the result may be generalized to a uniform heterogeneous lattice, provided the buckling is small; it gives the leakage as F (&%3
where Lk2 is one-half the mean square of the k component of the vector joining the points of birth and death (k = x, y, z), Bk2 is the component of ihe geometric buckling in the k-direction. If B2 is large; the heterogeneity produces buckling-dependent corrections to L2. The author’s new theory gives a closed expression for L2 in a lattice, based on exact one-group transport theory, in the limit of small B2. b is then defined as L2za, since the essence of any definition of D is that it should give the leakage correctly. The paper therefore begins with an account of this new theory. 2. THE
CALCULATION
OF L”
From now on, we restrict ourselves to the usual cylindrical lattice and to a single neutron group. It is infinite in all directions, being uniform in the z-direction and periodic in the x- and y-directions. Let P(r,r’) be the flux at a point r due to a unit source of neutrons at r’ [r is the vector (x, y, z), the origin being arbitrary]; it is therefore the Green’s function for the lattice. The rate of absorption at the point r EC,(r) being the macroscopic is just &(r)P(r,r’), absorption cross section. The mean square distance travelled by neutrons from this source is therefore sss m
&(r)P(r,r’)(r
- r’)2 d%
the integral being taken over the whole of the lattice,
D. C. LESLIE
2
that is, over the whole of space. To obtain R2, the overall mean square distance for the lattice as a whole, we have to weight this according to the distribution of sources q(f) in the lattice. Thus
becomes
lengthH d3r’C,(r)P(r,r’)q(r’)(x - x’)~ NW116
Q) +
unit
B. (2.6)
B is a contribution from the boundary of the region (N,H)introduced by the change in the limits of
where the second integral is taken over one cell in the x - y plane and over unit length in the z-direction. The source function q(d) is normalized to a unit rate of neutron production in this element of volume. For example, in a metal-fuelled lattice, q will be constant in the moderator and zero in the channel and in the fuel: it will be equal to l/V,, V, being the ‘volume’ (actually an area) of moderator in one cell. In a homogeneous medium, we would put 2 = 6L2, L2being the thermal diffusion area. In a lattice, there may be asymmetries in the diffusion area, and we must put i? = 2L,2 + 2L,2 -+-2L, where
(2.2)
integration. It can be made relatively as small as we wish by taking the region to be large enough. Jgnoring this contribution, equation (2.6) may be written length
LZ21Lr2=-
N
1
2NH sss
d%$Jr){$)(r)
- 2x@(r)
N cells
+
x24dr)I
(2.7)
where (2.8a)
J$) (r) =
x’q(r’)P(r,r’) d%' sss m
and
7-g)(r) =
(2.8b)
xr2q(r’)P(r,r’) d31.'.(2.8~) US
03
unit
Similar expressions may be derived for L,"and Lz2. The significance of the suffix 0 will be made clear in (2.3) the next section. Because of the definition of P(r,r’), #0(r) is simply the definition of Lv2and Lz2being similar. In lattices the flux due to a source distribution q(r’). Since this with square, hexagonal or trigonal symmetry is the actual source distribution in the lattice, &,0(r) must be the actual flux. &j(r) is a new type of flux, L,2= L,2= L,2 (2.4) due to a ‘tilted’ source distribution x’q(r’). Figure 1 shows a comparison of a tilted source with the normal the ‘radial’ diffusion area, while Ls2 is the axial source in a one-dimensional lattice. This type of diffusion area: in general, Lr2is not equal to Lz2. flux would appear to be new, although SELENGUT Equation (2.3) represents the fate of neutrons starting from one unit volume (unit length of one cell) Tilted Conventional source source ; and migrating over the whole of space. We could, if we wished, start the neutrons from a large region ccl-Nentlonor of length H in the z-direction and containing N cells. source We would then have
7L
-Moderator --)
the neutrons still being tracked over the whole of space. We now come to the vital step in this analysis. Instead of being start&d in the region (N,H) and tracked over the whole of space, the neutrons are started everywhere, a death only being counted if it occurs inside the region (NJ-Z). Equation (2.5) then
Fuel ------+
--Moderator-+
FIG. l.-Illustration of conventional and tilted sources in a circular lattice cell. (Seclion along the x-axis.)
The weighting of diffusion coefficients in cell calculations
I -
J Fuel
+ /
Moderotor i
-
Mrxierator-
problem was almost intractable if heterogeneity was taken properly into account, since the lattice had to be treated in its entirety. Now the calculation of leakage has been put on the same footing as the calculation of thermal utilization, in that both computations can be done in a single cell, the only restriction being that the geometric buckling B2 should be small: this restriction applies equally to the calculation off. 3.
FIG. 2.--Illustration of the fluxes % and tin a circular lattice cell. (Section along the x-axis.) b = 2a; D, = 40,; No absorption.
(1960) has noticed the relevance of flux tilt to leakage problems. It is shown in the next section that x$)(r) = .$&r) - Et)(r)
(2.9)
where 4,(r) is the ‘normal’ flux defined by (2.8a) and [f)(r) is another new flux which, like q&(r), is periodic in the lattice. It is antisymmetric in the s-direction, symmetric in they-direction and uniform in thez-direction. Its behaviour in a one-dimensional lattice is illustrated in Fig. 2. It is also shown in the next section that $‘)(r), which is the flux due to a ‘curved’ source .u’2q(r’),is given by 7:‘(r) = s*&(r) L 2s$Yr)
- $‘(r)
(2.10)
where v:‘(r) is yet another new flux, which is also periodic in the lattice cell. It is symmetric in both the x’ and y’ directions (though its dependence on the two co-ordinates is different) and uniform in the z-direction. Substituting from (2.9) and (2.10) into (2.7), we find
3
SIMPLIFICATION OF THE FORMULA
BASIC
We shall now show that the flux 1’can be eliminated from the basic formula (2.12), which can be expressed in terms of ,t only. This reduces considerably the labour of using (2.12), since (2.10) shows that it is necessary to know 5 before 1’can be calculated: one step in the work can now be eliminated. For this further step we need to consider the equation satisfied by 1~. We shall take this in its exact or transport form, so that the work of this section does not involve any further approximations. BENOIST (1959) prefers to use the transport equation in its integral form, but we find the integro-differential or Boltzmann form more convenient. We write this equation in a tensor notation, in which
is the co-ordinate vector, normally denoted by r or (s, ,v, r): thus x1 is to be identified with X, xq with _r and s3 with :. Similarly {~4~} is a unit vector giving the direction cosines of the neutron path. It is often denoted by Q, and since it is a unit vector, it must satisfy
From this it will be seen that we are using the summation convention : a repeated index is to be summed over 1, 2, 3. With this notation, the Boltzmann equation for a vector flux + may be written
(2.11) ui x
(-yj,uI;)
7
c,(xj)b(sj,uk)
- 1
Since t’?)(r) is periodic, this integral can be taken over one cell only, giving as the basic formula of the new theory unit rrr Lx2 = 3
JJJ
C,(r) v?‘(r) d3r.
(2.12)
rell
Similar formulae can be derived for L,2 and Lz2. The advantage of this new procedure is that it reduces the calculation of leakage to the computation of a function in a single lattice cell. Previously the
=
C,(.Uj;U,U,')~(Xj,Ui')
d*Um’+
q(XjvU&
(3.1)
ss
Here Ct(xj) is the total cross section at the point (xj‘; &,(.x~;u,u,‘) is the scattering cross section for deflections through an angle cos-‘(u,u,‘) at the point {_uj} q(xj,uk) is the source of the vector flux I$(x~,u,). This equation brings out explicitly the dependence of the vector flux on direction as well as position. The first term on the right-hand side of (3.1) represents
D. C. LESLIE
4
the contribution to a given directed flux at a particular point from neutrons moving in other directions which suffer a suitable collision at that point. SJ d2ui’ represents an integration over all possible directions of motion of the incident neutron, and is normalized to 1. In computing reaction rates, a full knowledge of the vector flux (x,, uJ is not required. These are determined by the ‘ordinary’ flux
which is the average of the vector flux over all directions of motion, and is a function of position only. The suffix 0 in the formulae of the previous section indicate that it is this average flux which is involved. In this paper we shall assume that the scattering cross section has a simple cosine anisotropy only, so that X,(X,; I(k1lk))= E,(Xj){l
+ 3E(Xj)Uj$~)a
(3.3)
Here & is the usual isotropic scattering cross section, and p is the average of the cosine of the scattering angle, which would be zero if the scattering were purely isotropic. In a more accurate treatment there would be further terms on the right-hand side of (3.3) involving P,(u,u,‘) and perhaps higher-order spherical harmonics. However (3.3) will suffice for illustration, and the conclusions of this section are not affected by the assumption. Substituting from (3.3) into (3.1) and suppressing implicit variables, we finb
equation (3.4). We then find that 5@)satisfies uj y
+
c,5 @)- C,(@’ + P&r,+)} = -ui+
X,
ax1
This follows since axi = 6ii, the Kronecker
(3.7) delta
symbol, which is 1 if i = 1 and 0 otherwise. (3.7) is a normal Boltzmann equation, showing that 5(x) is derived from the anisotropic source -ui$, Since 4 and the cross sections are periodic in the lattice, it follows at once from (3.7) that 5(x) is also periodic. In exactly the same way, we find that Y(X)satisfies
adz)
uiaxi +
Cp
- lqvp
+ puiYi(+)} =
-2up
(3.8) v@) and the BY) also being defined by equations similar to (3.2) and (3.5). Equation (3.8) shows that Y@)must also be periodic in the lattice. We now integrate (3.8) over all directions of neutron motion, and use the definitions of VP’, Y!*) I and El”): we find I
3
a+)
awYj
+
I;&)
=
2
r”’
(3.9)
3 l
where &, is of course equal to Zc, - Z,: the second term inside-the curly brackets in (3.8) gives no contribution since the average value of ui is zero. Thus from equation (3.9), equation (2.12) may be rewritten unit
unit
L,” =
-2
Z(r) sss ! sss p(r) d% -
CC11
6
cell
t
dyr
(3.10)
where (bOis defined by equation (3.2) and the vector +i, which is also a function of position only is defined by +ifxj)
=
3.f j”iNxj4
@%2
Using Green’s theorem, the second term is equal to
(3.5)
+i is proportional to the neutron currentj,, the exact relation being jd = &$, u’being the average neutron velocity in the single group. We shall take (3.4) as the standard form of the Boltzmann equation for a flux 4 due to a source q: a more extensive discussion of this topic is given by DAVISON(1957). Thus the gux x(~), which is derived from a source qx, must satisfy
# and the xi(*) being defined by equations similar to (3.2) and (3.5); in (3.6), we have represented x as x,. In this equation we write x(*) as x1+ + P, and use
This is proportional to the current out of the unit cell of neutrons having the distribution v@). Since v,,tr) is symmetric in the x- and y-directions and uniform in the z-direction, this current is zero. Therefore the second term of (3.10) vanishes, and we have unit’ L,2 = -
f
SIJP(r)
dar.
(3.11)
cell
This is the simplification which was promised at the beginning of this section. Similar equations may be deduced for L,"and LZ2: provided that B2 is small, these equations are all exact.
The weighting of diffusion coefficients in cell calculations 4. THE
CALCULATION OF THEORY
special ‘flux’
5 BY DIFFUSION
The calculation of I# by transport theory is a formidable problem even in the simplest geometries, and the calculation of 6 is even more difficult. It therefore seems profitable to start by calculating L2 from equation (3.11) using diffusion theory, and we shall now do this. It must be remembered that diffusion theory is not very accurate in real lattices, and that a better procedure is likely to be needed. However if 4 has been calculated by diffusion theory, it seems perfectly reasonable to compute E in the same way. Strictly, we are going to use the PI procedure rather than diffusion theory, since the former can always be carried through in a consistent manner. In this procedure, we put 4 = 40 +
(4.1)
%di
and E(I) = ,$’ + n.EW. z,
5
(4.2)
Here &, and +i are defined by (3.2) and (3.5) respectively, the definition of $1 and @) being similar. Substituting from (4.1) and (4.2) into (3.7) we find
(4.7) In fact, (4.6) may be written L,” = b,. z,
(4.8)
where unit d,
=
unit D(dUr)
us
d3r
L(r) II(Sss
Cdl
d+)
(4.9)
Cdl
and unit z,
=
S,(r)
d3r.
(4.10)
sss
Cdl
Here D* is a cell-averaged diffusion coefficient, the weighting function being t,(r). The form (4.8) suggests that we should investigate the quantity Z,. As an illustration, suppose that the lattice cell is square, the half-width being a. Since the ‘unit cell’ is of unit height, we have
The integral of the second term in the inner integral is just @‘(Xl = a) - @)(x, = -a).
- &{@’ + Ff4&@) = -z&
- urujl&. (4.3)
We now multiply this equation by u, and integrate over all u. Using the obvious results
But .$’ is antisymmetric in the x,-direction and periodic in the cell, so that it must be zero on the faces x1 = a and x1 = --a of the unit cell. It follows that the second term in (4.11) is zero, and that unit
I&
d2U +u,u,dpu
= 0,
(6, being the Kronecker viously) . we find g’
=
(.&
d% = g aij
(4.12)
delta symbol defined pre-
_‘,&(A+$!g.
But the Pr approximation is just
D=
]/Wj
(4.4)
to the diffusion coefficient 1
3(% - $3
*
(4.5)
Substituting from (4.4) and (4.5) into (3.11) we find
a very suggestive formula. This formula is not exact. since it depends on the use of the PI approximation. Equation (4.6) shows that Lz2 is obtained by weighting the local diffusion coefficient D(r) with the
From this equation, it follows that Z, = Z, = Z,. Now in the Introduction, we remarked that the average absorption cross section in the lattice is obtained by weighting with the normal flux &r), that is unit uuit c, = V%Mr) d3r )/(///M) @rJ. (4.13) (IIS cell cell The numerator of this expression is the mean rate of absorption in the cell. This must be equal to the mean rate of production, which is 1, from the definition of q(r). It follows at once that z, = l/E,
(4.14)
and therefore L,” = DJZ,.
(4.15)
This is, of course, just the sort of relation which we
D. C. LESLIE
6
had hoped to find. It shows that the averaging process defined by equation (4.9) is compatible with simple definitions of migration area and diffusion coefficient. In exactly the same way, we may show that L,2 = I@,
(4.16)
I), being obtained by weighting D(r) with
holes. If diffusion theory @ valid, the axial diffusion coefficient b, is obtained by weighting the local diffusion coefficient D(r) with the normal flux, and the radial diffusion coefficient b, by weighting with the special flux c*(r) defined by equation (4.7). The corresponding migration areas are obtained by dividing the appropriate diffusion coefficients by the cell-averaged absorption cross section. 5.
In a square cell we shall obviously find (4.17)
D, = D, = b,
the radial diffusion coefficient. We would expect that the weighting function for the axial diffusion coefficient would be
this following from the equation
which is analogous to (4.4). Integrating the equation for P) analogous to (4.3) over all U, we find
f.as:“~C,J(;‘ -;$3=Dz = ay 3
a.Yi
=D
2)
(4.19)
The right-hand side of (4.19) is zero, since the reactor is uniform in the z-direction, and #+, is therefore independent of z. Combining (4.18) and (4.19), we find that $1 satisfies
CALCULATION
OF THE ASYMMETRY
Thus the axial diffusion coefficient is obtained by weighting the local diffusion coefficient with the normal flux y$(r) unit
unit
B, =
(4.20) Cdl
f3 = 0 being one of the axes of symmetry of the cell (other terms in the Fourier series being excluded by symmetry) and then assuming that g,r is the dominant term. From the remarks in Section 2 following the definition of EF’, the boundary conditions on this flux are* ” on
ap
ay@) =O On
P _x=---
-2 +;.
From this it follows that the Fourier expansion of lt’ is (2n + 1>e. More terms are now allowed, since the diagonals are no longer axes of symmetry. If we take only the first term of this series, namely j,(r) cos 6
Cdl
the axial migration area being then given by L,2 = D,/&
AREA
Having discovered that d, is not identical with b,, we must see whether the discrepancy is large enough to be important in practice. We do this by calculating both quantities for a typical lattice cell. Since we are using diffusion theory, this cell must be free from voids. We take it to consist of a fuel rod of radius a, absorption cross section C,,,, and diffusion coefficient D,. embedded in a block of moderator in the form of a square of side P. The moderator diffusion coefficient is D,, and moderator absorption is negligibly small. It is conventional in normal flux calculations to replace the square by a circle of radius b = P/15. This is equivalent to expanding the flux in a series of the form d,(r) = &g,(r) ~0s 4~~9
Eg)(r) = 0
since D is also independent of z. From this, it follows that Et) is also independent of z, and therefore that
MIGRATION
(4.21)
Thus we may conclude: There is a migration area asymmetry in reactors with cylindrical symmetry, even in the absence of
the calculation of Eg’ can be cylindricalized in the same way as that of &,. The author (unpublished) has calculated 6 in a cell with square boundaries and has * These conditions follow at once from 5 being continuous, periodic in the lattice, anti-symmetric in the x-direction and symmetric in the y-direction.
The weighting of diffusion coefficients in cell calculations
thus been able to show that as long as the fuel diameter is not more than half the lattice pitch, this cylindricalization introduces very little error. It is used here to simplify the work. With this approximation, the boundary condition on j,(r) is j,(b) = 0
so that, using the boundary condition (5.3) L,2 = nD1 ob2rg(r) dr + (Do s
x {nab(a) - na2g(u)} + 27 l&r)
~&j(r) = x+,(r) + $?(r) = r cos 0$,(r) + @j(r)
(5.2)
b, = D,+ @h(a) (Do
-
4)
- m2g(a)} -I- Zn[rg(r)
x
rb
J0
(5.7a) Similarly, from (4.20), the axial diffusion coefficient is given by 2n ‘rg(r) dr s
D, = D, $ (Do - D3
Our task is therefore to calculate the functions g(r) and h(r) (we can drop the suffices since we are only keeping one term in each series) and to average them suitably. Now in the diffusion approximation unit B,
=
’Wr)Ur)
1, SSJ
.
r”b
h,(r) cos 6.
(5.3)
dr
27r -rg(r) dr
2x
h,(b) = bg&) [and h,(O) finite].
(5.6)
whence
rather than with @‘. Following the line of argument used above, we approximate to $)(r) by
From (5.1) and (5.2) the function h,(r) must satisfy the boundary conditions
dr]
[
(5.1)
b being, as before, the outer radius of the equivalent cylindricalized cell. The general theory has been formulated in terms of the flux 5, since this is periodic in the lattice while x is not. Now that we have reduced the problem to one of calculation in a single cell, this distinction is not important, and it is convenient to work with
Dl)
J0
(5.8a)
rg(r) dr
The migration area asymmetry is therefore due to the term in curly brackets in (5.7a). These formulae can be simplified by introducing two new symbols E: the ratio of the mean flux in the fuel to the flux on the interface r = a. S: the ratio of the mean flux in the moderator to the flux on the interface. (In the usual notation, E = l/G). Using these new symbols we find
d3r
b, = D, + (Do - D3 ;‘E”+“v”E;
cell
0
(5.7b)
1
from (4.6) and (4.15) where
axf) ag) L.(r)= +o(r)f (r> = (r>ax1 ax1
b,
i$O (4 1
x1 -
from (5.2). Now
sothat 5,(r)=
-=a axI
cose---a
sin 0
ar
r
$ (30s~e + ” sin2 13-
I
D(r) = Do =D,
44
X = --
a&)
‘OE VoE 4 v,S
1; Vo=rra2;
(5.9)
ae COG
if
0 < r
if
a
(5.8b)
Vl=7r(b2-a2)
a
0 2
(5.4)
Since D(r) in fact depends only on r, and C,(r) is independent of z, we have
Now
where
D1 C (Do - 01)
=
The calculation of E and S is quite standard, but the computation of X is in general rather tedious. However if the absorption in the lattice cell is negligible, this latter computation is quite easy: provided Do + D,, this case is not trivial. Since x@)is the flux produced by a source qx,, it follows that in the diffusion approximation the isotropic component &) satisfies -gi(D(r)z) Since we are putting
i-&&)=4X1. xe) = h(r) cos 8, we find that
D. c.
8
h(r) satisfies -- l&)(r):)
(5.10)
+ (X.+$+.
LESLIE
This will be zero if D,, = D,. Since the mean free path in the fuel is presumably always less than in the moderator the asymmetry will be greatest for D,, = 0, when it is given by
If the absorption is zero, the source must also vanish and in a region where D is constant, h(r) = Ar + ;“. We shall take this form for h in the moderator region a < r < b: in the fuel region 0 < r < a the second term must be excluded, so that h(r) = Cr, 0 < r < a.
The boundary conditions are that h(r) and D dh/dr should be continuous at the interface r = a, together with equation (5.3). Since there is no absorption, the ‘normal’ flux g(r) must be constant, and we can write this equation h(b) = b&a).
Eliminating the arbitrary find that
h(a) f&a)
-=x+l=(l_v)D,+(l+v)Dl
constants A, B and C we
24
=
’
(1 + VP, + (1 - W,
Since v may be as large as 1 (equal volumes of fuel and moderator) it would seem that asymmetries as large as 50 per cent are possible. It should be noted that the asymmetry is considerably reduced if Do is not actually zero. For example, (5.13) gives, for v = a, a = 20 per cent if Do = 0. If D, = ) D,, a is only 9.9 per cent: this reduction is larger than might be expected. Asymmetries of this size would have considerable of practical consequences for the interpretation exponential experiments and in the optimization of small power reactors. There is therefore a need for a more precise calculation in which allowance is made for absorption in the fuel. For simplicity, we shall suppose that the source is zero in the fuel and flat in the moderator: we shall also ignore moderator absorption. Using the diffusion approximation in the fuel as well as in the moderator, the flux rise parameters are given by the standard formulae E
D
(1 - v)D, + (1 + v)D, "
21,&d
=
(5.14)
Kd,( KU)
.y = 1 + ; (W7)”D”E D 1
3b2 - a2 b4 b (b2 - a2)2 In ii - 4(b2 - a2)I
(5.12)
This formula has already been deduced by SELENGUT(1960) using an entirely different method. Since the physical models used by SELENGUTand by the author are both reasonable (although quite different) it is highly satisfactory that they give the same answer. The author’s method would seem to be the more general, since it covers the case of finite absorption, it can be used with more exact methods than diffusion theory, and it has revealed an asymmetry in the migration area. For this case of zero absorption, equation (5.8b) gives b, = vD, + (1 - v)D, so that the migration area asymmetry is L2 a=-4_1=f_1 L,”
(5.13)
1+ 0
(5*11)
where v = u2/b2 is the volume fraction of the cell occupied by the fuel. We can now substitute back into (5.7b), remembering that E = S = 1 since g(r) is flat. We find D
L2
orZ~_llL. L,2
b D,
(1 - v)D, + (1 + v)D, (1 - v)D, + uD, _ I = (1 + u)D, + (1 - v)D, ’ D,
(5.15) where, as usual $=3X
a0Z tr.0.
We shall also need to know the total flux rise in the moderator, which is specified by T=
(5.16)
1
We can now calculate the ‘tilted flux’ h(r), which is the solution of (5.10). In the fuel, where there is absorption but no source, the solution is h(r) = CIl(Kr)
O
the second solution K,(Kr) being excluded because it is singular at r = 0. In the moderator there is a
source but no absorption and the appropriate solution is h(r) = Ar + 5 - &qr3
acr-cb.
1
The constants A, B and C are determined
by the
The weighting of diffusion coefficients in cell calculations
continuity and boundary conditions
Jf-
cr,(Ka) = Au +
&$
(1
Ql
Kc&‘(KU)
=
A
-
;
(5.17) 1 (5.18)
-
B bg(6) = bTg(a) = Ab + h - $
and
qb3
(5.19)
1 The total source in the cell must be equal to the total absorption, so that n(b2 -
a2)q =
?ra2.
= n-(~a)~D, . E&a)
&&(a)
9
this new theory is to the calculation of cell-averaged diffusion coefficients, this averaging being done with a new kind of flux. This flux should, strictly, be calculated using transport theory and such a calculation will clearly be lengthy and difficult. To obtain some preliminary results, the flux has been calculated by diffusion theory. Provided diffusion theory is valid, the cell-averaged axial diffusion coefficient is obtained by weighting the local diffusion coefficient with the normal flux, while the averaged radial diffusion coefficient is derived from a new flux due to a ‘tilted’ source. There is thus a migration area asymmetry, even in the absence of voids.
(5.20) Eliminating C between (5.17) and (5.18), we find A--$-
where
&qa2=Z
A+;;,C
1
B
1
1 5
1
4a2 \
Z=--
‘I
[from the recurrence relations satisfied by the function I,(z)1 or (5.21)
I )
j
Finally, eliminating (5.20), we find
A
B
and
and using the relation
0
2T-t =x+1=
2
E(Kcz)~(~ -
(l+t?;+z(l-U)
u)
(5.22)
where, as before, u = a2/b2. If K = 0 and E = 1, the above equation becomes identical with (5.1 l), as it should. The formulae given above enable d, and b, to be computed from equations (5.7b) and (5.8b). Mr. H. M. Brown has carried out such a computation for the special case v = 1, Do = tDl, and has deduced the variation of the migration area asymmetry ccwith KU: his results are shown in Fig. 3. It will be seen that realistic amounts of absorption do not have much effect on the asymmetry. When the absorption is large, the statistical weight of the fuel will become very small. When this occurs the asymmetry must vanish, but the figure shows that this only happens if the absorption is extremely large. A priori, one might expect that the asymmetry would decrease monotonically to zero : the figure shows that this is not so. 6. CONCLUSIONS
I
4
6
8
IO
12
14
16
18
20
K0
FIG.
h(a) a&)
2
AND
FURTHER
WORK
A new theory of neutron currents in a heterogeneous assembly has been formulated. One application of
3.-Variation
of migration area asymmetry with fuel 1 DO 1 absorption for a lattice cell with a = b
2’%=4’
The conclusions of the previous paragraph are only valid if the individual regions of the lattice cell are large enough for diffusion theory to be applicable: they certainly do not apply if the cell consists of large numbers of very small regions, or if there are voids. It is most desirable to remove this restriction by calculating the required fluxes more accurately. This programme of work has already been started. In Appendix 1, we show how Lr2 can be calculated for a cell consisting of a circular hole surrounded by moderator, using boundary conditions worked out by CARTER and JARVIS(1960). The analysis leads to a formula already deduced by CARTER(1960), who has concluded that it is in good agreement with available evidence. At present, the theory is based on BENOIST’lemma S which relates leakage out of a large region to the mean square distance travelled by neutrons inside the region. It is hoped to reformulate the work of Section 2 so that leakage can be calculated directly. With this improvement, it may be possible to study flux matching between different zones of a reactor,
D. C. LESLIE
10
and the dependence spectrum effects.
of diffusion
coefficients
on
Acknowledgments-The author is indebted to M. P. BENOI~Tof C.E.A., Saclay, for making his work on streaming available in advance of publication and for helpful comment, to DR. C. CARTER of A.E.R.E. for many interesting discussions, and to MR. H. M. BROWNof A.E.E. for help with computation.
arguments can be used about 7 which, unlike v, is a ‘real’ flux derived from a real source qx,“). It follows that $‘(a)
C.E.A. 1354. CARTERC. (1960) Streaming Due to Holes in a Reactor. AERE R-3367. CARPER C. and JARVIS R. J. (1960) Zntuence of a Cylindrical
Press. SCHAEFERG. W. and PARKYND. M. (1958) Proceedings of the Second International Conference on the Peaceful Uses of Atomic Energy, Geneva, A/CONF/l5/P.310. United Nations, N. Y. SELENGUTD. S. (1960) Trans. amer. Nucl. Sot. 3, No. 2, 398.
q = &#.
(A.11
Since Z,, is zero in the hole, the basic result (2.12) for the radial migration area may be written &Jr)+‘(r)
and this removes the flux v, which would be difficult to calculate. Now X’“’= X’$ + 5:“’ and since the +i are zero, ,iz’ must be equal to$‘. may be written
mod
the integral being taken over the moderator only. The diffusiontheory methods of Section 4 do not apply to cells containing holes, but CARTER(1960) has shown how the transport calculation of DAVISON (1957) may be combined with diffusion theory to confirm the Behrens value of L,” for this case. It would no doubt be possible to translate CARTER’Scalculation into the language of Section 3 of this paper. The transformation of Section 3 reduces the previous equation to sss
mod
x1$‘(a)
L,a = )
X:‘(r) d%.
dS, - +
i
Thus (A.3)
sss
+
@‘(a)
i
(A.4)
mod
h0le
We shall now follow the procedure of the main text, and will calculate x1”’by diffusion theory in the moderator: equation (A.4) has been put into a form in which it is not necessary to know x(Z) in the hole. We begin by cylindricalizing the problem, replacing the square outer boundary of the cell by a circle of radius b = Zy&. of the form
The isotropic component
of the flux is then (A.5)
so that, in the diffusion approximation $
2:“’ = -30
~09 8 + 2 sine e
r
,
(A.6)
The first term of equation (A.4) is simply - 1
2VxF’(a) . aa cos 0 d0 s0
where Xj”‘(r) = -30
dh G cos 0
64.7)
(in the diffusion approximation) is the current of xtz’ in the radial direction. Substituting into equation (A.4) from (A.6) and (A.7), we find
d3r
sss
-
(A.3)
mod
xb”’ = h(r) cos f3
We now want to consider a lattice cell of a type slightly different from that treated in Section 5. It consists of a square block of moderator of side P, containing a concentric circular hole of radius a. The moderator absorption cross section and diffusion coefficient are &, and D respectively, and themoderator contains a constant source of intensity q per unit area. The conventional flux in the cell is constant and is given by
l:“‘(r) &
t:.)(r) dJr
It
sss
1
Streaming in a cell containing a circular hole
L> = -)
dS* - )
hole
Channel on a Neutron Flux which Varies Linearly in a Direction Perpendicular to the Axis of the Channel. AERE R-3366. DAVISONB. (1957) Neutron Transport Theory, Oxford University
= 4
x&“‘(a)
I
Coeflcient de Diffusion darts un Reseau Comportant de Cavites.
Ip
hole
hole
L;L = f
ASHFORDJ. R. Unpublished material. BENOISTP. (1959) Formulation GPnt?rale et Calcul Pratique du
+ p:“‘(a)} dSi = 0.
f
But $ is constant and the #Q are therefore zero everywhere. follows that
REFERENCES
APPENDlX
{x?+<(a) + 2x&“(a)
d& =
I
dSi
(A.2)
hole
second integral being taken round the channel wall: it represents the current of the flux Y out of the moderator into the hole (the contribution from the cell boundary vanishes). Now the net current of the flux 7 into the hole is zero, since there is neither absorption nor production within the hole (such
L:=D~e~rdr~~Ms’B+~sin’Bj 2=al dh z (a)
+D
COS*e
de
s0 = SOD bh(b) - ah(a) + a’ z dh (a)
1.
(A.81
Now h(r) must satisfy equation (5.10), with constant C,, D and q. It is therefore given by h(r) = AI,
+ BK,(ur)
the
4 + - r c.
(‘4.9)
where K’= 3B,&, the A and B being determined by the boundary conditions. On the outer boundary we have simply h(b) = b+
(A.10)
The weighting of diffusion coefficients in cell calculations from equation (5.3). CARTER(1960) has remarked that the most general form of condition on the inner boundary is h(a) =
5;
af
(a)
(All)
and that all the information about streaming across the hole must be contained in the value of ,!J. CARTERand JARVIS(1960) have discussed the calculation of B and have concluded that the formula 1 + 2aC,, tA.12) P= 2(1 + a&,) should be quite accurate for all sizes of hole. Eliminating A and B from equations (A.9), (A.10) and (A.11) and substituting into (A@, we find L 1 L nDb*# (1 + B + 2Bw’)M + (1 - B + 2W)N r (1 + B)M -t (1 - P)N
(A.13)
where w’ = a*]bP is the fraction of the cell volume occupied by the hole, and M = Z,(~b)K&ra) - Zl(m)Kt(Kb) N = /u?Z,‘(txa)~,(!rb)- Z&b) .
KUK,‘(KU).
Now in Section 2 it was remarked that the source 4 must be normalized to unit rate of production in the cell, so that r(b* - a*)q = 1. Combining
this with equation (A.]), we have rrDbL$ =
D (1 -
w’)C, = (1 -
1 W’)2
=-
L” 1 -
w’
where L,B = l/G = D/C, is the migration area in pure moderator. The radial streaming ratio S, is defined by S, = Lr2/L,,*, so that our final result is S,_tl
+/9+2LIw?M+(l -p+2/%‘)N (1 - w’)((l + p)6)M+ (1 - B)N}
’
(A.14)
11
This result is new. However if the moderator very weak K tends to zero, and we find
absorption
M=
1 Y (1 - W) i O(KW log Kfl) 24/w’
N=
-$
(1 + w’) +
0(K2d
log
is
Ka).
Equation (A.14) then takes the very simple form (1 + BW’) Sr = (1 - w’)(l - Pw’)’
(A.15)
This result has previously been deduced by CARTER (1960). Indeed, combining (A.15) and (A.12) we derive an equation identical to CARTER’Sequation (3.15). This provides valuable confirmation of CARTER’Sanalysis. CARTERhas had difficulty in finding material to compare with equation (A.15). There are very few relevant experiments, and some of these are not of the highest accuracy. Also before comparing this formula with experiment it is necessary to verify that the model is adequate and that effects such as anisotropic scattering and finite pile size are not important. There are also a few Monte Carlo calculations, notably those of SCHAEFERand PARKYN (1958) and ASHFORD(unpublished). These have been done in lattices of low absorption (since the effect of absorption on S, was previously unknown). The lifetime of each individual neutron in these calculations is therefore long, the number of individuals which can be tracked rather small, and the statistical uncertainty large. The number of collisions which are recorded is large but since (as BENOIST(1959) and ASHFORD(unpublished) have emphasized) there is substantial correlation between paths, this does not help: it is necessary to track the neutron from birth to death. Equation (A.14) offers an escape from this difficulty. The absorption can now be increased by a large factor to shorten the neutron life, and the results compared with this equation. If the comparison is satisfactory, the equation can be used to extrapolate to real situations of low absorption. Some work on these lines has been started by ASHFORD.