- Email: [email protected]
On calculations for cylindrical reactors
100
Letters to the editors
On calculations for cylindrical reactors* (Received 24 December
1956)
IN calculations of neutron multiplication factors and critical sixes for reactors of the shape of spheres or infinite cylinders using the age-diffusion approximation, a useful procedure is to use the multigroup method, and to solve the group equations by the method of difference factorization developed by MARCHUK.“) In these cases the neutron flux nu depends only on a single spatial variable and the Laplacian V2 has the form
v&La
-
rG(ar
a %r ( 1
(CL= 1, 2).
If the reactor has the shape of a finite cylinder, and the flux depends on two variables, r and z, then the Laplacian has the form
and the method of difference factorization is not immediately applicable for solving the group equations. This is because this method involves the inversion of a matrix formed from the coefficients of the unknowns in the one-dimensional group equations and essential use is made of the fact that this matrix contains non-vanishing elements only in the principal diagonal and the two adjacent diagonals. For two-dimensional systems, the corresponding matrices also have elements in the next two diagonals, five diagonals in all, in which case they cannot be inverted by the method of difference factorization. MARCHUK’S method of reactor calculation consists of two distinct stages. Firstly the group constants are calculated by averaging cross-sections over lethargy groups and secondly the group equations are solved by difference factorization. In the case of two spatial variables the first stage of the calculation is unaltered, and we shall not discuss it here. The purpose of this paper is to draw attention to the possibility of using an approximate Laplacian enabling the entire procedure of the second stage of the method to be applied. We start with the following equation for the slowing down of neutrons in an intermediate reactor
--
1
3&
agz,nu V%u + C,nv+ = vS(r, au
2, u),
s co
S(r, 2, u) = F(u)
where
C,nu du
0
Here C,,, I;,, EE., F and & are experimentally known functions of the lethargy u. The required neutron flux must be a finite non-negative function of u throughout the reactor and must vanish at the extrapolated boundary. At the boundaries between physically distinct media both nv and the current density 1/(3&)V,no must be continuous, where V, is the component of the gradient operator normal to the boundary. Also, for u = 0 we must have nv = 0. For the numerical solution of(l), we form a lattice in space given by ri = ih, and .zj = jk, and for simplicity we put hr = h, = h, and also assume the same value of h in all regions of the reactor. We also refrain from discussing possible lattice points at region boundaries. Likewise, we divide the lethargy scale into steps of length 1. For brevity we introduce the following notation nu(ih, jh, kl) = nvif’ ; r,ha A;dm$ ha Arant$
=
= nu:F]_l - 2nz$’ + nojI;‘,l;
/A,fCF’no:f’ = AXE”+”nt$+” - &Z:” no::‘.
The indices i, j and k have the following ranges:
l
by N. KEMMERfrom Aromnaya Energiya 3. NO. 7, 53 (1957).
nv
is to be
101
Letters to the editors
computed. Finally, we shall attempt to make the statement of the problem clearer by repeating each equation in matrix notation. We now state two approximations to equation (1) which, in a sense, are opposed to each other. (I) The explicit scheme of equations
or more briefly nylk+ll = ,zjU~‘~v” +
co
Q'k'
(II) The “inverted” scheme of equations
or more briefly S’“+Unv’L+” = nvlL’ + Q'k' (3) Here A”’ and Pfl’ are matrices of order N, x N, each row of which has no more than five elements differing from zero: nP’ is the (N, x NJ-dimensional vector which we have to determine and Q”’ is a similar vector that is given if the sources are known. A comparison of the systems (2) and (3) shows that the former set of equations is considerably simpler than the latter because (2) requires multiplication of the matrices P’ by a vector, whereas (3) requires the inversion of the matrices B lk+l’. However, in order that the numerically greatest eigenvalues of AIL’should be less than 1, as is necessary for stability, the lethargy step I must be small in comparison with the spatial step h. On the other hand system (3) is stable for any 1 but it requires the inversion of a high order matrix for each k. Thus both Method 1 and Method 2 require a considerable number of arithmetical operations; Method 1 because the lethargy step has to be kept small, Method 2 because of the slow convergence of the interaction procedure for inverting the B(*’(non-iterative methods for solving (3) are of no use). Bearing in mind the particular, rectangular, shape of the spatial region we are concerned with, we propose that such reactor calculations be performed using another method, which is free from these disadvantages. (III) Hybrid scheme of equations
or, more briefly
Cy+l’nv, 12ktl)= ~‘1Bk’,#k’ + ~‘2k+2’nv;2k+2’ a
=
&2k+l',$k+l'
+
Q;““’ Qy+l,
(i = 1, 2, . . . N,), (j
=
1,2, . . . Nz)
Here Ciartl’ and Ciaktz’ and likewise 0:““’and. Dp’+” are respectively matrices of order N, and N, which have no rows with more than three non-vanishing elements; the unknown vectors nv!” I and nvy’ are of orders N, and N, respectively and so are the known vectors Qj”’ and Qy’. Like Scheme II, Scheme III is stable but requires the inversion of matrices, but in place of the inversion of one matrix of order (Nr x N,) one has here to invert N, matrices of order N, or N, matrices of order N,. The advantage of III over II lies in the fact that the matrices Ci” and Cl” have elements only in three diagonals. This makes it possible to invert them exactly, instead of by iteration, by the use of factorization. In the number of arithmetical operations involved, Method III is surpassed only by the method of reduction to Cauchy’s problem,‘*’ which is not stable and therefore unsuitable for an electronic computer. V. K. SAULEV REFERENCES 1. MARCHUKG. I., J. Nucl. Energy 3, 238 (1956). 2. KANTOROVICH L. V. and K~YLOVV. I., Approximate Methods of Higher Analysis p. 240. Gostekhizdat, Moscow (1949).
Letters to the editors
On calculations for cylindrical reactors* (Received 24 December
1956)
IN calculations of neutron multiplication factors and critical sixes for reactors of the shape of spheres or infinite cylinders using the age-diffusion approximation, a useful procedure is to use the multigroup method, and to solve the group equations by the method of difference factorization developed by MARCHUK.“) In these cases the neutron flux nu depends only on a single spatial variable and the Laplacian V2 has the form
v&La
-
rG(ar
a %r ( 1
(CL= 1, 2).
If the reactor has the shape of a finite cylinder, and the flux depends on two variables, r and z, then the Laplacian has the form
and the method of difference factorization is not immediately applicable for solving the group equations. This is because this method involves the inversion of a matrix formed from the coefficients of the unknowns in the one-dimensional group equations and essential use is made of the fact that this matrix contains non-vanishing elements only in the principal diagonal and the two adjacent diagonals. For two-dimensional systems, the corresponding matrices also have elements in the next two diagonals, five diagonals in all, in which case they cannot be inverted by the method of difference factorization. MARCHUK’S method of reactor calculation consists of two distinct stages. Firstly the group constants are calculated by averaging cross-sections over lethargy groups and secondly the group equations are solved by difference factorization. In the case of two spatial variables the first stage of the calculation is unaltered, and we shall not discuss it here. The purpose of this paper is to draw attention to the possibility of using an approximate Laplacian enabling the entire procedure of the second stage of the method to be applied. We start with the following equation for the slowing down of neutrons in an intermediate reactor
--
1
3&
agz,nu V%u + C,nv+ = vS(r, au
2, u),
s co
S(r, 2, u) = F(u)
where
C,nu du
0
Here C,,, I;,, EE., F and & are experimentally known functions of the lethargy u. The required neutron flux must be a finite non-negative function of u throughout the reactor and must vanish at the extrapolated boundary. At the boundaries between physically distinct media both nv and the current density 1/(3&)V,no must be continuous, where V, is the component of the gradient operator normal to the boundary. Also, for u = 0 we must have nv = 0. For the numerical solution of(l), we form a lattice in space given by ri = ih, and .zj = jk, and for simplicity we put hr = h, = h, and also assume the same value of h in all regions of the reactor. We also refrain from discussing possible lattice points at region boundaries. Likewise, we divide the lethargy scale into steps of length 1. For brevity we introduce the following notation nu(ih, jh, kl) = nvif’ ; r,ha A;dm$ ha Arant$
=
= nu:F]_l - 2nz$’ + nojI;‘,l;
/A,fCF’no:f’ = AXE”+”nt$+” - &Z:” no::‘.
The indices i, j and k have the following ranges:
l
by N. KEMMERfrom Aromnaya Energiya 3. NO. 7, 53 (1957).
nv
is to be
101
Letters to the editors
computed. Finally, we shall attempt to make the statement of the problem clearer by repeating each equation in matrix notation. We now state two approximations to equation (1) which, in a sense, are opposed to each other. (I) The explicit scheme of equations
or more briefly nylk+ll = ,zjU~‘~v” +
co
Q'k'
(II) The “inverted” scheme of equations
or more briefly S’“+Unv’L+” = nvlL’ + Q'k' (3) Here A”’ and Pfl’ are matrices of order N, x N, each row of which has no more than five elements differing from zero: nP’ is the (N, x NJ-dimensional vector which we have to determine and Q”’ is a similar vector that is given if the sources are known. A comparison of the systems (2) and (3) shows that the former set of equations is considerably simpler than the latter because (2) requires multiplication of the matrices P’ by a vector, whereas (3) requires the inversion of the matrices B lk+l’. However, in order that the numerically greatest eigenvalues of AIL’should be less than 1, as is necessary for stability, the lethargy step I must be small in comparison with the spatial step h. On the other hand system (3) is stable for any 1 but it requires the inversion of a high order matrix for each k. Thus both Method 1 and Method 2 require a considerable number of arithmetical operations; Method 1 because the lethargy step has to be kept small, Method 2 because of the slow convergence of the interaction procedure for inverting the B(*’(non-iterative methods for solving (3) are of no use). Bearing in mind the particular, rectangular, shape of the spatial region we are concerned with, we propose that such reactor calculations be performed using another method, which is free from these disadvantages. (III) Hybrid scheme of equations
or, more briefly
Cy+l’nv, 12ktl)= ~‘1Bk’,#k’ + ~‘2k+2’nv;2k+2’ a
=
&2k+l',$k+l'
+
Q;““’ Qy+l,
(i = 1, 2, . . . N,), (j
=
1,2, . . . Nz)
Here Ciartl’ and Ciaktz’ and likewise 0:““’and. Dp’+” are respectively matrices of order N, and N, which have no rows with more than three non-vanishing elements; the unknown vectors nv!” I and nvy’ are of orders N, and N, respectively and so are the known vectors Qj”’ and Qy’. Like Scheme II, Scheme III is stable but requires the inversion of matrices, but in place of the inversion of one matrix of order (Nr x N,) one has here to invert N, matrices of order N, or N, matrices of order N,. The advantage of III over II lies in the fact that the matrices Ci” and Cl” have elements only in three diagonals. This makes it possible to invert them exactly, instead of by iteration, by the use of factorization. In the number of arithmetical operations involved, Method III is surpassed only by the method of reduction to Cauchy’s problem,‘*’ which is not stable and therefore unsuitable for an electronic computer. V. K. SAULEV REFERENCES 1. MARCHUKG. I., J. Nucl. Energy 3, 238 (1956). 2. KANTOROVICH L. V. and K~YLOVV. I., Approximate Methods of Higher Analysis p. 240. Gostekhizdat, Moscow (1949).