- Email: [email protected]
A unified formalism for spatial discretization schemes of transport equations in slab geometry
Annals of Nuclear Energy, Printed in Great Britain
Vol. 8, pp. 671 to 681, 1981
0306-6549/81/I 10677-05$02.00/0 Pergamon Press Ltd
A UNIFIED FORMALISM FOR SPATIAL DISCRETIZATION SCHEMES OF TRANSPORT EQUATIONS IN SLAB GEOMETRY J. P. HENNART* Institute de Jnvestigaciones en MatemLticas Aplicadas y en Sistemas, Universidad National Aut6noma de Mkxico, Apartado Postal 20-726, Mkxico is shown in this paper that most of the spatial discretization schemes of transport equations in slab geometry which have been developed recently are particular applications of a general finite element oriented formalism developed by this author and his collaborators for the numerical integration of systems of
Abstract--It
stiff ordinary differential equations.
subject to the initial condition
1. INTRODUCTION
engineering literature, a In the recent nuclear considerable amount of work has been devoted to the development of improved spatial discretization schemes for discrete ordinates transport equations in slab geometry (Alcouffe et al., 1979 ; Gopinath et al., 1980; Larsen and Miller, 1980; Lee and Vaidyanathan, 1980). Independently, in previous works by this author and his collaborators (Hennart, 1977; Hennart and England, 1979; Hennart and Gourgeon, 1981). a general theory has been developed for the numerical integration of initial value ordinary differential equations by piecewise continuous methods providing a piecewise polynomial and (or) exponential behaviour in the independent variable. After recalling in Section 2 the main results of this finite element oriented general formalism, we show in Section 3 how most of the recently explored spatial discretization schemes can be seen as particular applications of this formalism, leading to some interesting reinterpretations of some well-known schemes. In Section 4 finally, some possible extensions to more complicated situations are briefly mentioned.
Q%%)= $0.
(2)
For the sake of simplicity, let n : xi = ih, i = 0,. , N where h = (xN- x,)/N, be a uniform mesh of mesh-size h, although a variable mesh could have been considered without affecting our arguments. Integration of equation (1) over [xi, xi+ 1] gives ti(xi+ I)-$(xi)
=
XI”‘f($,x) dx, s
which will be approximated
(3)
by
Xi+1 tii+l-tii=
x,
fWpp-4
(4)
dx,
1, where p, q are non-negative
integers with p + q = I + 1. Let Sk = Sh[xi, xi+ 1] be the linear space S*={x”exp(l,x)In=O ,..., N,--I; m= l,..., M; C, N, = p + q}, where the d,‘s are given real numbers, with Li # 1, for i # j, one of them being eventually zero. The function $p4 in equation (4) will be chosen in S” subject to the (p + q) interpolation conditions $:4(X{)
=
$y’,
r =
0,.
,p-
1,
S=O,...,q-l.
$~~(Xi+l)=$~?l,
(5a) (Sb)
2. FINITE ELEMENT ORIENTED DISCRETIZATION SCHEMES FOR INITIAL VALUE PROBLEMS
Let us consider the following non-linear differential equation
first-order
The existence and uniqueness of this generalized Hermite interpolant is guaranteed by Theorem 1 of Hennart and Gourgeon (1981). In virtue of equation (5), $,,(x) can be rewritten as 9-l
p-1 tip&4
=
c
~l”G4
I.=0
+
1
Hi,
v,(x),
(6)
S=O
where u,, v, E S* and verify *Also consultant to “Comisibn National de Seguridad Nuclear y Salvaguardias”.
@(Xi) = (srl, r,2=0 ,..., p-l, 677
(74
678
J. P.
u(fyx.+l) = 0, r = 0 I
I
,...,
p-l;
I=0
,..., q-l,
HENNART
(7b)
and
exist constants C, j = 0,. . . ,I, such that IIJIQ(x)-#$(x)ll ~oo[xj.xi+ 11G Cjh’+l-j.
#‘(Xi) = 0, s = 0 ,..., q-l;
I = 0,...) p-l,
u(l)(x.+1) = C&r, s,z = 0 ,..., q-l. s I
(7c) (74
Finally, the successive derivatives @, r = 1,. . . , p- 1 (if p>2) and I#~, s=l...,q-1 (if 422) are eliminated by collocation at x = xi and x = xi + 1using equation (I), namely $,p’= f”-‘)(Jli,xi),
r = 1,...,p-1
(withp > 2)
(8a)
and
(8b; assuming that f is smooth enough. As a consequence of the collocation conditions in equation (8), there remain at most two parameters in the expression (6) for JI, and hence in equation (4), namely tii z I(Iw(xi) and$,+ I = tiw(xi+ l).If$iappearsexplicitly,it iseither known from the initial condition (i = 0) or from integration over the previous interval [xi- i, xi], i=l ,...,N-l,while$i+l is obtained from equation (4). A computationally convenient form for equation (4) implies some numerical quadrature scheme (QS) over [xi,xi+ i] unless the integral may be evaluated analytically. As a result of the numerical or analytical integration, equation (4) becomes a non-linear equation in the only unknown Jli+ 1. Equation (4) thus defines a class of discrete one-step integration methods depending on which quadrature formula is used. In Hennart and Gourgeon( 1981), the existence of a unique solution to equation(4) is proved for sufficiently small h. Moreover, discrete and continuous error bounds are derived for the methods presented above in Theorems 2 and 3, respectively of Hennart and Gourgeon (198 1). These theorems are quoted here without demonstration :
(10)
The readers interested in detailed demonstrationsin the general case where S” has polynomial and exponential basis functions should consult Hennart and Gourgeon (1981). The simpler case where S* = P,, the linear space ofpolynomials ofdegree less than or equal to 1,is treated in Hennart (1977). In Sections 4 and 5 of Hennart and Gourgeon (1981), the above schemes are compared with existing schemes ofthe Runge-Kutta type(seee.g. Hennart and England, 1979) or with schemes exhibiting exponential fitting capabilities : it is shown in particular that as long as the quadrature rule used is exact for members of Sh, the schemes are exponentially fitted or order N, - 1 at A,,,, m = l,... , M (Bjurel, 1972). An automatic exponential fitting at I (leading to a scheme which is exact if the solution behaves as exp(h)) can in particular be obtained after the following manipulations. First multiply equation (1) by I#J= exp( - 1x) to get
or equivalently
Ifequation (12)instead ofequation (1) is integrated over [xi>xi+ 1] and JI is replaced as above by $, in the r.h.s. of the resulting expression, it is easy to realize that if $ behaves like exp@x), it will be integrated exactly, independently of Sh and of the particular quadrature scheme chosen as long as it is acceptable (i.e. exact for members ofS’). In general, q5can be any function ofx and
not only exp ( -Lx).
If the emphasis is on the order rather than on the exponential fitting capabilities of a given scheme, quadrature schemes exact for P, even when Sk f P, are (a) Discrete error bounds (Theorem 2 of Hennart and perfectly sufficient. More details and in particular Gourgeon (198 1)) general stability considerations can be found in Assume that f (J/,x) is of class CL, k 3 I, over Hennart and Gourgeon (1981). All the schemes proposed are definitely finite element oriented in the Rx[x,,x,], and let the discrete one-step integration method be defined as abovefor $I, E Sh. Then there exists sense that they provide a piecewise continuous a constant C such that behaviour of the approximation with respect to the independent variable, with comparable discrete and ]$(x~)-$~] c Ch’+‘, i = l,..., N, (9) continuous error bounds. In other words, no under the hypothesis that the numericalquadrature rule is superconvergence results at some discrete points are exact for members of S”. claimed, which would lead to finite difference oriented schemes. (b) Continuous error bounds (Theorem 3 of Hennart ana’ As a final comment, a given scheme is completely Gourgeon (1981)) specified if the pair (p,q) the space Sh, the weight With the hypotheses of the previous Theorem, there function 4 and the quadrature scheme QS are defined.
679
Spatial discretization schemes-a unified formalism
At this point, it is convenient to introduce. some notation : in the following Gn, Ln, Rn will denote Gauss, Lobatto and Radau right-hand n-point formulae.
each spatial cell Pil= p,
3xi = txi- 112+xi+ I/J27 l=O,l,...,
3. SPATIAL DISCRETIZATION TRANSPORT
EQUATIONS
SCHEMES OF
(16)
where P,(x) is the usual lth Legendre polynomial, then
IN SLAB GEOMETRY
If for the sake of simplicity the group and directional subscripts are dropped in the standard multigroup slab geometry discrete ordinates transport equations (Carlson and Lathrop, 1968), these read as
h Pij(X)Pik(X)dx = 6>12k+l’
(17)
and S(x) can be written as S(x) = f
SilPil(x), xE [Xi- 1/2,xi+ 1/J.
I=0
(13)
(18)
where Assuming that p > 0, particles flow from left to right and it is appropriate to solve the transport equation from left to right across a spatial cell [xi _ i,t, xi + I,21 : in that direction, the corresponding equation is actually stable. In general, we can assume that tii- 1,2 has been determined as the exiting flux from the previous cell Cxi-3/2,xi- &J and that S(x), which is in fact a function of $(x), is known from some previous iteration. Knowing $i- ,,2 and S(x), we want to determine the exiting flux tii+ 1,z and the cell average flux eti With p > 0, we can rewrite equation (13) as
4W = -.dX
-!y(x) +y,
(14)
which is in the form of equation (1) to which the formalism of Section 2 can be applied. In all the schemes developed up to now, a spatially constant cross-section set is associated to each interval so that a(x) is in fact replaced by some cell average gi. In the next section, we shall see how this restriction can easily be removed. Let us now survey most of the schemes recently explored (Alcouffe et al., 1979; Gopinath et al., 1980; Larsen and Miller, 1980; Lee and Vaidyanathan, 1980) and show that they are particular cases of our general formalism. (a) The Diamond LX&-ewe scheme (DD)
The DD scheme (Carlson and Lathrop, 1968) is a scheme with S” = P1 and 4 = 1. Assuming that 4~) = ci, the integration in equation (4) can be performed analytically leading to equation (13) of Larsen and Miller (1980), namely (p,q) = (1,l)
21+
1
S(X)PiXX)dx.
Sil = __ h
(19)
In particular 1
Xi+10
sio = h s xi_ *,*
S(x) dx.
(20)
Clearly, the DD scheme is a second-order scheme from equation (9) since p + q = l+ 1 = 2. (b) The Weighted Diamond D@rence
scheme (WDD)
The WDD scheme (Carlson, 1976; Lathrop, 1969) has the same characteristics as the DD scheme, the only difference being that Sh = P, = { 1, x> is replaced here by (1, exp@x)). From equation (14), we get $‘i+1/2-$‘i-l/2=
-zS$‘(x)dx+%Sio,
(21)
where using an L2 QS correct for Sh
s
~(X)dx=hC~i-,,~(1_8(ah))+~i+1/2B(Lh)l. (22)
In equation (22) e(z) = (1 + z - eZ)/(z(1 - e’)),
(23)
so that as 1 varies from 0 to -co, 0 varies from 0.5 (i.e. the DD scheme) to 1.0. In the classical presentation of the WDD scheme, f3 is a fixed number comprised between 0.5 and 1.0. Here, 0 appears as a function of h, 1 being fixed, and this reinterpretation of the WDD scheme leads to a truly second-order scheme (p+q = 2). (c) The Step Characteristic scheme (SC)
where zi = o,h/p while S, is the first coefficient in the Legendre polynomials expansion of S(x). Explicitly if we introduce the modified Legendre polynomials pi, in
The SC scheme (Lathrop, 1969) is simply a particular WDD scheme with 1 = -o& in the same way as the DD scheme can be seen as a particular WDD scheme with 1= 0. The SC scheme again is of second order.
J. P. HENNART
680 (d) The Linear Discontinuous scheme (LD)
The LD scheme is presented in Alcouffe et al. (1979) as a finite element method where $ and S are taken to be piecewiselinear with jumps ordiscontinuities at the cell edges. Within our general formalism, it can be reinterpreted as a (1, 2) scheme with S” = P, = { 1, x,x’} and Cp= x-xi_ 1,2. With a linear weight function, the source contribution retains not only the SiOterm but also the Sii one. The approximation of $ over a cell is parabolic, namely ICIN $12 = ~10’1,2~~(~)+~10?1,2~~(~) + til:‘l,Zr,(x)Y
(24)
where by collocation
Consistent with the fact that only the first two terms of the Legendre expansion (equation (18)) of S(x) are retained s(xi+ i/Z) N si0 + sil 7
4. EXTENSIONSOF THR PROPOSED FORMALISM
(26)
and a QS like R2 (exact for PZ) leads to equation (14) of Alcouffe et al. (1979), namely
.
amounts to truncating S(x) in equation (28) to its first term in the general Legendre expansion equation (18). If we keep two terms in the same expansion, we get the LC scheme. Three terms would lead to the quadratic method described in Larsen and Miller (1980). If we look at these truncations with the numerical quadrature point of view in mind, it is easy to realize that truncating to say n terms leads to a scheme which is of order 2n. This of course heavily depends on the fact that u(x) is replaced by a cell average value cti In this section, we believe to have shown that most space discretization schemes recently mentioned in the nuclear engineering literature are particular cases of a general formalism that we have developed for general non-linear ordinary differential equations (Hennart, 1977 ; Hennart and England, 1979 ; Hennart and Gourgeon, 1981). To be fair, we should say that some schemes like the exponential method (Barbucci and Di Pasquantonio, 1977) still remain outside of our formalism in its present form.
. (27)
As predicted by equation (9), the LD scheme is a thirdorder scheme. It is interesting to point out that the interpretation given here of the LD scheme gives rise to a continuous, piecewise parabolic approximation to $, and not to a discontinuous, piecewise linear one. We believe that our approach is more natural while consistent with the obtained orders.
As it should be clear from the previous section, the general formalism that we have developed is far from being completely exploited and many existing schemes (see Section 5 of Hennart and Gourgeon (1981)) have not yet even been tried with discrete ordinates transport equations. Our feeling is that maybe the most interesting developments consist in extending the applicability of existing schemes. Let us take a simple example to illustrate what we have in mind, namely consider the WDD scheme when a(x) is any function of x. In that case, equation (21) is replaced by tii+
112
-$i-,,z
= -;
o-(x)$(x) dx+;
Si,.
(29)
(e) The Linear Characteristic scheme (LC)
With the LC scheme (Alcouffe et al., 1979 ; Gopinath et al., 1980), another feature is encountered: here 4 = exp(oix/p) with the result that after integration over a cell, we obtain Xi+1,* $‘i+ i/2 = exp(-ci)i/‘-
112+ $
s X,-L,2
x exp - ;(xi+ {
I/* -x)
S(x) dx.
(28)
1
Different schemes can be obtained depending on how weintegrate thesourceterminther.h.s.ofequation(28). The step characteristic scheme presented above
From the results of Section 2 (Theorem 2 of Hennart and Gourgeon (1981)), any Gl or L2 QS correct for Sh = {1, exp 1x} is consistent with the second order we are expecting from such a scheme, so that general WDD-like schemes are easily derived where a(x) can be any space-dependent function. The same is true for any (p, q) pair and any space Sh. Such extensions as well as extensions to multidimensional situations will be considered in a forthcoming paper. Acknowledgement-The author would like to thank W. L. Filippone for some interesting discussions which constituted the initial motivation for this investigation.
Spatial discretization schemes-a
unified formalism
681
REFERENCES Alcouffe R. E., Larsen E. W., Miller W. F. Jr and Wienke B. R. (1979) Nucl. Sci. Engng 71, 111. Barbucci P. and Di Pasquantonio F. (1977) Nucl. Sci. Engng 63, 179.
Bjurel B. (1972) BIT 12, 142. Carlson B. G. (1976) Rucl. Sci. Engng 61,408. Carlson B. G. and Lathrop K. D. (1968)Computing Methods in Reactor Physics (Eds Greenspan H., Kelber C. N. and Okrent D.), pp. 167-266. Gordon & Breach, New York. Gopinath D. V., Natarajan A. and Sundararaman V. (1980) Nucl. Sci. Engng 75, 181.
Hennart J. P. (1977) Maths Comput. 31,24. Hennart J. P. and England R. (1979) Preliminary VersionWorking Papers.for the 1979 Signum Meeting on Numerical Ordinary fX@rential Equations (Ed. Skeel R. D.), pp. 33.1-
33.4. Department of Computer Science, University of Illinois, Urbana-Champaign. Hennart J. P. and Gourgeon H. (1981) Comun. TPc. Serie NA, No. 263, IIMAS-UNAM. Larsen E. W. and Miller W. F. Jr (1980)Nucl. Sci. Engng 73,76. Lathrop K. D. (1969) J. comput. Phys. 4,475. Lee S. M. and Vaidyanathan R. (1980) Nucl. Sci. Engng 76,l.
Vol. 8, pp. 671 to 681, 1981
0306-6549/81/I 10677-05$02.00/0 Pergamon Press Ltd
A UNIFIED FORMALISM FOR SPATIAL DISCRETIZATION SCHEMES OF TRANSPORT EQUATIONS IN SLAB GEOMETRY J. P. HENNART* Institute de Jnvestigaciones en MatemLticas Aplicadas y en Sistemas, Universidad National Aut6noma de Mkxico, Apartado Postal 20-726, Mkxico is shown in this paper that most of the spatial discretization schemes of transport equations in slab geometry which have been developed recently are particular applications of a general finite element oriented formalism developed by this author and his collaborators for the numerical integration of systems of
Abstract--It
stiff ordinary differential equations.
subject to the initial condition
1. INTRODUCTION
engineering literature, a In the recent nuclear considerable amount of work has been devoted to the development of improved spatial discretization schemes for discrete ordinates transport equations in slab geometry (Alcouffe et al., 1979 ; Gopinath et al., 1980; Larsen and Miller, 1980; Lee and Vaidyanathan, 1980). Independently, in previous works by this author and his collaborators (Hennart, 1977; Hennart and England, 1979; Hennart and Gourgeon, 1981). a general theory has been developed for the numerical integration of initial value ordinary differential equations by piecewise continuous methods providing a piecewise polynomial and (or) exponential behaviour in the independent variable. After recalling in Section 2 the main results of this finite element oriented general formalism, we show in Section 3 how most of the recently explored spatial discretization schemes can be seen as particular applications of this formalism, leading to some interesting reinterpretations of some well-known schemes. In Section 4 finally, some possible extensions to more complicated situations are briefly mentioned.
Q%%)= $0.
(2)
For the sake of simplicity, let n : xi = ih, i = 0,. , N where h = (xN- x,)/N, be a uniform mesh of mesh-size h, although a variable mesh could have been considered without affecting our arguments. Integration of equation (1) over [xi, xi+ 1] gives ti(xi+ I)-$(xi)
=
XI”‘f($,x) dx, s
which will be approximated
(3)
by
Xi+1 tii+l-tii=
x,
fWpp-4
(4)
dx,
1, where p, q are non-negative
integers with p + q = I + 1. Let Sk = Sh[xi, xi+ 1] be the linear space S*={x”exp(l,x)In=O ,..., N,--I; m= l,..., M; C, N, = p + q}, where the d,‘s are given real numbers, with Li # 1, for i # j, one of them being eventually zero. The function $p4 in equation (4) will be chosen in S” subject to the (p + q) interpolation conditions $:4(X{)
=
$y’,
r =
0,.
,p-
1,
S=O,...,q-l.
$~~(Xi+l)=$~?l,
(5a) (Sb)
2. FINITE ELEMENT ORIENTED DISCRETIZATION SCHEMES FOR INITIAL VALUE PROBLEMS
Let us consider the following non-linear differential equation
first-order
The existence and uniqueness of this generalized Hermite interpolant is guaranteed by Theorem 1 of Hennart and Gourgeon (1981). In virtue of equation (5), $,,(x) can be rewritten as 9-l
p-1 tip&4
=
c
~l”G4
I.=0
+
1
Hi,
v,(x),
(6)
S=O
where u,, v, E S* and verify *Also consultant to “Comisibn National de Seguridad Nuclear y Salvaguardias”.
@(Xi) = (srl, r,2=0 ,..., p-l, 677
(74
678
J. P.
u(fyx.+l) = 0, r = 0 I
I
,...,
p-l;
I=0
,..., q-l,
HENNART
(7b)
and
exist constants C, j = 0,. . . ,I, such that IIJIQ(x)-#$(x)ll ~oo[xj.xi+ 11G Cjh’+l-j.
#‘(Xi) = 0, s = 0 ,..., q-l;
I = 0,...) p-l,
u(l)(x.+1) = C&r, s,z = 0 ,..., q-l. s I
(7c) (74
Finally, the successive derivatives @, r = 1,. . . , p- 1 (if p>2) and I#~, s=l...,q-1 (if 422) are eliminated by collocation at x = xi and x = xi + 1using equation (I), namely $,p’= f”-‘)(Jli,xi),
r = 1,...,p-1
(withp > 2)
(8a)
and
(8b; assuming that f is smooth enough. As a consequence of the collocation conditions in equation (8), there remain at most two parameters in the expression (6) for JI, and hence in equation (4), namely tii z I(Iw(xi) and$,+ I = tiw(xi+ l).If$iappearsexplicitly,it iseither known from the initial condition (i = 0) or from integration over the previous interval [xi- i, xi], i=l ,...,N-l,while$i+l is obtained from equation (4). A computationally convenient form for equation (4) implies some numerical quadrature scheme (QS) over [xi,xi+ i] unless the integral may be evaluated analytically. As a result of the numerical or analytical integration, equation (4) becomes a non-linear equation in the only unknown Jli+ 1. Equation (4) thus defines a class of discrete one-step integration methods depending on which quadrature formula is used. In Hennart and Gourgeon( 1981), the existence of a unique solution to equation(4) is proved for sufficiently small h. Moreover, discrete and continuous error bounds are derived for the methods presented above in Theorems 2 and 3, respectively of Hennart and Gourgeon (198 1). These theorems are quoted here without demonstration :
(10)
The readers interested in detailed demonstrationsin the general case where S” has polynomial and exponential basis functions should consult Hennart and Gourgeon (1981). The simpler case where S* = P,, the linear space ofpolynomials ofdegree less than or equal to 1,is treated in Hennart (1977). In Sections 4 and 5 of Hennart and Gourgeon (1981), the above schemes are compared with existing schemes ofthe Runge-Kutta type(seee.g. Hennart and England, 1979) or with schemes exhibiting exponential fitting capabilities : it is shown in particular that as long as the quadrature rule used is exact for members of Sh, the schemes are exponentially fitted or order N, - 1 at A,,,, m = l,... , M (Bjurel, 1972). An automatic exponential fitting at I (leading to a scheme which is exact if the solution behaves as exp(h)) can in particular be obtained after the following manipulations. First multiply equation (1) by I#J= exp( - 1x) to get
or equivalently
Ifequation (12)instead ofequation (1) is integrated over [xi>xi+ 1] and JI is replaced as above by $, in the r.h.s. of the resulting expression, it is easy to realize that if $ behaves like exp@x), it will be integrated exactly, independently of Sh and of the particular quadrature scheme chosen as long as it is acceptable (i.e. exact for members ofS’). In general, q5can be any function ofx and
not only exp ( -Lx).
If the emphasis is on the order rather than on the exponential fitting capabilities of a given scheme, quadrature schemes exact for P, even when Sk f P, are (a) Discrete error bounds (Theorem 2 of Hennart and perfectly sufficient. More details and in particular Gourgeon (198 1)) general stability considerations can be found in Assume that f (J/,x) is of class CL, k 3 I, over Hennart and Gourgeon (1981). All the schemes proposed are definitely finite element oriented in the Rx[x,,x,], and let the discrete one-step integration method be defined as abovefor $I, E Sh. Then there exists sense that they provide a piecewise continuous a constant C such that behaviour of the approximation with respect to the independent variable, with comparable discrete and ]$(x~)-$~] c Ch’+‘, i = l,..., N, (9) continuous error bounds. In other words, no under the hypothesis that the numericalquadrature rule is superconvergence results at some discrete points are exact for members of S”. claimed, which would lead to finite difference oriented schemes. (b) Continuous error bounds (Theorem 3 of Hennart ana’ As a final comment, a given scheme is completely Gourgeon (1981)) specified if the pair (p,q) the space Sh, the weight With the hypotheses of the previous Theorem, there function 4 and the quadrature scheme QS are defined.
679
Spatial discretization schemes-a unified formalism
At this point, it is convenient to introduce. some notation : in the following Gn, Ln, Rn will denote Gauss, Lobatto and Radau right-hand n-point formulae.
each spatial cell Pil= p,
3xi = txi- 112+xi+ I/J27 l=O,l,...,
3. SPATIAL DISCRETIZATION TRANSPORT
EQUATIONS
SCHEMES OF
(16)
where P,(x) is the usual lth Legendre polynomial, then
IN SLAB GEOMETRY
If for the sake of simplicity the group and directional subscripts are dropped in the standard multigroup slab geometry discrete ordinates transport equations (Carlson and Lathrop, 1968), these read as
h Pij(X)Pik(X)dx = 6>12k+l’
(17)
and S(x) can be written as S(x) = f
SilPil(x), xE [Xi- 1/2,xi+ 1/J.
I=0
(13)
(18)
where Assuming that p > 0, particles flow from left to right and it is appropriate to solve the transport equation from left to right across a spatial cell [xi _ i,t, xi + I,21 : in that direction, the corresponding equation is actually stable. In general, we can assume that tii- 1,2 has been determined as the exiting flux from the previous cell Cxi-3/2,xi- &J and that S(x), which is in fact a function of $(x), is known from some previous iteration. Knowing $i- ,,2 and S(x), we want to determine the exiting flux tii+ 1,z and the cell average flux eti With p > 0, we can rewrite equation (13) as
4W = -.dX
-!y(x) +y,
(14)
which is in the form of equation (1) to which the formalism of Section 2 can be applied. In all the schemes developed up to now, a spatially constant cross-section set is associated to each interval so that a(x) is in fact replaced by some cell average gi. In the next section, we shall see how this restriction can easily be removed. Let us now survey most of the schemes recently explored (Alcouffe et al., 1979; Gopinath et al., 1980; Larsen and Miller, 1980; Lee and Vaidyanathan, 1980) and show that they are particular cases of our general formalism. (a) The Diamond LX&-ewe scheme (DD)
The DD scheme (Carlson and Lathrop, 1968) is a scheme with S” = P1 and 4 = 1. Assuming that 4~) = ci, the integration in equation (4) can be performed analytically leading to equation (13) of Larsen and Miller (1980), namely (p,q) = (1,l)
21+
1
S(X)PiXX)dx.
Sil = __ h
(19)
In particular 1
Xi+10
sio = h s xi_ *,*
S(x) dx.
(20)
Clearly, the DD scheme is a second-order scheme from equation (9) since p + q = l+ 1 = 2. (b) The Weighted Diamond D@rence
scheme (WDD)
The WDD scheme (Carlson, 1976; Lathrop, 1969) has the same characteristics as the DD scheme, the only difference being that Sh = P, = { 1, x> is replaced here by (1, exp@x)). From equation (14), we get $‘i+1/2-$‘i-l/2=
-zS$‘(x)dx+%Sio,
(21)
where using an L2 QS correct for Sh
s
~(X)dx=hC~i-,,~(1_8(ah))+~i+1/2B(Lh)l. (22)
In equation (22) e(z) = (1 + z - eZ)/(z(1 - e’)),
(23)
so that as 1 varies from 0 to -co, 0 varies from 0.5 (i.e. the DD scheme) to 1.0. In the classical presentation of the WDD scheme, f3 is a fixed number comprised between 0.5 and 1.0. Here, 0 appears as a function of h, 1 being fixed, and this reinterpretation of the WDD scheme leads to a truly second-order scheme (p+q = 2). (c) The Step Characteristic scheme (SC)
where zi = o,h/p while S, is the first coefficient in the Legendre polynomials expansion of S(x). Explicitly if we introduce the modified Legendre polynomials pi, in
The SC scheme (Lathrop, 1969) is simply a particular WDD scheme with 1 = -o& in the same way as the DD scheme can be seen as a particular WDD scheme with 1= 0. The SC scheme again is of second order.
J. P. HENNART
680 (d) The Linear Discontinuous scheme (LD)
The LD scheme is presented in Alcouffe et al. (1979) as a finite element method where $ and S are taken to be piecewiselinear with jumps ordiscontinuities at the cell edges. Within our general formalism, it can be reinterpreted as a (1, 2) scheme with S” = P, = { 1, x,x’} and Cp= x-xi_ 1,2. With a linear weight function, the source contribution retains not only the SiOterm but also the Sii one. The approximation of $ over a cell is parabolic, namely ICIN $12 = ~10’1,2~~(~)+~10?1,2~~(~) + til:‘l,Zr,(x)Y
(24)
where by collocation
Consistent with the fact that only the first two terms of the Legendre expansion (equation (18)) of S(x) are retained s(xi+ i/Z) N si0 + sil 7
4. EXTENSIONSOF THR PROPOSED FORMALISM
(26)
and a QS like R2 (exact for PZ) leads to equation (14) of Alcouffe et al. (1979), namely
.
amounts to truncating S(x) in equation (28) to its first term in the general Legendre expansion equation (18). If we keep two terms in the same expansion, we get the LC scheme. Three terms would lead to the quadratic method described in Larsen and Miller (1980). If we look at these truncations with the numerical quadrature point of view in mind, it is easy to realize that truncating to say n terms leads to a scheme which is of order 2n. This of course heavily depends on the fact that u(x) is replaced by a cell average value cti In this section, we believe to have shown that most space discretization schemes recently mentioned in the nuclear engineering literature are particular cases of a general formalism that we have developed for general non-linear ordinary differential equations (Hennart, 1977 ; Hennart and England, 1979 ; Hennart and Gourgeon, 1981). To be fair, we should say that some schemes like the exponential method (Barbucci and Di Pasquantonio, 1977) still remain outside of our formalism in its present form.
. (27)
As predicted by equation (9), the LD scheme is a thirdorder scheme. It is interesting to point out that the interpretation given here of the LD scheme gives rise to a continuous, piecewise parabolic approximation to $, and not to a discontinuous, piecewise linear one. We believe that our approach is more natural while consistent with the obtained orders.
As it should be clear from the previous section, the general formalism that we have developed is far from being completely exploited and many existing schemes (see Section 5 of Hennart and Gourgeon (1981)) have not yet even been tried with discrete ordinates transport equations. Our feeling is that maybe the most interesting developments consist in extending the applicability of existing schemes. Let us take a simple example to illustrate what we have in mind, namely consider the WDD scheme when a(x) is any function of x. In that case, equation (21) is replaced by tii+
112
-$i-,,z
= -;
o-(x)$(x) dx+;
Si,.
(29)
(e) The Linear Characteristic scheme (LC)
With the LC scheme (Alcouffe et al., 1979 ; Gopinath et al., 1980), another feature is encountered: here 4 = exp(oix/p) with the result that after integration over a cell, we obtain Xi+1,* $‘i+ i/2 = exp(-ci)i/‘-
112+ $
s X,-L,2
x exp - ;(xi+ {
I/* -x)
S(x) dx.
(28)
1
Different schemes can be obtained depending on how weintegrate thesourceterminther.h.s.ofequation(28). The step characteristic scheme presented above
From the results of Section 2 (Theorem 2 of Hennart and Gourgeon (1981)), any Gl or L2 QS correct for Sh = {1, exp 1x} is consistent with the second order we are expecting from such a scheme, so that general WDD-like schemes are easily derived where a(x) can be any space-dependent function. The same is true for any (p, q) pair and any space Sh. Such extensions as well as extensions to multidimensional situations will be considered in a forthcoming paper. Acknowledgement-The author would like to thank W. L. Filippone for some interesting discussions which constituted the initial motivation for this investigation.
Spatial discretization schemes-a
unified formalism
681
REFERENCES Alcouffe R. E., Larsen E. W., Miller W. F. Jr and Wienke B. R. (1979) Nucl. Sci. Engng 71, 111. Barbucci P. and Di Pasquantonio F. (1977) Nucl. Sci. Engng 63, 179.
Bjurel B. (1972) BIT 12, 142. Carlson B. G. (1976) Rucl. Sci. Engng 61,408. Carlson B. G. and Lathrop K. D. (1968)Computing Methods in Reactor Physics (Eds Greenspan H., Kelber C. N. and Okrent D.), pp. 167-266. Gordon & Breach, New York. Gopinath D. V., Natarajan A. and Sundararaman V. (1980) Nucl. Sci. Engng 75, 181.
Hennart J. P. (1977) Maths Comput. 31,24. Hennart J. P. and England R. (1979) Preliminary VersionWorking Papers.for the 1979 Signum Meeting on Numerical Ordinary fX@rential Equations (Ed. Skeel R. D.), pp. 33.1-
33.4. Department of Computer Science, University of Illinois, Urbana-Champaign. Hennart J. P. and Gourgeon H. (1981) Comun. TPc. Serie NA, No. 263, IIMAS-UNAM. Larsen E. W. and Miller W. F. Jr (1980)Nucl. Sci. Engng 73,76. Lathrop K. D. (1969) J. comput. Phys. 4,475. Lee S. M. and Vaidyanathan R. (1980) Nucl. Sci. Engng 76,l.