A nodal solution of the multigroup neutron transport equation using spherical harmonics

Ann. nucl. Energy, Vol. 15, No. 10/11,pp. 501-509, 1988

0306-4549/88$3.00+ 0.00 Copyright (C 1988PergamonPress pie

Printedin Great Britain.All rights reserved

A NODAL SOLUTION OF THE MULTIGROUP NEUTRON TRANSPORT EQUATION USING SPHERICAL HARMONICS F. |NANC and A. F. ROHACH Department of Nuclear Engineering, Iowa State University, 261 SweeneyHall, Ames, IA 5001l, U.S.A. (Received 2 June 1988)

Abstraet--A nodal method has been developed for solving the neutron transport equation using spherical harmonics in 2-D Cartesian coordinates for rectangular nodes. The first-order differential equations resulting from the expansion of the angular flux into spherical harmonics are reduced to second-order equations by eliminating odd-order moments. Albedo type of boundary conditions analogous to the diffusion equation albedo boundary condition are developed. The resulting system is solved by the leastsquares minimization scheme employing Legendre polynomial expansions for the spherical harmonic moments. Some results presented for P~ and P3 approximations compare to the reference solutions favorably.

1. I N T R O D U C T I O N

The nodal method based on the least-squares minimization has been applied to the neutron diffusion equation (Rohach, 1986a~c, 1987) and to the 1-D neutron transport equation (Feiz, 1986) successfully. In this project, that nodal method is developed further to solve the P, equations in 2-D geometry for the multigroup transport equation problems. The developed scheme uses rectangular nodes in Cartesian geometry. The first-order differential equations resulting from the spherical harmonic expansion of the angular flux are reduced to a set of second-order differential equations by eliminating the odd-order moments from the first-order equations. These second-order differential equations are decoupled from each other and solved by using an iterative scheme. The boundary and the interface conditions are manipulated to have the same forms as the partial currents in the diffusion equation calculations. The nodal method is built upon these equations which resemble the diffusion equation and the partial currents. In the least-squares method, the spherical harmonic moments are approximated by 2-D Legendre polynomials. The residuals which arise from substitution of these approximations into the governing equations are minimized. The conditions provided by this minimization scheme is not enough for determining all coefficients of the Legendre polynomials. The rest of the conditions are provided by the implementation of the boundary and the interface conditions in an integral sense. The method is tested by solving some fixed source and eigenvalue problems for which the P, approximation solutions are already

available. The results are in good agreement with previously published results.

2. DEVELOPMENTOF THE GOVERNINGEQUATIONS The expansion of the angular neutron flux is done by using unnormalized spherical harmonics in the following manner : ud,(r,fl) = Z ( 2 n + l ) i n-0

[~nmY(r) cOS ~'/~0

m--0

+7,,,g(r) sin m~p]P,,,(cos 0),

(1)

where 0 and ~p are axial and azimuthal angles. P,,~(cos 0) is the associated Legendre polynomial. The source term in the neutron transport equation is also expanded into unnormalized spherical harmonics as done with the angular flux while the anisotropic scattering cross section is approximated by addition formula of the Legendre polynomials. The application of a Galerkin scheme to the neutron transport equation results in a set of first-order partial differential equations. These first-order partial differential equations allow odd moments (n = odd) to be expressed in terms of the even moments and their derivatives. If these expressions for the odd moments are substituted throughout all equations, the resulting system is obtained as a set of second-order partial differential equations. This second-order form of the P3 approximation to the one-group neutron transport is given below for an isotropic source and a homogeneous 501

502

F. |NANCand A. F. ROHACH

region case where (2,, = Y.,- 2~') :

1

1

2

2( ~2

3z, V 2 * ° ° - ~ V *2o+~ &~ 4 02722

+ Z, Ox@ 1

2

20*00 =

/ 1 8 2 ' +7Y3"\

02722

2

(18Y, +4223`]

48Z +42Z~'~ 2 - ~l-~ i ~7 ) V *22-- 1022*22

02)

~y2 "22

,

-Soo(x,y),

.(2)

x ~2

(3621 +8423`] 02

(2421+2123~ 2 " 2 ~ 1 ~ - ) V 722--5722 --

02

*20, (8)

1 02*00 3Y~l (.3X~y

(3Z, +7Y3"] 02*20

+ \ ~K

-1-522*20 = O, (3)

1 02 32, (0-~x22 ~ 2 ) * 0 0 - ( 3 2 ' - + 723'] ( 02 \ 212,Z3 /]~kaX 2 (48E1 +422,']

+ \ ~i~

1 02.oo

02 072"/*20

-/v2.~2-1°~2.22 =0,

(4)

(242,

-1- \

~

3. BOUNDARY AND INTERFACE CONDITIONS

} OxOy -1 V21222--522)]22 = O. (5)

Although the coupled differential equations given above can be solved simultaneously, an iterative scheme is preferred for computer storage concerns. In that scheme, the first equation can be solved for *00 and then the second equation for *20. The third and the fourth equations can be solved for *22 and Y=, respectively. Other moments in these equations are replaced by the estimates from the previous iterations. This leads to a scheme analogous to a GaussSeidel iterative scheme. For such an iterative scheme, the partial differential equations can be rewritten in the following form :

1

l V2,00__~..0,00(X, 1') = -- $O0(X,y)- ~ ~ 1 V2*20 32, 4 ~2722

Z, OxOy 18Zj q-7Z3'~

2 (~ £, ~ x 2

~2) ~

*22,

(6)

.n< 0 n~P"m (c°s 0)cos mcpUd(S,1~) df~ = f(S),

(10)

n.QP,,, (cos 0) sin m~pUd(S,[~) dD = f(S),

(11 )

fl<0

n is even and n < N, 0 ~< m ~< n, S is the surface of the domain.

,

\ ~iz~G I bT~y+ 3~, v *oo

(182' +42Y~3`]( 02 02)

+\ 21~

All the boundary and the interface conditions used in this project are manipulated such that they all have the same form similar to the partial currents used in the neutron diffusion equation. The advantage of using boundary conditions analogous to the partial currents is that a spectrum of boundary conditions varying from vacuum conditions to the reflective boundary conditions can be easily expressed with the same expression. This also simplifies implementation of the numerical method to be used. The boundary conditions used in this work are based upon the variational principle developed by Davis (1966). Although this principle is given only for the vacuum boundary conditions, it will be generalized such that conditions other than the vacuum conditions can also be represented through this principle. The generalized form of that principle is given by the two following equations :

where

2

~ - 2 ~ ) v .20- 522*20 _ (36 ,+842,

(9)

If the source term given in equation (6) contains external sources, the problem becomes a fixed source problem and if it has only the fission source, the overall problem becomes an eigenvalue problem.

(3E, + 7 2 , ' ] c?2*2o

32, OxOy \ N 2 ~

I ~-y

] Ox2 ~y2 *22,

(7)

When the angular integrations in the equations (10) and (11) are carried out for the proper angular intervals, the resulting expressions contain both the even and the odd moments. At that point, the odd moments are replaced by the odd moment expressions in terms

A nodal solution of the multigroup neutron transport equation using spherical harmonics of the even moments and their derivatives. After that substitution, the boundary conditions consist of even moments and their derivatives in a coupled manner. The further simplification of the boundary conditions is done by dropping all moments and derivatives which do not have the same indices as the weight function used in the angular weighting process for that specific boundary condition. The general forms resulting from that process are given by equations (12-15) in natural coordinates:

[O~Onm~0.,.7 + f l ~I u u

Ju.._ * =f~(w),

(12)

e@,,m+ fl~ ~ww Ju, w=, = f2(u),

(13)

1

[~O.,.-3~ aO.~] 0u Jo,.= , =f3(w), 1 a~.,.]

~O.m--~U;rw J~,w= , =f0(~),

Table 1. The coefficientsof the boundaryconditions Weight function

•

fl

1 3 cos20 1 sin 2 0 cos 2o

1 250 ~ 50

2/3~ (128/3E 0 + (768/7E~) (384/3Zj) + ( 1024/7Z0 (24/3E j)+ (64/7E3)

,,in:Oqn~,:~

2~

for the adjacent node in determining the right side. These functions have neutronic parameters of the adjacent nodes and the spatial coordinates are different from the left sides. The general form of these equations are given below. The constants and the indices in these equations are the same as the ones for the left sides of the boundary conditions :

(14)

J'l(W) = [O~,,,,, q'-fl l~~630n"7 -j/,j+

(15)

fz(u)

where 2~ is width and 2t/is height of the node ij shown in Fig. 1. The constants in the equations (12-15) are given in the Table 1 for all weight functions. If the fls in this table are studied, it is seen that these constants are exact multiples of the constants multiplying the Laplace operator in the corresponding differential equations. The right-hand-sides of the boundary conditions have two different forms for the interface and the boundary conditions. For the interface case, the same expressions as the equations (10) and (11) are used

503

/-

/

,.. ~ , ,

F 1 O~.~.]

(17) ....

F

1 a%.7

r

1 aO.,~7

(16)

A(w)= L~*nm--3~-U--J,j ,...,

,.

(18) (19)

The right sides for the external boundaries are developed by setting a proportionality between the inbound and the outbound weighted angular fluxes. The proportionality factor is the albedo 0c). The mathematical expressions for these albedo boundary conditions are given by the following two expressions : fnn< 0 n~'~P"m(c°s 0) cos rmp~(S. ~~)d~

i-ld-1

i-ld

i-l,j~-I

= ~cf

jm > 0

nD.Pnm(cosO)cosm~o~(S,~)d_Q, (20)

f.n< 0 n.QP.,. (cos 0) sin moW(S, ~) d~ i,j-1

i,j

i,j+t

= ~cf

n.QP,~(cosO)sinmq)W(S,~)d~, (21) >0

i+ld-1

i+ld

i+ld+l

Fig, 1. The nodal geometry. ANE 15:10/ll-C

If these conditions are processed as described in the beginning of this section, the right sides for the external boundary conditions are obtained in the following forms : ~f~~w~

KR [

7100.m]

(22)

504

F. INANCand A. F. ROHACH

F

f3(w) =

1

Ol~)nmG "

I-

1

....

,'

(24)

Oi~)nm-]

f4(u)=',.L~<*o.+n~Wl,.,=,

The coefficients el.is of the first polynomial in the residual function are actually functions of the coefficients of the original Legendre approximation to the spherical moments. These new coefficients are given as below :

Coo =

1

1

~2 (3a20 + 10a40) + ~ 1

All the as and fls in these equations are the same as in the Table 1 except the ones in the fourth equation where the weight function is (sin20 sin2go). The constants for that type are : c¢= - 25, +

64)

(3ao2+

10a04), (29)

i

Col = if: 3a21 + r/2 15a03, 1

(30)

1

C02 = ~2 3a22 q- r/2 35a04, 1

1

1

1

(31)

c,0 = ~ 1 5 a 3 0 + ~ 3 a , 2 ,

.

(32)

ell = ~ 15a31 q-- ~ 15al3,

C2o =

4. FORMALISM FOR NODAL METHOD

The nodal method developed in this project is based upon the least squares minimization and it uses fourthorder Legendre polynomials for approximating the spherical harmonic moments in the spatial domain. Since all the equations in the second section have been manipulated into the same form, the nodal formalism will be developed for the following generalized form in the natural coordinates : l- 1 02 A • +

1 02 -'1 w) + o r ( u , w) = g(u, w),

(26) where u = x/~ and w = y/rl. The right-hand-side of the differential equation will include both the source terms and the other spherical moments as given by equations (6-9). All spherical harmonic moments are expanded into fourth-order 2-D Legendre polynomials in the following manner. The advantage of using Legendre polynomials lies with its orthogonality property which results in simple final expressions for the coefficients : 4

1

~3a22.

(34)

The conventional application of the least-squares minimization is done by using the derivative of the residues with respect to the coefficients of the polynomials as the weight functions as shown below :

~ai,j

rea

3 Jarea

Oai,j

where i = 0 . . . . , 4 and j = 0 . . . . . 4 - i . This approach has difficulties in introducing the boundary conditions into the scheme. The application of the least-squares minimization in this project deviates from the conventional approach in the choice of the weight functions. In contrast to the common implementation, the derivatives of the residues are taken with respect to the coefficients of the first polynomial representing the Laplace operator in the original differential equations. ,

, ~

R(u, w) du dw = 0,

(36)

4 i

(l)(u, w) = ~ ~ a,jPi(u)Pj(w).

(27)

i=Oj-O

When the polynomial approximation for the moment is substituted into equation (26), the equality can not be satisfied and a residual function is obtained in the following form : 2

1

~735a40 +

(33)

2-i

wherei=0, 1,2andj=0,...,2-£ When these integrations are carried out, the following six conditions are obtained :

A[~(3azo+lOa,o)+~(3ao2+lOao4)] 1

+Baoo = ~Goo,

R(u,w) = A ~ Z cijPi(u)P/(w)

(37)

i-0/=0 4

4 i

+B ~ ~ a,jPi(u)Pj(w)--g(u,w). i-Oj=O

(28)

AI~23a21+~15ao3]+Baol

3 = ~Gol,

(38)

A nodal solution of the multigroup neutron transport equation using spherical harmonics

505

otaoo-- (~ + fl)alo + (a + 3fl)a2o-- (~ + 6fl)a3o A A

[;

3a22+

35ao4

1] 15a~o+r/~3a,~

+Bao~ = ~Go2,

3 +Ba,o=jG~o,

9 ,,, A I 1~ 1 5 a 3 , + ~ - 115 a , 3 ] + B a I , = ~ G

A[~35a4o+l~3a:z]

(39) + ( ~ + 10fl)a4o + ~[~ao, - ( ~ + f l ) a , , (40)

+ (~ + 3fl)a2, -- (~ + 6fl)a3,] - ~[~ao3 - (a + fl)a, 3]

+Ba2o= :G2o,

(42)

;o

--

(41)

f3(w) dw,

(47)

c~aoo- (a + fl)a,o + (a + 3fl)a2o - (c~+ 6fl)a3o + (~+ 10fl)a4o - l[eao, - (a+ fl)a,,

where :

j.o

+ (~ + 3 fl)a2, - (~ q- 6fl)a3, ] q- ~[o~a03- (~ -k-fl)al 3] Gi/ =

[g(u, w)P i (u)P~ (w)] du dw. I

I

=

The other conditions to be used in determining the coefficients of the polynomials are obtained from the implementation of the boundary conditions. The implementation of the boundary conditions are done in an integral sense. The boundary conditions given in the third section are integrated at the interfaces from 1 to 0 and then from 0 to - 1. This scheme gives two conditions for each nodal interface which total up to eight conditions. The expressions obtained from this scheme are listed in the equations (43-50). The as and the fls in these equations are given in the Table 1.

aaoo + (c~+ fl)at o + (~ + 3fl)a20 + (~ + 6fl)a30

f3(w)dw,

C~aoo- (e + fl)ao, + (e + 3fl)ao2 - (~ + 6fl)ao3 + (~ + 10fl)a04 + {[~a,0-- (a + fl)a,,

+ (~ + 3fl)a~2 - (ct + 6fl)a, 3] - ~ [aa 3o - (a + fl)a3,] =

(49)

j,o

+ (~t + 3fl)a12--(et +6fl)a13]+ ~[ota3o--(~ + fl)a3~]

j"

fT(w) dw,

(43)

0

:moo + (c~+ fl)a~ o + (a + 3fl)a2o + (c~+ 6fl)a3o +(c~+ 10fl)a4o- ~[~aol + ( ~ + f l ) a l l

=

fa(u)du,

+ (~ + 10fl)ao4 -- ½[aa,o -- (~ + fl)a,,

=

+ (~ + 3fl)a:~ + (~ + 6fl)a3 ~]-- 8~[~ao3 + (~ + fl)a~ 3]

+ (a + 3fl)a2, + (~ + 6fl)a3 i1 + ~[~ao~ + (~ +

fo

eaoo - (e + fl)ao, + (~ + 3fl)ao 2 - (~ + 6fl)ao3

+ ( ~ + 10fl)a4o + ½[~ao, + (~-I-fl)a, 1

=

(48)

I

f4(u)du.

(50)

I

As it can be seen, the minimization and the boundary conditions total up to fourteen which is one short for determining all fifteen coefficients. This problem is bypassed by setting the coefficient az2 equal to zero. Reason for choosing this specific coefficient is that it is the least used coefficient in the fourteen conditions.

B)a, 3]

f~(w) dw.

(44)

5. RESULTS

I

aaoo + (~ + fi)aoT + (~ + 3fl)a0~ + (~ + 6fl)a03 + (~-b 10fl)a04 + ½[~a,0 + (~+ fl)a,,

+ (~ + 3fl)a,~ + (~t + 6fl)a, 3 ] - ~[~a3o + (a + fl)a3,] =

f~(u) du,

(45)

aaoo + (a + fi)ao, + (~ + 3fl)a0~ + (~ + 6fl)a03 + ( e + 10fl)a04- ½[aa~o + (a+fl)a,,

j.o

+ (a + 3fl)at2 + (~ + 6fl)a, 3] + ~ [aa30 + (a + fl)a~] =

f2(u)du,

1

(46)

A computer program has been developed which implements the nodal equations described in Section 4. The method has been tested by solving four problems for which either P, approximation solutions or exact solutions are already available. The first problem solved has been developed by Fletcher (1977, 1983, 1986) and the problem has an exact solution. The second problem has been presented by Nastelson (1971) and the solution to that problem is taken from Kobayashi et al. (1986). The third problem is a twogroup source problem and the benchmark solution for that problem has been provided by Gelbard and Crawford (1977), The fourth problem has been presented by Kobayashi et al. (1986) and it is a fourgroup criticality problem for a fast reactor.

506

F. INANC a n d A. F. ROHACH

The results for the both approximations are tabulated in Table 2 and the results are also presented in Fig. 3. The P~ approximation underestimates the exact solution in a significant amount while the P3 approximation enhances accuracy considerably over the Pj approximation. The comparison of both P, approximation solutions to the solutions given by Fletcher show that the nodal method results overshoot the Fletcher's results near the reflective boundary while they underestimate the results near the vacuum boundary. The second problem is also a one-group source problem. As it can be seen from Fig. 4, the Nastelson's

Y 4.0

1.2

Source It

0.0

1.2

4.0

5.0

X

Fig. 2. The g e o m e t r y for Fletcher's problem.

Fletcher's problem is a one group source problem for which the geometry is shown in Fig. 2. As it can be seen from this figure, the problem is a square pure absorber region with sidelengths of 4.0 cm. It has a constant source in a 1.2 cm 2 region adjacent to the lower left corner. The neutronic parameters of the pure absorber are E, = 1 cm-1 and Zs = 0 cm J and a total unit source strength. The boundary conditions for the sides adjacent to the source are reflective, while the boundary conditions for other sides are vacuum conditions. This problem has been solved by partitioning the domain into a 6 × 6 node configuration.

~X-- PIapprox.

2.5

2.0

~x~ x ~n~ "X~x~

~

P3 approx.

~

Exact solution

I

"x

14.11.5

" X~'xx X ~ Xx x x

u.

1.0

%X~x, ~ X~X~.x.. "

0.,5

J

0

1

I

I

2 Distance

5

4

(cm)

Fig. 3. The total fluxes at y = 3.9 c m for Fletcher's problem.

Table 2. The total fluxes at y = 3.9 cm for Fletcher's problem Distance

-I

(cm)

P- I

P- 1 (Fletcher)

P-3

(Fletcher)

P-3

Exact

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5 1.7 1.9 2.1 2.3 2.5 2.7 2.9 3.l 3.3 3.5 3.7 3.9

2.134E - 03 2.095 2.025 1.922 1.790 1.639 1.475 1.305 1.136 9 . 7 3 4 E - 04 8.244 6.877 5.636 4.584 3.694 2.930 2.286 1.792 1.381 1.040

1.947E - 03 1.914 1.851 1.759 1.644 1.511 1.366 1.215 1.063 9 . 1 7 4 E - 04 7.797 6.537 5.412 4.428 3.583 2.871 2.279 1.794 1.401 1.087

2.589E - 03 2.547 2.477 2.363 2.226 2.071 1.886 1.708 1.531 1.351 I. 187 1.036 8.909E-- 04 7.660 6.553 5.520 4.628 3.871 3.147 2.460

2.532E - 03 2.500 2.438 2.349 2.234 2.100 1.950 1.791 1.627 1.464 1.305 1.153 1.010 8.787E - 0 4 7.588 6.510 5.552 4.709 3.974 3.338

2.603E - 03 2.596 2.502 2.408 2.286 2.143 1.985 1.818 1.646 1.475 1.307 1.149 1.000 8.652E-04 7.426 6,335 5.369 4.528 3.802 3.179

A n o d a l s o l u t i o n o f the m u l t i g r o u p n e u t r o n t r a n s p o r t e q u a t i o n u s i n g s p h e r i c a l h a r m o n i c s Y 3.0

1.0 Source

0.0

1.0

3.0

X

507

Although the results from Kobayashi are reliable up to only one significant digit due to the digitization from a plot, it can still be said that the differences between the nodal method results and the published results usually lie within a range of 2%. The third problem is a two-group source problem for which the geometry is shown in Fig. 6. This problem has been solved with reflective boundary conditions all around the domain. The domain has been divided into a 14x 13 array. The solutions for the problem are plotted in Fig. 7 and P 3 results are also compared to the exact solution and the solution given by Kobayashi in Table 4. As it can be seen from Fig. 7, the P, results are very poor. Similar behaviour is also reported by Kobayashi. As shown by the comparison in Table 4, the nodal method results are in good agreement with Kobayashi's results. As in the

Fig. 4. T h e g e o m e t r y f o r N a s t e l s o n ' s p r o b l e m .

problem resembles the first problem but this problem has two different material regions. In addition to that, all boundaries are reflective boundaries. The neutronic parameters of the source region are different from the outer region neutronic parameters. The E, in the source region is 1 c m - ' and Es is 0.25 cm '. The scattering cross section in the outer region is Es = 0.5 cm- ~while the total cross section is the same as that in the source region. The source density to be used is S = 1. This problem also has been partitioned into a 6 x 6 node configuration. The P, approximation solutions published by Kobayashi et al. (1986) have been used as the reference solutions for this problem. Since the data were in graphical form, the graphics were digitized for obtaining the numerical results in Table 3. As it can be seen from Table 3, the results obtained from the nodal method agree with Kobayashi's result in a very close manner.

4.0

X~x

--x--P 1 approx. I'O--" P3 approx'

~X

i,

\x

1.5

~.0

~'X~'x~x

0.5

I 1

I 2 Distance (cm)

Fig. 5. T h e t o t a l fluxes at y = 3.0 c m f o r N a s t e l s o n ' s p r o b l e m .

Table 3. The total fluxes at y = 3.0 cm for Nastelson's problem Distance (cm)

P- I

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00

4.06E-02 3.99 3.82 3.54 3.19 2.80 2.41 2.04 1.71 1.45 1.26 1.14 I. 10

P- 1 (Kobayashi) 4.05E-02 4.00 3.83 3.54 3.24 2.84 2.43 2.03 1.74 1.46 1.27 1.16 1.10

i 3

P-3 3.65E-02 3.59 3.48 3.27 3.02 2.74 2.45 2.19 1.94 1.75 1.60 1.51 1.48

P-3 (Kobayashi) 3.66E-03 3.60 3.50 3.30 3.04 2.75 2.46 2.19 1.94 1.77 1.62 1.54 1.49

508

F. |NANC a n d A. F. ROHACH

Table 4. The first group total fluxes at y = 139.0 cm for Gelbard's problem Distance 140.O

60,0

Source

P 0.0

133.0

65,0

X

Fig. 6. The geometry for Gelbard's problem.

second problem, the nodal method results lie within 2% of Kobayashi's results. The fourth problem, which is shown in Fig. 8, is a four group eigenvalue problem and the multiplication factors and the total fluxes given by Kobayashi have been used as the reference solutions. The fast reactor consists of a core and a blanket surrounding the core. The outer boundaries of this problem are vacuum boundaries while the boundaries adjacent to the core are reflective boundaries. The cross sections for the 233U-Th core and the Th blanket are given by Kobayashi and Nishihara (1967). This problem

P-3

0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0 50.0 55.0 60.0 65.0

5.8474E 5.8345 5.8003 5.7385 5.648l 5.5388 5.3829 5.1860 4.9331 4.6416 4.2946 3.9219 3.5196 3.1020

(Kobayashi) 06

5.8843E-06 5.87t6 6.0534 5.7657 5.6674 5.5341 5.3618 5.1472 4.8884 4.5857 4.2421 3.8641 3.4610 3.0446

Exact 6.1185E-06 6.1046 6.0624 5.9901 5.8848 5.7432 5.5618 5.3379 5.0701 4.7591 4.4083 4.0242 3.6158 3.1947

has been divided into a 10x 10 configuration and 10 -5 convergence criteria has been used for the outer iteration of the source iteration. The results obtained from the nodal method are compared to Kobayashi's results in Table 5. The difference between the kc~ for P~ approximations is 0.08% while that difference is 0.05% for the P3 approximation. The total fluxes also compare favorably. The maximum difference between the total fluxes provided by the nodal method and Kobayashi's is 1.6%. Although the nodal method seems to be able to handle P, equations, as can be seen from the four problems which have been solved, it still requires examination of the affects of the node sizes on the results and the further work is underway for increasing the execution speed and reducing the dynamic memory requirements.

Y

7E-006 I l 6E-0061 ~

/

100.(

--X~ PI approx. .--o-- P3 approx. "~:~_.--E)--- Exact solution

r Blanket

4E-006/ ~

Blanket Region III

I

\

2E-006 ~

i

_

3E-oo6

50.0

m

1E - O06x-x-x-x-x-,x-x-x-x-x.--Wx-r~-x..~~ ~~ .. ±- .

OE+O00~)

P-3

(cm)

. . . . . . . e 5I. . . . . . .

5 0"i~ ' ' ' y ~ x - x ' r ~7l ,5

Core Re,on I

'-"t~:l,~

Blanket Region II

" l O d " ~ -"". X

Distance (cm) Fig. 7. The first group total fluxes at y = 139.0 cm for

Gelbard's problem.

0.0

50.(

100.0

Fig. 8. The geometry for Kobayashi's problem.

x

A n o d a l s o l u t i o n of the m u l t i g r o u p n e u t r o n t r a n s p o r t e q u a t i o n using spherical h a r m o n i c s

509

Table 5. The results for Kobayashi's problem

P-I Multiplication factor

P- 1 (Kobayashi)

! .04536

1.04449

P-3 1.05002

P-3 (Kobayashi) 1.0495 l

Integrated flux Group l

Region l

1.000

1.000

1.000

1.000

II lll

5.424E- 02 3.998E-03

5.512E - 02 4.040E - 0 3

5.047E - 02 4,241 E - 0 3

5.141 E - 02 4.263E-03

2

I II lII

4.038 6.224E - 01 1.153E-01

4.046 6.234E- 01 1.149E-01

4,000 6.022E - 01 1.140E-01

4.031 6.034E- 01 1.138E-01

3

1 II I11

3.269 7.130E-01 1.666E - 01

3,275 7.120E-01 1.656E - 01

3.247 6.915E-01 1.642E - 01

3.258 6.957E-01 1,632E - 01

4

1 lI Ill

1.087 3.207E-01 8.726E-01

1.084 3.186E-01 8.681E-02

1.073 3.141E-01 8.647E-02

1.077 3.128E-01 8.553E-02

REFERENCES

D a v i s J. (1966) Nucl. Sci. Engng 25, 189. Feiz M. (1986) U n p u b l i s h e d Ph.D. Dissertation, I o w a State University. Fletcher J. K. (1977) Ann. nucl. Energy 4, 401. Fletcher J. K. (1983) Nucl. Sci. Engng 84, 33. Fletcher J. K. (1986) Transport Theory Statist. Phys 15, 157. G e l b a r d A. M. a n d C r a w f o r d B. (1977) A N L - 7 4 1 6 (Suppl. 2).

K o b a y a s h i K., O i g a w a H. and Y a m a g a t a H. (1986) Ann. nucl. Energy 13, 663. K o b a y a s h i K. and N i s h i h a r a H. (1967) Nucl. Sci. Engng 28, 93. N a s t e l s o n M. (1971) Nucl. Sci. Engng 43, 131. R o h a c h A. (1986a) Ann. nucl. Energy 13, 57. R o h a c h A. (1986b) Trans. Am. Nucl. Soc. 52, 426. R o h a c h A. (1986c) Ann. nacl. Energy 13, 549. R o h a c h A. (1987) Ann. nucl. Energy 14, 653.