- Email: [email protected]
A hexagonal geometry nodal expansion method for fast reactor calculations
0149~1970/86 $0.00+,50 Copyright ~ 1986 Pergamon Journals Ltd.
Progress in Nuclear Energy, Vol. 18, No. 1/2, pp. 113 121, 1986. Printed in Great Britain, All rights reserved.
A H E X A G O N A L G E O M E T R Y N O D A L E X P A N S I O N M E T H O D FOR FAST REACTOR C A L C U L A T I O N S J. M. PUTNEY* CEGB, Central Electricity Research Laboratories, Kelvin Avenue, Leatherhead, Surrey, U.K.
Abstract--A hexagonal geometry extension of the nodal expansion diffusion theory method is presented together with solution and acceleration techniques which have proved effective for two-dimensional (2D) and threedimensional (3D) fast reactor calculations. The method uses a radial nodal flux expansion which has hexagonal symmetry. A one-dimensional (1D) expansion is used for the axial direction. Numerical results are presented for four-group LMFBR benchmark problems and a 37-group PFR problem. A zeroth order (quadratic degree) radial expansion is found to be sufficient to permit calculations to be performed using one mesh per subassembly, whereas higher order expansions are more efficient in the axial direction because a coarser mesh can usually be used in this direction. For comparable accuracy the method is found to be a factor of 2 faster than the finite difference method for two-dimensional calculations and an order of magnitude faster for three-dimensional calculations•
1. INTRODUCTION In recent years the advent of nodal diffusion methods (Doming, 1979) has enabled computing times for LWR analysis to be significantly reduced. The nodal expansion method (NEM) (Finnemann et al., 1977) has proved particularly effective. Its basis is the solution of equations for the average flux on a coarse spatial mesh, combined with a high order local polynomial flux expansion, which is used to determine a relationship between the mean interface partial currents and the mesh average fluxes. Being an LWR method however, the N E M is formulated in rectangular geometry. For fast reactor analysis, extension to hexagonal geometry is required, In addition, suitable solution and acceleration techniques must be developed for general fast reactor multigroup problems. This paper describes a hexagonal geometry N E M and associated solution techniques which have proved effective for two-dimensional and three-dimensional fast reactor calculations. Numerical results are included which demonstrate the efficiency of the method over the finite difference method.
energy group g integrated over the volume of (homogeneous) node FI~ (the nodal balance equation) can be written in the form:
~w Aw's m m_ E ~ s:;,oT, ~ ~i÷m ~,.ws-j-m~ ~w~,+z.,
where 0~' is the average node flux (nodal flux) and Jowl +m and j g ~ are the average outgoing and incoming partial currents (interface currents) at the left (s=l) or right (s = r) node surface in the w-direction (surface F~). The interface currents at surface F ~ are related to the local flux gradient by Fick's law
Jo,.., -Y0w~ = - c~
•
.
2.1. General background The multigroup neutron continuity equation for
fl, s=r
as=].
(2)
1, s=l
In rectangular geometry, by reversing the order of integration and differentiation, equation (2) may be written as
dw
r,%
(3)
where ~,~w(W)is the one-dimensional nodal flux formed by averaging the flux in node FIm over the direction transverse to the w-direction. This flux satisfies a onedimensional nodal diffusion equation 2
* This work was undertaken while the author was a Research Student at Imperial College London, and was carried out at AEE Winfrith.
dFm,
Aw~J~,
Jg+~'~-J°-w;= - f "D~ 2. TWO-DIMENSIONAL NODAL EQUATIONS
i1)
g'=l
-D~
m
d ~baw
L"
G
dw 2 +E,mo~lgmw-I- gw = E Sg~'~am'w, (4) g'=l
where L~w(w) is the transverse leakage. 113
114
J. M. PUTNEY
The essence of the NEM is the determination of a polynomial expansion for ~bg%. In the zeroth order variant, a quadratic expansion is assumed and fitted to the average value of ~"~ in the node (¢b~")and its values m m m m at Fwt and F~, (denoted by ffgwt and ~b0~,). Introducing the Fick's corollary
Oo.~- 2(Jo.~ +3a~ ) m
--
• +m
• -m
(5)
allows the expansion coefficients to be expressed solely in terms of O~' and (jo+~+J0-w~)' Substituting the result into equation (3) then leads to a set of interface current equations which, together with the nodal balance equation (1) and the continuity and boundary conditions, is sufficient to determine all the nodal fluxes and interface currents. In the nth order variant (n_> 1), an n + 2th degree expansion is constructed by adding on higher order polynomial terms. The coefficients of these are determined from local weighted residual solutions of the one-dimensional diffusion equations, a process which requires the introduction of a polynomial approximation to the transverse leakage. The principal difficulty in applying the NEM to a hexagonal geometry, is that the process of changing the order of integration and differentiation in equation (3) introduces the flux values at the corner points of F~. However, the one-dimensional rectangular geometry expansions are equivalent to a general twodimensional nodal flux expansion. In the NEM the cross-terms in this expansion either drop out completely or appear only in weighted integrals of the transverse leakage• The remaining terms are simply the one-dimensional nodal flux expansions. With a two-dimensional expansion the interface current equations can be derived directly from equation (2). In hexagonal geometry, this avoids the need to estimate corner point fluxes. The problem is thus to find a suitable expansion for the flux in a hexagonal node.
GG.~..L)=Ag'ho+
Y { ,whgG)+b,~hdG)} a
m
m
w=x,u,v
+ coh l(~x)h 1(~u)h l(~v)'
(6)
where hi(~w) is a polynomial of degree i in Cwand {ho, hi (~0, h2 (¢~)} is orthogonal on integration over the volume of the node. A suitable set of expansion functions is h o = 1, h(~) = ~w, h2 (~w)= ~
5 72.
(7)
Introducing these into equation (6) and transforming into local rectangular cartesian variables, leads to a nodal expansion having seven independent terms. The coefficients are obtained by fitting to the node and surface average fluxes. Using the nodal expansion to evaluate equation (2) for each surface and eliminating the surface fluxes using equation (5), leads to six relations of the form: (D~o~+n7"].+m
T/ox.
7 +m
(.
o
=
)'29 +,
+
29._ m
Jow,;
7
T d/ . ; . 11._,.)
- 10j.w,
ll.+m~
,vO0s.w.+ w=u,v 2 36
5-* 7.
(8)
Inverting the resulting 6 × 6 system to express the outgoing currents in terms of the incoming currents, yields a set of interface current equations of the form: • +m
laxr
m
=
CffI
* -m
rn
. --m
Jffxr + CozJoxl
m m. t + C ~.ms , , . 3Jgwr+Co4Jffwl
+ E
(9)
w=u,v m Here, the nodal coefficients Co~ - C o ms are formally rational functions of D~"/H. In practice, they can be determined efficiently using a numerical inversion procedure based on Cholesky's decomposition• Equations (9) can be used to eliminate the outgoing currents from the nodal balance equation (1):
2.2. A zeroth order hexagonal geometry formulation It is convenient to introduce a set of'hexagonal axes' (x, u, v) where axes u and v are at - 120° and + 120° to the x axis. In this system the co-ordinates of a point P1 are (xl, ul, v0, where wl(=xl, ul, vl) is the distance along the w-axis to the perpendicular from P~. It is also convenient to define the local dimensionless variables, ~ = w / H , where H is the across flats width of the hexagon. For the interface current equations to have hexagonal symmetry the expansion should be symmetric in the variables x, u, v. The zeroth order method presented here is based on the following choice:
0'=1
2
+ 3H{ 1 - C ~ - C ; 3 - 2 c o m 3 - 2 C ; ~ }
Z w=x,u,v s=l,r
Jowl" --
(10)
Equations (1)/(10) and (9) represent zeroth order hexagonal geometry NEM equations, expressed in a form analogous to their application in rectangular geometry (Finnemann et al., 1977). Together with the continuity of interface currents and boundary conditions they are sufficient to determine the nodal fluxes
115
A hexagonal geometry nodal expansion method and interface currents. In contrast to rectangular geometry though, all coefficients in the expansion are determined and the problem has not been decoupled into a series of one-dimensional problems. Although the method can be generalized to higher orders of expansion (Putney, 1984) this has been found to be unnecessary for the hexagonal mesh required to represent fast reactor subassemblies. Other hexagonal geometry nodal expansion methods have been developed, for example by Duracz (1981), Lawrence (1983) and Arkuszewski (1982). Duracz and Lawrence treat the one-dimensional transversely average flux in each of the three hexagonal geometry directions, x, u, v. Half-range quadratic polynomial expansions are used to fit the averaged flux and the two surface fluxes in a way which treats the variation in transverse width of the node. Arkuszewski expresses the flux in the node in terms of six fundamental solutions of the diffusion equation inside the node. Using these the outgoing currents are expressed in terms of the incoming currents.
The zeroth order expansion coefficients are determined by fitting to the node and surface average fluxes: m
m
4
~.(~. ~, ~0=4,~(~. ~)+ Y~ d.~,g,(~O (11) i=1
where ~, = z/h " is a dimensionless axial variable for the node. Here ~b~2~ is the two-dimensional hexagonal expansion and ,q~(~=)is a polynomial of degree i in ~= such that
m
--
m
m
m
To determine the higher order expansion coefficients a local weighted residual equation is formed by multiplying the diffusion equation by a weighting function W k ( z ) ( k = l , 2) and integrating over the volume of the node. Two weighting schemes have been considered: Galerkin ( W k = gk + 2) and moments ( W k = g0. After integrating over the radial directions, the diffusion equation in the residual reduces to the onedimensional equation (4) (with w--z). The transverse leakage is expanded as m m m L,=(~=)= L,~ o + Lax i O,(~=)+ L,~2 92(~=)
(14)
and fitted to the average transverse leakages (functions of the radial interface currents) in the node and its two axial neighbours (Lawrence and Dorning, 1980). For axial boundary nodes the flux boundary condition is applied to Lg=(~z).The weighted residual equations for m m the higher order coefficients, dg~3 and do= 4 , are
}
3. THREE-DIMENSIONAL NODAL EQUATIONS The derivation of a three-dimensional hexagonal geometry expansion requires the combination of the two-dimensional expansion with a one-dimensional rectangular geometry expansion in the axial direction. Cross terms involving hexagonal and rectangular expansion functions need not be included explicitly, as due to the orthogonality of the directions they behave exactly as for the rectangular geometry method. Since the axial variations in composition in fast reactors generally permit a relatively coarse axial mesh to be used, it is appropriate to consider higher order axial expansions. The method presented here is based on second order (i.e. fourth degree) axial accuracy. Following Finnemann et al. (1977) for the axial treatment, the following nodal expansion is assumed:
m
dozl = @,a=,- d/g=, dog 2 - 3(q%=,- 2Og + q/a~,)- (13)
J'D•' "do:, )h,.2 + Ak Y~,g
°
2
= 0'=1
" (1)0'k m
Sgg,
q- O k zrmgdgmk - BkL,S~k
(15)
where k = i - 2 and the weighted flux (I),"k = (A k do==~ Bkdo~k). For Galerkin weighting A I = ~ , A2=~o, Bm= ~, B2 = 1½. For moments weighting A I = ~o, A 2 _-
1 ~16 ,
BI =16, B 2 = 11.
(16)
For the nodal expansion (11) the radial interface current equations are identical to their 2D counterparts. At the axial surfaces, the equations take the form .+ra m m - m ra . - " m m Jg~l = Ag= q)gm + Bg~jg~t + Cg=jg=r + fitEg= dg=3
(17)
+ F~m~dgmz4
where A o ~ - F ~ are rational functions of D ~ ' / h " obtained analytically by inverting the 2 × 2 matrix. Eliminating the outgoing currents from equation (1) gives:
c,"; + ~ -
+ z 7; * ;" = Z s.";,• ~; g'=l
{1, gl(~=), g2(~=)} and {l, g3(~z), ga(~z)} are orthogonal on integration over the volume of node l-I', and g3(~) and g4(~,) integrate to zero over its axial surfaces. A suitable set of expansion functions is •q l ( ~ z )
= ~z'
g2(~z) 2
1
=~2-
= ¢ z ( ~ z -- 4),
---m
w=x,u,v s=l,r
Jgws
1
12'1 g 3 ( ~ z )
2 1 2 O4(~z) = (¢z -- 4)(~z -- ~)"
2 + 3H {1-Co~ - C Z - 2Co"3 - 2C0=4 }
+ h" {(1 - B ~ - Co"=) (Jo-~' +Jo-=r") -- 2F~=~dg~,, } (12)
(18)
116
J.M. PUTNEY 4. SOLUTION TECHNIQUES
The nodal equations presented above have been implemented in a computer code HEXNEC. The solution procedure reduces the nodal balance equation (1) to a system of nonlinear finite difference type equations by defining leakage coefficients +m
• +m
ra
~gw~=Jgw~/~g
(19)
which are used to eliminate the interface currents in terms of nodal fluxes: • +m
+m
m
Jgws = ~:gw~~ g
(20)
with the incoming currents being expressed in terms of the leakage coefficients and average fluxes for the neighbouring nodes. The resulting 'pseudo finite difference equations' are solved using a conventional fission source iteration scheme in parallel with an (outer) iterative update of the leakage coefficients based on the other nodal equations. The following update procedure is used: Using one Gauss-Seidel sweep, for each node (i) For the three-dimensional problems with higher order axial treatment compute the transverse leakage moments from the latest values of the interface currents and solve the weighted residual equations (15) for the higher order expansion coefficients. (ii) Use the outgoing current eliminated nodal balance equation (10)/(18) to update the nodal flUX.
(iii) Use the interface current equations (9) and (17) to update the outgoing currents and apply the continuity and boundary conditions to update the corresponding incoming currents. When solving the weighted residual equations (15) the weighted sources are computed in the same way as the normal sources. The zeroth order expansion coefficients are calculated using the latest values of the nodal fluxes and the interface currents. As well as being used to update the nodal fluxes and interface currents, the resulting higher order expansion coefficients are combined with the zeroth order coefficients to form new weighted fluxes which are then stored for future use in calculating weighted sources. This differs from the approach outlined by Finnemann et al. (1977) for the rectangular geometry method, in which the expansion coefficients are re-calculated when needed. The adopted approach involves less arithmetic and the extra storage required is acceptable as the weighted fluxes are only present in one direction. The outer iterations in HEXNEC are accelerated by asymptotic extrapolation (AX) (Wagner, 1968) of the
fission source, followed by renormalization of the nodal fluxes and currents. This reduces the number of outer iterations by 25-40%. Coarse mesh rebalancing may be a more efficient acceleration technique, but this is dependent on choosing a coarse-mesh and rebalancing strategy appropriate for the application. The inner iterations use a line Gauss-Seidel method (Varga, 1962) accelerated by AX. For two-dimensional problems, the number of unaccelerated iterations required is not much greater than the 'threshold' number required before AX is permitted (and sometimes less). Consequently, the reduction in the number of inner iterations per outer is around 10% or less. On threedimensional problems however, a reduction of around 25% can be achieved. Over-relaxation was found to be less effective than AX. One some problems HEXNEC has been found to produce some negative nodal fluxes and/or interface currents. Although these generally only occur at unimportant locations, they can lead to divergence of the inner iterations. Consequently, the definition of the leakage coefficients was modified to ensure that they are always non-negative. If an interface current becomes negative then, in calculating the leakage coefficients, it is interpreted as positive current in the opposite direction. If a nodal flux becomes negative, it is multiplied by - 1. The inner iterations then yield the negative of the flux which is multiplied by - 1 before proceeding to the update stage.
5. NUMERICAL STUDIES HEXNEC calculations for three four-group LMFBR benchmark problems (Buckel et al., 1977) and a 37-group (7-downscatter) PFR problem have been compared with results obtained used the finite difference (FD) codes (TIGAR) (Matthews, 1972) and SNAP (McCallien, 1975). Both these F D codes employ mesh-centred differencing, produce virtually identical results and have been extensively validated (e.g. Butland et al., 1980). TIGAR however, required less computing time on the problems considered and so these times have been used to measure the speed-up achieved by HEXNEC. The LMFBR benchmark problems model the MARK I core design of the SNR300 (a 300 MW(E) sodium cooled fast reactor) at the start of life. Problem 1 is a two-dimensional plan model of the upper core (UC) and contains a ring of inserted control rods. Problem 2 is similar, being the lower core (LC) with rods removed. Problem 3 is a three-dimensional model. The PFR problem is a twodimensional centre plane model of the UKAEA Prototype Fast Reactor core at equilibrium burn-up. Tables 1 and 2 compare the results of HEXNEC and
A hexagonal geometry nodal expansion method
7 .o O
.o
II o ~ .
t"q
_Eo
t~
£ .i z'~ () O
,--4
[-
r-~ t-q O
O
t~--..E o II ~
117
J. M. PUTNEY
118
Table 2. Calculations on the PFR problem
Code
No. of nodes per
CPU time
S/A a
(s)b
HEXNEC TIGAR TIGAR
1 1 6
Error in keff
10.2 5.2 19.6
0.997185 1.007956 0.998354
Max-mean errors in average S/A values (%)d
kef f (%)¢
Powers
Damage rates
0.11 1.19 0.23
2.50 0.30 4.26 1.38 0.76 0.28
2.17 0.70 16.89 5.32 2.40 0.89
a 1 node per S/A--hexagonal geometry (across fiats = 14.59 cm), 6 nodes per S/A--triangular geometry (height = 7.30 cm). b As for Table 1. c Reference keff (=0.996077) obtained by extrapolation of SNAP kcff values to zero mesh interval. d Reference fluxes obtained by SNAP using a mesh of 96 triangles per S/A (height= 1.82 cm).
TIGAR calculations for the two-dimensional problems. In the case of the group fluxes, the errors shown are deviations expressed as a percentage of the average flux in the group (rather than the local value of the flux). A similar definition applies to the power and damage rate errors. Typical accuracy requirements for fuel management studies are around 0.25% in keff and 4% in subassembly (S/A) powers and damage rates. Evidently, for the problems considered, accuracies of this order are achieved by HEXNEC (zeroth order NEM) on the S/A mesh and TIGAR on a mesh of six triangles per S/A. The HEXNEC calculations however, require about half as much computing time. In Table 3, k~ff values computed for the benchmark models have been used to calculate the reactivity worth of the ring of control rods in the SNR 300 model. In this case, the accuracy of the worth computed by HEXNEC is comparable to that computed by the FD method using a mesh of 24 triangles per S/A. It is estimated that the speed-up achieved by HEXNEC over TIGAR for a calculation of control rod reactivity worth to within 1% is around an order of magnitude. Table 4 presents the results of HEXNEC calculations made on the three-dimensional benchmark Table 3. Control rod worth calculations on the benchmark problem Code
No. of nodes per S/A a
Worth
Error in worth (%)b
HEXNEC SNAP SNAP SNAP
1 1 6 24
0.072684 0.062919 0.071461 0.072721
0.7 14.0 2.3 0.6
a 1 and 6 nodes per S/A--as for Table 1. 24 nodes per S/A--triangular geometry (height= 2.8 cm). Reference worth =0.073177 (see note (c) of Table 1).
problem using six axial mesh intervals, or planes (one for each breeder region and four for the core region). The various methods considered correspond to different forms of axial treatment: W ° denotes a zeroth order expansion, and WnL k denotes an nth (= 1, 2) order expansion based on Galerkin (W= G)/moments (W= M) weighting and a kth ( = 0, 2) degree transverse leakage approximation. It should be noted that these results include the radial (hexagonal) mesh error. Although the radial and axial errors are not strictly separable, the radial error for keff has been estimated to be of the order of 0.13% by finding the minimum error achievable by refining the axial mesh. For this problem Table 4 shows that the first order expansion gives a significant improvement over the zeroth order method, while the further improvement on going to second order is smaller. The quadratic leakage method is much more accurate than the flat leakage approximation, and the use of moments weighting is marginally better than Galerkin weighting. Table 5 shows the fastest HEXNEC calculation for each variant of NEM which yielded a solution of a specified accuracy, i.e. 0.3% for keff and 3% in S/A powers and damage rates. (The accuracy requirement on S/A powers and damage rates has been reduced to account for the rather coarse reference solution.) Also shown is the fastest FD calculation which produced such a solution. According to these results, the first order methods with the quadratic expansion for transverse leakage achieve the greatest speed-up although the second order methods are more accurate and are only slightly more time-consuming (even though the axial mesh is the same in this case). 6. CONCLUSIONS A hexagonal geometry extension of the nodal expansion method (Finnemann et al., 1977) has been
A hexagonal geometry nodal expansion method
119
.E tr~ ¢-q
t-
O
z.
.o d~
t-~
E E em
9
.o
z.i 2~ i o
.o
9 I. ¢-q
z
.-2
k~
-r ~f r~
o
c-i r...) ' ~ , ~
()
120
J . M . PUTNEY
0
I"
)
0 m
¢ < a
e-
£
¢e £
z~ •" 0
o
7. ¢N C
o -.=
e" e~ 0
LO
-~t:o '=t~
.8
0~
¢qr~t'qcqcq~CqCqO0~
C
0
e~
t..
¢~
'= o - ~
o •. 0
LL ,~
o
o
[..
tr~
t'~
0
~[-[-
0
0.0 0
.o
A hexagonal geometry nodal expansion method obtained by viewing the method as the determination of a multidimensional flux expansion in each node, and constructing an expansion with hexagonal symmetry. The axial direction is treated using a conventional one-dimensional expansion. Effective solution and acceleration techniques have been developed for multigroup fast reactor calculations and implemented in a computer code HEXNEC. The efficiency of the code has been demonstrated for four-group twodimensional and three-dimensional LMFBR benchmark problems and a 37-group two-dimensional PFR problem. A plan mesh of one node per subassembly is used. On the two-dimensional problems, the zeroth order variant of HEXNEC (quadratic flux expansion) produced solutions within the chosen target accuracy. Similar accuracies were provided by the finite difference code TIGAR using 6 triangles per subassembly, but the calculations required about twice as much computing time. The fastest TIGAR calculations which provided a solution of sufficient accuracy to the three-dimensional benchmark problem needed a radial mesh of 6 triangles per subassembly and 36 axial mesh intervals. Similar accuracies were obtained by HEXNEC with 12 axial mesh intervals when a zeroth order axial expansion was used and with six mesh intervals when a first order (cubic degree) axial expansion was used. These calculations were faster than the TIGAR calculation by factors of 8.7 and 10.6 respectively. In the latter case, the calculation was only a factor of 7 slower than the TIGAR calculation on the two-dimensional problem. The insertion of absorber rods at different depths could force the adoption of a finer axial mesh at some locations. The performance of HEXNEC in such cases could be improved by allowing a different order of axial expansion in each
121
mesh plane (or even each node). This would be a relatively simple extension.
REFERENCES
Arkuszewski J. J. (1982) SIXTUS-2 A two-dimensional multigroup diffusion theory code in hexagonal geometry. EIR-470. Buckel G., Kfifner K. and Stehle B. (1977) Nucl. Sci. Engng 64, 75-89. Butland A. T. D., Putney J. M. and Sweet D. W. (1980) Comparisons of neutron diffusiontheory codes in two and three space dimensions using a sodium cooled fast reactor benchmark. AEEW R1346, UKAEA. Doming J. (1979) CONF-790402, 3-1-3-31. Duracz T. (1981) Proceedings of lnternational Topical Meeting on Advances in Mathematical Methods for the Solution of Nuclear Engineering Problems, Vol. 1, pp. 423~,29, Munich. Finnemann H., Bennewitz F. and Wagner M. R. (1977) Atomkernenergie 30, 123 128. Lawrence R. D. and Doming J. J. (1980)Nucl. Sci., Engng 76, 218 231. Lawrence R. D. (1983) The DIF3D nodal neutronics option for 2- and 3-dimensional diffusion theory calculations in hexagonal geometry. ANL-83-1, Argonne National Laboratory. Matthews J. D. (1972) A speed and storage optimised 3D diffusion code operating with the COSMOS standard interface. UKAEA Internal Document. McCallien C. W. J. (1975) SNAP-3D: A three-dimensional neutron diffusion code. TRG Report 2677 (R), UKAEA. Putney J. M. (1984) Nodal methods for solving the diffusion equation for fast reactor analysis. Ph.D. Thesis, University of London. Varga R. S. (1962) Matrix lterative Analysis. Prentice-Hall, Englewood Cliffs, NJ. Wagner M. R. (1968) GAUGE--A two-dimensional few group neutron diffusion-depletionprogram for a uniform triangular mesh. GA 8307, Gulf General Atomic Company.
Progress in Nuclear Energy, Vol. 18, No. 1/2, pp. 113 121, 1986. Printed in Great Britain, All rights reserved.
A H E X A G O N A L G E O M E T R Y N O D A L E X P A N S I O N M E T H O D FOR FAST REACTOR C A L C U L A T I O N S J. M. PUTNEY* CEGB, Central Electricity Research Laboratories, Kelvin Avenue, Leatherhead, Surrey, U.K.
Abstract--A hexagonal geometry extension of the nodal expansion diffusion theory method is presented together with solution and acceleration techniques which have proved effective for two-dimensional (2D) and threedimensional (3D) fast reactor calculations. The method uses a radial nodal flux expansion which has hexagonal symmetry. A one-dimensional (1D) expansion is used for the axial direction. Numerical results are presented for four-group LMFBR benchmark problems and a 37-group PFR problem. A zeroth order (quadratic degree) radial expansion is found to be sufficient to permit calculations to be performed using one mesh per subassembly, whereas higher order expansions are more efficient in the axial direction because a coarser mesh can usually be used in this direction. For comparable accuracy the method is found to be a factor of 2 faster than the finite difference method for two-dimensional calculations and an order of magnitude faster for three-dimensional calculations•
1. INTRODUCTION In recent years the advent of nodal diffusion methods (Doming, 1979) has enabled computing times for LWR analysis to be significantly reduced. The nodal expansion method (NEM) (Finnemann et al., 1977) has proved particularly effective. Its basis is the solution of equations for the average flux on a coarse spatial mesh, combined with a high order local polynomial flux expansion, which is used to determine a relationship between the mean interface partial currents and the mesh average fluxes. Being an LWR method however, the N E M is formulated in rectangular geometry. For fast reactor analysis, extension to hexagonal geometry is required, In addition, suitable solution and acceleration techniques must be developed for general fast reactor multigroup problems. This paper describes a hexagonal geometry N E M and associated solution techniques which have proved effective for two-dimensional and three-dimensional fast reactor calculations. Numerical results are included which demonstrate the efficiency of the method over the finite difference method.
energy group g integrated over the volume of (homogeneous) node FI~ (the nodal balance equation) can be written in the form:
~w Aw's m m_ E ~ s:;,oT, ~ ~i÷m ~,.ws-j-m~ ~w~,+z.,
where 0~' is the average node flux (nodal flux) and Jowl +m and j g ~ are the average outgoing and incoming partial currents (interface currents) at the left (s=l) or right (s = r) node surface in the w-direction (surface F~). The interface currents at surface F ~ are related to the local flux gradient by Fick's law
Jo,.., -Y0w~ = - c~
•
.
2.1. General background The multigroup neutron continuity equation for
fl, s=r
as=].
(2)
1, s=l
In rectangular geometry, by reversing the order of integration and differentiation, equation (2) may be written as
dw
r,%
(3)
where ~,~w(W)is the one-dimensional nodal flux formed by averaging the flux in node FIm over the direction transverse to the w-direction. This flux satisfies a onedimensional nodal diffusion equation 2
* This work was undertaken while the author was a Research Student at Imperial College London, and was carried out at AEE Winfrith.
dFm,
Aw~J~,
Jg+~'~-J°-w;= - f "D~ 2. TWO-DIMENSIONAL NODAL EQUATIONS
i1)
g'=l
-D~
m
d ~baw
L"
G
dw 2 +E,mo~lgmw-I- gw = E Sg~'~am'w, (4) g'=l
where L~w(w) is the transverse leakage. 113
114
J. M. PUTNEY
The essence of the NEM is the determination of a polynomial expansion for ~bg%. In the zeroth order variant, a quadratic expansion is assumed and fitted to the average value of ~"~ in the node (¢b~")and its values m m m m at Fwt and F~, (denoted by ffgwt and ~b0~,). Introducing the Fick's corollary
Oo.~- 2(Jo.~ +3a~ ) m
--
• +m
• -m
(5)
allows the expansion coefficients to be expressed solely in terms of O~' and (jo+~+J0-w~)' Substituting the result into equation (3) then leads to a set of interface current equations which, together with the nodal balance equation (1) and the continuity and boundary conditions, is sufficient to determine all the nodal fluxes and interface currents. In the nth order variant (n_> 1), an n + 2th degree expansion is constructed by adding on higher order polynomial terms. The coefficients of these are determined from local weighted residual solutions of the one-dimensional diffusion equations, a process which requires the introduction of a polynomial approximation to the transverse leakage. The principal difficulty in applying the NEM to a hexagonal geometry, is that the process of changing the order of integration and differentiation in equation (3) introduces the flux values at the corner points of F~. However, the one-dimensional rectangular geometry expansions are equivalent to a general twodimensional nodal flux expansion. In the NEM the cross-terms in this expansion either drop out completely or appear only in weighted integrals of the transverse leakage• The remaining terms are simply the one-dimensional nodal flux expansions. With a two-dimensional expansion the interface current equations can be derived directly from equation (2). In hexagonal geometry, this avoids the need to estimate corner point fluxes. The problem is thus to find a suitable expansion for the flux in a hexagonal node.
GG.~..L)=Ag'ho+
Y { ,whgG)+b,~hdG)} a
m
m
w=x,u,v
+ coh l(~x)h 1(~u)h l(~v)'
(6)
where hi(~w) is a polynomial of degree i in Cwand {ho, hi (~0, h2 (¢~)} is orthogonal on integration over the volume of the node. A suitable set of expansion functions is h o = 1, h(~) = ~w, h2 (~w)= ~
5 72.
(7)
Introducing these into equation (6) and transforming into local rectangular cartesian variables, leads to a nodal expansion having seven independent terms. The coefficients are obtained by fitting to the node and surface average fluxes. Using the nodal expansion to evaluate equation (2) for each surface and eliminating the surface fluxes using equation (5), leads to six relations of the form: (D~o~+n7"].+m
T/ox.
7 +m
(.
o
=
)'29 +,
+
29._ m
Jow,;
7
T d/ . ; . 11._,.)
- 10j.w,
ll.+m~
,vO0s.w.+ w=u,v 2 36
5-* 7.
(8)
Inverting the resulting 6 × 6 system to express the outgoing currents in terms of the incoming currents, yields a set of interface current equations of the form: • +m
laxr
m
=
CffI
* -m
rn
. --m
Jffxr + CozJoxl
m m. t + C ~.ms , , . 3Jgwr+Co4Jffwl
+ E
(9)
w=u,v m Here, the nodal coefficients Co~ - C o ms are formally rational functions of D~"/H. In practice, they can be determined efficiently using a numerical inversion procedure based on Cholesky's decomposition• Equations (9) can be used to eliminate the outgoing currents from the nodal balance equation (1):
2.2. A zeroth order hexagonal geometry formulation It is convenient to introduce a set of'hexagonal axes' (x, u, v) where axes u and v are at - 120° and + 120° to the x axis. In this system the co-ordinates of a point P1 are (xl, ul, v0, where wl(=xl, ul, vl) is the distance along the w-axis to the perpendicular from P~. It is also convenient to define the local dimensionless variables, ~ = w / H , where H is the across flats width of the hexagon. For the interface current equations to have hexagonal symmetry the expansion should be symmetric in the variables x, u, v. The zeroth order method presented here is based on the following choice:
0'=1
2
+ 3H{ 1 - C ~ - C ; 3 - 2 c o m 3 - 2 C ; ~ }
Z w=x,u,v s=l,r
Jowl" --
(10)
Equations (1)/(10) and (9) represent zeroth order hexagonal geometry NEM equations, expressed in a form analogous to their application in rectangular geometry (Finnemann et al., 1977). Together with the continuity of interface currents and boundary conditions they are sufficient to determine the nodal fluxes
115
A hexagonal geometry nodal expansion method and interface currents. In contrast to rectangular geometry though, all coefficients in the expansion are determined and the problem has not been decoupled into a series of one-dimensional problems. Although the method can be generalized to higher orders of expansion (Putney, 1984) this has been found to be unnecessary for the hexagonal mesh required to represent fast reactor subassemblies. Other hexagonal geometry nodal expansion methods have been developed, for example by Duracz (1981), Lawrence (1983) and Arkuszewski (1982). Duracz and Lawrence treat the one-dimensional transversely average flux in each of the three hexagonal geometry directions, x, u, v. Half-range quadratic polynomial expansions are used to fit the averaged flux and the two surface fluxes in a way which treats the variation in transverse width of the node. Arkuszewski expresses the flux in the node in terms of six fundamental solutions of the diffusion equation inside the node. Using these the outgoing currents are expressed in terms of the incoming currents.
The zeroth order expansion coefficients are determined by fitting to the node and surface average fluxes: m
m
4
~.(~. ~, ~0=4,~(~. ~)+ Y~ d.~,g,(~O (11) i=1
where ~, = z/h " is a dimensionless axial variable for the node. Here ~b~2~ is the two-dimensional hexagonal expansion and ,q~(~=)is a polynomial of degree i in ~= such that
m
--
m
m
m
To determine the higher order expansion coefficients a local weighted residual equation is formed by multiplying the diffusion equation by a weighting function W k ( z ) ( k = l , 2) and integrating over the volume of the node. Two weighting schemes have been considered: Galerkin ( W k = gk + 2) and moments ( W k = g0. After integrating over the radial directions, the diffusion equation in the residual reduces to the onedimensional equation (4) (with w--z). The transverse leakage is expanded as m m m L,=(~=)= L,~ o + Lax i O,(~=)+ L,~2 92(~=)
(14)
and fitted to the average transverse leakages (functions of the radial interface currents) in the node and its two axial neighbours (Lawrence and Dorning, 1980). For axial boundary nodes the flux boundary condition is applied to Lg=(~z).The weighted residual equations for m m the higher order coefficients, dg~3 and do= 4 , are
}
3. THREE-DIMENSIONAL NODAL EQUATIONS The derivation of a three-dimensional hexagonal geometry expansion requires the combination of the two-dimensional expansion with a one-dimensional rectangular geometry expansion in the axial direction. Cross terms involving hexagonal and rectangular expansion functions need not be included explicitly, as due to the orthogonality of the directions they behave exactly as for the rectangular geometry method. Since the axial variations in composition in fast reactors generally permit a relatively coarse axial mesh to be used, it is appropriate to consider higher order axial expansions. The method presented here is based on second order (i.e. fourth degree) axial accuracy. Following Finnemann et al. (1977) for the axial treatment, the following nodal expansion is assumed:
m
dozl = @,a=,- d/g=, dog 2 - 3(q%=,- 2Og + q/a~,)- (13)
J'D•' "do:, )h,.2 + Ak Y~,g
°
2
= 0'=1
" (1)0'k m
Sgg,
q- O k zrmgdgmk - BkL,S~k
(15)
where k = i - 2 and the weighted flux (I),"k = (A k do==~ Bkdo~k). For Galerkin weighting A I = ~ , A2=~o, Bm= ~, B2 = 1½. For moments weighting A I = ~o, A 2 _-
1 ~16 ,
BI =16, B 2 = 11.
(16)
For the nodal expansion (11) the radial interface current equations are identical to their 2D counterparts. At the axial surfaces, the equations take the form .+ra m m - m ra . - " m m Jg~l = Ag= q)gm + Bg~jg~t + Cg=jg=r + fitEg= dg=3
(17)
+ F~m~dgmz4
where A o ~ - F ~ are rational functions of D ~ ' / h " obtained analytically by inverting the 2 × 2 matrix. Eliminating the outgoing currents from equation (1) gives:
c,"; + ~ -
+ z 7; * ;" = Z s.";,• ~; g'=l
{1, gl(~=), g2(~=)} and {l, g3(~z), ga(~z)} are orthogonal on integration over the volume of node l-I', and g3(~) and g4(~,) integrate to zero over its axial surfaces. A suitable set of expansion functions is •q l ( ~ z )
= ~z'
g2(~z) 2
1
=~2-
= ¢ z ( ~ z -- 4),
---m
w=x,u,v s=l,r
Jgws
1
12'1 g 3 ( ~ z )
2 1 2 O4(~z) = (¢z -- 4)(~z -- ~)"
2 + 3H {1-Co~ - C Z - 2Co"3 - 2C0=4 }
+ h" {(1 - B ~ - Co"=) (Jo-~' +Jo-=r") -- 2F~=~dg~,, } (12)
(18)
116
J.M. PUTNEY 4. SOLUTION TECHNIQUES
The nodal equations presented above have been implemented in a computer code HEXNEC. The solution procedure reduces the nodal balance equation (1) to a system of nonlinear finite difference type equations by defining leakage coefficients +m
• +m
ra
~gw~=Jgw~/~g
(19)
which are used to eliminate the interface currents in terms of nodal fluxes: • +m
+m
m
Jgws = ~:gw~~ g
(20)
with the incoming currents being expressed in terms of the leakage coefficients and average fluxes for the neighbouring nodes. The resulting 'pseudo finite difference equations' are solved using a conventional fission source iteration scheme in parallel with an (outer) iterative update of the leakage coefficients based on the other nodal equations. The following update procedure is used: Using one Gauss-Seidel sweep, for each node (i) For the three-dimensional problems with higher order axial treatment compute the transverse leakage moments from the latest values of the interface currents and solve the weighted residual equations (15) for the higher order expansion coefficients. (ii) Use the outgoing current eliminated nodal balance equation (10)/(18) to update the nodal flUX.
(iii) Use the interface current equations (9) and (17) to update the outgoing currents and apply the continuity and boundary conditions to update the corresponding incoming currents. When solving the weighted residual equations (15) the weighted sources are computed in the same way as the normal sources. The zeroth order expansion coefficients are calculated using the latest values of the nodal fluxes and the interface currents. As well as being used to update the nodal fluxes and interface currents, the resulting higher order expansion coefficients are combined with the zeroth order coefficients to form new weighted fluxes which are then stored for future use in calculating weighted sources. This differs from the approach outlined by Finnemann et al. (1977) for the rectangular geometry method, in which the expansion coefficients are re-calculated when needed. The adopted approach involves less arithmetic and the extra storage required is acceptable as the weighted fluxes are only present in one direction. The outer iterations in HEXNEC are accelerated by asymptotic extrapolation (AX) (Wagner, 1968) of the
fission source, followed by renormalization of the nodal fluxes and currents. This reduces the number of outer iterations by 25-40%. Coarse mesh rebalancing may be a more efficient acceleration technique, but this is dependent on choosing a coarse-mesh and rebalancing strategy appropriate for the application. The inner iterations use a line Gauss-Seidel method (Varga, 1962) accelerated by AX. For two-dimensional problems, the number of unaccelerated iterations required is not much greater than the 'threshold' number required before AX is permitted (and sometimes less). Consequently, the reduction in the number of inner iterations per outer is around 10% or less. On threedimensional problems however, a reduction of around 25% can be achieved. Over-relaxation was found to be less effective than AX. One some problems HEXNEC has been found to produce some negative nodal fluxes and/or interface currents. Although these generally only occur at unimportant locations, they can lead to divergence of the inner iterations. Consequently, the definition of the leakage coefficients was modified to ensure that they are always non-negative. If an interface current becomes negative then, in calculating the leakage coefficients, it is interpreted as positive current in the opposite direction. If a nodal flux becomes negative, it is multiplied by - 1. The inner iterations then yield the negative of the flux which is multiplied by - 1 before proceeding to the update stage.
5. NUMERICAL STUDIES HEXNEC calculations for three four-group LMFBR benchmark problems (Buckel et al., 1977) and a 37-group (7-downscatter) PFR problem have been compared with results obtained used the finite difference (FD) codes (TIGAR) (Matthews, 1972) and SNAP (McCallien, 1975). Both these F D codes employ mesh-centred differencing, produce virtually identical results and have been extensively validated (e.g. Butland et al., 1980). TIGAR however, required less computing time on the problems considered and so these times have been used to measure the speed-up achieved by HEXNEC. The LMFBR benchmark problems model the MARK I core design of the SNR300 (a 300 MW(E) sodium cooled fast reactor) at the start of life. Problem 1 is a two-dimensional plan model of the upper core (UC) and contains a ring of inserted control rods. Problem 2 is similar, being the lower core (LC) with rods removed. Problem 3 is a three-dimensional model. The PFR problem is a twodimensional centre plane model of the UKAEA Prototype Fast Reactor core at equilibrium burn-up. Tables 1 and 2 compare the results of HEXNEC and
A hexagonal geometry nodal expansion method
7 .o O
.o
II o ~ .
t"q
_Eo
t~
£ .i z'~ () O
,--4
[-
r-~ t-q O
O
t~--..E o II ~
117
J. M. PUTNEY
118
Table 2. Calculations on the PFR problem
Code
No. of nodes per
CPU time
S/A a
(s)b
HEXNEC TIGAR TIGAR
1 1 6
Error in keff
10.2 5.2 19.6
0.997185 1.007956 0.998354
Max-mean errors in average S/A values (%)d
kef f (%)¢
Powers
Damage rates
0.11 1.19 0.23
2.50 0.30 4.26 1.38 0.76 0.28
2.17 0.70 16.89 5.32 2.40 0.89
a 1 node per S/A--hexagonal geometry (across fiats = 14.59 cm), 6 nodes per S/A--triangular geometry (height = 7.30 cm). b As for Table 1. c Reference keff (=0.996077) obtained by extrapolation of SNAP kcff values to zero mesh interval. d Reference fluxes obtained by SNAP using a mesh of 96 triangles per S/A (height= 1.82 cm).
TIGAR calculations for the two-dimensional problems. In the case of the group fluxes, the errors shown are deviations expressed as a percentage of the average flux in the group (rather than the local value of the flux). A similar definition applies to the power and damage rate errors. Typical accuracy requirements for fuel management studies are around 0.25% in keff and 4% in subassembly (S/A) powers and damage rates. Evidently, for the problems considered, accuracies of this order are achieved by HEXNEC (zeroth order NEM) on the S/A mesh and TIGAR on a mesh of six triangles per S/A. The HEXNEC calculations however, require about half as much computing time. In Table 3, k~ff values computed for the benchmark models have been used to calculate the reactivity worth of the ring of control rods in the SNR 300 model. In this case, the accuracy of the worth computed by HEXNEC is comparable to that computed by the FD method using a mesh of 24 triangles per S/A. It is estimated that the speed-up achieved by HEXNEC over TIGAR for a calculation of control rod reactivity worth to within 1% is around an order of magnitude. Table 4 presents the results of HEXNEC calculations made on the three-dimensional benchmark Table 3. Control rod worth calculations on the benchmark problem Code
No. of nodes per S/A a
Worth
Error in worth (%)b
HEXNEC SNAP SNAP SNAP
1 1 6 24
0.072684 0.062919 0.071461 0.072721
0.7 14.0 2.3 0.6
a 1 and 6 nodes per S/A--as for Table 1. 24 nodes per S/A--triangular geometry (height= 2.8 cm). Reference worth =0.073177 (see note (c) of Table 1).
problem using six axial mesh intervals, or planes (one for each breeder region and four for the core region). The various methods considered correspond to different forms of axial treatment: W ° denotes a zeroth order expansion, and WnL k denotes an nth (= 1, 2) order expansion based on Galerkin (W= G)/moments (W= M) weighting and a kth ( = 0, 2) degree transverse leakage approximation. It should be noted that these results include the radial (hexagonal) mesh error. Although the radial and axial errors are not strictly separable, the radial error for keff has been estimated to be of the order of 0.13% by finding the minimum error achievable by refining the axial mesh. For this problem Table 4 shows that the first order expansion gives a significant improvement over the zeroth order method, while the further improvement on going to second order is smaller. The quadratic leakage method is much more accurate than the flat leakage approximation, and the use of moments weighting is marginally better than Galerkin weighting. Table 5 shows the fastest HEXNEC calculation for each variant of NEM which yielded a solution of a specified accuracy, i.e. 0.3% for keff and 3% in S/A powers and damage rates. (The accuracy requirement on S/A powers and damage rates has been reduced to account for the rather coarse reference solution.) Also shown is the fastest FD calculation which produced such a solution. According to these results, the first order methods with the quadratic expansion for transverse leakage achieve the greatest speed-up although the second order methods are more accurate and are only slightly more time-consuming (even though the axial mesh is the same in this case). 6. CONCLUSIONS A hexagonal geometry extension of the nodal expansion method (Finnemann et al., 1977) has been
A hexagonal geometry nodal expansion method
119
.E tr~ ¢-q
t-
O
z.
.o d~
t-~
E E em
9
.o
z.i 2~ i o
.o
9 I. ¢-q
z
.-2
k~
-r ~f r~
o
c-i r...) ' ~ , ~
()
120
J . M . PUTNEY
0
I"
)
0 m
¢ < a
e-
£
¢e £
z~ •" 0
o
7. ¢N C
o -.=
e" e~ 0
LO
-~t:o '=t~
.8
0~
¢qr~t'qcqcq~CqCqO0~
C
0
e~
t..
¢~
'= o - ~
o •. 0
LL ,~
o
o
[..
tr~
t'~
0
~[-[-
0
0.0 0
.o
A hexagonal geometry nodal expansion method obtained by viewing the method as the determination of a multidimensional flux expansion in each node, and constructing an expansion with hexagonal symmetry. The axial direction is treated using a conventional one-dimensional expansion. Effective solution and acceleration techniques have been developed for multigroup fast reactor calculations and implemented in a computer code HEXNEC. The efficiency of the code has been demonstrated for four-group twodimensional and three-dimensional LMFBR benchmark problems and a 37-group two-dimensional PFR problem. A plan mesh of one node per subassembly is used. On the two-dimensional problems, the zeroth order variant of HEXNEC (quadratic flux expansion) produced solutions within the chosen target accuracy. Similar accuracies were provided by the finite difference code TIGAR using 6 triangles per subassembly, but the calculations required about twice as much computing time. The fastest TIGAR calculations which provided a solution of sufficient accuracy to the three-dimensional benchmark problem needed a radial mesh of 6 triangles per subassembly and 36 axial mesh intervals. Similar accuracies were obtained by HEXNEC with 12 axial mesh intervals when a zeroth order axial expansion was used and with six mesh intervals when a first order (cubic degree) axial expansion was used. These calculations were faster than the TIGAR calculation by factors of 8.7 and 10.6 respectively. In the latter case, the calculation was only a factor of 7 slower than the TIGAR calculation on the two-dimensional problem. The insertion of absorber rods at different depths could force the adoption of a finer axial mesh at some locations. The performance of HEXNEC in such cases could be improved by allowing a different order of axial expansion in each
121
mesh plane (or even each node). This would be a relatively simple extension.
REFERENCES
Arkuszewski J. J. (1982) SIXTUS-2 A two-dimensional multigroup diffusion theory code in hexagonal geometry. EIR-470. Buckel G., Kfifner K. and Stehle B. (1977) Nucl. Sci. Engng 64, 75-89. Butland A. T. D., Putney J. M. and Sweet D. W. (1980) Comparisons of neutron diffusiontheory codes in two and three space dimensions using a sodium cooled fast reactor benchmark. AEEW R1346, UKAEA. Doming J. (1979) CONF-790402, 3-1-3-31. Duracz T. (1981) Proceedings of lnternational Topical Meeting on Advances in Mathematical Methods for the Solution of Nuclear Engineering Problems, Vol. 1, pp. 423~,29, Munich. Finnemann H., Bennewitz F. and Wagner M. R. (1977) Atomkernenergie 30, 123 128. Lawrence R. D. and Doming J. J. (1980)Nucl. Sci., Engng 76, 218 231. Lawrence R. D. (1983) The DIF3D nodal neutronics option for 2- and 3-dimensional diffusion theory calculations in hexagonal geometry. ANL-83-1, Argonne National Laboratory. Matthews J. D. (1972) A speed and storage optimised 3D diffusion code operating with the COSMOS standard interface. UKAEA Internal Document. McCallien C. W. J. (1975) SNAP-3D: A three-dimensional neutron diffusion code. TRG Report 2677 (R), UKAEA. Putney J. M. (1984) Nodal methods for solving the diffusion equation for fast reactor analysis. Ph.D. Thesis, University of London. Varga R. S. (1962) Matrix lterative Analysis. Prentice-Hall, Englewood Cliffs, NJ. Wagner M. R. (1968) GAUGE--A two-dimensional few group neutron diffusion-depletionprogram for a uniform triangular mesh. GA 8307, Gulf General Atomic Company.