A finite element method for multigroup diffusion-transport problems in two dimensions

Annals of Nuclcw Ener,qy, Printed in Great Bntnm

Vol. 8, pp. 689 to 698, 1981

0306-6549/81/i 10689.10%02.00/0 Pergamon Press Lki

A FINITE ELEMENT METHOD FOR MULTIGROUP DIFFUSION-TRANSPORT PROBLEMS IN TWO DIMENSIONS A. K. ZIVERand A. J. H.

GODDARD

Imperial College of Science and Technology, Nuclear Power Section, Department of Mechanical Engineering, Exhibition Road, London SW7 2BX, U.K. Abstract-The method uses linear triangular and rectangular elements for the spatial dependence of the angular flux and spherical harmonic expansions for the angular dependence. Multigroup solutions are based onamaximumprinciplefor theone-groupsecond-ordereven-paritytransportequation.Anovelfeatureofthe method is an enumeration scheme for the spherical harmonic moments at the nodes of the finite element network which minimizes the bandwidth of the stiffness matrix. For coarse-mesh solutions the finite element method employs a direct elimination scheme and to proceed to fine-mesh solutions an iterative scheme is employed. The method is illustrated with the aid of four studies ; a group of ray-effect problems, a streaming problem, a PWR cell problem and a bulk-shield problem.

1. INTRODUCTION

The finite element method is a versatile numerical technique used in many branches of engineering and

sciences. Applications of the method to radiation transport is complicated by the fact that, in phase-space, three positional and two directional coordinates have to be considered. The energy dependence of the cross sections and the anisotropy of the scattering can cause further complications and can increase computational effort significantly. To date, trials of finite element methods for neutron transport have been concerned with elucidating particular features of the general problem. In this work one-dimensional formulations of Ackroyd et al. (1980) have been extended to treat problems with two space dimensions. The even-parity second-order form of the transport equation has been considered in conjunction with an orthogonal set of spherical harmonic expansions. Linear Lagrangian elements have been used for the spatial representation. Algorithms have been implemented in the twodimensional, group-dependent X-Y geometry diffusion-transport code FEED2 (Ziver, 1981). Results have been compared with published results obtained from discrete ordinate finite difference codes TWOTRAN (Lathrop and Brinkley, 1970) and DOT (Rhoades and Mynatt, 1973)and with published results obtained using the discrete-angle finite element codes DAFE (Lillie and Robinson, 1976) and FENT (Briggs et&., 1975). Ray-effect benchmark problems, in one and

two energy groups, and a one-group streaming problem have been considered. These are supplemented by a PWR cell problem and a two-group benchmark bulk-shield problem. Comparisons have been made to show the coarse-mesh superiority of the present finite element results over finite difference and discrete angle finite element predictions. In the next section a brief introduction to the formulation is given and this is followed by the results for benchmark problems. 2. DEVELOPMENT

OF THE NUMERICAL

MODEL

2.1. Even-Parity Transport and Dtjiision Functionals

In this work, the finite element method is applied to a functional giving the second-order even-parity transport equation as an Euler equation. A Rayleigh-Ritz approach is thus used to construct the discretized finite element model. The quadratic functional has the following form in the case of isotropic scattering and isotropic neutron sources in x-y geometry. Transport Functional

+2sinwcoso

84

a4

( >(> >I> ax

ay

+sin* w 3 ’ ( 8Y 689

- [(a, @ > -

%,b#m

690

A. K. ZIVW and A. J. H. GODDARD

2.3. Trial Function for Spatial Representation

+2(&Y+) - w)I$~) dS. In the case of the diffusion theory approximation in rectangular geometry the functional reduces to the following form :

Linear triangular and rectangular elements have been used to represent the spatial dependence of the even-parity angular flux. Shape functions are taken from Zienkiewin (1977) and have the form : G(x, Y) = NLfi+ Njfi + NJ,

+ NJ;,

(5)

where,

D@iision Functional

Ni=(l-;)(l-$Nm=g

for rectangular shaped elements. +2

S+(x, Y)&,Y) dx dy. IS

(2)

R

2.2. Trial Functions for Angular Representation

In this work spherical harmonic expansions are used for the angular dependence of the even-parity flux. For one- and two-moment approximations the trial functions take the following form : One Moment

2.4. Construction of the Finite Element Equations (The Global Stijiiss Matrix) The trial functions, equations (4) and (5), are used to construct the set oflinear equations for both eigenvalue and fixed-source problems. An algorithm has been developed to construct the global stiffness matrix for diffusion and transport problems. This procedure can be summarized as follows : (i) Store coordinates ofelements using the numbering scheme in Figs la and lb which gives the minimum half bandwidth. (ii) Create the material library associated with each index of the network given in Fig. la. (iii) Calculate the coefficients ofthe matrix associated with a particular node, as required, by using an algorithm.

+

& J-

sin 2oPt$z

(3)

L Moment

Each moment in the series expansion, in equation (4), gives an additional unknown and a better representation of the even-parity angular flux. However, the increase in the number of unknowns emphasizes the substantial storage requirement. This problem has been considered by Ziver et al. (1981) where core storage has been reduced to a certain extent by the application of an iterative method to solve the finite element equations. We believe that more work should be done to reduce the core-storage requirements arising from high-order angular approximations in the evenparity transport equation.

It is found that the algorithm has two parts, one for boundary elements and the other for mid-elements. As a result, the element equations are formed as required in the iterative procedure. The nodes of the grid of Fig. la, with (E - 1) x (H- 1) meshes are given a single coordinate or index by numbering the nodes on a line from left to right, starting with the bottom line and moving upwards one line at a time. Thus the index for the node with coordinates (3,2) is 2E + 3. I

I

I

I

I

I

HE

3E

2E

I

2

3

4

5

6

E

Fig. la. Global numbering of elements in two dimensions.

FE method for multigroup diffusion-transport problems

I

hlfN+P

1

1

2N+I

3Ntl

4

N+I

&+IlEtll+,:!

I I

1

1

1

(E-IlN+I

I WE-IIN+2

691

been found that the nodes on the left and right boundary, on the finite element mesh shown in Fig. 1b, give two distinct submatrices from which two recurrence formulae are obtained. Likewise, another relationship is found for the nodes that lie between the two boundaries. The derivation of the formulas used in the iterative approach are as follows : The functional K(4) is minimized for the elements shown in the figures for the relevant range of m, by moving from left to the right. Ke, Ke’, IC”etc. denote the contributions to K from elements e, e’, e” etc. for left, right and mid-indices.

(a) Nodes on the left boundary

2

N+2

2N+2

3#+2

8K”’ + 3K’“” __ PC af, af,

LFIw+2

o

-_B

1
(6)

I

I

W

e

m (Mo~~ENTNI

N

2N

3N

4N

EN

Fig. lb. The numbering of moments in the network for one-, two- and N-moment angular approximations.

-_-

e

-_-

(b) Nodes in between two boundaries + aK(“? aK(““) = o -aKce) + __ ~ aK”” 8fk af, af, + ah (N-t 1) < m < N(HE-

The coefficients of the moments in the expansion of the angular flux are numbered m = 1,2,. . . , N starting with the coefficient of the moment which is independent of direction, as shown in Fig. lb. The following scheme in Fig. lb is used to give an index to each moment at each node. In giving an index at each node for the first moment the grid is scanned in the same way as that used for assigning an index to a node but instead of increasing the index by one at each step along a scanned line it is increased by N. The first moment at the node with index 1 is given the index 1. Thus the index associated with the first moment for the node with coordinates (3,2) of node index 2E + 3 is N(2E + 2) + 2. The indices for the second moment of the nodes are assigned in a similar way. Thus for the node 2E + 3 the index of the second moment is N(2E + 2) + 2. In this way an index can be assigned to each moment at each node. The moments arefJm = 1,. . . HEN). 2.5. Algorithms for Element Equations The element equations are derived from the submatrices in the main global system matrix. It has

--

-1)

(7)


(8)

I ,

I

H u,

e

e"

m

e

--

e’ -I

I

(c) Nodes on the right boundary aK’“1 + aK”” af, af, =* [(E-l)N+l]

--n e”

_-

--

III m

et

Equations (6), (7) and (8) can be further put into the matrixformThe halfbandwidthofthesystemshownin

A. K. ZIVW and A. J. H.

692

Fig. lb can be calculated as : S,,, = [N((H-- l)E++((H--2)EN+

I)]+ 1

or

(9) S1,2 = (2fE)N.

Note that in this calculation N <
remained. Figure 4 shows finite element one- and twomoment solutions obtained in the present work; no oscillations attributable to the ray effect have been observed. Comparison with the Natelson D 16 solution shows how satisfactory is the few moment finite element solution when compared with a high-order discrete ordinate finite difference solution. Table 1. Cross-section and source data for Natelson’s onegroup ray-effect problem (Case A), illustrated in Fig. 1 Cross sections (cm- ‘)

Region 1

Region 2

z,

1.0 0.5 1.0

1.0 0.25 0.0

RESULTS

3.1. One-Group Ray-E$ect Problems

GODDARD

z, Source density (cm-“)

Case A. Natelson’s problem This problem has been studied in order to compare FEED2 finite element solutions with discrete ordinate results of Natelson. Figure 2 shows the problem geometry ; a one mean-free-path square source in a three mean-free-path shield. The source is uniform in space and isotropic in direction and source and cross sections are given in Table 1. Natelson (1971) applied the variational method to discrete ordinate equations with the aim of mitigating the ray effect. Results of discrete ordinate approximations, denoted by S,, VS,, D,, together with Natelson’s modified discrete ordinate solutions MD,, are shown in Fig. 3. It can be seen that flux shapes plotted along the x-axis (y = 0) showed anomalies for both discrete ordinate and modified discrete ordinate approximations. However, application of the variational method to discrete ordinate equations shown as MD,, suppressed the major oscillations, but minor oscillations still

0.045 ,

I

/

1

0.030 1

0.025

a005

VS

I s4 0

I

3

2

x in mean- free- paths

Fig. 3. Case A one group ray-effect problem. Flux along y = 0 axis of Fig. 2 as calculated in Natelson’s quadrature test.

Y

3.0

Region

I

2.0

Region 2

J

I 00

kfect

3.0

-X (m.f

1.0

0.0

p.l

reflector boundary conditions all sides

Fig. 2. Geometry of one-group ray-effect problem of Natelson (Case A).

x

2.0

3.0

in mean -free - paths

Fig. 4. Case A one-group ray-effect problem. Flux along y = 0 axis from FEED 2 one- and two-moment calculations compared to Natelson D,, calculations.

FE method for multigroup diffusion-transport problems Table 2. Calculated absorption rates in the source region 1 of the Case A one-group ray-effectproblem

D,, (Natelson) FEED2 one-moment 15 x 15 mesh FEED2 two-moments FEED2 two-moments VS,, (Natelson) D, (Natelson)

Table 3. Cross-section and source data for the second onegroup ray-ekct problem (Case B), illustrated in Fig. 5

Absorption rates

Method

6 x 6 mesh 15 x 15 mesh

0.4660 0.3905 0.4192 0.4464 0.4699 0.4719

In Table 2 absorption rates in the souxe region for 6 x 6 and 15 x 15 finite element mesh calculations are compared with Natelson’s discrete ordinate D16 finemesh result. With a two-moment 15 x 15 calculation the agreement is satisfactory hearing in mind the limitations of the DIG calculation. Case B. Lathrop-Lewis-Lillie-Reed

Problem

This second ‘ray effect’problem has been analysed by Lathrop (1970), Lewis (Briggs et al., 1975), Lillie (Lillie and Robinson, 1976) and Reed (Miller and Reed, 1977). Geometry and the one-group constants are given in Fig. 5 and Table 3. Lathrop (1970) considered the first-order transport equation with regard to ray effects in two-dimensional discrete ordinate equations. He concluded that the partially effective remedy of using specially selected quadrature sets is practical for most applications, but that to provide reference calculations for certain problems (benchmark), where the ray effect is potentially very severe, spherical harmonic like formulations should be used. Lewis (Briggs et al., 1975) used phase-space finite elements in the second-order form of the one-group two-dimensional transport equation. They showed

Cross sections (cm-‘)

Region 1

Region 2

z, z, Source density (cmW3)

0.75 0.25 1.0

0.75 0.25 0.0

that piecewise constant approximations in angle substantially mitigate ray effects. Calculated flux along the line x = 2 cm are shown in Fig. 6, together with the results obtained from FEED2 in one- and two-moment modes. It can be seen that the Lewis et al. piecewise constant and piecewise bilinear fluxes for coarse and fine finite element structure show no sign of the ray effect,decreasingmonotonically, as yincreases. There is good agreement between the coarse-mesh FEED2 results and the Lewis et al. fine-mesh solutions. Our conclusion is that the coarse-mesh finite element results based upon the present work are as accurate as results obtained by using fine-mesh phase-space finite elements. The advantage of using continuous representations of the angular dependence, as in the present work, is that it guarantees theelimination of ray effects. In Fig. 7, FEED2 results for edge (x = 2 cm) flux for the same problem are compared with the predictions of Miller and Reed (1977). Reed used the first-order discrete ordinate transport equation with projection ,

v P

0.30

P

IF

&! u- 0.20

Vacuum

$ s v) 0.15

Perfect reflector

l-----l

QIO

I

0.0

ccastant

I

I

1.0

-x

klld

Fvrfect reflector Fig. 5. Geometry

piewise

FEED2 2-rOma)cwrsa

Vacuum

Region

FENT

FENT piwwise bitiir FEED2 I-mommt

Region 2

1.0

I

0.35

Q25

z,oi

693

of Case B one-group

ray-effect problem.

I ().O

0.5

IO

1.5

2.0

Y. cm Fig. 6. Case B one-group ray-effect problem. Flux along axis x = 2 cm from coarse mesh one- and two-moment FEED2 calculations compared with coarse and fine mesh calculations of Lewis.

694

A. K. ZIWR and A. J. H.

operators to obtain spherical harmonic like equations, with the aim of mitigating ray effects. In the conversion procedure from SNto P,_ 1 equations, fictitious sources were defined as ray-effect mitigating devices and the resulting equations were solved using an iterative algorithm. Figure 7 shows the best Reed solution, the discrete ordinate unmodified S, solution, and the FEED2 moment coarse-mesh solution. The agreement between Reed’s best solution and the FEED2 calculation is satisfactory. Finally, the results of Lillie (Lillie and Robinson, 1976) for this problem are compared with results of the present work, in Fig. 8. Lillie has also used the secondorder even-parity equation with a linear triangular finite element spatial and discrete-angle representations. Ray effects were reduced in the E, calculations shown in Fig. 8 but low-order Ez calculations yielded incorrect results. FEED2 one- and two-moment solutions are replotted in Fig. 8 to permit comparison between the two codes FEED2 and DAFE. DOT S, and Ss results of Lillie are also shown and exhibit oscillations. From the above comparison it has been shown that the present finite element solutions are not subject to ray effects even for low-order angular approximations. This is due to the continuous representation ofthe evenparity elliptic leakage operator which is invariant under angular rotation. 3.2. Two-Group Ray-Efect Problem This problem has the same geometry and boundary conditions as the second one-group ray-effect problem

1

GODDARD

0.15

0.10 05

00

1.0

1.5

ye Fig. 8. Case B one-group ray-effect problem. Flux along axis x = 1.875 cm from DAFE and DOT calculations of Lillie compared with one- and two-moment FEED2 calculations.

(Case B). The two-group constants shown in Table 4 represent a uniform isotropic source in the first energy group with strong removal of neutrons from group 1 to group 2. Due to the isolated neutron source and the large removal cross section in the first group, ray-effect distortions are to be expected in discrete ordinate finite difference results. Miller and Reed (1977) examined this problem by using an S4 -+ P, approximation to mitigate the ray effect. Their results for the group-l flux along the axis x together with the present one- and two-moment approximation iinite element results, are shown in Fig. 9. The Sr6 discrete-ordinate results of Reed is also shown to demonstrate the ray effect oscillations. Good agreement can be seen between modified ordinate solutions and the finite element results. 3.3. Streaming Problem This problem has been studied to investigate finite element modelling of streaming, where the dimensions of the void are of the same order as the mean-free-path in the surrounding medium. Geometry of the problem

0.15 .

Table 4. Cross-section and source data for two-group rayeffect problem with geometry illustrated in Fig. 5

00TS2

?? $(a43

/v=39) Ret Miller and Read (1977) O.lO- . FEED2I-munml

I

I 0.0

0.5

1.0

I.5

2.0

Y,cm

Fig. 7. Case B one-group raycffect problem. Flux along axis x = 2 cm from mod&d and unmodi6ed (DOT S,) calculations of Reed, compared one-moment REED2 _ _ with . cakulattons.

Cross sections km-9 z, z g-P r 1 Source densrty r-v.+ (cn-“)

Group 1 region 1 and 2 0.75 0.40 0.10 1.0 (reg.l)

Group 2 region 1 and 2 0.75 0.14 0.10 0.0

FE method for multigroup dilfusion-transport

695

problems

*FEE02 I-nK4mnt OTvmlwSE

0.24

oFEED2

2-malMnt

1.25-

9 c

Region I

Pwf6ct RfbCtOr

0.16

Vacuum

075-

0.16

0.14

0.12 00

2.0km)

05

Perfeet reflector

0.10

Fig 10. Geometry of one-group streaming problem. 0.06

a06 0.0

1.0

0.5

I.5

2.0

Y,cm Fig. 9. Two-group ray-effect problem. Group 1 flux along axis x = 2 cm of Fig 5 from one.- and two-moment FEED2 calculations compared to Reed’s S,, and S, -P P3 results. is shown in Fig. 10 and cross section and source data in

Table 5 ; absorption and total mean-free-paths are 2 and 1 cm, respectively, in regions 1 and 2. Lillie and Robinson (1976) analysed this problem using the discrete ordinate finite difference code DOT and the discrete angle finite element code DAFE. His results for flux along the line x = 1.875 cm are shown in Fig. 11, together with coarse-mesh FEED2 one- and two-moment results. Dot Sz and S, calculations show incorrect flux shapes in the void region. The onemoment FEED2 calculation shows a flux in the void which is flat and higher than the DOT S12 and DAFE E., results, while the latter two calculations show peaks in the void that Lillie has attributed to the ray effect. In the two-moment FEED2 calculation, the void flux is again flat, and is close to the DOT SIz flux. The differences between two-moment and S, z results, away Table 5. Cross-section and source data for the one-group streaming problem illustrated in Fig. 10 Cross sections (cm-‘)

z, R, Source density tcm-3)

Region 1

Region 2

1.0

0.5

1.0 0.5

1.0

0.0

Region 3 O.lE-0.8 0.0 0.0

0.0

0.5

1.0

Y,

1.5

2.0

cm

Fig 11. One-group streaming problem. Flux along line x = 1.875 of Fig. 10 from Lillie’s DOT and DAFE calculations compared with FEED2 flux.

from the void, may indicate the need for an improved angular representation. 3.4. P WR Cell Problem To demonstrate the ability of present finite element method to cope with coarse-mesh multiregion source problems, a PWR cell problem has been anaiysed. Geometry of the problem is given in Fig. 12 and represents an idealized PWR cell, for which group constants and isotropic source strengths are given in Table 6.

696

A. K. ZIVERand A. J. H. GODDARD

0 9066

7.7761km)

I.6602

77761

0.9086

Fig. 12. Group PWR cell problem showing regions and coarse finite element mesh. Natelson (1971) applied a variational method to discrete ordinate equations to solve this problem using a 44 x 46 finite difference mesh, while in the present calculations a coarse-mesh (16 x 17) finite element model has been used. Region absorption rate results from both methods are compared in Table 7. For this problem absorption rates for the discrete ordinate D,, solution are regarded as the reference solution. A comparison, in Table 8, of the percentage deviations from the reference solution, for Sr4 fine mesh and for the two-moment finite element coarse mesh demonstrates the coarse-mesh potential of the present finite element method in that roughly comparable accuracies are achieved. 3.5. Two-Group Bulk-Shield Problem A shield with a comer source has been considered to investigate the performance of the finite element Table 6. Cross-section data for one-group PWR problem shown in Fig. 12 Cross sections (a-‘)

method in deep penetration problems. Both iterative and direct Gauss elimination techniques have been applied and the results are compared with S, discrete ordinate finite difference solution (Argonne, 1972). This particular problem has been considered by Lathrop (Lathrop and Brinkley, 1970), Blomquist and Lewis (1980) and Asaoka et al. (1978), to test finite Table 7. One-group PWR cell problem region absorption rates (arbitrary units)

Regions

Two-moment finite element (coarse-mesh)

: 3 4 5

0.2642 0.2265 0.1772 0.1767 0.1488

(Ref. Nat%&

1971)

0.2409 0.2600 0.1743 0.1740 0.1498

Table 8. One-group PWR cell problem percentage deviation from reference.

Source densities (cm-“)

Regions

z,

x,

S,

Regions

Two-moment finite element (coarse mesh)

1 2 3 4 5

0.148934 0.151461 0.159109 0.159109 0.159109

0.049591 0.062158 0.064256 0.064256 0.064256

0.0 0.01048083 0.00217729 0.00214231 0.00215024

1 2 3 4 5

5.98 1.60 1.66 1.55 0.64

S4 (fint mesh) (NateLeon) 1.87 0.92 0.40 0.63 0.27

FE method for multigroup diffusion-transport

difference discrete ordinate transport and phase-space finite element codes against Pi9 reference solutions (Argonne, 1972). The geometry of the problem is given in Fig. 13 and the energy group constants in Table 9. In comparison with the reference solution, an oscillatory behaviour of the neutron distribution function (ray effect) was observed (Blomquist and Lewis, 1980)in the S, approximation. Further work in converting S, equations to P,_l equations was carried out (Blomquist and Lewis, 1980) to eliminate the ray effect. The present authors believe that the presence of the ray effect is due not only to the discrete representation of the angular dependence of the directional flux but also to the presence of the first-order form of the leakage operator in the transport equation. The results from the present finite element code show no ray effect even for P, (low order) angular approximation. When comparing the finite element coarse-mesh solutions against Ss discrete ordinate finemesh results, very good agreement has been obtained both for group 1 and 2. The comparison is shown in Fig. 14 together with the P, diffusion theory fine-mesh iterative solution. It can be seen that diffusion theory

Grwp

problems

697

results agree with the S, reference solution in the source region but not near the side of the shield, as expected. The use of the P, transport approximation reduces the discrepancy at the right boundary significantly (-4% for group 1 and 13% for group 2). Y

140 I

Perfect reflector

Perfect reflector

Perfect reflector

Shield

0

65

133

(cm)

Perfect reflector

Fig. 13. Geometry of two-group bulk-shield problem.

I

A .

60 x 60 FEM/P,

A

16 x 17 M/P3

,

61 x 64 Grwp

FD/SglReference) 2

A 0

6Ox6OFEMq

A

16 x I7 EM/~

0

61 x 64

FDE$(Refs~ca)

I

0

20

I

40

60

60

I

loo

Ix)

I

I40

x,cm

Fig. 14. Calculated neutron flux in the two energy groups, along the line y = 80 cm of Fig. 13.

698

A. K. ZIVER and A. J. H. &BDDARD

Table 9. Two-group constants for two-group bulk-shield problem Cross sections (cm- ‘)

Group 1

Group 2

z,

0.092104 0.006947 0.006546

0.100877 0.004850 0.023434 0.017701

z B-B z e-s+ 1 Source density (cmT3)

4. CONCLUSIONS

Linear finite elements have been used for spatial dependence, in conjunction with spherical harmonic expansions for directional representation, to solve the second-order even-parity transport functional in x-y geometry. It is shown from the results obtained with the codes FEED1 (Ziver, 1979) and FEED2 that this approach has the following potential relative to finite difference methods.

-ACCURACY: For coarse mesh models it has been shown that the finite element method gives greater accuracy. -RAY-EFFECT : Ray effects have not been observed in the present finite element solutions. -CONVERGENCE : For eigenvalue problems where the multigroup source iteration technique is employed, faster convergence has been achieved compared to discrete ordinate finite difference methods. -SPEED: It has been shown (Ackroyd et al., 1980) that present finite element solutions can be obtained at least as fast as finite difference solutions. Bearing in mind also the accurate treatment of irregular boundaries that is possible with finite elements, we believe that finite element methods have

future role in reactor radiation transport calculations.

an important

physics and

Acknowledgement-This work has been carried out with the aid of Science Research Council Grant GRA/A/8815.6.

REFERENCES

Ackroyd R. T., Ziver A. K. and Goddard A. J. H. (1980) Ann. nucl. Energy 7,335-849. Argonne National Laboratory Benchmark Problem Book (1972) ANL-7416. Asaoka et al. (1978) J. nucl. Sci. Tech&. 15,56-71. Blomquist R. N. and Lewis E. E. (1980) Nucl. Sci. Engng 73, 125-139. Briggs L. L., Miller W. F. and Lewis E. E. (1975) Nucl. Sci. iigrg 57,20.5-217. Lathroo K. D. 11970)Nucl. Sci. Enma 45.255-268. Lathrob K. D. aid B&kley F. W.( l%)T&ory and use of the general geometry TWOTRAN program. Report LA-9432. Lillie R. A. and Robinson J. C. (1976) A linear triangle finite element formulation for multigroup neutron transport analysis with anisotropic scattering. ORNL/TM-5281. Miller W. F. Jr and Reed W. H.(1977) Nucl. Sci. Engv62,391411. Natelson M. (1971) Nucl. Sci. Engng 43,131-144. Rhoades W. A. and Mynatt F. R. (1973) The DOT II twodimensional discrete ordinates transport code. ORIVLTM-4280. Zienkiewin 0. C. (1977) The Finite Element Method in Structural and Continuum Mechanics. McGraw-Hill, London. Ziver A. K. (1979) A finite element program for two-group diffusion and transport problems and its application to various benchmarks. Imperial College, Nuclear Power Section Report, SAF 017. Ziver A. K. (1981) A finite element program for two dimensional diffusion transport problems. Imperial College, Nuclear Power Section Report SAF 034. Ziver A. K., Goddard A. J. H. and Ackroyd R. T. (1981) Inz. J. numer. Meth. Engng. Submitted for publication.