Comparison of finite differences and finite elements in the case of a large fast power reactor

Annab o/-Nuclear Energy, Printed in Great Britain

0306-6549/81/l 10609-12$02.00/O Pergamon Press Ltd

Vol. 8, pp. 609 to 620, 1981

COMPARISON OF FINITE DIFFERENCES AND FINITE ELEMENTS IN THE CASE OF A LARGE FAST POWER REACTOR J. C. ESTIOT,D. HONDE,*G. PALMIO~TIand M. SALVATORES DRNR/SEDC/SPNR/LPR, Centre d’Etudes Nucleaires de Cadarache, B.P. No 1, F. 13115 Saint-Paul-lez-Durance, France Abstract-A large number of test calculations in two-dimensional hexagonal geometry for different configurations of SUPER-PHEIUIX 1 type Fast Reactor have been performed to compare finite ditference theory vs finite element theory performances. At present, no definitive advantages were found for the application of the finite element method to two-dimensional hexagonal design calculations.

1. INTRODUCTlON A large number of diffusion calculations are required at different stages of a Fast Power Reactor (FPR) design. Parametric studies are necessary at an early stage to investigate the relative advantages and compatibility of different technological solutions. Typical in this respect is the problem of optimizing the control-rod system, or, in the case of a heterogeneous core concept, the problem ofinternal fertile location. Usually at this stage of a design, a reduced multigroup scheme is used, but a full geometrical description is necessary, often in a twodimensional schematization and, much less frequently, due to computing costs, in actual three dimensions. At later stages of the design, validation calculations are necessary to provide the needed detail and accuracy on the final performance parameters of the system. It is then evident that there is a strong incentive to optimize the computing performances of the diffusion codes used. Many improvements have been done in recent years on the performances of diffusion codes based on the finite difference (FD) algorithms (Tucson Conf., 1977). A new possibility however could be presented by the use of algorithms other than the finite difference, as the finite element algorithms, used with success for reactor applications, other than fast reactors. In the present paper, we have investigated the possibility of using the finite element algorithms for typical fast reactor configurations. Our analysis will be focused on computer parameter performances, mainly the computing cost, as far as

possible comparing different algorithms at the same level of required accuracy. First, we will introduce the computing tools that we have used, which were developed in the framework of the CEA Reactor Programs. We will successively define the geometrical configurations chosen for the comparison, based on plane hexagonal representations of SUPER-PHENIX type reactor configurations. Finally, we will present and discuss several of the most interesting results so far obtained in the comparison between the finite element and the finite difference algorithms.

2. CALCULATIONMETHODSAND CODES

To compare finite element and finite difference algorithm performances we used two codes very recently developed in France: the BILAN code (Kavenoky and Lautard, 1977; Lautard, 1981; Dumas and Lautard, to be published) and the HEXAGONE code (Ravier et al., 1980). Both of them are based on the finite element and the finite difference methods. (The finite element method is not yet fully operational for HEXAGONE.) 2.1. Iteration strategy Below, we will briefly describe ticularities of these two codes :

the main

par-

A. the HEXAGONE code accelerates the fission sources for each outer iteration by means of the Tchebycheff polynomial method and uses for the inner iterations a method of resolution by alternating even * Present address: CISI (Compagnie Intemationale des and uneven lines ; B. the BILAN code uses for external iterations either Services en Informatique), 35, Boulevard Brune, 75680 Paris, Cedex 14, France. the Tchebycheff acceleration, applied to the flux, or the 609

J. C.

610

coarse-mesh rebalancing acceleration. iteration, the code uses -iteration by punctual overrelaxation, -Cholesky direct matrix inversion.

For

hTlOT

inner

built-in criteria), to obtain a stable (i.e. ‘well converged’) solution, both for the eigenvalue and for the fluxes.

2.2. Convergence criteria

3. GEOMETRY SPECIFICATION

The same convergence criteria imposed on the eigenvalue have been imposed on both codes. For the HEXAGONE code, these criteria are applied on the fission source defined as (1)

with

and Z, = t#~,= VP= n= g =

macroscopic fission cross section ; flux at point p ; mesh volume ; iteration number ; energy-group number.

No point-flux convergence criteria are imposed and the minimal and maximal punctual values of the fission source are provided : max min or@i I

et a/.

e

-1

(2)

I

The BILAN code has the following criteria : EPSI (eigenvalue) = l(K:f: ’- &)I

< 10e4

OF TEST PROBLEMS

Calculations were performed for a fast reactor of SUPER-PHENIX type (see Fig. 1). This includes two enrichment zones for the core, one blanket zone and lateral-shield zone. Its control system includes 21 main control rods divided in two rings, and three shutdown complementary control rods. The calculations are performed in two-dimensional hexagonal geometry with cross sections in a six-group energy structure (Nisan et al., 1979). The size of the configuration is limited to 18 assembly rings for a total of 919 assemblies. We considered three particular configurations with different numbers ofcontrol rods inserted in the reactor (see Table 1). The first one represents a startup configuration without control rods inserted. The geometrical description is limited to a third of the reactor, exploiting the periodicity of the configuration. The second one differs from the previous one by insertion of all control rods except three symmetric ones. The last configuration is a complete geometrical description ofthe reactor, with all control rods inserted except one. This configuration was intended to test an asymmetrical situation with a large number of points involved.

(3)

and 4. NUMERICAL CALCULATION OPTIONS

EPSI (flux) = max

4: I,,+,

< 10-z. I

(4)

These are built-in criteria. However, we were often obliged to increase the desired precision (i.e. beyond the

The two codes have several options which can be used either in finite dilference, for the HEXAGONE code, or in finite element for the BILAN code; these are resumed below (see Tables 2 and 3).

Table 1 No. of

No. of

Periodicity

calculation assemblies

non-inserted assemblies

Case 1: reference, no control rods inserted

2x13

307

l-24

Case 2 : all control rods inserted except 3

21113

307

18, 13,23

Case 3 : all control rods inserted except 1

Without (2@

919

7

Configuration

Finite dilkrences and tinhe elements in a large FPR

Case

I:

Control Case

assemblies

Control

I-

24

not i

location

for flux

2:

Control Case

611

assemblies

18,23,13

not inserted

3: assembly

7 not inserted

Fig. 1

Code HEXAGONE

The calculations are carried out in finite differences on a triangular basis with variable mesh sixe(see Fig. 2) : (a) FD7: 10.4 cm/mesh with 7 points/assembly (hexagonal assembly); (b) FD19: 5.2 cm/mesh with 19 points/assembly; (c) FD37 : 3.47 cm/mesh with 37 points/assembly. Code BILAN

The following options were used (see Fig. 2) : (d) for a hexagonal basis of Wachspress type (Haipern, 1979)-six nodes/assembly-the basis func-

tions are rational fractions which correspond to the ratio between a fourth-degree polynomial and seconddegree polynomial (FER6); (e) ‘non conforming’ elements (FENClZ), (Crouseix and Raviart, 1973; Lautard, 1981). They are linear polynomials applied on a triangular basis (with 12 nodes/assembly). Generally, there is no continuity between the six planes defined in the hexagonal assembly; (f) linear polynomials (FEL7) applying on a triangular basis to build a hexagon with 7 nodes; (g) same as option (f) but with a lumping technique which consists of diagonalixing the matrices of mass.

J. c.

612

ESI-IOT et d.

6 nodes / assemblyFER6

option (BILAN

/assembly FD7 option option

for

for the

rational fraction tinite elements code)

-Finite

difference

the HEXAGONE BILAN

code

12 nodes/assembly - FENCIZ option for (non -conformlng element 1

I9 nodes / assembly-EFPIS or points BILAN

code

and

FEL7

code

and FDl9

37 points / assembly FD37 option for the

the BILAN

and EFLIS option

options

for the

HEXAGONE

code

for the

HEXAGONE

code

code

Fig. 2. Geometrical representation of calculation options (see Tables 2 and 3 for option specifications).

This option is equivalent to the finite difference method (FEL7 with lumping) ; (h) linear polynomials (FEL19) applied on a triangular basis to form a hexagon with 19 nodes ; (i) parabolic polynomials (FEP19) applied on a triangular basis to form a hexagon with 19 nodes. 5. RESULTSAND DISCUSION We will analyse the results in terms of finite difference vs finite element performances, The chosen performance parameters are the calculation precision and some computer parameters such as: the CPU time; the central memory occupation; and the related cost, which also takes into account the I/O operations.

Even if the cost is the most interesting parameter, it should be used carefully, since it is strongly dependent on each particular computer environment. In our case, we used the IBM 370/168 of the CEN/CADARACHE. 5.1. Case 1 In Table 4, the eigenvalue results for Case 1 are shown. We recall that this configuration, a startup situation for a SUPER-PHENIX 1 fast-reactor type, is a relatively ‘simple’ case since no sharp gradients are present (see.Fig. 3). Therefore no large-mesh-size elk&s are involved. Nevertheless some interesting results arise from the examination of Table 4. It should be noted that for all the cases the FD19 and FD37 HEXAGONE

Finite dilTerencesand finite elements in a large FPR Table 2. Number of calculation points or nodes for all the

options No. of points or nodes in the calculated configuration Used options FD7 FD19 FD37 Rational fractions FER6 ‘Non conforming’ elements FENC12 Linear polynomial FEL7 Linear polynomial FEL7 with lumping Linear polynomial FEL19 Parabolic polynomial FEP19

calculations respectively,

2~/3

2n

981 3799 8455

2863 11238 -

674 3684 981

2863

981

2863

3799 3799

are initiated using the flux obtained, from the FD7 and FD19 calculations, this

explains the small number of outer iterations needed to converge. Unfortunately, this option is still not available in the BILAN code. A further remark regards the coarse-mesh rebalancing acceleration used in the BILAN code, which gives some advantage on the total number of outer iterations, but a greater effort, i.e. calculation time, is needed for each outer iteration. A comparison on the FEL7 option using Tchebycheff acceleration has shown a substantially equal cost. On the contrary, if the configuration involved sharper gradients as in Case 2, or a large number of Table 3

Method FD7 FD19 FD37 FEL7 FEL19 FER6 FENCl2 FEP19

Finite difference Finite difference Finite ditkence Finite element with linear polynomials Finite element with linear polynomials Finite element with rational fractions Finite element with ‘non conforming’ elements Finite element with secondorder degree polynomials

No. of points or nodes for each hexagonal assembly 7 19 37 7 19 6 12 19

613

calculation points, some instability was found when the coarse-mesh rebalancing acceleration was used. Therefore it was decided for all future cases to use the Tchebycheff acceleration. First of all, as we expected, the performances of the two FEL7 calculations with and without lumping transformations (i.e. equivalent to FD7) are almost equivalent, and this is due to the same number of unknowns being involved in the calculation. Therefore the only advantage that we can expect from a finite element calculation of this type is an improvement in the precision of the calculated eigenvalue and fluxes. Actually, since the HEXAGONE and BILAN codes use the point-centered flux scheme, for this type of situation the calculated eigenvalue, related to large mesh size, is underestimated with respect to the eigenvalue calculated at a mesh size which approaches zero. To obtain the eigenvalue we must apply the Kato Law (Kato et al., 1976) which allows a linear extrapolation with an inverse dependence on the surface associated with each calculated point. This law is an immediate result of the finite difference approximation when a regular mesh-size grid is used, since third or higher-order terms are neglected in the Taylor expansion of the gradient operator in the diffusion equation. Therefore, if we want to calculate the KFEm, we must use the following formula

This represents the usual formula applied in the SUPER-PHENIX 1 design which takes account of the mesh-size effects on the eigenvalue. Since for design purposes we need to take account of effects of this type, at least two calculations of the same configuration with two different mesh sizes are always necessary. Clearly, at this time, an important problem arises : how can we calculate the eigenvalue at a mesh size approaching zero in the finite element theory using different orders of polynomial approximation? This is theoretically unclear. The only possibility is represented by increasing the number of nodes without changing the polynomial approximation(i.e. FEL7 and FEL19 calculations). We will examine this type of calculation for the next case. Theextrapolatedeigenvaluefor Case 1 usingequation (5) is 1.04247 ; therefore the FEL7 calculation without lumping shows an effective improvement of the eigenvalue of 0.023% on a total difference of 0.117% between the FD7 calculation and Kzm. If we look at the relative discrepancies in the total flux calculated at the points shown in Fig. 1 (seeTable 5 and Fig. 4), where the FD37 calculation was taken as reference and the

614

J.

c.

Esn0.r et al.

@

: Case I

0:

Case 2

@

: Care 3

R, cm

Fig. 3. Point total flux distribution for the three cases in a selected direction (A --*F of Fig. 1). fluxesarenormalizedat thecentreofthereactor, wecan observe a relevant discrepancy starting from the blanket zone down to the lateral-shield zone, namely where the flux gradient begins to increase. At this point, we have to make two remarks -first, the FD37 eigenvalue is slightly out of the linear law of extrapolation (see Fig. 5), and we will talk more extensively about thii discrepancy for Case 2 ; -second, the performances of FD7 calculation with the HEXAGONE code are better than those of the BILAN code, and this is essentially due to the exploitation in HEXAGONE of an extension of the FORTRAN language called ESOPE (Mougin, 1979) which allows a better data structuralization and automatic management of auxiliary memories. If we go back to Table 4, we observe a very poor result from the FER6 option, which is confirmed by the

examination ofthe total flux discrepancies distribution. This makes the FER6 option unsuitable for our Purposes. An interesting result was obtained applying the direct inversion calculation (Cholesky method) to this option, since this case has a limited number of calculation points and therefore it is the most suitable method. The corresponding result in Table 4 shows an increased total amount of CPU time, and consequently a higher cost. Therefore it was decided to eliminate this method for future cases. The FENCl2 calculation shows an overestimation of the eigenvalue with regard to Kgm, whereas the FEP19 is better than the FD19 value. Moreover, the total flux discrepancies are very pronounced for the FENCl2 and some problems also appear for the FEP19 (see Fig. 4). An essential point is the prohibitive cost of these two

615

Finite di6erenaa and finite elements in a large FPR Table 4. Case 1: Relative pCiLli0n EipenVdd

Method (code) FD37 FD19 FD7 FELT FEL7d FER6 FER6C FENCl2 FEP19

(K,)

(BILW (BILAN) (BILAN) (BILAN)

RclatiVe pre&ion

Total number

mQnorY

CPU

cost

oampation

time

(cost

(s)

unities)

71 36 24 47 47 33 46 199 520

172 63 35 118 111 85 115 350 loo0

readdon

reachedon

of outer

C!igenVdUe

flUXa

iterationsb

2OE-03 24E-03 2OE-03 3.3E-03 4.OE- 03 6.OE- 03 8.OE-03 3.6E-03 l.OE-03

7 7 19 14 14 11 10 39 47

1.6E-05 1.4E-05 8.6E-06 24E-05 4.2E - 05 5.5E - 07 26E-05 l.lE-05 l.OE-05

1.04238 1.04218 1.04130 1.04153 1.04140 1.04080 1.04080 1.04257 1.04228

(HESAGONE) (HEXAGONE) (HEX&ZONE) (BILAN) (BILAN)

results Central

(WW

1800 loo0 500 z 450 450 loo0 1000

.For the HEXAGONE axle the precision is reached on fission sources, see Section 2.2. b If not specSed, five inner iterations for each outer arc performed. ‘Coarse-mesh rebalancing aaderation applied on outer iterations. d With lumping transformation, quivdent to FD7 (see text).

‘Direct invemion (Cholesky method), no inner iterations peifonncd. r Zero meshsizeextrapolated eigenvalue Kp = 1.04247. ‘1.6E-05

= 1.6 x lo-‘.

* With lumping transformation (see text 1

FER 6 -2o-

Core

I

1 Core 2

/Blanket 1 Loteral

Fig. 4. Case 1. Relative discrepancy in the total flux.

shield

616

J. C.

&TOT

et d.

Table 5. Case 1: relative discrepancies in percentage of total 5ux distribution’ No.

of assembly

Method (code)

2b

B”

CC

Dd

E’

F’

FD7 (HEXAGONE) FD19 (HEXAGONE) FER6 (BILAN) FER6 Cholesky (BILAN) FEL7 with lumping (BILAN) FEL7 (BILAN) FENC12 (BILAN) FEP19 (BILAN)

-0.5 -0.1 0.6 0.6 0.4 0.05 3.4 0.8

0.05 0.01 1.4 1.4 0.01 0.1 4.8 2.0

0.1 0.02 2.0 2.0 23.0 0.4 1.6 3.7

0.1 0.1 0.7 0.2 0.5 0.5 9.9 4.5

-0.5 -0.5 12.3 8.1 -0.6 -2.9 14.3 4.5

2.7 0.21 - 17.0 15.3 2.2 -15.9 13.7 0.15

“The 5ux reference distribution is the FD37 calculation. b Main control rod. ‘Inner fuel. *External fuel. ‘Blanket. f Lateral shield.

options, due to the fact that it is impossible at the moment to initiate the calculation exploiting the flux obtained from a lower-order approximation calculation.

This seems due (Rowland et al., 1979) to the technique used to calculate the eigenvalue. We recall that in the HEXAGONE code the Keff is calculated as

5.2. Cases 2 and 3 If we examine Case 2, where sharp gradients due to rod insertion cause large-mesh-size effects. The discrepancies found in the previous case are underlined in this situation. If we look at Table 6, we can observe a very relevant overestimation of the eigenvalue of the FENC12 calculation with respect to the KFim. The FEP19 also shows an overestimation, but the FD calculation has similar problems. In fact, if we use for instance the FD7 and FD37 eigenvalues to extrapolate the Kz@ eigenvalue, we would find a discrepancy of 0.037% in K,, (see Fig. 6).

Ksff

=

n = iteration index p = point index.

Npp=’ c pti.-lp*” p=1

This formulation may cause rounding errors to appear due to the fact that the absorption term is some order ofmagnitude lower than the term of the removal. The calculation ofK,, as the balance ofproduction and loss in the diffusion equation, should eliminate this inconvenience. For Case 2, we also ran the FEL19 calculation to allow an extrapolation of Keff for the finite element

Table 6. Case 2 : results

Method (code) FD37 FDL19 FD7 FEL7 FELT FEL19 FENC12 FEP19

(HEXAGONE) (HEXAGONE) (HEXAGONE) (BILAN) (BILAN) (BILAN) (BILAN) (BILAN)

Eigenvalue* U&r)

Relative precision reached on eigenvalue

Relative precision reached on 5ux’

Total number of outer iterationsb

Central memory occupation (kbytes)

CPU time (s)

Cost (cost unities)

0.95170 0.95022 0.94404 0.94612 0.94404 0.95117 0.95482 0.95286

2.4E-06 3.4E-06 3.8E-07 4.5E-06 3.OE- 07 5.OE-06 8.OE-06 l.OE-05

8.6E-04 2.4E-03 2.2E-04 l.OE-02 l.OE-03 1.3E-03 1.6E-03 l.OE-03

19 18 24 18 23 45 49 53

2ooo 1000 475 450 450 loo0 loo0 1000

195 83 29 43 57 305 207 547

450 158 47 106 127 725 401 1097

‘For the HEXAGONE code, the precision is reached on fission wurces, see equation (2). b If not specifid 5ve inner iterations for each outer are performed. ‘With lumping transformation equivalent to FD7 (see text). *Zero mesh size extrapoled eigenvalue Kz:” = 0.95228.

617

Finite differencesand finite elements in a large FPR

1.043cn

4

+

FD 37

\

\

FD7 +\

I so

X_Assemblyawfwa

I 100

L

x'

Assembly

surface

Na.af paints/assembly

No.ofpdnfs/asssrnbly

Fig. 6. Case 2. Eigenvalue extrapolation.

Fig. 5. Case 1. Eigenvalue extrapolation.

calculations also. The extrapolated eigenvalue in this case is Kz* = 0.95285, thismeans a difference of 0.057% with respect to the Kzm. Thiscangenerate an ambiguity in the interpretation of the results obtained for this case. Actually, some suspicions exist for the

finite element calculations,

since we found negative fluxes, in particular in the control-rod regions. This explains the large discrepancies shown in Fig. 7 and Table 7. We recall that for design requirements, the point

Table 7. Case 2 : relative discrepancies in percentage of total flux distribution No. of assembly Method (code)

2b

BG

FD7 (HEXAGONE) FD19 (HEXAGONE) FEL7 I:::$ FEP19 FENCl2 (BILAN) FEL7 with lumping (BILAN) FEL19 (BILAN)

0.6 -1.6 -15.4 -7.5 -17.3

7.2 -1.4 13.0 -7.6 - 10.0 9.0 -5.5

-1.1 -6.8

CC

Dd

7.0 -2.1 14.3 -9.2

11.4 -1.4 16.8

12.7 - 1.2 14.7

- 10.1 - 14.0 15.4 -7.2

- 10.2 -0.3 15.0 -7.8

- 12.5 9.3 -6.9

‘The flux distribution reference is the ED37 calculation. bMain control rod. cInternal fuel. *External fuel. ‘Blanket. ‘Lateral shield.

E’

F’ 18.2 -1.9 1.6 - 12.7 -11.0 20.8 - 12.4

618

J. C. Esrtrrr et CJ.

/

Core

1 Core2

I

]Btonkrr

Fig. 7. Case 2. Relativediscrepancy power distribution

must be known with an uncertainty

N + 3%, and that this must include the uncertainties due to the cross-section data. In this case, the discrepancies between the three finite

1 Lateral

:FD7

shield

in the total flux.

difference calculations, seem to be more regular than those observed among all the finite element calculations. Finally, we examine the Case 3 results (see Tables 8

Table 8. Case 3 : results

Method (code) FD19 FD7 FEL7 FEL7’

(HEXAGONE) (HEXAGGNE) (BILAN) (BILAN)

Eigenvalued KU)

Relative precision reached on eigenvalue

Relative precision reached on flux

Total number of outer iterationsb

central memory occupation &bytes)

CPU time (s)

Cost (cost unities)

0.92498 0.91540 0.91858 0.91496

4.1E-06 5.3E-06 l.OE-05 l.OE-05

LlE-03 4SE-03 2.OE- 02 2.OE- 02

16 38 23 23

u)oo 900 1000 loo0

198 109 161 169

550 173 325 334

’For the HEXAGONE code, the precision is reached on fission sources, see equation (2). b If not specSed, five inner iterations for each outer one are performed. GWith lumping transformation equivalent to FD7 (see text). *Zero mesh size extrapoled eigenvalue KrtIrn = 0.92814.

619

Finite differences and finite ekments in a large FPR

FEL 7 ,+ (Lumping) / / / /

Core I

Laterol

Blanket

Core 2

shield w

Fig. 8. Case 3. Relative discrepancy in the total flux.

Table 9. Case 3 : relative discrapancies in percentage total flux distribution No. of assembly

Method (code)

Zb

B’

(75

Dd

E’

FD7 (HEXAGONE) FEW (BILAN) FEL7 with lumping (BILAN)

- 3.3 -20.0b’ -2.6

-0.2 1.4 2.2

2.7

4.0

9.3

1.5

::: 3.3

5.1

3.1 9.9

-7.9 16.1

- 14.1 -13.1

“The flux reference distribution is the FD19 calculation. ’Main control rod. b* Negative flux. CInternal fuel. d External fuel. ’Blanket. ’Lateral shield. gAuxiliary control rod.

8.8

F’

78

620

J. C. Esrmr et al.

and 9), where only a few of the possible optional calculations were used due to the cost related to the large number of points involved. As in the previous cases, the FEL7 calculation shows an improvement of the same order of magnitude on the eigenvalue. In Fig. 8, we show the relative discrepancies in the total flux. In this case, the FD19 was taken as reference. It seems that also the FEL7 case with lumping has some problem because of the large differences found in the flux distribution, with regard to the finite difference calculation; this also explains the difference found between the two eigenvalues. 6. CONCLUSIONS

A large number of test calculations in hexagonal geometry for different configurations of a SUPERPHENIX type fast reactor have been performed to compare finite difference vs finite element approximation performances. No definitive advantages were found for the use of the finite element method, since in final design stages eigenvalues corrected for mesh-size effects are required and, consequently, at least two calculations with different mesh sixes are always required. This means that when the same mesh size is used, the improvement found for the finite element eigenvalue is not very relevant, since the finite difference and finite element calculations have a substantially equal cost, due to the substantial, equivalent, number of unknowns involved. Moreover, the same discrepancies were found when the total fluxes were compared. For historical reasons, namely the use of finite differences as a reference in reactor physics design for several years, we still consider that the finite difference calculation is a reference. Therefore we request further

investigations to explain the discrepancies found between the finite difference and finite element results. Some interesting possibilities still remain for the finite element method in threedimensional geometries, in particular Hexagonal, 2. In this case, the reduction of planes in the Z-direction could give a consistent gain in calculation time possibly without loss of accuracy, when the finite element method is used.

REFERENCES

Crouseix M. and Raviart P. A. (1973)Conforming and nonconforming Finite Element methods for solving the stationary Stokes equation RAIRO Analyse Numerique. Dumas M. and Lautard J. J. Le code BILAN. SERMA/CEA. To be published. Halpern P. (1979) Approximation d’un problZme elliptique d’ordre 2 sur un ouvert compose d’hexagone par &ment finis rationnels de type Wachspress. Rapport inteme SERMA/CEA. Kato Y., Takeda T. and Takeda S. (1976) Nucl. Sci. Engng 61, 127. KavenokyA.andLautardJ. J.(1977)Nucl.Sci.Engng64,563575. LautardJ.J.(1981)NewFiniteElementrepresentationfor3-D reactor calculation. Advances in mathematical methods for the solution of nuclear engineering problems. Munich Conference, April 27-29. Mougin J. (1979) Systeme RAMSES. Rapport CISI (Compagnie Interttationale des Services Informatiques). Nisan S., Hammer P., Lautard J. J. and Mus M. (1979) Colloque international sur la Physique des Rkacteurs Rapides, Aix-en-Provence, France.,September 24-28.

Ravier M., Santamarina C. and Le. Cardinal G. (1980) So&&cation des modules neutroniaues. Ratmort CISI No. 79.018 (Compagnie Internatidnale &s Services Infotmatiques). Rowland J. L., Sweet D. W., Franklin B. M. and Sunderland R. (1979) Collogue International sw la Physique des Rlacteurs Rapides, Aix-en-Provewe, September 24-28, IAEA-SM244136. Tucson Conference (1977) Nucl. Sci. Engng 64,2.