CONF-870601--2

DE87 004250

A NODAL INTEGRAL METHOD FOR THE NEUTRON DIFFUSION EQUATION

IN CYLINDRICAL GEOMETRY

Y. Y. Azmy
Oak Ridge National Laboratory
Oak Ridge, Tennessee 37831

DISCLAIMER

This report was prepared as an account of work sponsored by an agency of the United States
Government. Neither the United States Government nor any agency thereof, nor any of their
employees, makes any warranty, express or implied, or assumes any legal liability or responsi-
bility for the accuracy, completeness, or usefulness of any information, apparatus, product, or
process disclosed, or represents that its use would not infringe privately owned rights. Refer-
ence herein to any specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recom-
mendation, or favoring by the United States Government or any agency thereof. The views
and opinions of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

Submitted for presentation at the 1987 Annual Meeting of the American

Nuclear Society, Dallas, Texas, June 7-11, 1987.

Unsponsored research performed in part at the University of Virginia,

Charlottesville, VA, and, in part at the Oak Ridge National Laboratory,
Oak Ridge, Tennessee.

"The submitted manuscript has

authored bv > con,r,cto, of the

~ " r s ' ° * s° to us G
— ' a ] S T B , B U T ,Q}j OF THIS DOCUMENT IS UNLIMITED

O p e r a t e d by M a r t i n Marietta Enorr-y S y s t e m s , I n c . u n d e r Contract No. DE-

AC05-84OR21400 w i t h the U.S. D e p a r t m e n t of. E n e r g y .
 A NODAL INTEGRAL. METHOD FOR THE NEUTRON DIFFUSION EQUATION
IN CYLINDRICAL GEOMETRY
Y. Y. Azmy
Oak Ridge National Laboratory
Oak Ridge, Tennessee 37831

The nodal methodology is based on retaining a higher degree of

analyticity in the process of deriving the discrete-variable equations
compared to conventional numerical methods [1-3]. As a result,
extensive numerical testing of nodal methods developed for a wide
variety of partial differential equations [1-3] and comparison of the
results to conventional methods, have established the superior accuracy
of nodal methods on coarse-meshes. Moreover, these tests have shown
that nodal methods are more computationally efficient than
finite-difference and finite-element methods in the sense that they
require shorter CPU times to achieve comparable accuracy in the
solutions. However, nodal formalisms and the final discrete-variable
equations they produce are, in general, more complicated than their
conventional counterparts. This, together with anticipated difficulties;
n applying the transverse-averaging procedure in curvilinear
coordinates, has limited the application;; of nodal methods, so far, to
Cartesian geometry , and with additional approximations to hexagonal geometry.
In this Summary we report recent progress in deriving and
numerically implementing a nodal integral method (NIM) for solving the
neutron diffusion equation in cylindrical, r-z, geometry. Also, we
present comparisons of our numerical solutions to two test problems
with those obtained by the code EXTERMINATOR-2 [5], which indicate the
superior accuracy of the nodal integral method solutions on much coarser
meshes.. The derivation briefly described here is an extension of the
 -2-

general nodal integral approach [6] previously developed for deriving

nodal methods for linear partial differential equations in Cartesian
geometry. The main difference is the weight function, r, applied in
performing the transverse-integration in the radial direction. That is,
the r-transverse averaging operator is defined by
r
?

r
2 1 1
while the z-transverse averaging operator, Z, is defined exactly as in
Cartesian geometry; e.g., see Ref. 5. In Eq. (1), r, and r ? are the
inner and outer radii of a typical node, or computational cell, in the
r-z mesh covering the domain of the problem. The r- [z-] averaged, z-
[r-] dependent flux is defined, as usual by $ r (z) = R <f> (r,z) [IfJz (r) E
Z * (r,z)].

The nodal method discrete-variable equations are derived as

follows. First, we operate on the continuous-variable neutron diffusion
equation by R, then by Z. This results in two second order ordinary
differential equations in the transverse-averaged fluxes, 4 and <ji
with differential operators identical to the one-dimensional radial, and
the one-dimensional axial diffusion operators. The RHS's of these ODE's
include the transverse-averaged external source (which also contains the
in-scattering and fission sources in the commonly used source iteration
technique) as well as the "leakage" terms. The latter arise from the
r- (z-) transverse-averaging of the second order derivative with respect
to r (z), and are responsible for the only approximation in the NIM [6].
This approximation is necessary in order to express the leakage terms in
terms of the discrete-variables in the method as explained shortly.
 -3-

The transverse-integrated ODE's can be solved exactly; the solution

is written as the sum of the complementary function (which includes tv/o
arbitrary constants) and a particular integral arising from the terms on
the RHS. The complementary function for the z-averaged equation
involves modified Bessel functions of zeroth order, and for the
r-averaged equation involves hyperbolic functions. Various possible
choices for the discrete-variables have been discussed before [6]; here,
we use the arbitrary constants of the complementary functions. To
determine the particular integral, however, we need to know the
functional form of the leakage terms, which as yet are unavailable.
Hence, we expand them locally in Legendre series, truncated (in the
lowest order approximation discussed here) at the constant terms. The
expansion coefficients are determined by imposing two constraints on the
solution [L6]: the uniqueness of the nodal average of the flux, and the
balance of the neutron population over a node. These conditions are
sufficient to express the leakage in terms of the discrete-variables
[6].

The final discrete-variable equations are essentially local

boundary conditions that force the solution to satisfy the continuity of
the flux, and the continuity (or known discontinuity) of the current
across node boundaries. At external boundaries of the system, explicit
boundary conditions characteristic of the problem being solved are
assigned. The resulting linear algebraic problem can be solved exactly
or iteratively by any of the numerous available techniques.
We have coded the discrete-variable equations described above into
a computer program, and used an exact, sparse matrix solver routine to
solve them. In order to evaluate the performance of our new numerical
method, it is applied to two test problems and the results are
 -4-

compared with those of a finite-difference code, EXTERMINATOR-2. The

first test problem, a five-region infinite cylinder problem, has an
analytic solution; hence, it is used as a basic check on the derivation
and implementation of the method. The NIM yields the exact analytic
solution on various meshes, the most coarse mesh being one node, or
computational cell, per region. The region-averaged fluxes, and the
leakage on the outer surface of the cylinder are calculated by
EXTERMINATOR-2 for the same problem. The error in these quantities is
of the order of 8% for a node size approximately half the diffusion
length.
The second test problem, shown in Fig. 1, is a four-region finite-
cylinder problem. This problem is solved on 4x4, 8x8, and 12x12 NIM
meshes (6 = .25, .125, and .083, respectively) and on 14x14, 22x22, and
42x42 EXTERMINATOR-2 meshes (6 = .1, .05, and .025, respectively). Six
quantities are calculated from the solution on each mesh and are used to
determine the accuracy of the solution on this mesh with respect to values
obtained by h2-extrapolation of the 22x22 and 42x42 EXTERMINATOR-2 results.
For this problem, the six quantities are the four region-averaged fluxes,
the average radial leakage on the outer surface of the cylinder, and the
average axial leakage on the upper surface of the cylinder. Comparison of
the results of this problem are presented in Table 1. In this problem also,
the NIM achieves a reasonably good accuracy on very coarse meshes as can be
seen from the 4x4 (6 = .25) mesh. The NIM 8x8 (6 = .125) results are observed
to be more accurate than the EXTERMINATOR-2 42x42 (6 = .025)results. Further-
more, the latter are about two to five times less accurate than the much
coarser mesh, NIM 12x12 (6 = .083) results.
 / -5-

We have derived and implemented a nodal method for the numerical

solution of the neutron diffusion equation in r-z cylindrical coordinates.
When applied to two test problems it results in solutions which are far
more accurate than those obtained on fine meshes by a finite-difference
code. This demonstrates that it is possible to develop nodal methods for
partial differential equations in cylindrical geometry which have many
attractive features similar to those displayed by their Cartesian geometry
counterparts.
 -6-

REFERENCES

1. J. Doming, "Modern Coarse-Mesh Methods--A Development of the

70's," in Computational Methods of Nuclear Engineering, p. 3-1,
Williamsburg, VA (1979).
2. J. Doming, "A Review of Local Green's Function Coarse-Mesh and
Nodal Methods for the Numerical Solution of Fluid Flow Problems,"
in Advances in Math. Methods for the Solution of Nuclear
Engineering Problems, Vol. 1, p. 437, Fachinformationszentrurn
Energie, Physik, Mathematik GmbH, Karlsruhe, West Germany (1981).
3. R. D. Lawrence, "Progress in Nodal Methods for the Solution of the
Neutron Diffusion and Transport Equations," to appear in Prog.
Nucl. Energy.

4. R. D. Lawrence, "A Nodal Method for Three-Dimensional Fast

Reactor Calculations in Hexagonal Geometry," in Advances in
Reactor Computations, p. 1030, Am. Nucl. S o c , LaGrange Park,
IL (1983).
5
- T. B. Fowler, M. L. Tobias, D. R. Vondy, "EXTERMINATOR-2: A
FORTRAN IV Code for Solving Multigroup Neutron Diffusion Equations
in Two Dimensions," ORNL-4078.
5
- Y. Y. Azmy and J. J. Doming, "A Nodal Integral Approach to the
Numerical Solution of Partial Differential Equations," in Advances
in Reactor Computations, p. 893, Am. Nucl. S o c , LaGrange Park, IL
"0983).
 Table 1. Region-averaged fluxes, and radial and axial leakages for the four-region, finite-cylinder problem
(see Fig. 1) as calculated by the NIM on 4x4, 8x8, and 12x12 meshes, and by EXTERMINATOR-? on
14x14, 22x22, and 42x42 meshes using a 1.0E-8 pointwise, relative convergence criterion. The
percentage error is calculated relative to the extrapolated solution which is obtained by h -
extrapolation of the 22x22 and 42x42 meshes EXTERMINATOR-2 results.
Mesh NIM NIM NIM EXT.-2 EXT.-? EXT.-2 Extrapolated
4x4 8x8 12x12 14x14 22x22 42x42 value
cell size .25 .125 .083 -.1 -.05 -.025
(6)

.200016E-2 .199144E-2 .198842E-2 .193274E-2 .195056E-2 .197591E-2 .198436E-?

err) (.80) (.36) (.20) (-2.60) (-1.70) (-.43)

<t>9 .259792E-3 .263875E-3 .265420E-3 .292119E-3 .284623E-3 .271826E-3 .267560E-3

(% err) (-2.90) (-1.38) (-.80) (9.18) (6.38) (1.59)

*, .193276E-4 .201016E-4 .202513E-4 .240390E-4 .220640E-4 .208181E-4 .204028E-4

(% err) (-5.27) (-1.48) (-.74) (17.82) (8.14) (2.04)

.236527E-5 .240976E-5 .241952E-5 .286088E-5 .263681E-5 .248272E-5 .243136E-5

err) (-2.72) (-.89) (-.49) (17.67) (8.45) (2.11)

LR .996853E-4 J01318E-3 .101684E-3 .123044E-3 .111294E-3 .104422E-3 .101382E-3

(% err) (-2.39) (-.80) (-.44) (20.48) (8.97) (2.24)

LT .314715E-4 .329100E-4 .331387E-4 .405846E-4 .364596E-4 .341418E-4 .333692E-4

i err) (-5.69) (-1.38) (-.69) (21.62) (9.26) (2.32)
 v a c u u m
1.0

.75

V
a
.50 c
u
u
m

Z=o
///////////77///7777///////////777/
r=0 .25 .50 .75 1.0

C = T. C -C =C =f)
1 ' ? ^ 4

D D
r 2=D3=D4=1 •
a-j=400;

Fig. I. 1w -The second test problem:, the four-region,

finite-cylinder problem.