Available online at www.sciencedirect.com annals of
NUCLEAR ENERGY Annals of Nuclear Energy 34 (2007) 958–966 www.elsevier.com/locate/anucene
The spectral Green’s function linear-nodal method for one-speed X, Y-geometry discrete ordinates deep penetration problems Dany S. Domı´nguez 1, Ricardo C. Barros
*
Departamento de Modelagem Computacional, Instituto Polite´cnico, Universidade do Estado do Rio de Janeiro, Rua Alberto Rangel s/n, 28630-050 Nova Friburgo, RJ, Brazil Received 18 January 2007; received in revised form 22 April 2007; accepted 26 April 2007 Available online 2 July 2007
Abstract A new spectral nodal method is developed for the solution of one-speed discrete ordinates (SN) problems with isotropic scattering in X, Y geometry. In this method, the spectral Green’s function (SGF) scheme, originally developed for solving SN problems in slab-geometry with no spatial truncation error, is generalized to solve the one-dimensional transverse-integrated SN linear-nodal equations with linear polynomial approximation for the transverse leakage terms. The resulting SGF-linear-nodal (SGF-LN) equations are solved with the full-node inversion (FNI) iterative scheme, which uses the best available estimates for the node-entering quantities to evaluate the node angular quantities in all the exiting directions as the equations are swept across the system. We give numerical results that illustrate the accuracy of the SGF-LN method for coarse-mesh calculations. 2007 Elsevier Ltd. All rights reserved.
1. Introduction Neutron radiation protection relies on radiation shielding. A large mass of hydrogen-rich material, such as water or concrete, serves to slow down the neutrons by elastic scattering, so they can then be absorbed by nuclear reactions. A topic of considerable interest for more than two decades is the development of coarse-mesh nodal methods for the efficient numerical solution of the neutron transport equation (Wagner, 1979; Lawrence and Dorning, 1980; Walters and O’Dell, 1981; Badruzzaman et al., 1984; Walters, 1986; Lawrence, 1986; Azmy, 1988a,b; Barros and Larsen, 1992; Mello and Barros, 2002). By efficient numerical solution we mean achieving accurate results in shorter computational running time, as the reduction in the number of discrete variable unknowns is assumed to offset the computational penalty arising from the relatively more *
Corresponding author. Tel.: +55 21 99969271; fax: +55 22 25288536. E-mail address:
[email protected] (R.C. Barros). 1 Current address: Departamento de Cieˆncias Exatas e Tecnolo´gicas, Universidade Estadual de Santa Cruz, Rodovia Ilhe´us/Itabuna km 16, Salobrinho, 45662-000 Ilhe´us, BA, Brazil. 0306-4549/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2007.04.015
complicated nature of the discretized algebraic nodal equations (Azmy, 1988a). In this paper, we develop a new nodal method for onespeed discrete ordinates (SN) deep penetration problems in X, Y-geometry by following the essence of the spectral nodal methods, that solve without spatial truncation errors the one-dimensional transverse-integrated SN nodal equations, only approximating the transverse leakage terms; the scattering source terms are treated exactly. This is in contrast to the conventional weighted diamond-difference form of nodal methods that are based on polynomial approximations for both the transverse leakage and the scattering source terms (Azmy, 1988a). To solve the transverse-integrated SN nodal equations, in our method, we follow the steps taken in the spectral Green’s functionconstant nodal (SGF-CN) method (Barros and Larsen, 1992) and in the SGF-exponential nodal (SGF-ExpN) method (Mello and Barros, 2002) for one-speed X, Ygeometry SN problems. We approximate the transverse leakages through the edges of each spatial node by linear functions, so we call our method the SGF-linear-nodal (SGF-LN) method. Therefore, in contrast to the previous SGF-CN and SGF-ExpN methods, the offered SGF-LN
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966
method is based on two node transverse-integrations of the SN equations for each spatial variable. This procedure yields the transverse-integrated SN nodal equations for the zeroth-order spatial moments of the angular flux, and the transverse-integrated SN nodal equations for the firstorder spatial moments of the angular flux (Walters, 1986). To obtain iterative solutions of the discretized SGF-LN equations, we use the full-node inversion (FNI) scheme. This iterative scheme uses the most recent available estimates for the incoming node-edge spatial moments of the angular fluxes to calculate such quantities in all exiting directions as the equations are swept across the spatial grid set up on the rectangular domain (Dominguez and Barros, 2007). In this paper, Section 2 contains the derivation of the transverse-integrated SN linear-nodal equations for onespeed deep penetration neutron transport problems in X, Y-geometry. In Section 3 we describe the SGF-LN method by deriving the SGF-LN equations for an arbitrary node and describing the FNI iterative scheme to solve the SGF-LN equations in a system consisting of many nodes. We present numerical results for two classical deep penetration SN model problems in Section 4, and conclude with a discussion in Section 5. 2. The transverse-integrated SN linear-nodal equations Let us consider the SN equations in a rectangular domain D of width X and height Y with isotropic scattering lm o gm o wm ðx; yÞ þ w ðx; yÞ þ wm ðx; yÞ rT ðx; yÞ ox rT ðx; yÞ oy m M cðx; yÞ X Qðx; yÞ ¼ ; wn ðx; yÞxn þ 4 rT ðx; yÞ n¼1 where we have M = N(N + 2)/2,
defined
(x, y) 2 D,
m = 1:M
ð1Þ with
N = order of the angular quadrature set, c(x, y) = rs(x, y)/rT(x, y) = scattering ratio, rT(x, y) = total macroscopic cross section, rS(x, y) = isotropic scattering macroscopic cross section, Q(x, y) = isotropic neutron source in D, {lm, gm} = discrete angular directions, m = 1:M, xm = angular quadrature weights, m = 1:M. We assume that these material quantities and the neutron source are piecewise constant functions in D and we consider a rectangular spatial grid on D, where each spatial node is denoted as Di,j of width hi and height kj, as shown in Fig. 1. Therefore, each node Di,j has constant cross sections rTi,j, rSi,j, (ci,j = rsi,j/rsi,j), and a constant interior source Qi,j, i = 1:I and j = 1:J. To obtain the one-dimensional transverse-integrated SN nodal equations, we follow the standard procedure (Walters, 1986). That is, we define the generalized transverse-integration operator
959
yJ+1/2 yJ-1/2 yj+1/2 yj-1/2 y5/2
. . .
kj
Di,j
. . .
y3/2 hi
. . .
y1/2 x1/2
x3/2
x5/2
xi-1/2
. . . xi+1/2
xI-1/2
xI+1/2
Fig. 1. Rectangular spatial grid on D.
Lð‘Þ u ¼
ð2‘ þ 1Þ rs
Z
usþ1=2
P‘ us1=2
2ðu us Þ du; rs
ð2Þ
where u = x or y, s = i or j, r = h or k, ‘ = 0 or 1, and the midpoint usþ1=2 us1=2 ; ð3Þ us ¼ 2 with P‘ being defined as the ‘th degree Legendre polynomial. First we apply Lyð0Þ and Lxð0Þ to Eq. (1) and obtain the onedimensional transverse-integrated SN nodal equations for the zeroth-order spatial moments of the angular flux for the x direction and the y direction, respectively. That is lm d ~ ~ m;j ðxÞ wm;j ðxÞ þ w rTi;j dx M ci;j X ~ n;j ðxÞxn þ Qi;j y1 ½w ðx;y ¼ w jþ1=2 Þ wm ðx;y j1=2 Þ 4 n¼1 rTi;j am;i;j m ð4Þ and gm d ^ ^ m;i ðyÞ wm;i ðyÞ þ w rTi;j dy M ci;j X ^ n;i ðyÞxn þ Qi;j 1 ½w ðxiþ1=2 ;yÞ w ðxi1=2 ;yÞ; ¼ w m 4 n¼1 rTi;j axm;i;j m ð5Þ where we have defined hi rTi;j axm;i;j ¼ ; ð6aÞ lm k j rTi;j aym;i;j ¼ ð6bÞ gm and the averages of the angular fluxes over each spatial dimension within node Di,j (zeroth-order spatial moments) Z y jþ1=2 ~ m;j ðxÞ ¼ 1 w w ðx; yÞ dy; ð7aÞ k j y j1=2 m Z 1 xiþ1=2 ^ wm;i ðyÞ ¼ w ðx; yÞ dx: ð7bÞ hi xi1=2 m
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966
960
Now we apply Lyð1Þ and Lxð1Þ to Eq. (1) and obtain the onedimensional transverse-integrated SN nodal equations for the first-order spatial moments of the angular flux for the x direction and the y direction, respectively lm d ~ m;j ðxÞ ~ m;j ðxÞ þ u u rTi;j dx M ci;j X 3 ~ n;j ðxÞxn y ½wm ðx; y jþ1=2 Þ ¼ u am;i;j 4 n¼1 ~ m;j ðxÞ þ wm ðx; y j1=2 Þ 2w
ð8Þ
and
and ~ m;j ðxi1=2 Þ; ~ m;i1=2;j ¼ U U
ð13bÞ
with U = w or u. At this point, we substitute Eq. (11) into Eqs. (4) and (8), and Eq. (12) into Eqs. (5) and (9). The results for the x coordinate direction are lm d ~ ~ m;j ðxÞ wm;j ðxÞ þ w rTi;j dx M ci;j X ~ n;j ðxÞxn þ a0 2ðx xi Þ þ b0 ; ¼ w ð14Þ m;x m;x hi 4 n¼1 where we have used the definitions
gm d ^ m;i ðyÞ ^ m;i ðyÞ þ u u rTi;j dy M ci;j X 3 ^ n;i ðyÞxn x ½wm ðxiþ1=2 ; yÞ ¼ u am;i;j 4 n¼1 ^ m;i ðyÞ; þ wm ðxi1=2 ; yÞ 2w
ð9Þ
where we have defined the first-order spatial moments of the angular flux Z y jþ1=2 2ðy y j Þ 3 ~ m;j ðxÞ ¼ u wm ðx; yÞ dy; ð10aÞ k j y j1=2 kj Z 1 xiþ1=2 2ðx xi Þ ^ m;i ðyÞ ¼ u wm ðx; yÞ dx: ð10bÞ hi xi1=2 hi At this point we remark that in the conventional linear– linear-nodal (LLN) method and its companion, the linear-nodal (LN) method, both the scattering source terms and the transverse leakage terms in Eqs. (4), (5), (8) and (9) are approximated by linear polynomial functions (Walters, 1986; Azmy, 1988b), and then the resulting system of ordinary differential equations (ODE’s) is solved analytically within each spatial node. This is in contrast to the present SGF-LN method, wherein we use the linear polynomial approximations only for the transverse leakage terms within each spatial node. Therefore, the offered SGF-LN method yields exact numerical solution for problems that are one-dimensional in any of the coordinate directions. On the other hand, this is not true for the conventional LLN method unless the domain is a purely absorbing medium (c = 0). Therefore, we assume that the transverse leakage terms are linear polynomial functions along the edges of each node Di,j. This yields the linear approximations ^ m;i;j1=2 þ 2ðx xi Þ u ^ m;i;j1=2 ; wm ðx; y j1=2 Þ ¼ w hi
ð11Þ
and ~ m;i1=2;j þ 2ðy y j Þ u ~ m;i1=2;j : wm ðxi1=2 ; yÞ ¼ w kj
ð12Þ
ð15aÞ
b0m;x
ð15bÞ
and lm d ~ m;j ðxÞ ~ m;j ðxÞ þ u u rTi;j dx M ci;j X 2ðx xi Þ 6 ~ ~ n;j ðxÞxn þ a1m;x ¼ þ b1m;x þ y w u m;j ðxÞ; hi am;i;j 4 n¼1 ð16Þ with similar definitions ^ m;i;j1=2 Þ 3ð^ um;i;jþ1=2 þ u ; aym;i;j ^ m;i;jþ1=2 þ w ^ m;i;j1=2 Þ 3ðw ¼ : y am;i;j
a1m;x ¼
ð17aÞ
b1m;x
ð17bÞ
We obtain similar results for the y coordinate direction. As we see, Eqs. (14) and (16) are coupled for each angular direction lm, as the solution of Eq. (14) appears on the right-hand side of Eq. (16). In order to simplify the tedious algebraic problem arising from this coupling, we consider a further decoupling approximation. That is, in Eq. (16) we assume that ~ ~ ~ m;j ðxÞ ¼ wm;iþ1=2;j þ wm;i1=2;j þ ðx xi Þ ðw ~ m;iþ1=2;j w ~ m;i1=2;j Þ: w hi 2 ð18Þ A similar decoupling procedure is considered for the y direction. Therefore, Eq. (16) appears as lm d ~ m;j ðxÞ ~ m;j ðxÞ þ u u rTi;j dx M ci;j X 2ðx xi Þ 1 ~ n;j ðxÞxn þ a1m;x ¼ þ bm;x ; u hi 4 n¼1
ð19Þ
where
Here we define ^ m;i ðy j1=2 Þ ^ m;i;j1=2 ¼ U U
^ m;i;j1=2 ^ m;i;jþ1=2 u u aym;i;j ^ m;i;jþ1=2 w ^ m;i;j1=2 Qi;j w ¼ y rTi;j am;i;j
a0m;x ¼
ð13aÞ
a1m;x ¼
3 aym;i;j
~ m;iþ1=2;j w ~ m;i1=2;j u ^ m;i;jþ1=2 u ^ m;i;j1=2 Þ ðw ð20aÞ
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966
3 ~ ~ ^ ^ b1m;x ¼ y ðw m;iþ1=2;j þ wm;i1=2;j wm;i;jþ1=2 wm;i;j1=2 Þ: am;i;j ð20bÞ An analogous result is obtained for the y direction. At this point, we remark that we expect that this decoupling procedure leads to non-significant loss of accuracy of the numerical results. We believe this is a reasonable expectation since in the conventional linear–linear-nodal (LLN) method, this decoupling procedure was also considered in order to derive the companion linear-nodal (LN) method, with non-significant loss of accuracy (Walters, 1986). Therefore, Eqs. (14) and (19) for the x coordinate direction, together with the corresponding equations for the y direction form the transverse-integrated SN linear-nodal equations that we write in the following generalized form Km d ðkÞ F ðuÞ þ F ðkÞ m ðuÞ rTi;j du m M ci;j X 2ðu us Þ ¼ F ðkÞ ðuÞxn þ aðkÞ þ bðkÞ m;u : m;u rs 4 n¼1 n
M X
ð21Þ
C ‘ F hm;‘ ðuÞ þ F p;ðkÞ m ðuÞ:
ð22Þ
‘¼1
Here, C‘ is an arbitrary constant, F hm;‘ ðuÞ is an elementary solution of the homogeneous ODE’s associated with Eq. (21) and F p;ðkÞ m ðuÞ is the particular solution. We remark that the homogeneous solution is independent of k = 0 or 1. However, the particular solution depends upon k, as the coefficients a and b in Eq. (21) do depend upon the value Table 1 Definition of the generalized quantities in Eq. (21) u=x
k=0
k=1
u=y
Km = l m
Km = gm
~ F ð0Þ m ðuÞ ¼ wm;j ðxÞ
^ F ð0Þ m ðuÞ ¼ wm;i ðyÞ
0 að0Þ m;u ¼ am;x ; cf. Eq. (15a)
0 að0Þ m;u ¼ am;y
0 bð0Þ m;u ¼ bm;x ; cf. Eq. (15b)
0 bð0Þ m;u ¼ bm;y
~ m;j ðxÞ F ð1Þ m ðuÞ ¼ u
^ m;i ðyÞ F ð1Þ m ðuÞ ¼ u
að1Þ m;u
cf. Eq. (20a)
að1Þ a1m;y m;u ¼
1 bð1Þ m;u ¼ bm;x ; cf. Eq. (20b)
1 bð1Þ m;u ¼ bm;y
¼
a1m;x ;
of k. Therefore, in order to determine the homogeneous solution, we use the spectral analysis we describe in (Barros and Larsen, 1992). In addition, in order to determine the expression for the particular solution of Eq. (21), we assume that F p;ðkÞ m ðuÞ in Eq. (22) is a linear polynomial function of u, substitute it into Eq. (21) and then we determine the coefficients. Following this procedure, we are able to solve Eq. (21) analytically inside each node Di,j of the spatial grid set up on the domain. In the next section, we derive a numerical method that preserves the local general solution of the transverse-integrated SN nodal equations inside each node Di,j. Using the boundary conditions and the continuity conditions, we generalize the resulting method to solve numerically SN problems on arbitrary rectangular grids. This is the essence of the offered SGF-LN method. 3. The spectral Green’s function linear-nodal method
Here the generalized quantities are defined in Table 1. Our ensuing analysis differs from the conventional LLN method and its simplified companion LN method in that we solve Eq. (21) analytically within each spatial node, whereas in previous work, linear approximations to the scattering source terms in Eq. (21) are introduced, and then, the resulting equations are solved analytically within each spatial node. Moreover, as the transverse-integrated SN linear-nodal equations for the zeroth-order and first-order spatial moments of the angular flux in the x and y variables have the same shape, cf. Eq. (21), we may write the generalized form of the local general solution F ðkÞ m ðuÞ ¼
961
Now we derive the discretized SGF-LN equations for one spatial node and describe the FNI iterative scheme for solving these equations on systems consisting of many nodes. To begin, we consider the operator ð2‘1 þ 1Þð2‘2 þ 1Þ hi k j Z y jþ1=2 Z xiþ1=2 2ðy y j Þ 2ðx xi Þ P ‘1 P ‘2 dx dy; hi kj y j1=2 xi1=2 ð23Þ where P‘ is the ‘th-degree Legendre polynomial and xi and yj are the midpoints given in Eq. (3) for (u, s) = (x, i) or (y, j). In order to obtain the zeroth-order spatial moments discretized SN balance equations, we set ‘1 = ‘2 = 0 in Eq. (23) and apply the resulting operator to the SN Eq. (1). The result is ~ m;iþ1=2;j w ~ m;i1=2;j Þ ðw ^ m;i;jþ1=2 w ^ m;i;j1=2 Þ ðw m;i;j þ þw y x am;i;j am;i;j ¼
M ci;j X n;i;j xn þ Qi;j ; w 4 n¼1 rTi;j
m ¼ 1 : M;
ð24Þ
where we have defined the node-average angular flux Z y jþ1=2 Z xiþ1=2 1 wm;i;j ¼ w ðx; yÞ dx dy: ð25Þ hi k j y j1=2 xi1=2 m Furthermore, we set ‘1 = 1 and ‘2 = 0 in Eq. (23) and apply the resulting operator to the SN Eq. (1). The result is the discretized SN balance equations for the first-order x-moments of the angular flux averaged in the y direction ~ m;iþ1=2;j þ w ~ m;i1=2;j 2w m;i;j Þ ð^ ^ m;i;j1=2 Þ 3ðw um;i;jþ1=2 u þ axm;i;j aym;i;j M X x xn ; m ¼ 1 : M: x ¼ ci;j w þw ð26Þ m;i;j 4 n¼1 n;i;j
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966
962
Here we define Z y jþ1=2 x ¼ 1 ^ m:i ðyÞ dy: w u m;i;j k j y j1=2
ð27Þ
Now we set ‘1 = 0 and ‘2 = 1 in Eq. (23) and apply the resulting operator to the SN Eq. (1). As a result we obtain the discretized SN balance equations for the first-order ymoments of the angular flux averaged in the x direction ^ m;i;jþ1=2 þ w ^ m;i;j1=2 2w m;i;j Þ ~ m;i1=2;j Þ 3ðw ð~ um;iþ1=2;j u þ y x am;i;j am;i;j M X y xn ; m ¼ 1 : M; y ¼ ci;j w þw ð28Þ m;i;j 4 n¼1 n;i;j where y ¼ 1 w m;i;j hi
Z
Therefore, we find that the discretized SN balance equations (24), (26) and (28) have more unknowns than equations, thus leading to an underdetermined system with infinitely many solutions. As it is, additional equations are needed to uniquely specify the solution. In the present SGF-LN method, we follow our earlier work (Barros and Larsen, 1992; Barros et al., 1999; Mello and Barros, 2002) and use the 4M SGF auxiliary equations X X m;i;j ¼ ~ n;i1=2;j þ ~ n;iþ1=2;j þ Hð0Þ ; w hm;n w hm;n w ð30Þ m ln >0
m;i;j ¼ w
X
ln <0
^ n;i;j1=2 þ cm;n w
gn >0
y ¼ w m;i;j ~ m;j ðxÞ dx: u
ð29Þ
xi1=2
x ¼ w m;i;j
X
X
~ n;i1=2;j þ hm;n u ^ n;i;j1=2 þ cm;n u
X
~ n;iþ1=2;j þ Hmð1Þ ; hm;n u
X
^ n;i;jþ1=2 þ Cð1Þ cm;n u m ; m ¼ 1 : M: ð33Þ
Since the transverse-integrated SN linear-nodal equation (21) has the same general form for u = x or y and k = 0 or 1 in Table 1, we use the SGF auxiliary equations (30) and (32) for the x coordinate direction with k = 0 and k = 1, respectively. Moreover, for the y direction, we use the SGF auxiliary equations (31) and (33) for k = 0 and k = 1, respectively. Therefore, in the SGF auxiliary equations (30)–(33), hm,n ð‘Þ and cm,n as well as Hð‘Þ m and Cm , ‘ = 0, 1, are determined by requiring that the general solution of Eq. (21), given by Eq. (22), has node-average and node-edge average quantities that for all constants C‘ satisfy the appropriate SGF auxiliary equation. In other words, evaluating the node-edge average and node-average quantities that correspond to Eq. (22), introducing these into the SGF auxiliary equations, and requiring that the resulting equations hold for all C‘, we obtain, after tedious algebra, expressions for 2 2 ð‘Þ Hð‘Þ m and Cm , ‘ = 0,1, and M linear equations for the M constants hm,n and cm,n, respectively. ψˆ m.i. j+1/ 2 ϕˆ m.i. j+1/ 2
ηm > 0
~ ψ m .i −1 / 2. j ~ ϕ m .i −1 / 2. j
ηm > 0
D i. j
μm < 0
ψ m.i. j m = 1: M ψˆ m.i. j−1/ 2 ηm < 0
~ ψ m .i +1 / 2. j ~ ϕ m .i +1 / 2. j
μm > 0
μm < 0
~ ψ m .i −1 / 2. j
ð32Þ
gn <0
ψˆ m.i. j+1/ 2
D i. j
ð31Þ
ln <0
gn >0
For the nodal methods that use only the zeroth-order spatial moments of the angular flux, such as the constant–constant nodal (CCN) method (Azmy, 1988a), the SGF-CN method (Barros and Larsen, 1992) and the SGF-ExpN method (Mello and Barros, 2002), Eq. (24) together with the boundary conditions and the continuity conditions constitute a system of M algebraic linear equations in 3M unknowns: M outgoing node-edge average angular fluxes in the x coordinate direction, M outgoing node-edge average angular fluxes in the y direction, plus M node-average angular fluxes; viz Fig. 2a. Moreover, a simple count indicates that for the nodal methods that use both the zerothand the first-order spatial moments of the angular flux, such as the LLN method, the LN method (Walters, 1986) and the present SGF-LN method, Eqs. (24), (26) and (28) form a system of 3M algebraic linear equations in 7M unknowns: 2M zeroth- and first-order outgoing node-edge average quantities [Eqs. (7a) and (10a)] in the x coordinate direction, 2M outgoing node-edge average quantities [Eq. (7b) and (10b)] in the y direction, plus 3M node-average quantities given by Eqs. (25), (27) and (29); viz Fig. 2b.
^ n;i;jþ1=2 þ Cð0Þ ; cm;n w m
gn <0
ln >0
xiþ1=2
X
~ ψ m .i +1 / 2. j
ψ m.i. j ψ mx ,i , j ψ my ,i , j
μm > 0
m = 1: M ψˆ m.i. j−1/ 2 ϕˆ m.i. j−1/ 2 ηm < 0
Fig. 2. Count of unknowns in node Di,j. (a) Zeroth-order nodal methods. (b) First-order nodal methods.
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966
Once we have determined Hmð‘Þ and Cð‘Þ m , ‘ = 0, 1, and the constants hm,n and cm,n for each node Di,j, Eqs. (30)–(33) constitute the SGF auxiliary equations, which, with the discretized SN balance equations (24), (26) and (28), are exactly satisfied by the general solution of each one-dimensional transverse-integrated SN linear-nodal equation (21) in the generalized form, bearing in mind Table 1. In fact, the balance equations (24), (26), (28) combined with the auxiliary equations (30)–(33) and appropriate boundary and continuity conditions form the discretized SGF-LN equations. Next, we describe the FNI scheme for iteratively solving the SGF-LN equations. Before presenting the details of this numerical iterative scheme, we remark that the auxiliary equations in any spectral nodal method couple the node-average quantity in any direction to the node-edge average quantities for all the incoming directions. Therefore, spectral nodal methods are not amenable to the conventional source iteration (SI) scheme (Lewis and Miller,
1993), which treats each angular direction independently as the discretized equations are swept across the system. The FNI scheme iterates on the node-edge average angular quantities by performing ‘‘full-node inversions’’. That is, the FNI scheme uses the most recent available estimates for the incoming node-edge average angular quantities, i.e., the zeroth- and the first-order spatial moments of the angular flux (thinner arrows in Fig. 3), to calculate the exiting node-edge average quantities in all the angular directions (thicker arrows in Fig. 3). Each arrow in Fig. 3 represents N(N + 2)/8 directions in each quadrant. As we see, we need to iterate only on the node-edge average angular quantities. Therefore, we substitute the SGF auxiliary equations (30) and (31) into the balance equation (24), the SGF auxiliary equation (32) into the balance equation (28), and the SGF auxiliary equation (33) into the balance equation (26) to remove the node-average angular quantities from the SN balance equations (24), (26) and (28). After convergence has been achieved, we use Eq. (30) or Eq. (31) to determine the node-average angular fluxes, in case we need them. In the following section, we present numerical results illustrating the performance of the SGF-LN method with the FNI iterative scheme.
ψˆ m.i. j+1/ 2
ηm > 0
ϕˆ m.i. j+1/ 2
963
4. Numerical results We consider in this section two numerical experiments. To model the two classical test problems, we used the level-symmetric LQN quadrature set with N = 4 (model problem no. 1) and N = 6 (model problem no. 2) (Lewis and Miller, 1993), and a convergence criterion requiring the discrete maximum norm of the relative deviation between two estimates of the node-edge average scalar fluxes to be less than or equal to 105. The first model problem consists of a uniform isotropic neutron source (Q = 1) surrounded by a shielding material (Q = 0); viz Fig. 4 that represents one-fourth of the whole
ηm < 0 ~ ψ m .i −1 / 2. j μm < 0 ~ ϕm.i −1/ 2. j
~ ψ m .i +1 / 2. j μm > 0 ~ ϕ
μm > 0 μ <0 D i. j m
m .i +1 / 2. j
ηm > 0
ψˆ m.i. j−1/ 2
ηm < 0
ϕˆ m.i. j−1/ 2
Fig. 3. FNI iterative scheme illustration.
10 Q=1
III
σ Τ=1 σ S = 0.5
IV
6
2
5
Q=0
σ T=2 σ S = 0.1 II
I
Reflective boundary Vacuum boundary
1
0
6
0
5
10
Fig. 4. Model Problem No. 1.
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966
964
Table 2 Numerical results for model problem no. 1 Spatial grid
Region I
Region II
1.676d (0.00) 1.676 (0.00) 1.676 (0.00) 1.676
0.4170E1a (0.26)c 0.4160E1 (0.02) 0.4159E1 (0.00) 0.4159E1
0.1986E2 (0.30) 0.1990E2 (0.00) 0.1992E2 (0.00) 0.1992E2
0.2151 (0.33) 0.2145 (0.00) 0.2144 (0.00) 0.2144
0.4148E1 (3.62) 0.4006E1 (0.07) 0.4003E1 (0.03) 0.4004E1
SGF-CN method-(Mello and Barros, 2002) 10 · 10 1.676 (0.00) 20 · 20 1.676 (0.00) 40 · 40 1.676 (0.00)
0.4290E1 (3.14) 0.4165E1 (0.14) 0.4163E1 (0.09)
0.2850E2 (43.07) 0.2019E2 (1.33) 0.2000E2 (0.40)
0.2371 (10.58) 0.2146 (0.00) 0.2145 (0.00)
0.4648E1 (16.11) 0.4013E1 (0.24) 0.4006E1 (0.07)
SGF-ExpN method-(Mello and Barros, 2002) 10 · 10 1.676 (0.00) 20 · 20 1.676 (0.00) 40 · 40 1.676 (0.00)
0.4169E1 (0.24) 0.4161E1 (0.04) 0.4161E1 (0.04)
0.2000E2 (0.40) 0.1995E2 (0.00) 0.1993E2 (0.00)
0.2149 (0.23) 0.2145 (0.00) 0.2144 (0.00)
0.4140E1 (3.42) 0.4007E1 (0.09) 0.4004E1 (0.02)
0.4163E1 (0.10) 0.4160E1 (0.02) 0.4159E1 (0.00) 0.4159E1 (0.00)
0.1991E2 (0.05) 0.1992E2 (0.00) 0.1992E2 (0.00) 0.1992E2 (0.00)
0.2144 (0.00) 0.2144 (0.00) 0.2144 (0.00) 0.2144 (0.00)
0.4021E1 (0.42) 0.4005E1 (0.02) 0.4004E1 (0.00) 0.4004E1 (0.00)
LN method-(Azmy, 1988b) 10 · 10b 20 · 20 40 · 40 100 · 100
SGF-LN method 10 · 10 20 · 20 40 · 40 100 · 100 a b c d
1.676 (0.00) 1.676 (0.00) 1.676 (0.00) 1.676 (0.00)
Region IV
Region 1
Region 2
Read as 0.4170 · 101. Number of nodes in the x direction · number of nodes in the y direction. Relative deviation (%) with respect to the fine-grid (100 · 100) LN results. Bold face for the relative deviations less than 1%.
shielding structure. This shielding problem was considered in (Azmy, 1988b; Mello and Barros, 2002). Our numerical experiment consists in calculating the average scalar flux in regions I, II and IV, as well as in the more localized regions 1 and 2 in Fig. 4. Table 2 shows the numerical results generated by the conventional LN method (Azmy, 1988b) and three distinct spectral nodal methods, namely the SGF-CN method (Barros and Larsen, 1992), the SGF-ExpN method (Mello and Barros, 2002) and the present SGF-LN method on various spatial grids. The numerical results for the average scalar fluxes, listed in bold face in Table 2, are referred to as valid results in this numerical experiment, as the relative deviations with respect to the fine-grid (100 · 100) LN results are less than 1%. As we see, the present SGF-LN method generated very accurate results, and all the results are valid (< 1%) for the spatial grid composed of 10 nodes in each spatial direction. This is in contrast to the LN and the SGF-ExpN methods that generated valid numerical results on finer spatial grid (20 · 20). Moreover, in this numerical experiment, the
SGF-CN method generated valid results on even finer spatial grid (40 · 40). We remark that for the 10 · 10 SGF-LN run, the CPU execution time was roughly 42% shorter than the CPU execution time of the 20 · 20 LN run, that generated valid results in this numerical experiment, as we see in Table 4. The second model problem is a monoenergetic version of the oil well-logging problem (Azmy, 1988b; Mello and Barros, 2002). The geometry and nuclear data for this test problem are shown in Fig. 5. Table 3 displays the numerical results generated for the average scalar fluxes in regions D1 and D2, that represent the locations of the detectors traditionally used in this configuration. These results were generated by the same nodal methods used to solve the previous test problem on various spatial grids. In this second numerical experiment, the bold faced results in Table 3 are referred to as valid, since the relative deviations with respect to the fine-grid (112 · 128) LN results are less than 1%. We note that the SGF-LN method generated very accurate numerical results
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966 VACUUM 24 32
8
Table 3 Numerical results for model problem no. 2
56
64
Spatial grid
Average scalar flux in detector D1
STEEL 56 D2 48
L I M E S T O N E
V A C U U M
WATER
STEEL
LIMESTONE
965
V A C U U M
24 D1 16 STEEL 8 Q=1
REFLECTIVE BOUNDARY
Fig. 5. Model problem No. 2. Limestone: rT(rs) = 0.330263 (0.314419); water: rT(rs) = 0.694676 (0.634883) and steel: rT(ra) = 0.499122 (0.494460), width = 56 cm and height = 64 cm.
LN method-(Azmy, 1988b) 1.647 7 · 8a (4.07)b 14 · 16 1.732c (0.87) 28 · 32 1.718 (0.05) 56 · 64 1.717 (0.00) 112 · 128 1.717
0.1374 (10.18) 0.1260 (1.03) 0.1247 (0.00) 0.1247 (0.00) 0.1247
SGF-CN method-(Mello and Barros, 2002) 7·8 1.527 (11.06) 14 · 16 1.697 (1.17) 28 · 32 1.705 (0.70)
0.1376 (10.34) 0.1260 (1.03) 0.1250 (0.24)
SGF-ExpN method-(Mello and Barros, 2002) 7·8 1.638 (4.60) 14 · 16 1.699 (1.05) 28 · 32 1.709 (0.46) SGF-LN method 7·8
1.692 (1.46) 1.719 (0.12) 1.717 (0.00) 1.717 (0.00) 1.717 (0.00)
14 · 16 28 · 32 56 · 64
for this problem. As we see, D1 and D2 results are valid (< 1%) for the spatial grid composed of 14 nodes in the x direction and 16 nodes in the y direction. This is in contrast to the LN, SGF-CN and SGF-ExpN runs, that generated valid results for the average scalar fluxes in D1 and D2 in finer spatial grid (28 · 32). In addition, we remark that the CPU execution time of the 14 · 16 SGF-LN run was 80% shorter than the CPU time of the 28 · 32 LN run to generate valid results in the sense we have defined in this numerical experiment, viz Table 4. 5. Discussion We described a spectral nodal method, that we refer to as the SGF-LN method for monoenergetic SN deep penetration problems in X, Y-geometry. In the class of spectral nodal methods, the only approximation involved is in the transverse leakage terms of the transverse-integrated SN nodal equations. The source terms, that include the scattering events and the external source, are treated analytically. Therefore, in slab-geometry SN problems, the spectral nodal methods generate numerical solutions that are completely free from spatial truncation errors (Barros and Larsen, 1990). Moreover, based on the results generated by the present SGF-LN method for the two model problems con-
Average scalar flux in detector D2 (·10)
112 · 128
0.1359 (8.98) 0.1249 (0.16) 0.1247 (0.00) 0.1289 (3.37) 0.1248 (0.08) 0.1247 (0.00) 0.1247 (0.00) 0.1247 (0.00)
a
Number of nodes in the x direction · number of nodes in the y direction. b Relative deviation (%) with respect to the fine-grid (112 · 128) LN results. c Bold face for the relative deviations less than 1%.
sidered in this paper, we conclude that the offered SGF-LN method is more accurate than the conventional LN method for coarse-mesh calculations. In addition the SGF-LN coarse-grid results turned out to be more accurate than
Table 4 Iteration counts and CPU time for the valid LN versus SGF-LN runs Model problem
Numerical method
Spatial grid
Iteration counts
Model problem no. 1
LN SGF-LN
20 · 20 10 · 10
21 SIa 9 FNIb
Model problem no. 2
LN SGF-LN
28 · 32 14 · 16
337 SI 53 FNI
a b
Source iterations. FNI iterations.
CPU time (s) 0.350 0.201 13 2.65
966
D.S. Domı´nguez, R.C. Barros / Annals of Nuclear Energy 34 (2007) 958–966
the numerical results generated by the companion SGF-CN and the SGF-ExpN methods. A negative feature of the SGF-LN method is that it requires more storage than the LN method. This is due to the fact that the FNI iterative scheme requires the storage of the node-edge average angular quantities. However, because the SGF-LN method is more accurate than the LN method, this extra storage requirement can be compensated by using coarser spatial grids. This certainly helps to alleviate the amount of storage and the quantity of floating point operations necessary in the iterative scheme. As a result, the execution time decreases to the point that the SGFLN method generates acceptably accurate solutions much faster than the standard unaccelerated LN method. As future work, we suggest the investigation of an efficient technique to accelerate the FNI iterative scheme used in the SGF-LN method. We also suggest the extension of the present method to multigroup SN problems, even considering three-dimensional Cartesian geometry. One serious drawback of the nodal methods, as a whole, is the reduced resolution, i.e., accurate point values of the neutron flux, or other related quantities, e.g., absorption rate, are somewhat difficult to obtain from the converged coarse-mesh solution (Azmy, 1988a). In this context, we intend to develop approximate nodal reconstruction schemes in order to recover the spatial profile of the neutron flux averaged over each spatial dimension within each spatial node. In conclusion, we remark that the tedious algebra involved in the derivation of the SGF-LN equations and in setting them in the appropriate form for the computational implementation of the FNI iterative scheme has been overwhelmed by the efficiency of the SGF-LN code that generates accurate coarse-mesh results in shorter CPU time.
Acknowledgements The work by R.C. Barros was supported by Conselho Nacional de Desenvolvimento Cientı´fico e Tecnolo´gico (CNPq, Brazil), and the work by D.S. Dominguez was supported by Coordenac¸a˜o de Aperfeic¸oamento de Pessoal de Nı´vel Superior (CAPES, Brazil). References Azmy, Y.Y., 1988a. Nuclear Science and Engineering 98, 29. Azmy, Y.Y., 1988b. Nuclear Science and Engineering 100, 190. Badruzzaman, A., Xie, Z., Dorning, J.J., Ullo, J.J., 1984. In: Proc. Topl. Mtg. Reactor Physics and Shielding, Chicago, Illinois, 1984. American Nuclear Society, p. 170, September 17–19. Barros, R.C., Larsen, E.W., 1990. Nuclear Science and Engineering 104, 199. Barros, R.C., Larsen, E.W., 1992. Nuclear Science and Engineering 111, 34. Barros, R.C., da Silva, F.C., Alves Filho, H., 1999. Progress in Nuclear Energy 35 (3–4), 293. Dominguez, D.S., Barros, R.C., 2007. Proc. Mathematics and Computations and Supercomputing in Nuclear Applications, Monterey, California, American Nuclear Society, April 15–19, Session 8A, CD/ISBN: 0-89448-059-6. Lawrence, R.D., 1986. Progress in Nuclear Energy 17, 271. Lawrence, R.D., Dorning, J.J., 1980. In: Proc. Topl. Mtg. Advances in Reactor Physics and Shielding, Sun Valley, Idaho, vol. 2. American Nuclear Society, p. 840, September 14–17. Lewis, E.E., Miller, W.F., 1993. Computational Methods of Neutron Transport, American Nuclear Society. La Grange Park, Illinois, USA. Mello, J.A.M., Barros, R.C., 2002. Annals of Nuclear Energy 29, 1855. Wagner, M.R., 1979. In: Proc. Topl. Mtg. Computational Methods in Nuclear Engineering, Williamsburgh, VA, vol. 2. American Nuclear Society, p. 4117. April 23–25. Walters, W.F., 1986. Progress in Nuclear Energy 18, 21. Walters, W.F., O’Dell, R.D., 1981. In: Proc. Int. Topl. Mtg. Advances in Mathematical Methods for the Solution of Nuclear Engineering Problems, Munich, FRG, Kernforschungszentrum Karlsruhe, April 27–29, vol. 1, p. 115.