- Email: [email protected]
A well-balanced coarse-mesh flux expansion method
Ann. Nucl. Eneryy, Vol. 21, No. 1, pp. 45-53, 1994 Printed in Great Britain. All rights reserved
A WELL-BALANCED
0306-4549/94 $6.00+0.00 Copyright © 1993 Pergamon Press Ltd
COARSE-MESH
FLUX EXPANSION
METHOD
B. MONTAGNINft, P. SORAPERRA*, C . TRENTAVIZI* Dipartimento di Costruzioni Meccaniche e Nucleari, Universit~ di Pisa, via Diotisalvi 2, 1-56016 Pisa, Italy M. StIMINI, D.M. ZARDINI* Laboratori di Ingegneria Nucleate di Montecuccolino, Universit~ di Bologna, via dei Colli 16, 1-40136 Bologna, Italy (Received 3 July 1993) A b s t r a c t - An efficient low-order coarse-mesh method for LWR criticality problems, which allows for fast routine computations as well as for more accurate, but still unexpensive, reference solutions, is presented. The efficiency of the method rests upon a very simple nodal flux expansion formula and an improved scheme for outer iterations, coupled with a fast ADI technique.
1. I N T R O D U C T I O N
NOMENCLATURE a
#
=
¢g(r)
Neutron diffusion coarse-mesh (CM) methods are usually grouped into several classes, according to the nature of the underlying approximation technique. A more pragmatic classification, consisting of only two classes, refers to their domain of application. The first class concerns the rather rough (nodal power max error < 7%0, mean error < 2.5%, say) but very fast CM techniques used in routine computations. The second class is constituted by the much more accurate (max error < 1.5%, or so) and theoretically refined CM methods which are used for reference purposes, mostly to assess the reliability of the methods of the first class. In this paper we wish to illustrate a wellbalanced CM method for steady-state (criticality) problems, which allows for both kinds of calculations (with a proper choice of the m e s h size), within approximately the error limits indicated above. No attempt has been made in order to lower such limits, since they seem to be acceptable (for the appropriate class of problems), if allowance is made for the uncertainties of the physical parameters. Care has been taken, instead, to achieve a good computational efficiency.The present work has, however, a further, and perhaps deeper, motivation. The success of the CM methods based on 'analytical' formulations or, more precisely, on the representation of the 'transverse-integrated' flux by means of 1-D solutions of the diffusion equation (Shober et hi., 1977, Bonalumi et aL, 1978, Lawrence and Dorning, 1980, Ougouag and Rajic, 1988, Rajic and Ougouag, 1989, and others) has a little obscured the family of the CM methods based on polynomial approximations. The latter methods have, however, their own merits. They do not require the lengthy computation of many transcendental functions
= Total number of neutron energy
=
groups Neutron flux in group g (at the point r) (n/cm 2 s)
Dg = Dg(r) = Diffusion coefficient in group g (cm) 2 ; / = 2~rg(r) = Macroscopic removal cross-section in g r o u p g (cm -1)
xsg'->g =
27sg'->g(r) = Macroscopic transfer cross-section from g r o u p g ' to group g (cm -l)
v~g = v ~ g ( r )
= Macroscopic fission cross-section
times the average number of secondaries, in group g (cm -I)
z g = Total fission spectrum fraction in group g k = Effective multiplication factor t To whom all the correspondence should be addressed. * W o r k p e r f o r m e d as a partial fulfilment for a M.Sc.Thesis.
45
B. MONTAGNINI et al.
46
(this may be an important advantage, e.g. if the scheme is a candidate for an extension to time-dependent problems) and their structure is simple, at least if the degree of the polynomials is not too high.The relative importmlce of the various terms of the approximation formula can be easily investigated and it may be possible, by suppressing some of them, to obtain a consistent simplification of the scheme without compromising the accuracy in an essential way. One of our aims was just to show that, for a definite class of diffusion problems as those concerning the LWR cores, a quite simple polynomial formulation can be more efficient and not much less accurate than other schemes (analytical or not). The method is based on a low-order expansion of the nodal flux, namely a sum of cubic polynomials in each one of the space variables x,y,z, as in the CUBBOX method ( ' s u m ' expansion formula, Langenbuch et al., 1977a).The solution procedure has undergone, however, several modifications with respect to CUBBOX.
2. T H E FLUX E X P A N S I O N F O R M U L A The time-independent multigroup neutron diffusion equations are (see the Nomenclature): G
V.D~V¢
~-
- Z~ + Z,,)~ +kl-F~ = 0
(V.DV
Having m mind a solution procedure based, as customary,
on some outer/inner iterative process, we first write down the outer iteration scheme (i.e. the iterative process involving the fission source): -(V.DV +Z'
o~ F i ~ ( " + l ) = I - ~ F ~ ( " ) ' k(")
(3a)
k(")
where E = - T-r+ ~'~s (thus,
fee'
=
.
~,rg~gg ' + .~$g">g
with ~gg" = 0 Of 1 according to g a g ' o r g = g ' ), n is the iteration index and tx, 0_
0~I//(n+l)+ (I - OQI//(n) (g = l ..... ~ )
(I) These equations are to be solved subject to the following homogeneous conditions at the reactor boundary OV :
(2)
where b/Ou denotes the outward normal derivative and (r g, z g are real constants. It is convenient to rewrite eqs.(1) and (2) in matrix form. Let D - - c o l [4) l ..... ¢G] and 1
G
'~r = diag [Zr ..... Zr ] T = diag ['r I ,... z G]
0
~)~I
...
~g~l
Z~ I~2
0
...
Z ~2
where
g=l V being the reactor volume.The denominator at the r.h.s. of eq.(3b) thus coincides with the total fission source (depending on a linear combination of both fluxes ¢I~n) and
(g = 1..... G)
OU
S = diag [o"1..... O'G],
I
k ('+I) = k (')
g'=l G + ~ - z ~ Y .-- v z : % ~ ' = 0 , g'=l
D = diag [D 1..... D G],
(23
s o + T ~---~= 0
X:¢g+~ r.:'-~ ¢~;
(Tgog+r gong =0,
(I')
Zs =
¢b("+t)) which sustains the flux in eq.(3a). Since the r.h.s. of eq.(3a), i.e. st") = I - t~ F~(n) k (")
(4)
is known (when the n-th step is completed) we have actually to solve an inhomogeneous system of partial differential equations. We assume that the reactor domain cau be represented as a regular array of rectangular parallelepipedes VO.t ('cells', or 'nodes'), with centres (xi,Yi, Zt) and edges of length hxi, h W, hz/, and that all the nuclear parameters are spatially cohstant over each Vqt. [laving introduced local dimensionless variables
= (x - xi)/hxi, 77 = (y - yj)/hyj, ( = (z - Zl)/hzl we approximate the g-group flux inside the node V//t by means of the following polynomial formula
...
0
~g(~,)?,Q = ~gl + axg~ + bxg(~ 2- -~2) + cxg(~- 4~ 3) F
+ aye 0 + byg(o 2 - 1-~2) + cyg(71 - 4r/3)
=
G
1
+ a z g ( + bzg(( 2- 1~) + czg( ( - 4(3) (these matrices depend, in general, on r). Then eqs.(1) and (2) can be written as follows:
(5) which coincides essentially with the 'sum' formula used in CUBBOX (Langenbuch et al., 1977a). It is written,
A well-balanced coarse-mesh flux expansion m e t h o d however, in such a way that its coefficients may be directly expressed in terms of the average node flux and the average interface fluxes (see below), rather than the corresponding centrepoint fluxes. We have, in fact, t
[ I ¢#(x,y,z) d V hxihyjhzt Jv =
d~'~)g(~,r/,~ ") :
d /2
where A = V.V and we = wg(x,y,z) are weight functions (eq.(4) has also been taken into account). Let us first take wg = 1 for all g. If the approximate expression for ¢g. eq.(5), is substituted into the integral over V/fl.written in terms of the variables ~,r/,(, we get, by performing the partial derivatives, integrating and taking account of eqs.(7) and (8):
$ l
6D~0-t\
t2
so that ¢g/fl actually represents the (approximate) average node flux. As it concerns the interface fluxes, let us define the following 'transversally averaged' flux:
fYJ'~Jdy Ft*~-hadz ~)g(x,y,z)
= ~
47
g g g " 2~/~1 +, ¢h~+tn,t 4" ~i:,/-112 - 2¢~/~1+ ~/1,1+112~ 4- ¢i~-I/2,1
hy2 -
ZG
h 2z, g,
g'
Qgd.Ol(POI= S~Ot ,
(g = 1..... G)
g'=t
= I,= dr/
= d~'¢g(~,r/,O
Then, using eq.(5), t~gx(~) = O~t + axg~ + bxg(~ 2- i=~2) + Cxg(~- 4~ 3) (6) Similar equations hold for
tpyg(r/) mad (~zg(O. The fluxes
0x~(:t:l/2) are, in particular, the average interface fluxes at the right (left) face of the cell and will be denoted henceforth by ~g~l/2dl. Putting ~ = :t: 1/2 into eq.(6), then adding and subtracting, we g e t , as in CUBBOX, the desired expression of the axg, bxg coefficients in terms of the average (volume or interface) fluxes:
ag =
bg
¢~+t/2jt- Cgi-m,it
(7)
= 3(~gi+ll2d/- 2~)g/fl+ Ogi-i/2jt)
(8)
Similar formulas hold for ayg, byg and azg, b.g. Let us postpone, for a moment, the problem of deterrr~ining the r e m a i n i n g c o e f f i c i e n t s Cxg, Cyg, Czg and try to approximately solve, by the weighted residual method, the inbomogeneous diffusion equations, eq.(3a). We denote by D q ff, ,Sr,Otg,..the constant values of Dg, ,Stg .... over the node Vijt. We also put
Qa, ql = ~"qt + ~Fqt
(9)
According to the weighted residual procedure, the i n t e g r a t e d f o r m o f eq.(3a), written in terms of components, or elements, of the involved vectors and matrices and suppressing the index n, is
o~.ql
] wg dV
sgwgdV,
(g = ! ..... G)
Jv#t
(10)
)
(11)
where sgOl is the mean value of s g over Vqt. This formula, which represents the conservation of the neutron balance over the node, is a first relationship between the average nodal flux and the (average) fluxes at the six interfaces with the neighbouring nodes. Let us now return to the c-coefficients. In order to minimize the numerical effort we set t~ = 1, which implies sg = 0 (this will be done only for the present calculation: a generalised use of this recipe would, in fact, destroy the convergence of the method). Then we take w e = (2p+3)22p+2~ 2p+l,p a positive integer. Eq.(10) gives, after a few calculations and writing Qijlgg' for QI.Vtss' = Zijffg' + (Ilk)FitS g',
-241~Otc~g+ ~ QO/g'[a:'+(l-~)c:']---0 g'=l
h2xi
(g= I,...,G),
a GxG linear system which allows for determining the coefficients cxg. A n y numerical work implied by the solution of this system can be avoided by simply letting p - > ,0, since the equations of the s y s t e m become uncoupled and we get, recalling eq.(7). 2 h~ cxg = 24
yG
g. g,
D~ ,'~'=1Q~qt(~ i+1124l - ~g]l/2j/) q g=
(12)
Similar formulas are obtained for Cy g and Czg.This procedure could be considered as a kind of collocation technique, based on the weight functions wg = 8(~-1/2) t$(~+1/2). T h i s also m e a n s t h a t , differently from C U B B O X , only the points at the node boundaries contribute to the determination of the c-coefficients. Expressions (7)(8)(12) for the coefficients ax g, b:d, cxg and their analogues for the y,z directions are now substituted into eq.(5). T h e resulting p o l y n o m i a l representation of the nodal flux is only dependent on the average (volume) node flux of the g-group, #gt, and the six average interface fluxes (of all groups) ea'i±l/2jt, ~g'i.j±ll2.1, ~g'ij,13:ll2. It is remarkably simple and almost as manageable as the corresponding Q U A B O X - i i k e expression, obtained by setting all the c ' s equal to zero. Now we exploit the continuity conditions at the interface between neighbouring cells, Vijt and Vt+ljt say. Flux
B. MONTAGNINI et al.
48
continuity is automatical, since we assume that ~gi+ll2jl -~g(i+l)-ll2jl. The continuity of the (average) normal current,
-
t~ij,'-~'-Ix,+~a=, d,~
d~
= - D~i+lj,~
Ix.,, - ~ , ÷ 1
leads to the following relation (dldx = (l/h.~i) d]d~ ) G
g'=l
i-ll2dl
12 t
+ (_4_ 1~. 8 ~' ~ nxi +
v~ 4.
hx,i+l + (
c~'
- 12 ~/J/
Dg/+ 1 "1 ~gg'- h x ' i * l ,1 ~
2-...~ r~g+ ldl ,ggg'..L hx,i+ 1 ~ -
g'
/)l"g' "w~g' . 1 ~.i+l,j/l~ • 1+312j/ J
12
Fig =
g'
Q~i+gl,jl)~ i+l/2jI
6D~+ i,j/ - r ~ q )
hxi
"x.i+ 1
i+l,jl
(13a)
If the node Vijt is the first (last) of a line of nodes, so that its left (right) face coincides with a piece of the boundary, the continuity condition must be replaced by the boundary condition, eq.(2). This gives, at the left boundary, " g'=l
,
ogg'
{[-(Tilgl~gg''F ~'g ,__4.__xgg' l~xil i,~l,l~,g' ~itjl~.hxil v
,
,xit
~,gl
,~ DSS
ihjl
,
- 12 Fig
~il~/
)J~il-ll2jl
} :--6-
h.~i, lid/
(13b) where il is the first value of the index i in the line j/. A similar formula holds at the right boundary. Eqs.(13a,b) and the similar ones for the other extreme and the other directions, together with eq.(11), determine completely the solution of the problem. If no upscattering is present these equations can be solved group-by-group, at each outer iteration, a procedure which is convenient if the number of energy groups, G, is relatively large. If, however, G is small (G < 3, say), as is the case for most LWR calculations, block inversion procedures are easy. We shall be concerned henceforth with the latter procedures only. R e m a r k 1. The expansion formula based on the average volume and interface fluxes (eq.(5)) yields generally more accurate results than the corresponding formula based on the centrepoint fluxes. R e r m a r k 2. When our first calculations, by an earlier version of this method, were shortly presented (Botta et al., 1987), we became aware of the systematic analysis of a class of nodal schemes, including the QUABOXCUBBOX scheme as well as the one here presented, undertaken by Heunart and coworkers (Hennart,1984,1986, Del Valle et a/.,1985, 1986, Hennart et a/.,1988). The family of nodal schemes of increasing order r "climbing up correctly in order" exhibited by these authors seems to give, in particular, a definitive rational basis to any attempt to increase the accuracy of the QUABOX method by adding higher order
terms to its simple and intuitive basic expansion formula. Careful attention is given to whether rectangular terms should be included or not, a widely known problem which, however, has been dealt with insofar (if exception is made for some other high-order schemes, H6bert, 1987, Rohach, 1987, Ougouag ad Rajic, op.cit.)) by means of such a recipe as the 'quadratic transverse leakage fit' over three adjacent nodes (Bennewitz et al.,1975, Langenbuch et al., 1977b, Lawrence et al., op. cir.) It has been shown for instance that, while a QUABOX (or QUABOX-like) sum ot three quadratic polynomials in x,y,z without rectangular terms (r = 0) is consistent with the use of only 'zero order moments' at the boundary, i.e. the familiar continuity conditions for the average fluxes and currents, this is not the case for a CUBBOX-Iike scheme, which involves cubic polynomials: according to Hennart, in order to get an O(h 2) convergence in the H l norm (h = max(hxi,hyj,hzl)) one should introduce rectangular terms and, consequently, impose further continuity conditions at the node boundaries, involving 'first order moments' ( obtained by first order Lagrange polynomial weighting of the flux along each side of the node; this consistent scheme is, of course, considerably more complicate than the CUBBOX 'sum' scheme). Otherwise, the order of approximation is only O(h), as for QUABOX, no matter how high the degree of the polynomials for the flux expansion may be. This explains, by the way, the rather surprising results of some calculations concerning typical benchmark problems (Trentavizi, 1985), namely that, in absence of higher order boundary moments, no accuracy improvement could be observed, if the cubic polynomials were replaced by sums of I-D analytic solutions of the diffusion equations in each of the variables x,y,z, a scheme which can be considered as a r = =, extension of the CUBBOX method, along the line of the 'analytical' approach (see also Hennart, 1988). The above error estimates have, however, only an asymptotic character.There is therefore no wonder if, for the typical values of the physical parameters and the mesh sizes conmlonly used in LWR core calculations, the cubic formula without rectangular terms and only zero order moments at the boundary gives definitely better results than the quadratic one. The effect of the non-optimal o'der of convergence of the scheme is sensible only if the mesh size becomes exceedingly small. We have even found that, due to the simplicity of our expansion formula, the accuracy improvement obtained by including mixed terms (although in a slightly different form with respect to Hennart's) was more than balanced by the increased complexity of the computations: the best way to get a higher accuracy, within the present calculational method, is the trivial one of reducing (e.g. halving) the mesh size. Thus, we decided not to include any rectangular terms in our cubic scheme, nor higher order moments at the node boundary.We even decided, in the same spirit, to avoid the use of the 'transver~ leakage fit'.
4. S O L U T I O N EQUATIONS
OF
THE
DISCRETIZED
Eq.( 1 I), written in matrix form, for the 2-D case, reads: ( tl)i.ll2 J - 2 ~ , j + ~i+tr2j
6 D~ k
h2~
+ ~ij+1/2 - 2~ij + tI)ij+l/2 ) - Qa,ij dPV = s U
(14)
A well-balanced coarse-mesh flux expansion method Interpreting this equation a~ a classical five point formula, we may write down the following ADI scheme -
49
where
~,~(n,) hxi ~ ij = h~iDij + " ~ i J
_ . ~ . ar.(m+l/2) 9dh(m+ I/2) = t i n + l / 2 ) 2 ~,"~'i.1125 " ~ i j +qJi+ll2j )
hxi
N(..O ,j = h~xixiDij _ ~
+ tox(m) d~(m+l/2) - D'i) Qa,ij d~('m+l/2)
..,3__Dij U-m1 ij - hxi
QiJ . h-~xi DiJ Um,q -1
and ..(1.- r,~ (rn) ,~,~(m) . ,~(m) ~ .... era) ,~(m) . r~- 1 2 I"~'ij'll2 - L"t'iJ ar ,,..k-i,]+l/2 ) I- tojt "~'ij "1- Jtaij Sij hyj
=
(15a) - ~ t,ar~(m+l) ___(re+l) =(re+l) . h2" I."~'ij-112 - 2¢Pij + q)ij+l/2 ) .-yj
_ Di] ~~ , t j . ~(m+l) + toO,,) y ~(m+l) =ij ~tj = ~
t,,r,(m+l/2) 9d~(m+l/2) ~(m+l/2) 2 ['t"i'll2'] " ~ i j + tPi+l/2,j )
hxi
C0(yr")*~7+"~)+ Di] s,j
+
(15b) where m is the 'inner' iteration index and COx0"), a,~y0") are the ADI acceleration parameters. However, we shall not make a direct use of this scheme. Let us denote the r.h.s. of eq.(15a) by r/jr"). Then, having introduced the matrix h2
eq.(15a) is rewritten as follows ¢'I'r~(re+l/2)
"V'( tn+il2)"
h2i r(.?O
+','i+mj )+ ~
v
and we get. upon inversion of U.,,9. (an easy matter if G <3), 2 a~(~+~rz) I -I .~(,n+l12) ~(m+I/2). + hxi iri. I .. r(,n ) ~U ---- ~ Um,ij (°)i-ll2j +q'~i+l/2j ) 12 ~m,o q (16) We see therefore that the average node flux, d~q ('+m), is only dependent on the arithmetic mean of the values of the interface fluxes along the chosen direction, plus a term which is known from the preceding semi-iteration. Let us now recall eq.(13a), also written, in matrix form, for the 2-D case and the semi-iteration of index m + 1/2:
(~xi Dij + --~ hxi -- ~¢ra+l/2) {,~ij)cPi. lt2~
+(
2
**x,i+1
(re+l/2)
DiJdPq
~_
+ t.
t;x.i+ I
(m+l121 Di+ljf~i+lj
it\T(tn)
,,~F~(m+l12)
m) + ~ ; i + l j ) "~'i+ll2j l~m)
+
((1-ff)/k) FdPij(") and Qa, ij = Qa, ij('t') = ~ij + (Odk(m)) Fij where
k(,,) = k(,,-u2)
~v(m)
t x ~ ' ) + (1 - O0VO ' ' m )
(18a)
,a-,.On+l/2)
Oij tin+u2) and sij, Qa;iy are modified correspondingly, with ~¢(m+I/2) a W <"+m) + (1 - o 0 ¢ " )
Using eq.(16), after a few manipulations we arrive at the following block three-diagonal system with only interface fluxes as unknowns: ~ i j m ) .a-(m+ I/2) ± 'qt'i'l/2,J ~ (
1. W h e n the m e s h size is small, a wise convergence strategy may require a relatively large number of inner iterations per each outer iteration. If the node size is comparable to that of a LWR fuel assembly, however, it is convenient to update the fission source and the eigenvalue as often as possible, for example at every inner iteration or even, as in the present method, at every inner semi-iteration. This means that, in our case, the outer and inner (semi-) iterations coincide. Thus, Oij (") is to be identified with Oij (') in eq.(15a), so that sij = sij(m) =
k(m+l/2) = k (m)
Di+Ij + hx,i+l r~ . .~(nl+l/2) ~i+1 j)'~'i+312j
=
l--OT-~T~- l
At the next semi-iteration, O//,0 is to be identified with r~ ,J)"tl)(m+ 1/21i+1/2j Di+l.j - ~hx'i+l ~i+1
+(h_.~/Dij _ ~hxiQij + ~
EQ.(I7) holds for internal nodes; obvious changes (eq.(13b)) are to be made for the first node (the leftmost node of the line) and, similarly, for the last node. Once the three-diagonal equations have been solved, by means of a direct inversion (e.g. Gaussian elimination) subroutine and the interface fluxes of each line of nodes computed, the average volume fluxes can be immediately updated, via eq.(16) ( see fig.l) and the scheme is ready to be applied to the columns, along the y-direction. No auxiliary mesh sweeps (e.g. for updating coupling factors) are needed.
Remark
1 ¢,"1~(re+l/2)
= ~-.-i-mj
Di+l,i LI~,I+i,/~r~j
~),j U-lmij r~T) + ~
Fig.1 A line of nodes: o old interface flux @ new interface flux ® volume flux
u,.,~ = ~ + ~ (~") ~ - D}) Q,,,0)
U,.,ij~v
p(m) = ~
(m)
i+~,/'-'i+3t24 = Pij
(17)
Remark
(18b)
2. The procedure for the choice of the to acceleration parameters, although classical (Varga, 1962; Wachspress, 1966) requires some words of comment. An 'equivalent model reactor', i.e. a one-group, bare, homogeneous, rectangular model of the actual (2-D) reactor is first constructed and a rectangular grid of equally spaced points superimposed to it. The eigenvalues and eigenvectors of the two commuting matrices representing the (finite difference) central difference terms with respect to x and y being known (Varga, op.cit., p.213 ff.), a standard minimax procedure allows to determine the
50
B. MONTAGNINI et al.
'optimal' values of cox and toy, which are taken as 'good' ones also for the actual problem.It should he remarked that the spectral intervals to be explored (in order to eliminate the higher order eigenvectors) must obviously exclude the fundamental eigenvalue of both matrices (see the Appendix).
$.
3-D
P R O B L E M S
Sharp estimates of the convergence ratio of the ADI schemes have been obtained, as well known, only for certain simple 2-D finite difference problems. For 3-D CM problems the situation is certainly worse. In order to extend our method to 3-D problems we have therefore adopted a quite prudential, almost totally implicit scheme. Eq.(15a) is replaced by the following one 6 t,,-~ip+I13) - -ip+ll3 ip+ll3) x - ~'~ U't'i-l/'2jI - 2q~ijl + {~i+ll2jl) hxi _~/,~ip+l/3) - 2,~ip4-1/3 + .~ip+l/3)., h2j u't'ijll2"l ,,t.ifl ".t.,i~/+ll2.1)
#'~ip) ,~(~,+1/3) + COip)tI)(~''+l/3)l}l- T-~-I l.,~ijl ~a,ijl "~ijt = 6
ip) ip) ip) ,,,gp)m~.1+ r~:!. s~;) 2 (tYPij't'll2 - 2¢y~ijl + ¢~ij,t+l/2) + -- Y q l ~lii ijr hzl (19) where COO,) is a new acceleration parameter (>_ 0; experience has shown that if tx is not IOO high, a < 0.8 say, then td,") can be pul equal to zero) and Qtp) _ + ...R_ Fi/i a,qt - Zqt ktp) sip) i/t -- ~ )-F i j t
the Gaussian inversion of each line of the layer a r e s t o r e d , to be used for all the remaining iterations of the m-cycle (which imply therefore only a multiplication of a matrix by a vector; k , cox, coy are kept fixed until the cycle is completed). This speeds up the computation of each layer, at the expense of a negligible additional memory requirement. R e m a r k 2. It has been observed that a I-D analytic treatment of the top and bottom reflectors (for details s e e Camiciola et al., 1986) may lead to some improvement of both accuracy and convergence.
6. N U M E R I C A L R E S U L T S Some numerical results concerning the classical 3-D PWR IAEA Benchmark Problem (Wagner, 1975) are shown in fig.2, which reports the errors (%) of channel power (the nodal power integrated along z ) with respect to the reference calculation (Wagner et al., 1977) for a fullassembly calculation (9x9x17 nodes, each one of 20x20x20 cm, with analytic treatment of top and bottom reflectors) as well as a quarter-of-assembly calculation ( 17x 17x 17 nodes, 10x 10x20 cm each). These calculations correspond to the two classes of accuracy referred to in t h e Introduction. Max errors are 6.5% and 1.4%, respectively (they correspond to the underlined figures), whereas the mean (r.m.s.) errors amount to 2.1% and 0.8%, respectively. Such errors are within the error bounds mentioned in sect. 1. With a convergence criterion of 10 -4 on the pointwise fluxes and 3 m-iterations for each cycle, the two computations require 15 and 20 ADI cycles respectively, and a CPU time (1/4 core symmetry) of 24 and 111 seconds, respectively, for the VAX 6400 machine. For the Macintosh IIfx PC the computing time may be 10 times longer. The CPU time of a test run on the CRAY YMP8/464 computer (used as a sequential computer) has been, for the two cases, of 1.85 and 8.28 seconds, respectively.
".vijl"~ip) k = I. 02972 k = I. 02918 k = I. 02903
9x9 17x17 Refer. - 2.9 6.5 - 2.4 - 0.8 0.9 - 0.8 0.476 0,700 0.610 0.0 1.9 2.7 - 0.3 0.1 0.4 0.0 - 0.6 1.177 0,971 0.923 0.866 - 2.0 - 0.6 0.6 0.0 0.9 - 4.3 - 0.1 0.2 0.4 0.0 0.0 - 1.3 1.368 1.309 1.180 1.088 0.999 0,709
k~) is updated according to the following relation kip+l/S) _-- kip)
IP(p+l/3) t~q/(p+I/3) + (I - ~)q.rip)
and the ADI scheme is completed by two similar equations (with a cyclic permutation of the indexes i j , l ), for ~ip+2/3) and ~ip+l) and the corresponding eigenvalues. As eq.(191 represents, for fixed l, a 2-D (discretised) diffusion problem in the xy plane, we only need apply the method of sect.4. We obtain by this way, after a moderate number of m-iterations (eqs.(15a,b)), reasonably converged flux values for the l-th node layer or, otherwise stated, a 'reasonably good' inversion for the 2-D diffusion problem for that layer (the whole 3-D process becomes therefore similar to the classical 2-D scheme, where 'exact' 1-D inversions are alternately performed along the two directions). When all the /-layers have been processed, a second cycle of iterations is started, the nodes being now regrouped according 1o layers parallel to the yz plane. The third cycle, parallel to the xz plane, concludes the scheme. R e m a r k 1. A variant of this scheme is the following one. After the fast m-iteration, the block triangular matrices for
- 0.9
- 3.2 - 0.6 0.729
- 0.9
- 0.6
0.4
- 0.I
0.6
- 0.4
0.1 0.2 0.0 0.4 - 0.1 - 0.2 - 0.5 1.396 1.430 1.290 1.071 1.054 0.975 0.756 - 0.1 0.2 0.6 - 3.6 1.1 1.3 - 1.4 0.5 0.5 0.5 - 0.7 0.1 0.0 - 0.7 1.281 1.421 1.193 0.610 0.954 0.959 0.777 Fig.2.3-D IAEA PWR Benchmark Problem
7. C O N C L U S I O N S A non-conforming and even non-consistent, but nevertheless quite efficient coarse-mesh flux expansion method has been worked out, which can be used as a fast
A well-balanced coarse-mesh flux expansion method running LWR core calculational tool as well as a reasonably accurate method for reference problems.
Acknowledgements- We are indebted with F.Oriolo for many discussions on the numerical part of the method. This work was performed under the auspices of the National Group for Mathematical Physics of the Consiglio Nazionale delle Ricerche, with the help of a contribution of the Ministero della Pubblica Istruzione. An early version of the computer program was supported by the Ente Nazionale Energie Alternative, Italy.
APPENDIX
The itemtive fonnula for the eigenvalue k Let L = V . D V + ~ . The outer iterative process represented by eq.(3a), i.e. _ (L+...II-F)tl~(n+l) = ! - a F t b (n)
k (n)
Bennewitz F., Finnemann H. and Wagner M.R. (1975) Trans.Am.Nucl.Soc. 22, 250 Bonalumi R.A., Rouben B., Dastur A.R., Dondale C.S. and Li H.Y.H. (1978) Proc. of A.N.S. Top. Meeting CONF-78.0401, Gatlinburg, p. 169 Botta E., Camiciola P., Cundari D., Montagnini B., Raffaelli P., Soraperra P. (1987) ANS-ENS
lnt.Top. Meeting on Adv.in Reactor Phys., Math. and Computation, Paris, p.1247 Camiciola P., Cundari D. and Moutagnini B. (1986) Ann. nucl.Energy 13,629 Del Valle E., Hennart J.P. and Meade D. (1985) Proc Int. Meet. on Adv.Nucl.Eng. Comput.Methods, Knoxville, Tenn., p.473 Del Valle E.,Hennart J.P. and Meade D. (1986) Nucl. Sci.Engng. 92, 204 H6bert A. (1987) Ann. nucl. Energy 14, 527 Hennart J.P. (1984) Nodal Methods for the Numerical Solution of Partial Differential Equations, in Lect. Notes in Mathematics, Springer-Verlag, p. 1230 (Proc. Conf. on Numerical Analysis, Guanajuato, Mexico) Hennart J.P. (1986) SIAM J. Sci. Stat. Comp. 2,918 Hennart J.P. Numer.Methods for Partial Diff.Eqs. (1988) 4. 233 Hennart J.P., Jaffre J. and Roberts E. (1988) Numer. Math. :~., 701 Langenbuch S., Maurer W. and Werner W. (1977a) Nucl. Sci. Engng. 63,437 Langenbuch S., Maurer W. and Werner W. (1977b) Nucl. Sci. Engng. 64, 508 Lawrence R.D. and Doming J.J. (1980), Nucl. Sci. Engng. 76, 218 Ougouag A.M. and Rajic H.L. (1988) Nucl. Sci. Engng. 100, 332 Rajic H.L. and Ougouag A.M. (1989), Nucl. Sci. Engng. 103,392 Rohach A.F. (1987)Ann. nucl. Energy 14,653 Shober R,A., Sims R.N. and Henry A.F. (1977) Nucl. Sci. Engng. 64, 582 Trentavizi C. (1985) M. Sci. Thesis, Univ.of Pisa, Italy Varga R.S. (1962) Matrix lterative Analysis, PrenticeHall, Englewood Cliffs, N.J. Wachspress E.L. (1966) lterative Solution of Elliptic Systems, Prentice-Hall, Englewood Cliffs, N.J. Wagner M.R.(1975) Proc. Conf. on Comp. Meth. in Nucl. Eng. CONF-75.0413, Charleston,J., 1 Wagner M.R., Finnernann H., Koebke K. and Winter H.J. (1977) Atomkernenergie 30, 129
k (n)
(A 1)
can be given the following equivalent form: ¢b( n + l ) = - ( I +
REFERENCES
51
t~ L ' I F ) ' I 1 - a L - 1 F ~ ( n ) k (n) k (n)
(A2)
Let kh, ~h, be the eigenvalues, resp. the eigenvectors of the operator K = - L "1 F , so that I ~ h = k h ~ h . When referring to the fundamental eigenvalue and eigenvector we shall omit the index: k = k l , C p = ~ 1 . As the fundamental eigenvalue is simple, we can assume that 0 < kh < k (h -> 2).Then, by expanding the initial guess flux ~(o) in the eigenvectors ~h , i.e. ~(o) = Zh ChdPh , eq.(A2) gives, for n = 1, ~O) = _ (I +-fl[-LlF)
"l 1 - ot L - I F ~ ( ° )
k (o)
k (0)
= ~'h Ch (1 - ot)(kh/k(0)) ¢Ph I - ~(k~Ik (°)) where k (°) = i is the initialguess for k. For the second
iterate we have
dP(2) --" ~h Ch ( 1 - ot)(khlk(°)) ( I - ot)(kh/k(1)) ~h 1 - ot(kh/kt°)) 1 - ct(kh/k ~)) and, for the n-th iterate, ~(") = ~ h ch (1 - ot)(kh/k(°)) .. (1 - ot)(kh/k(nt)) ¢p~ 1 - a(kh/k(O))
I - a(kh/k (n'D)
= p(O)"p(n'l)[Cl,+h~2£h ('h/k)n{Y(O'..~(ha-1) ,h] (A3) where we have introduced the following abbreviations:
p(m) = p(ra)(Of) = (1 - oD(k/k (m))
1 - ot(klk (m)) O.(hm)= o.~m)(o[) = 1 - ~'(]Gtk (m>) 1-
of(kh/k(m))
We observe that the inequalities 0 < kh < k (h > 2), together with the inequality 0 < o~(k/k(n)) < I
(A4)
0 < crh(n) < 1
(A5)
imply that Note that a too high value of a could result into a violation of this inequality. Let, for instance, k = 1.02, k2 = 1.01 and take de = 0.99. With k(°) = 1 we get o~(k/k(°)) = 1.0098 and o'2(°)(0.99) = -98, with disastrous
B. MONTAGN1NI et al.
52
consequences.We shall assume in the following that (A4) is always fulfilled. Then we have, for any vector w with non negative components, using (A4) and the relations K~b
= k¢~, K~h = kh % (h > 2),
where use has been made of eq.(A1). Note that the last quotient represents the appropriate production-to-loss ratio at the (n+l)-th iteration. Now LTe = (V.DV + Z T) • = zTe = (-Zr + •T) e-- - Za
(K~(n),w) (in fact, - Zr~ * Z ~ = t Z ~ - ~ ' = - Za~ ).Then we have
(~(n),w)
Cl(¢P,w) + ~
ch (kh/k) n+l cr~°) ''
crh("'l)(~,w)
(- L O ("+l), e) = (O (n+t), - LTe) = (~(n+t), Za)
h>2
=k
c , ( ~ , w ) + Y~ ~
( F(b (n+t), e) = ( L L - t F ~ ("+l), e)
(kak)",,~o, ..o~,-,)(,~,w)
h>2
= (- L-tFO (n+t), - Lre) = ( KO ("+1), Za )
Here
(u,v)=
f(~ ugvg)dV Jv g = l
(when the discrete form of the equations is considered, the integral is replaced by a s u m over the nodes). Then, on account o f (A5), the inequalities 0 < kt, < k (h >2) and (~,w) > 0 (where the latter inequality follows from the positivity of the fundamental eigenvector ~ and the nonnegativity of w) we conclude (KO(n)'w)
*k
(A6)
(~("),w)
k ("+t) = k (n) (~(n+l), w) (¢'("), w) (w > 0), then we can easily show, using (A3), that k(n+l) k (")p(") = k (1 - t~)/[1 - ot(klk (n))] or, in the limit as n -> oo 1 - ot(klk(")) = (1 - t~)(k/k("+l)), a relationship which can be written, putting Mk(m) = 1 + ~,n), -
as n -> oo. Let now v2;f, Za be two G-vectors with c o m p o n e n t s v2;/l(r) ..... v2;fe(r) and 2~al(r),...,2;ae(r), respectively, ,~ag being the absorption cross section in group g . W e show that the r.h.s, of the following recursive formula, which is clearly equivalent to eq.(3b) of sect.2,
k ("+t) = k(")
and substitution into eq.(A8) leads to the conclusion. Eq.(A6), with w = Za, then shows that /¢(n+l).> k as n -> oo. Note that other recursive formulas, although asymptotically valid, m a y introduce u n n e c e s s a r y restrictions on the parameter ~t. If we choose, for instance, the following formula
~(n+l) __ _ 1 ~O~ ?(")
(A9)
O f course, we would like to have 7(") -> 0. But (A9) shows that this is not true, in general, if f f ~ 1/2.
(°("+l)'vZf)
&'((I ~("+I),VZf) + (I - £Z)(dP(n),VZf) (A7) can be written in the form (Kq~(n+l),w)/(~("+l),w) ( i.e. in the same form as the l.h.s, of eq.(A6)), with w = Za. Let e be the G-vector the components of which are all equal to one, over the whole reactor domain. It is immediately seen that FTe = vet, since the sum of the Xg is the unity. Then eq.(A7) can be given the following form
(dp("+1), •Te)
k (n+l) = k (")
C ~ ("+1), FTe) + (1 - C~)(~(n), FTe) ( F ~ (n+l), e) ( a F~("+i) + l - f f F c b ( n ) , e ) k(") k (") ( F ~ ("+l), e) (- L(b ("+1), e)
(AS)
Choice o f the parameters ~
toy
Let us assume that the approximate eigenvalues k(n) have already converged to a sufficiently accurate value, k*.The iterative process described by eqs.(15a, b) o f sect.4 requires a slight modification of the classical procedure of the cyclic variation of the co-parameters. The modification is necessitated by our use of only one inner iteration pex outer iteration, which implies that we are dealing with the acceleration of the power method for an eigenvalue problem rather than the acceleration of the inner iteration cycle for an inhomogeneous linear system. As previously said, a rectangular, bare, homogeneous equivalent 'model reactor' is first constructed and a rectangular grid x~ =/hx, yi=jhy (i = 0,1,...,nx, j = 0,1,...,ny), with approximately the same m e s h width as the actual problem is then superimposed over the domain {x,y; O.~x
A well-balanced coarse-mesh flux expansion method _4.- r a,(n+l)
where use has been made of eq.(Al2) and, similarly,
2 O ( n + l ) + a,(n+l) x
h~ ~v,i-1:2,i-
53
v,i+v2jj
+A ~yy ._(n+l) (~)i,]-tl2 - ~ - _i(jn + l ) + --(n+l) . + ...~ ~i,]+I/2) (,7.~--~
~,)~(n+l)
~)
+ ~-
a(F'/k*) + ~.~+/~
-
(A14)
(o~n) + Z - a(Y/k*) + / ~
~!~)
- ~ -
k*
-,',1
(A10)
with the condition of zero flux at the boundary Using vector notations, the ADI procedure for solving each step of this iterative process (Varga, 1962 or Wachspress, 1966) can be given the following form:
one sees that In has eigenvalues 0p4(n) = 0,1u(n). 0~,~(n). For a oJ-cycle with n* values of the ~a (~),coy(n) parameters, we may define a global iterationmatrix
n* with eigenvalues
J* = H In n=l
O;q -- fi 8p,q (n) n=l
--[-.
From (AI3)(AI4) w e o b t a i n 01A (n) = 0&l (n) . 0~,1 (n) -- l [- X +
~
~
(n+l/2)
where ~(n) has components ~i~(n)q(suitably ordered) and X~(n) y~(n) represent the central difference terms with respect to x and y, respectively. It is known (Varga, op.cit,, p.214) that the vectors ¢(P.q) with components ¢iff "q) = C O''q) sin (pax ilia) sin (qT~'yj Ily) = C(P'q ) sin (pgilnx) sin (qajlny) are eigenvectors of both X and Y, with eigenvalues 37,= 2[1..cos(prclnz) ]lha 2 and #q = 2[ l-cos(qftlny)]/hy 2, respectively. In particular, {~(1,1) is the fundamental eigenvector of X and Y and therefore of X + Y and represents the solution of the model problem. Substitution of ~O.~) into eq.(A10), as 'converged flux', leads immediately to the criticality equation:
~.l + I.q + Z - (o~/k*)F = (1 - o:)~'/k*
(AI2)
Now, the ADI iterative process, eqs.(A10), can be summarised l~y the single relationship
(1)(n+l)= ~'n(I}(n) where the iterationmatrix. Jn, has the following form
• [ x+ (o: +
+
Then, setting
o;(,)
o6") + (1- a)(P/k*) - X~ ~ " ) + x- a(P/~*) +
(0~n) + T. - a('ff/k*) +
),~
(A13)
for all n and therefore 0* m -- 1. If n* is large enough and we were able to determine the co-parameters in such a way that, for any p,q ~ 2, Ùx~(") = 0 or 0a,q(n)= 0 for some n, then the process would give the fundamental eigensolution ¢~O,t)exactly, in a finite number of steps. This result is, however, of no interest, since ¢~OA)is explicitly known. Rather, we try to solve, by the classical techniques (Varga, op.cit., Wachspress, op.cit.) the following minimax problem
min ~.,~
"" c~') + T. - a(':/k*) + ~h + I~ " x max 1-'[ a,~ " = ~ c4 n) + ~. - a('f/k*) + x eof")+ Z - o~(~/k*) + 3,i + gl - Y o V ) + r- -
a(':/k*)+ y
with x.y ranging over an interval containing all the ~, resp. #q, with p.q > 2, but excluding 3,1 and #l. The optimal fo's so determined are then used as 'good' acceleration parameters for the actual problem to be solved.
A WELL-BALANCED
0306-4549/94 $6.00+0.00 Copyright © 1993 Pergamon Press Ltd
COARSE-MESH
FLUX EXPANSION
METHOD
B. MONTAGNINft, P. SORAPERRA*, C . TRENTAVIZI* Dipartimento di Costruzioni Meccaniche e Nucleari, Universit~ di Pisa, via Diotisalvi 2, 1-56016 Pisa, Italy M. StIMINI, D.M. ZARDINI* Laboratori di Ingegneria Nucleate di Montecuccolino, Universit~ di Bologna, via dei Colli 16, 1-40136 Bologna, Italy (Received 3 July 1993) A b s t r a c t - An efficient low-order coarse-mesh method for LWR criticality problems, which allows for fast routine computations as well as for more accurate, but still unexpensive, reference solutions, is presented. The efficiency of the method rests upon a very simple nodal flux expansion formula and an improved scheme for outer iterations, coupled with a fast ADI technique.
1. I N T R O D U C T I O N
NOMENCLATURE a
#
=
¢g(r)
Neutron diffusion coarse-mesh (CM) methods are usually grouped into several classes, according to the nature of the underlying approximation technique. A more pragmatic classification, consisting of only two classes, refers to their domain of application. The first class concerns the rather rough (nodal power max error < 7%0, mean error < 2.5%, say) but very fast CM techniques used in routine computations. The second class is constituted by the much more accurate (max error < 1.5%, or so) and theoretically refined CM methods which are used for reference purposes, mostly to assess the reliability of the methods of the first class. In this paper we wish to illustrate a wellbalanced CM method for steady-state (criticality) problems, which allows for both kinds of calculations (with a proper choice of the m e s h size), within approximately the error limits indicated above. No attempt has been made in order to lower such limits, since they seem to be acceptable (for the appropriate class of problems), if allowance is made for the uncertainties of the physical parameters. Care has been taken, instead, to achieve a good computational efficiency.The present work has, however, a further, and perhaps deeper, motivation. The success of the CM methods based on 'analytical' formulations or, more precisely, on the representation of the 'transverse-integrated' flux by means of 1-D solutions of the diffusion equation (Shober et hi., 1977, Bonalumi et aL, 1978, Lawrence and Dorning, 1980, Ougouag and Rajic, 1988, Rajic and Ougouag, 1989, and others) has a little obscured the family of the CM methods based on polynomial approximations. The latter methods have, however, their own merits. They do not require the lengthy computation of many transcendental functions
= Total number of neutron energy
=
groups Neutron flux in group g (at the point r) (n/cm 2 s)
Dg = Dg(r) = Diffusion coefficient in group g (cm) 2 ; / = 2~rg(r) = Macroscopic removal cross-section in g r o u p g (cm -1)
xsg'->g =
27sg'->g(r) = Macroscopic transfer cross-section from g r o u p g ' to group g (cm -l)
v~g = v ~ g ( r )
= Macroscopic fission cross-section
times the average number of secondaries, in group g (cm -I)
z g = Total fission spectrum fraction in group g k = Effective multiplication factor t To whom all the correspondence should be addressed. * W o r k p e r f o r m e d as a partial fulfilment for a M.Sc.Thesis.
45
B. MONTAGNINI et al.
46
(this may be an important advantage, e.g. if the scheme is a candidate for an extension to time-dependent problems) and their structure is simple, at least if the degree of the polynomials is not too high.The relative importmlce of the various terms of the approximation formula can be easily investigated and it may be possible, by suppressing some of them, to obtain a consistent simplification of the scheme without compromising the accuracy in an essential way. One of our aims was just to show that, for a definite class of diffusion problems as those concerning the LWR cores, a quite simple polynomial formulation can be more efficient and not much less accurate than other schemes (analytical or not). The method is based on a low-order expansion of the nodal flux, namely a sum of cubic polynomials in each one of the space variables x,y,z, as in the CUBBOX method ( ' s u m ' expansion formula, Langenbuch et al., 1977a).The solution procedure has undergone, however, several modifications with respect to CUBBOX.
2. T H E FLUX E X P A N S I O N F O R M U L A The time-independent multigroup neutron diffusion equations are (see the Nomenclature): G
V.D~V¢
~-
- Z~ + Z,,)~ +kl-F~ = 0
(V.DV
Having m mind a solution procedure based, as customary,
on some outer/inner iterative process, we first write down the outer iteration scheme (i.e. the iterative process involving the fission source): -(V.DV +Z'
o~ F i ~ ( " + l ) = I - ~ F ~ ( " ) ' k(")
(3a)
k(")
where E = - T-r+ ~'~s (thus,
fee'
=
.
~,rg~gg ' + .~$g">g
with ~gg" = 0 Of 1 according to g a g ' o r g = g ' ), n is the iteration index and tx, 0_
0~I//(n+l)+ (I - OQI//(n) (g = l ..... ~ )
(I) These equations are to be solved subject to the following homogeneous conditions at the reactor boundary OV :
(2)
where b/Ou denotes the outward normal derivative and (r g, z g are real constants. It is convenient to rewrite eqs.(1) and (2) in matrix form. Let D - - c o l [4) l ..... ¢G] and 1
G
'~r = diag [Zr ..... Zr ] T = diag ['r I ,... z G]
0
~)~I
...
~g~l
Z~ I~2
0
...
Z ~2
where
g=l V being the reactor volume.The denominator at the r.h.s. of eq.(3b) thus coincides with the total fission source (depending on a linear combination of both fluxes ¢I~n) and
(g = 1..... G)
OU
S = diag [o"1..... O'G],
I
k ('+I) = k (')
g'=l G + ~ - z ~ Y .-- v z : % ~ ' = 0 , g'=l
D = diag [D 1..... D G],
(23
s o + T ~---~= 0
X:¢g+~ r.:'-~ ¢~;
(Tgog+r gong =0,
(I')
Zs =
¢b("+t)) which sustains the flux in eq.(3a). Since the r.h.s. of eq.(3a), i.e. st") = I - t~ F~(n) k (")
(4)
is known (when the n-th step is completed) we have actually to solve an inhomogeneous system of partial differential equations. We assume that the reactor domain cau be represented as a regular array of rectangular parallelepipedes VO.t ('cells', or 'nodes'), with centres (xi,Yi, Zt) and edges of length hxi, h W, hz/, and that all the nuclear parameters are spatially cohstant over each Vqt. [laving introduced local dimensionless variables
= (x - xi)/hxi, 77 = (y - yj)/hyj, ( = (z - Zl)/hzl we approximate the g-group flux inside the node V//t by means of the following polynomial formula
...
0
~g(~,)?,Q = ~gl + axg~ + bxg(~ 2- -~2) + cxg(~- 4~ 3) F
+ aye 0 + byg(o 2 - 1-~2) + cyg(71 - 4r/3)
=
G
1
+ a z g ( + bzg(( 2- 1~) + czg( ( - 4(3) (these matrices depend, in general, on r). Then eqs.(1) and (2) can be written as follows:
(5) which coincides essentially with the 'sum' formula used in CUBBOX (Langenbuch et al., 1977a). It is written,
A well-balanced coarse-mesh flux expansion m e t h o d however, in such a way that its coefficients may be directly expressed in terms of the average node flux and the average interface fluxes (see below), rather than the corresponding centrepoint fluxes. We have, in fact, t
[ I ¢#(x,y,z) d V hxihyjhzt Jv =
d~'~)g(~,r/,~ ") :
d /2
where A = V.V and we = wg(x,y,z) are weight functions (eq.(4) has also been taken into account). Let us first take wg = 1 for all g. If the approximate expression for ¢g. eq.(5), is substituted into the integral over V/fl.written in terms of the variables ~,r/,(, we get, by performing the partial derivatives, integrating and taking account of eqs.(7) and (8):
$ l
6D~0-t\
t2
so that ¢g/fl actually represents the (approximate) average node flux. As it concerns the interface fluxes, let us define the following 'transversally averaged' flux:
fYJ'~Jdy Ft*~-hadz ~)g(x,y,z)
= ~
47
g g g " 2~/~1 +, ¢h~+tn,t 4" ~i:,/-112 - 2¢~/~1+ ~/1,1+112~ 4- ¢i~-I/2,1
hy2 -
ZG
h 2z, g,
g'
Qgd.Ol(POI= S~Ot ,
(g = 1..... G)
g'=t
= I,= dr/
= d~'¢g(~,r/,O
Then, using eq.(5), t~gx(~) = O~t + axg~ + bxg(~ 2- i=~2) + Cxg(~- 4~ 3) (6) Similar equations hold for
tpyg(r/) mad (~zg(O. The fluxes
0x~(:t:l/2) are, in particular, the average interface fluxes at the right (left) face of the cell and will be denoted henceforth by ~g~l/2dl. Putting ~ = :t: 1/2 into eq.(6), then adding and subtracting, we g e t , as in CUBBOX, the desired expression of the axg, bxg coefficients in terms of the average (volume or interface) fluxes:
ag =
bg
¢~+t/2jt- Cgi-m,it
(7)
= 3(~gi+ll2d/- 2~)g/fl+ Ogi-i/2jt)
(8)
Similar formulas hold for ayg, byg and azg, b.g. Let us postpone, for a moment, the problem of deterrr~ining the r e m a i n i n g c o e f f i c i e n t s Cxg, Cyg, Czg and try to approximately solve, by the weighted residual method, the inbomogeneous diffusion equations, eq.(3a). We denote by D q ff, ,Sr,Otg,..the constant values of Dg, ,Stg .... over the node Vijt. We also put
Qa, ql = ~"qt + ~Fqt
(9)
According to the weighted residual procedure, the i n t e g r a t e d f o r m o f eq.(3a), written in terms of components, or elements, of the involved vectors and matrices and suppressing the index n, is
o~.ql
] wg dV
sgwgdV,
(g = ! ..... G)
Jv#t
(10)
)
(11)
where sgOl is the mean value of s g over Vqt. This formula, which represents the conservation of the neutron balance over the node, is a first relationship between the average nodal flux and the (average) fluxes at the six interfaces with the neighbouring nodes. Let us now return to the c-coefficients. In order to minimize the numerical effort we set t~ = 1, which implies sg = 0 (this will be done only for the present calculation: a generalised use of this recipe would, in fact, destroy the convergence of the method). Then we take w e = (2p+3)22p+2~ 2p+l,p a positive integer. Eq.(10) gives, after a few calculations and writing Qijlgg' for QI.Vtss' = Zijffg' + (Ilk)FitS g',
-241~Otc~g+ ~ QO/g'[a:'+(l-~)c:']---0 g'=l
h2xi
(g= I,...,G),
a GxG linear system which allows for determining the coefficients cxg. A n y numerical work implied by the solution of this system can be avoided by simply letting p - > ,0, since the equations of the s y s t e m become uncoupled and we get, recalling eq.(7). 2 h~ cxg = 24
yG
g. g,
D~ ,'~'=1Q~qt(~ i+1124l - ~g]l/2j/) q g=
(12)
Similar formulas are obtained for Cy g and Czg.This procedure could be considered as a kind of collocation technique, based on the weight functions wg = 8(~-1/2) t$(~+1/2). T h i s also m e a n s t h a t , differently from C U B B O X , only the points at the node boundaries contribute to the determination of the c-coefficients. Expressions (7)(8)(12) for the coefficients ax g, b:d, cxg and their analogues for the y,z directions are now substituted into eq.(5). T h e resulting p o l y n o m i a l representation of the nodal flux is only dependent on the average (volume) node flux of the g-group, #gt, and the six average interface fluxes (of all groups) ea'i±l/2jt, ~g'i.j±ll2.1, ~g'ij,13:ll2. It is remarkably simple and almost as manageable as the corresponding Q U A B O X - i i k e expression, obtained by setting all the c ' s equal to zero. Now we exploit the continuity conditions at the interface between neighbouring cells, Vijt and Vt+ljt say. Flux
B. MONTAGNINI et al.
48
continuity is automatical, since we assume that ~gi+ll2jl -~g(i+l)-ll2jl. The continuity of the (average) normal current,
-
t~ij,'-~'-Ix,+~a=, d,~
d~
= - D~i+lj,~
Ix.,, - ~ , ÷ 1
leads to the following relation (dldx = (l/h.~i) d]d~ ) G
g'=l
i-ll2dl
12 t
+ (_4_ 1~. 8 ~' ~ nxi +
v~ 4.
hx,i+l + (
c~'
- 12 ~/J/
Dg/+ 1 "1 ~gg'- h x ' i * l ,1 ~
2-...~ r~g+ ldl ,ggg'..L hx,i+ 1 ~ -
g'
/)l"g' "w~g' . 1 ~.i+l,j/l~ • 1+312j/ J
12
Fig =
g'
Q~i+gl,jl)~ i+l/2jI
6D~+ i,j/ - r ~ q )
hxi
"x.i+ 1
i+l,jl
(13a)
If the node Vijt is the first (last) of a line of nodes, so that its left (right) face coincides with a piece of the boundary, the continuity condition must be replaced by the boundary condition, eq.(2). This gives, at the left boundary, " g'=l
,
ogg'
{[-(Tilgl~gg''F ~'g ,__4.__xgg' l~xil i,~l,l~,g' ~itjl~.hxil v
,
,xit
~,gl
,~ DSS
ihjl
,
- 12 Fig
~il~/
)J~il-ll2jl
} :--6-
h.~i, lid/
(13b) where il is the first value of the index i in the line j/. A similar formula holds at the right boundary. Eqs.(13a,b) and the similar ones for the other extreme and the other directions, together with eq.(11), determine completely the solution of the problem. If no upscattering is present these equations can be solved group-by-group, at each outer iteration, a procedure which is convenient if the number of energy groups, G, is relatively large. If, however, G is small (G < 3, say), as is the case for most LWR calculations, block inversion procedures are easy. We shall be concerned henceforth with the latter procedures only. R e m a r k 1. The expansion formula based on the average volume and interface fluxes (eq.(5)) yields generally more accurate results than the corresponding formula based on the centrepoint fluxes. R e r m a r k 2. When our first calculations, by an earlier version of this method, were shortly presented (Botta et al., 1987), we became aware of the systematic analysis of a class of nodal schemes, including the QUABOXCUBBOX scheme as well as the one here presented, undertaken by Heunart and coworkers (Hennart,1984,1986, Del Valle et a/.,1985, 1986, Hennart et a/.,1988). The family of nodal schemes of increasing order r "climbing up correctly in order" exhibited by these authors seems to give, in particular, a definitive rational basis to any attempt to increase the accuracy of the QUABOX method by adding higher order
terms to its simple and intuitive basic expansion formula. Careful attention is given to whether rectangular terms should be included or not, a widely known problem which, however, has been dealt with insofar (if exception is made for some other high-order schemes, H6bert, 1987, Rohach, 1987, Ougouag ad Rajic, op.cit.)) by means of such a recipe as the 'quadratic transverse leakage fit' over three adjacent nodes (Bennewitz et al.,1975, Langenbuch et al., 1977b, Lawrence et al., op. cir.) It has been shown for instance that, while a QUABOX (or QUABOX-like) sum ot three quadratic polynomials in x,y,z without rectangular terms (r = 0) is consistent with the use of only 'zero order moments' at the boundary, i.e. the familiar continuity conditions for the average fluxes and currents, this is not the case for a CUBBOX-Iike scheme, which involves cubic polynomials: according to Hennart, in order to get an O(h 2) convergence in the H l norm (h = max(hxi,hyj,hzl)) one should introduce rectangular terms and, consequently, impose further continuity conditions at the node boundaries, involving 'first order moments' ( obtained by first order Lagrange polynomial weighting of the flux along each side of the node; this consistent scheme is, of course, considerably more complicate than the CUBBOX 'sum' scheme). Otherwise, the order of approximation is only O(h), as for QUABOX, no matter how high the degree of the polynomials for the flux expansion may be. This explains, by the way, the rather surprising results of some calculations concerning typical benchmark problems (Trentavizi, 1985), namely that, in absence of higher order boundary moments, no accuracy improvement could be observed, if the cubic polynomials were replaced by sums of I-D analytic solutions of the diffusion equations in each of the variables x,y,z, a scheme which can be considered as a r = =, extension of the CUBBOX method, along the line of the 'analytical' approach (see also Hennart, 1988). The above error estimates have, however, only an asymptotic character.There is therefore no wonder if, for the typical values of the physical parameters and the mesh sizes conmlonly used in LWR core calculations, the cubic formula without rectangular terms and only zero order moments at the boundary gives definitely better results than the quadratic one. The effect of the non-optimal o'der of convergence of the scheme is sensible only if the mesh size becomes exceedingly small. We have even found that, due to the simplicity of our expansion formula, the accuracy improvement obtained by including mixed terms (although in a slightly different form with respect to Hennart's) was more than balanced by the increased complexity of the computations: the best way to get a higher accuracy, within the present calculational method, is the trivial one of reducing (e.g. halving) the mesh size. Thus, we decided not to include any rectangular terms in our cubic scheme, nor higher order moments at the node boundary.We even decided, in the same spirit, to avoid the use of the 'transver~ leakage fit'.
4. S O L U T I O N EQUATIONS
OF
THE
DISCRETIZED
Eq.( 1 I), written in matrix form, for the 2-D case, reads: ( tl)i.ll2 J - 2 ~ , j + ~i+tr2j
6 D~ k
h2~
+ ~ij+1/2 - 2~ij + tI)ij+l/2 ) - Qa,ij dPV = s U
(14)
A well-balanced coarse-mesh flux expansion method Interpreting this equation a~ a classical five point formula, we may write down the following ADI scheme -
49
where
~,~(n,) hxi ~ ij = h~iDij + " ~ i J
_ . ~ . ar.(m+l/2) 9dh(m+ I/2) = t i n + l / 2 ) 2 ~,"~'i.1125 " ~ i j +qJi+ll2j )
hxi
N(..O ,j = h~xixiDij _ ~
+ tox(m) d~(m+l/2) - D'i) Qa,ij d~('m+l/2)
..,3__Dij U-m1 ij - hxi
QiJ . h-~xi DiJ Um,q -1
and ..(1.- r,~ (rn) ,~,~(m) . ,~(m) ~ .... era) ,~(m) . r~- 1 2 I"~'ij'll2 - L"t'iJ ar ,,..k-i,]+l/2 ) I- tojt "~'ij "1- Jtaij Sij hyj
=
(15a) - ~ t,ar~(m+l) ___(re+l) =(re+l) . h2" I."~'ij-112 - 2¢Pij + q)ij+l/2 ) .-yj
_ Di] ~~ , t j . ~(m+l) + toO,,) y ~(m+l) =ij ~tj = ~
t,,r,(m+l/2) 9d~(m+l/2) ~(m+l/2) 2 ['t"i'll2'] " ~ i j + tPi+l/2,j )
hxi
C0(yr")*~7+"~)+ Di] s,j
+
(15b) where m is the 'inner' iteration index and COx0"), a,~y0") are the ADI acceleration parameters. However, we shall not make a direct use of this scheme. Let us denote the r.h.s. of eq.(15a) by r/jr"). Then, having introduced the matrix h2
eq.(15a) is rewritten as follows ¢'I'r~(re+l/2)
"V'( tn+il2)"
h2i r(.?O
+','i+mj )+ ~
v
and we get. upon inversion of U.,,9. (an easy matter if G <3), 2 a~(~+~rz) I -I .~(,n+l12) ~(m+I/2). + hxi iri. I .. r(,n ) ~U ---- ~ Um,ij (°)i-ll2j +q'~i+l/2j ) 12 ~m,o q (16) We see therefore that the average node flux, d~q ('+m), is only dependent on the arithmetic mean of the values of the interface fluxes along the chosen direction, plus a term which is known from the preceding semi-iteration. Let us now recall eq.(13a), also written, in matrix form, for the 2-D case and the semi-iteration of index m + 1/2:
(~xi Dij + --~ hxi -- ~¢ra+l/2) {,~ij)cPi. lt2~
+(
2
**x,i+1
(re+l/2)
DiJdPq
~_
+ t.
t;x.i+ I
(m+l121 Di+ljf~i+lj
it\T(tn)
,,~F~(m+l12)
m) + ~ ; i + l j ) "~'i+ll2j l~m)
+
((1-ff)/k) FdPij(") and Qa, ij = Qa, ij('t') = ~ij + (Odk(m)) Fij where
k(,,) = k(,,-u2)
~v(m)
t x ~ ' ) + (1 - O0VO ' ' m )
(18a)
,a-,.On+l/2)
Oij tin+u2) and sij, Qa;iy are modified correspondingly, with ~¢(m+I/2) a W <"+m) + (1 - o 0 ¢ " )
Using eq.(16), after a few manipulations we arrive at the following block three-diagonal system with only interface fluxes as unknowns: ~ i j m ) .a-(m+ I/2) ± 'qt'i'l/2,J ~ (
1. W h e n the m e s h size is small, a wise convergence strategy may require a relatively large number of inner iterations per each outer iteration. If the node size is comparable to that of a LWR fuel assembly, however, it is convenient to update the fission source and the eigenvalue as often as possible, for example at every inner iteration or even, as in the present method, at every inner semi-iteration. This means that, in our case, the outer and inner (semi-) iterations coincide. Thus, Oij (") is to be identified with Oij (') in eq.(15a), so that sij = sij(m) =
k(m+l/2) = k (m)
Di+Ij + hx,i+l r~ . .~(nl+l/2) ~i+1 j)'~'i+312j
=
l--OT-~T~- l
At the next semi-iteration, O//,0 is to be identified with r~ ,J)"tl)(m+ 1/21i+1/2j Di+l.j - ~hx'i+l ~i+1
+(h_.~/Dij _ ~hxiQij + ~
EQ.(I7) holds for internal nodes; obvious changes (eq.(13b)) are to be made for the first node (the leftmost node of the line) and, similarly, for the last node. Once the three-diagonal equations have been solved, by means of a direct inversion (e.g. Gaussian elimination) subroutine and the interface fluxes of each line of nodes computed, the average volume fluxes can be immediately updated, via eq.(16) ( see fig.l) and the scheme is ready to be applied to the columns, along the y-direction. No auxiliary mesh sweeps (e.g. for updating coupling factors) are needed.
Remark
1 ¢,"1~(re+l/2)
= ~-.-i-mj
Di+l,i LI~,I+i,/~r~j
~),j U-lmij r~T) + ~
Fig.1 A line of nodes: o old interface flux @ new interface flux ® volume flux
u,.,~ = ~ + ~ (~") ~ - D}) Q,,,0)
U,.,ij~v
p(m) = ~
(m)
i+~,/'-'i+3t24 = Pij
(17)
Remark
(18b)
2. The procedure for the choice of the to acceleration parameters, although classical (Varga, 1962; Wachspress, 1966) requires some words of comment. An 'equivalent model reactor', i.e. a one-group, bare, homogeneous, rectangular model of the actual (2-D) reactor is first constructed and a rectangular grid of equally spaced points superimposed to it. The eigenvalues and eigenvectors of the two commuting matrices representing the (finite difference) central difference terms with respect to x and y being known (Varga, op.cit., p.213 ff.), a standard minimax procedure allows to determine the
50
B. MONTAGNINI et al.
'optimal' values of cox and toy, which are taken as 'good' ones also for the actual problem.It should he remarked that the spectral intervals to be explored (in order to eliminate the higher order eigenvectors) must obviously exclude the fundamental eigenvalue of both matrices (see the Appendix).
$.
3-D
P R O B L E M S
Sharp estimates of the convergence ratio of the ADI schemes have been obtained, as well known, only for certain simple 2-D finite difference problems. For 3-D CM problems the situation is certainly worse. In order to extend our method to 3-D problems we have therefore adopted a quite prudential, almost totally implicit scheme. Eq.(15a) is replaced by the following one 6 t,,-~ip+I13) - -ip+ll3 ip+ll3) x - ~'~ U't'i-l/'2jI - 2q~ijl + {~i+ll2jl) hxi _~/,~ip+l/3) - 2,~ip4-1/3 + .~ip+l/3)., h2j u't'ijll2"l ,,t.ifl ".t.,i~/+ll2.1)
#'~ip) ,~(~,+1/3) + COip)tI)(~''+l/3)l}l- T-~-I l.,~ijl ~a,ijl "~ijt = 6
ip) ip) ip) ,,,gp)m~.1+ r~:!. s~;) 2 (tYPij't'll2 - 2¢y~ijl + ¢~ij,t+l/2) + -- Y q l ~lii ijr hzl (19) where COO,) is a new acceleration parameter (>_ 0; experience has shown that if tx is not IOO high, a < 0.8 say, then td,") can be pul equal to zero) and Qtp) _ + ...R_ Fi/i a,qt - Zqt ktp) sip) i/t -- ~ )-F i j t
the Gaussian inversion of each line of the layer a r e s t o r e d , to be used for all the remaining iterations of the m-cycle (which imply therefore only a multiplication of a matrix by a vector; k , cox, coy are kept fixed until the cycle is completed). This speeds up the computation of each layer, at the expense of a negligible additional memory requirement. R e m a r k 2. It has been observed that a I-D analytic treatment of the top and bottom reflectors (for details s e e Camiciola et al., 1986) may lead to some improvement of both accuracy and convergence.
6. N U M E R I C A L R E S U L T S Some numerical results concerning the classical 3-D PWR IAEA Benchmark Problem (Wagner, 1975) are shown in fig.2, which reports the errors (%) of channel power (the nodal power integrated along z ) with respect to the reference calculation (Wagner et al., 1977) for a fullassembly calculation (9x9x17 nodes, each one of 20x20x20 cm, with analytic treatment of top and bottom reflectors) as well as a quarter-of-assembly calculation ( 17x 17x 17 nodes, 10x 10x20 cm each). These calculations correspond to the two classes of accuracy referred to in t h e Introduction. Max errors are 6.5% and 1.4%, respectively (they correspond to the underlined figures), whereas the mean (r.m.s.) errors amount to 2.1% and 0.8%, respectively. Such errors are within the error bounds mentioned in sect. 1. With a convergence criterion of 10 -4 on the pointwise fluxes and 3 m-iterations for each cycle, the two computations require 15 and 20 ADI cycles respectively, and a CPU time (1/4 core symmetry) of 24 and 111 seconds, respectively, for the VAX 6400 machine. For the Macintosh IIfx PC the computing time may be 10 times longer. The CPU time of a test run on the CRAY YMP8/464 computer (used as a sequential computer) has been, for the two cases, of 1.85 and 8.28 seconds, respectively.
".vijl"~ip) k = I. 02972 k = I. 02918 k = I. 02903
9x9 17x17 Refer. - 2.9 6.5 - 2.4 - 0.8 0.9 - 0.8 0.476 0,700 0.610 0.0 1.9 2.7 - 0.3 0.1 0.4 0.0 - 0.6 1.177 0,971 0.923 0.866 - 2.0 - 0.6 0.6 0.0 0.9 - 4.3 - 0.1 0.2 0.4 0.0 0.0 - 1.3 1.368 1.309 1.180 1.088 0.999 0,709
k~) is updated according to the following relation kip+l/S) _-- kip)
IP(p+l/3) t~q/(p+I/3) + (I - ~)q.rip)
and the ADI scheme is completed by two similar equations (with a cyclic permutation of the indexes i j , l ), for ~ip+2/3) and ~ip+l) and the corresponding eigenvalues. As eq.(191 represents, for fixed l, a 2-D (discretised) diffusion problem in the xy plane, we only need apply the method of sect.4. We obtain by this way, after a moderate number of m-iterations (eqs.(15a,b)), reasonably converged flux values for the l-th node layer or, otherwise stated, a 'reasonably good' inversion for the 2-D diffusion problem for that layer (the whole 3-D process becomes therefore similar to the classical 2-D scheme, where 'exact' 1-D inversions are alternately performed along the two directions). When all the /-layers have been processed, a second cycle of iterations is started, the nodes being now regrouped according 1o layers parallel to the yz plane. The third cycle, parallel to the xz plane, concludes the scheme. R e m a r k 1. A variant of this scheme is the following one. After the fast m-iteration, the block triangular matrices for
- 0.9
- 3.2 - 0.6 0.729
- 0.9
- 0.6
0.4
- 0.I
0.6
- 0.4
0.1 0.2 0.0 0.4 - 0.1 - 0.2 - 0.5 1.396 1.430 1.290 1.071 1.054 0.975 0.756 - 0.1 0.2 0.6 - 3.6 1.1 1.3 - 1.4 0.5 0.5 0.5 - 0.7 0.1 0.0 - 0.7 1.281 1.421 1.193 0.610 0.954 0.959 0.777 Fig.2.3-D IAEA PWR Benchmark Problem
7. C O N C L U S I O N S A non-conforming and even non-consistent, but nevertheless quite efficient coarse-mesh flux expansion method has been worked out, which can be used as a fast
A well-balanced coarse-mesh flux expansion method running LWR core calculational tool as well as a reasonably accurate method for reference problems.
Acknowledgements- We are indebted with F.Oriolo for many discussions on the numerical part of the method. This work was performed under the auspices of the National Group for Mathematical Physics of the Consiglio Nazionale delle Ricerche, with the help of a contribution of the Ministero della Pubblica Istruzione. An early version of the computer program was supported by the Ente Nazionale Energie Alternative, Italy.
APPENDIX
The itemtive fonnula for the eigenvalue k Let L = V . D V + ~ . The outer iterative process represented by eq.(3a), i.e. _ (L+...II-F)tl~(n+l) = ! - a F t b (n)
k (n)
Bennewitz F., Finnemann H. and Wagner M.R. (1975) Trans.Am.Nucl.Soc. 22, 250 Bonalumi R.A., Rouben B., Dastur A.R., Dondale C.S. and Li H.Y.H. (1978) Proc. of A.N.S. Top. Meeting CONF-78.0401, Gatlinburg, p. 169 Botta E., Camiciola P., Cundari D., Montagnini B., Raffaelli P., Soraperra P. (1987) ANS-ENS
lnt.Top. Meeting on Adv.in Reactor Phys., Math. and Computation, Paris, p.1247 Camiciola P., Cundari D. and Moutagnini B. (1986) Ann. nucl.Energy 13,629 Del Valle E., Hennart J.P. and Meade D. (1985) Proc Int. Meet. on Adv.Nucl.Eng. Comput.Methods, Knoxville, Tenn., p.473 Del Valle E.,Hennart J.P. and Meade D. (1986) Nucl. Sci.Engng. 92, 204 H6bert A. (1987) Ann. nucl. Energy 14, 527 Hennart J.P. (1984) Nodal Methods for the Numerical Solution of Partial Differential Equations, in Lect. Notes in Mathematics, Springer-Verlag, p. 1230 (Proc. Conf. on Numerical Analysis, Guanajuato, Mexico) Hennart J.P. (1986) SIAM J. Sci. Stat. Comp. 2,918 Hennart J.P. Numer.Methods for Partial Diff.Eqs. (1988) 4. 233 Hennart J.P., Jaffre J. and Roberts E. (1988) Numer. Math. :~., 701 Langenbuch S., Maurer W. and Werner W. (1977a) Nucl. Sci. Engng. 63,437 Langenbuch S., Maurer W. and Werner W. (1977b) Nucl. Sci. Engng. 64, 508 Lawrence R.D. and Doming J.J. (1980), Nucl. Sci. Engng. 76, 218 Ougouag A.M. and Rajic H.L. (1988) Nucl. Sci. Engng. 100, 332 Rajic H.L. and Ougouag A.M. (1989), Nucl. Sci. Engng. 103,392 Rohach A.F. (1987)Ann. nucl. Energy 14,653 Shober R,A., Sims R.N. and Henry A.F. (1977) Nucl. Sci. Engng. 64, 582 Trentavizi C. (1985) M. Sci. Thesis, Univ.of Pisa, Italy Varga R.S. (1962) Matrix lterative Analysis, PrenticeHall, Englewood Cliffs, N.J. Wachspress E.L. (1966) lterative Solution of Elliptic Systems, Prentice-Hall, Englewood Cliffs, N.J. Wagner M.R.(1975) Proc. Conf. on Comp. Meth. in Nucl. Eng. CONF-75.0413, Charleston,J., 1 Wagner M.R., Finnernann H., Koebke K. and Winter H.J. (1977) Atomkernenergie 30, 129
k (n)
(A 1)
can be given the following equivalent form: ¢b( n + l ) = - ( I +
REFERENCES
51
t~ L ' I F ) ' I 1 - a L - 1 F ~ ( n ) k (n) k (n)
(A2)
Let kh, ~h, be the eigenvalues, resp. the eigenvectors of the operator K = - L "1 F , so that I ~ h = k h ~ h . When referring to the fundamental eigenvalue and eigenvector we shall omit the index: k = k l , C p = ~ 1 . As the fundamental eigenvalue is simple, we can assume that 0 < kh < k (h -> 2).Then, by expanding the initial guess flux ~(o) in the eigenvectors ~h , i.e. ~(o) = Zh ChdPh , eq.(A2) gives, for n = 1, ~O) = _ (I +-fl[-LlF)
"l 1 - ot L - I F ~ ( ° )
k (o)
k (0)
= ~'h Ch (1 - ot)(kh/k(0)) ¢Ph I - ~(k~Ik (°)) where k (°) = i is the initialguess for k. For the second
iterate we have
dP(2) --" ~h Ch ( 1 - ot)(khlk(°)) ( I - ot)(kh/k(1)) ~h 1 - ot(kh/kt°)) 1 - ct(kh/k ~)) and, for the n-th iterate, ~(") = ~ h ch (1 - ot)(kh/k(°)) .. (1 - ot)(kh/k(nt)) ¢p~ 1 - a(kh/k(O))
I - a(kh/k (n'D)
= p(O)"p(n'l)[Cl,+h~2£h ('h/k)n{Y(O'..~(ha-1) ,h] (A3) where we have introduced the following abbreviations:
p(m) = p(ra)(Of) = (1 - oD(k/k (m))
1 - ot(klk (m)) O.(hm)= o.~m)(o[) = 1 - ~'(]Gtk (m>) 1-
of(kh/k(m))
We observe that the inequalities 0 < kh < k (h > 2), together with the inequality 0 < o~(k/k(n)) < I
(A4)
0 < crh(n) < 1
(A5)
imply that Note that a too high value of a could result into a violation of this inequality. Let, for instance, k = 1.02, k2 = 1.01 and take de = 0.99. With k(°) = 1 we get o~(k/k(°)) = 1.0098 and o'2(°)(0.99) = -98, with disastrous
B. MONTAGN1NI et al.
52
consequences.We shall assume in the following that (A4) is always fulfilled. Then we have, for any vector w with non negative components, using (A4) and the relations K~b
= k¢~, K~h = kh % (h > 2),
where use has been made of eq.(A1). Note that the last quotient represents the appropriate production-to-loss ratio at the (n+l)-th iteration. Now LTe = (V.DV + Z T) • = zTe = (-Zr + •T) e-- - Za
(K~(n),w) (in fact, - Zr~ * Z ~ = t Z ~ - ~ ' = - Za~ ).Then we have
(~(n),w)
Cl(¢P,w) + ~
ch (kh/k) n+l cr~°) ''
crh("'l)(~,w)
(- L O ("+l), e) = (O (n+t), - LTe) = (~(n+t), Za)
h>2
=k
c , ( ~ , w ) + Y~ ~
( F(b (n+t), e) = ( L L - t F ~ ("+l), e)
(kak)",,~o, ..o~,-,)(,~,w)
h>2
= (- L-tFO (n+t), - Lre) = ( KO ("+1), Za )
Here
(u,v)=
f(~ ugvg)dV Jv g = l
(when the discrete form of the equations is considered, the integral is replaced by a s u m over the nodes). Then, on account o f (A5), the inequalities 0 < kt, < k (h >2) and (~,w) > 0 (where the latter inequality follows from the positivity of the fundamental eigenvector ~ and the nonnegativity of w) we conclude (KO(n)'w)
*k
(A6)
(~("),w)
k ("+t) = k (n) (~(n+l), w) (¢'("), w) (w > 0), then we can easily show, using (A3), that k(n+l) k (")p(") = k (1 - t~)/[1 - ot(klk (n))] or, in the limit as n -> oo 1 - ot(klk(")) = (1 - t~)(k/k("+l)), a relationship which can be written, putting Mk(m) = 1 + ~,n), -
as n -> oo. Let now v2;f, Za be two G-vectors with c o m p o n e n t s v2;/l(r) ..... v2;fe(r) and 2~al(r),...,2;ae(r), respectively, ,~ag being the absorption cross section in group g . W e show that the r.h.s, of the following recursive formula, which is clearly equivalent to eq.(3b) of sect.2,
k ("+t) = k(")
and substitution into eq.(A8) leads to the conclusion. Eq.(A6), with w = Za, then shows that /¢(n+l).> k as n -> oo. Note that other recursive formulas, although asymptotically valid, m a y introduce u n n e c e s s a r y restrictions on the parameter ~t. If we choose, for instance, the following formula
~(n+l) __ _ 1 ~O~ ?(")
(A9)
O f course, we would like to have 7(") -> 0. But (A9) shows that this is not true, in general, if f f ~ 1/2.
(°("+l)'vZf)
&'((I ~("+I),VZf) + (I - £Z)(dP(n),VZf) (A7) can be written in the form (Kq~(n+l),w)/(~("+l),w) ( i.e. in the same form as the l.h.s, of eq.(A6)), with w = Za. Let e be the G-vector the components of which are all equal to one, over the whole reactor domain. It is immediately seen that FTe = vet, since the sum of the Xg is the unity. Then eq.(A7) can be given the following form
(dp("+1), •Te)
k (n+l) = k (")
C ~ ("+1), FTe) + (1 - C~)(~(n), FTe) ( F ~ (n+l), e) ( a F~("+i) + l - f f F c b ( n ) , e ) k(") k (") ( F ~ ("+l), e) (- L(b ("+1), e)
(AS)
Choice o f the parameters ~
toy
Let us assume that the approximate eigenvalues k(n) have already converged to a sufficiently accurate value, k*.The iterative process described by eqs.(15a, b) o f sect.4 requires a slight modification of the classical procedure of the cyclic variation of the co-parameters. The modification is necessitated by our use of only one inner iteration pex outer iteration, which implies that we are dealing with the acceleration of the power method for an eigenvalue problem rather than the acceleration of the inner iteration cycle for an inhomogeneous linear system. As previously said, a rectangular, bare, homogeneous equivalent 'model reactor' is first constructed and a rectangular grid x~ =/hx, yi=jhy (i = 0,1,...,nx, j = 0,1,...,ny), with approximately the same m e s h width as the actual problem is then superimposed over the domain {x,y; O.~x
A well-balanced coarse-mesh flux expansion method _4.- r a,(n+l)
where use has been made of eq.(Al2) and, similarly,
2 O ( n + l ) + a,(n+l) x
h~ ~v,i-1:2,i-
53
v,i+v2jj
+A ~yy ._(n+l) (~)i,]-tl2 - ~ - _i(jn + l ) + --(n+l) . + ...~ ~i,]+I/2) (,7.~--~
~,)~(n+l)
~)
+ ~-
a(F'/k*) + ~.~+/~
-
(A14)
(o~n) + Z - a(Y/k*) + / ~
~!~)
- ~ -
k*
-,',1
(A10)
with the condition of zero flux at the boundary Using vector notations, the ADI procedure for solving each step of this iterative process (Varga, 1962 or Wachspress, 1966) can be given the following form:
one sees that In has eigenvalues 0p4(n) = 0,1u(n). 0~,~(n). For a oJ-cycle with n* values of the ~a (~),coy(n) parameters, we may define a global iterationmatrix
n* with eigenvalues
J* = H In n=l
O;q -- fi 8p,q (n) n=l
--[-.
From (AI3)(AI4) w e o b t a i n 01A (n) = 0&l (n) . 0~,1 (n) -- l [- X +
~
~
(n+l/2)
where ~(n) has components ~i~(n)q(suitably ordered) and X~(n) y~(n) represent the central difference terms with respect to x and y, respectively. It is known (Varga, op.cit,, p.214) that the vectors ¢(P.q) with components ¢iff "q) = C O''q) sin (pax ilia) sin (qT~'yj Ily) = C(P'q ) sin (pgilnx) sin (qajlny) are eigenvectors of both X and Y, with eigenvalues 37,= 2[1..cos(prclnz) ]lha 2 and #q = 2[ l-cos(qftlny)]/hy 2, respectively. In particular, {~(1,1) is the fundamental eigenvector of X and Y and therefore of X + Y and represents the solution of the model problem. Substitution of ~O.~) into eq.(A10), as 'converged flux', leads immediately to the criticality equation:
~.l + I.q + Z - (o~/k*)F = (1 - o:)~'/k*
(AI2)
Now, the ADI iterative process, eqs.(A10), can be summarised l~y the single relationship
(1)(n+l)= ~'n(I}(n) where the iterationmatrix. Jn, has the following form
• [ x+ (o: +
+
Then, setting
o;(,)
o6") + (1- a)(P/k*) - X~ ~ " ) + x- a(P/~*) +
(0~n) + T. - a('ff/k*) +
),~
(A13)
for all n and therefore 0* m -- 1. If n* is large enough and we were able to determine the co-parameters in such a way that, for any p,q ~ 2, Ùx~(") = 0 or 0a,q(n)= 0 for some n, then the process would give the fundamental eigensolution ¢~O,t)exactly, in a finite number of steps. This result is, however, of no interest, since ¢~OA)is explicitly known. Rather, we try to solve, by the classical techniques (Varga, op.cit., Wachspress, op.cit.) the following minimax problem
min ~.,~
"" c~') + T. - a(':/k*) + ~h + I~ " x max 1-'[ a,~ " = ~ c4 n) + ~. - a('f/k*) + x eof")+ Z - o~(~/k*) + 3,i + gl - Y o V ) + r- -
a(':/k*)+ y
with x.y ranging over an interval containing all the ~, resp. #q, with p.q > 2, but excluding 3,1 and #l. The optimal fo's so determined are then used as 'good' acceleration parameters for the actual problem to be solved.