- Email: [email protected]
Solution of one-speed neutron transport equation for strongly anisotropic scattering by TN approximation: Slab criticality problem
annals of
NUCLEAR ENERGY Annals of Nuclear Energy 34 (2007) 743–751 www.elsevier.com/locate/anucene
Solution of one-speed neutron transport equation for strongly anisotropic scattering by TN approximation: Slab criticality problem Ayhan Yilmazer
*
Department of Nuclear Engineering, Hacettepe University, Beytepe, Ankara 06800, Turkey Received 27 March 2006; received in revised form 17 March 2007; accepted 18 March 2007 Available online 9 May 2007
Abstract In this study, a recently proposed version of Chebyshev polynomial approximation which was used in spectrum and criticality calculations by one-speed neutron transport equation for slabs with isotropic scattering is further developed to slab criticality problems for strongly anisotropic scattering. Backward–forward-isotropic model is employed for the scattering kernel which is a combination of linearly anisotropic and strongly backward–forward kernels. Further to that, the common approaches of using the same functional form for scattering and fission kernels or embedding fission kernel into the scattering kernel even in strongly anisotropic scattering is questioned for TN approximation via taking an isotropic fission kernel in the transport equation. As a starting point, eigenvalue spectrum of one-speed neutron transport equation for a multiplying slab with different degrees of anisotropy in scattering and for different cross-section parameters is obtained using Chebyshev method. Later on, the spectra obtained for different degree of anisotropies and cross-section parameters are made use of in criticality problem of bare homogeneous slab with strongly anisotropic scattering. Calculated critical thicknesses by Chebysev method are almost in complete agreement with literature data except for some limiting cases. More importantly, it is observed that using a different kernel (isotropic) for fission rather than assuming it equal to the scattering kernel which is a more realistic physical approach yields in deviations in critical sizes in comparison with the values presented in literature. This separate kernel approach also eliminates the slow convergency and/or non-convergent behavior of high-order approximations arising from unphysical eigenspectrum calculations. Ó 2007 Elsevier Ltd. All rights reserved.
1. Introduction A large collection of works exists on the one-speed neutron transport equation with isotropic, linearly anisotropic, and strongly anisotropic scattering for spheres and slabs either bare or reflected. Among the most widely used methods in those works is spherical harmonics or PN method whose earlier applications could be found in the monograph of Davison (1958). The critical slab problem in one-speed neutron transport theory was studied by Yildiz (1998) using an admixture of linearly anisotropic and strongly peaked backward and/or forward scattering kernel, and it was concluded that spherical harmonics method *
Tel.: +90 312 2977300; fax: +90 312 2992122. E-mail address: [email protected]
0306-4549/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2007.03.010
is very suitable for the numerical solution of neutron transport equation with linearly anisotropic scattering. As seen from the data of critical thicknesses given in this study, the effect of existence of strongly peaked backward and forward scattering components together in the scattering kernel resulted in missing data when obtaining critical thicknesses for large backward and forward scattering probabilities and collision parameter. Similar behavior at higher order PN approximation was reported in the study of Sahni et al. (2004) in which the criticality problem of the one-speed transport equation with isotropic scattering for reflected spheres was studied using PN approximation. In this study, it was reported that the method diverges as one increases the order PN approximation, but, the method produces results that agree with those of the integral equation. Williams (2002) developed an integral form of the
744
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
transport equation for bare spheres incorporating fission plus and arbitrary proportion of isotropic, backward and forward scattering. In this study, resulting integral equation was solved numerically and obtained critical radii were compared with those obtained by Yildiz and Alcan (1995) who accounted for strong anisotropy using a synthetic kernel. William’s work should be particularly noted since it uses a separate and isotropic kernel for fission instead of taking its functional form equal to the scattering kernel or embedding it into the scattering kernel in the transport equation. We will discuss later that this physically more convenient choice could also be used in Chebyshev polynomial-TN approximation or in any proper polynomial approximation as it was used in the integral transport equation formulation of William’s study. Among other various methods employed in solving onespeed neutron transport equation is singular eigenfunction or Case’s normal-mode expansion technique which was successfully applied reflected slab and sphere criticality problems with linearly anisotropic scattering by Atalay (1997). Atalay, in his study, has shown that normal-mode expansion technique serves to evaluate the complete spectrum of a critical system. Recently, a new version of Chebyshev polynomial expansion for the angular flux has been proposed to obtain eigenspectrum of one-speed neutron transport operator with isotropic scattering in a slab (Yasßa et al., 2006). Then, the so-called TN method was applied to bare and reflected critical slab problems with isotropic scattering (Anli et al., 2006a,b). One of our main objectives in this study is to extend Chebyshev method in slab criticality problem with anisotropic scattering which is an admixture of linearly anisotropic and strong forward–backward parts (backward–forward-isotropic model). When doing this, to compare the critical thicknesses computed by TN method with literature values we will initially use the common approach of using same kernel for fission with anisotropic kernel. Later on, to observe the effect of fission kernel on spectrum or equivalently slab critical size, we will use an isotropic kernel for fission and the backward–forward-isotropic model for the scattering kernel. It was reported in studies (Anli et al., 2006a,b) that TN approximation in criticality of bare and reflected slab problems with isotropic scattering applies equally well as the PN approximation. One presumably can predict that TN approximation could equally be applicable in criticality of bare and reflected slab problems with anisotropic scattering as well as the PN approximation since both Legendre and Chebyshev polynomials belong to a general class of polynomials known as Jacobi polynomials. However, it is still worthwhile to adopt Chebyshev polynomials expansion for the angular part of the neutron angular flux in the solution of the transport equation because of its comparatively simple algebraic manipulation and it is appealing to investigate whether it could alleviate the slowly convergent or divergent cases encountered in PN approximation. Of course, a more rigorous study would be to employ the most general form of the Jacobi polynomials expansion for the
angular flux and to discuss qualitatively on purely mathematical basis for which Jacobi polynomial best convergent solution is obtainable for a particular problem. 2. Basic equations To obtain an accurate solution to transport equation when anisotropy is taken into account is more difficult and involved and effect of anisotropy was studied by many authors. The one-speed neutron transport equation could written as Xrwðr;XÞþRT wðr;XÞ Z Z ¼ RS fs ðXX0 Þwðr;X0 ÞdX0 þmRF ff ðXX0 Þwðr;X0 ÞdX0 4p
4p
ð1Þ where w(r, X) is the angular flux, RT is the total cross-section, RS is the scattering cross section, RF is the fission cross-section, m is the mean number of neutrons per fission, fs and ff stand for distribution functions (kernel) for scattered and fission neutrons, respectively. Since fission kernel is an isotropic event in the laboratory reference frame, it is reasonable to assume it in the form ff(X Æ X 0 ) = 1/4p. It is commonly assumed that scattering function fs(X 0 Æ X) which is representative of anisotropy in the interactions of neutrons with matter could be developed into a Legendre polynomial series (Davison, 1958) of the argument X 0 Æ X as follows: fs ðX0 XÞ ¼
N 1 X ð2n þ 1Þbn P n ðX0 XÞ 4p n¼0
ð2Þ
where X 0 and X stand for the direction of the neutron velocity before and after scattering, respectively, with the normalization condition, and b1 is the mean cosine of scattering angle. For the linearly anisotropic scattering, it is enough to take only the first two terms of series expansion. However, we will use backward–forward-isotropic model given as fs ðX0 XÞ ¼
1ab a ð1 þ 3b1 X0 XÞ þ dðX0 X 1Þ 4p 2p b ð3Þ þ dðX0 X þ 1Þ 2p
for scattering kernel where a and b are forward and backward scattering probabilities in a collision, respectively. Such a kind of anisotropic scattering has been used by many investigators. The first term on the right hand side of Eq. (3) represents linearly anisotropic scattering, second and third terms represent forward and backward scattering, respectively. The backward–forward-isotropic model that we have used in this study was proposed by Fermi for inclusion of anisotropy into the transport equation and has been used by many other investigators; for example, Yildiz (1998) used this kernel to calculate critical slab thickness for various degrees of anisotropy and cross section parameters in PN. Using scattering kernel given by
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
Eq. (3), we will first calculate eigenspectrum of one-speed neutron transport equation for a homogeneous slab with anisotropic scattering and use resulting spectrum in the slab criticality problems. For one-dimensional case Eq. (3) takes the following form: 1 fs ðll0 Þ ¼ ð1 a bÞð1 þ 3b1 ll0 Þ þ adðl l0 Þ 2 þ bdðl þ l0 Þ
ð4Þ
and the one-speed neutron transport equation given by Eq. (1) becomes l
o ½wðx;lÞ þ RT wðx;lÞ ox Z 1 Z 1 fs ðll0 Þwðx;l0 Þdl0 þmRF ff ðll0 Þwðr;l0 Þdl0 ¼ RS 1
ð5Þ
1
where ff(ll 0 ) is the one-dimensional fission kernel for a slab and is of the form ff(ll 0 ) = 1/2. If we scale the distance in units of mean free path and define collision and fission ratios as cs ¼
Rs Rt
and
cf ¼
mRf Rt
ð6Þ
the transport equation given by Eq. (5) takes the following form: Z 1 o l ½wðx; lÞ þ wðx; lÞ ¼ cs fs ðll0 Þwðx; l0 Þ dl0 ox 1 Z 1 þ cf ff ðll0 Þwðr; l0 Þ dl0 ð7Þ 1
If we use the conventional approaches of either using the same kernel for scattering and fission or embedding isotropic fission kernel into the scattering kernel, or the scattering event is isotropic Eq. (7) will become Z 1 o l ½wðx; lÞ þ wðx; lÞ ¼ c fs ðll0 Þwðx; l0 Þ dl0 ð8Þ ox 1 where total number secondary neutrons generated per collision c is related to cs and cf as c ¼ cs þ cf
ð9Þ
3. Chebyshev polynomial (TN) approximation As we have already touched on that Legendre and Chebyshev polynomials belong to a more general class of polynomials known as Jacobi or hypergeometric polynomials which are solutions to the Jacobi differential equation. Jacobi polynomials are represented by P ðg;jÞ ðlÞ and form n a complete set of orthogonal system in the interval [1, 1] with respect to the weighting function (1 l)g(1 + l)j. For g = j = 0, P nð0;0Þ ðlÞ reduces to Legendre polynomial, and for g = j = 1/2, P nð1=2;1=2Þ ðlÞ reduces to the Chebyshev polynomials of the first kind Tn(l). The orthogonality relation of Jacobi polynomials are of the following form:
Z
745
1 g
1
j
ð1 xÞ ð1 þ xÞ P mðg;jÞ ðlÞP nðg;gÞ ðlÞ dl
¼
2gþjþ1 Cðm þ g þ 1ÞCðm þ j þ 1Þ dmn m!Cðm þ g þ j þ 1Þ 2m þ g þ j þ 1
ð10Þ
Following previous information, it could be inferred that angular flux could be expanded in any Jacobi polynomial series since any function defined in the interval [1, 1] satisfying certain requirements could be expressed in series expansion of an orthogonal set of functions which is complete. Apart from some mathematical concerns such as convergency and possible singularities over the domain, it seems choice of Jacobi polynomials are completely arbitrary for the representation of the angular flux on purely mathematical basis. But, in addition to these mathematical annoyances, we will indicate some physical inconsistencies which could arise, etc. assigning Tn(l) for series representation of angular flux results in a difference form of diffusion coefficient or Fick’s Law obtained using Pn(l) expansion. Notwithstanding the mathematical and physical considerations, we could expand angular flux in terms of Chebyshev polynomials as proposed by Anli et al. (2006a) as 1 X 1 wðx; lÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi ð2 dn0 ÞT n ðlÞ/n ðxÞ ð11Þ p 1 l2 n¼0 where Tn(l) is nth order Chebyshev polynomial of the first kind. The mth Chebyshev moment of the angular flux can be defined as follows: Z 1 1 T m ðlÞ X pffiffiffiffiffiffiffiffiffiffiffiffiffi /m ðxÞ ¼ ð2 dn0 ÞT n ðlÞ/n ðxÞdl ð12Þ 1 l2 n¼0 1 p If we introduce the expansions given by Eq. (11) into Eq. (5) and use isotropic kernel for fission and the kernel given by Eq. (4) for scattering kernel transport equation takes the following form: 1 X ð2 dn0 ÞT n ðlÞ o pffiffiffiffiffiffiffiffiffiffiffiffiffi l /n ðxÞ þ ð1 cS ða þ ð1Þn bÞ/n ðxÞ ox p 1 l2 n¼0 cs cf ¼ ð1 a bÞ½/0 ðxÞ þ 3b1 l/1 ðxÞ þ /0 ðxÞ ð13Þ 2 2 Eq. (13) may be arranged if we use the recurrence relation for Chebyshev polynomials given below 1 lT n ðlÞ ¼ ðð1 þ dn;0 ÞT nþ1 ðlÞ þ T n1 ðlÞÞ 2 as " 1 X ð2 dn0 Þ½ð1 þ dn0 ÞT nþ1 ðlÞ þ T n1 ðlÞ o pffiffiffiffiffiffiffiffiffiffiffiffiffi ½/ ðxÞ ox n 2p 1 l2 n¼0 ð2 dn0 ÞT n ðlÞ n pffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 cs ða þ ð1Þ bÞ/n ðxÞ þ p 1 l2 cs cf ¼ ð1 a bÞ½/0 ðxÞ þ 3b1 l/1 ðxÞ þ /0 ðxÞ 2 2
ð14Þ
ð15Þ
After multiplying Eq. (15) by Tm(l) and integrating over l 2 (1, 1) and using the orthogonality relation of Chebyshev polynomials given below
746
Z
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
1
1
ð2 dn0 ÞT m ðlÞT n ðlÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi dl ¼ dnm p 1 l2
ð16Þ
we obtain 1d m ½ð1þdm0 Þ/mþ1 ðxÞþ/m1 ðxÞþð1cs ðaþð1Þ bÞ/m ðxÞ 2 dx ( " #) ð1þð1Þm Þ ð1þð1Þmþ1 Þ ¼cs ð1abÞ / ðxÞ3b1 /1 ðxÞ 2ðm2 1Þ 0 2ðm2 4Þ cf
ð1þð1Þm Þ / ðxÞ 2ðm2 1Þ 0
ð17Þ
The quantities /m(x) (m = 0, 1, 2 . . .) could be referred as the Chebyshev moments of the angular flux similar to spherical harmonics moments in the PN approximation. It is obvious that the first two moments /0(x) and /1(x) are identical with the scalar flux and the net current. A discrepancy merely appears in the definition or restatement of the Fick’s Law as will be shown later. This could be interpreted as if a new definition for diffusion coefficient were introduced in TN approximation. Although the exact solution of the differential equations system in Eq. (17) is impossible, an approximate solution can be found by assuming that /N+1(x), say, is negligibly small. The resulting equations are called the TN approximation: 1 d ½ð1 þ dm0 Þ/mþ1 ðxÞ þ /m1 ðxÞ 2 dx þ ð1 cs ða þ ð1Þm bÞ/m ðxÞ ¼ cs ð1 a bÞ " # ð1 þ ð1Þm Þ ð1 þ ð1Þmþ1 Þ / ðxÞ 3b1 /1 ðxÞ 2ðm2 1Þ 0 2ðm2 4Þ cf
ð1 þ ð1Þm Þ / ðxÞ; 2ðm2 1Þ 0
m ¼ 0; 1; 2; . . . ; N 1
m¼N
ð18Þ
ð19Þ
If we use the common approach of embedding fission kernel into the scattering kernel, or use same kernels for both events as in the isotropic scattering, c = cs + cf would be introduced as a lumped parameter for secondary neutrons, and the coupled set of equations given by Eqs. (18) and (19) take the form 1 d ð1 þ dm0 Þ/mþ1 ðxÞ þ /m1 ðxÞ 2 dx þ ð1 cða þ ð1Þm bÞ/m ðxÞ ¼ cð1 a bÞ " # ð1 þ ð1Þm Þ ð1 þ ð1Þmþ1 Þ / ðxÞ 3b1 /1 ðxÞ ; 2ðm2 1Þ 0 2ðm2 4Þ m ¼ 0; 1; 2; . . . ; N 1
It should be reminded that the approach just indicated before obtaining Eqs. (20) and (21) has also been used in PN approximation somehow in spectrum and criticality calculations of slab and sphere for anisotropic scattering problems. T1, namely taking N = 1, approximation results in two coupled differential equations for scalar flux /0(x) and net current /1(x) as following: d ½/ ðxÞ þ ð1 cÞ/0 ðxÞ ¼ 0 ð22Þ dx 1 1 d ½/ ðxÞ þ ð1 cfð1 a bÞb1 þ ða bÞgÞ/1 ðxÞ ¼ 0 2 dx 0 ð23Þ Eq. (23) could be rearranged for isotropic scattering as /1 ðxÞ ¼
1 d ½/ ðxÞ 2ð1 cb1 gÞ dx 0
ð24Þ
introducing another version of Fick’s Law stating that diffusion coefficient D is equal to 1 2ð1 cb1 Þ
ð25Þ
In order to obtain eigenvalue spectrum, we employ the well known procedure of seeking a solution of the form
N
ð1 þ ð1Þ Þ /0 ðxÞ; 2ðN 2 1Þ
ð21Þ
D¼
1 d N ½/ ðxÞ þ ð1 cs ða þ ð1Þ bÞ/N ðxÞ ¼ cs ð1 a bÞ 2 dx " N 1 # ð1 þ ð1ÞN Þ ð1 þ ð1ÞN þ1 Þ /0 ðxÞ 3b1 /1 ðxÞ 2ðN 2 1Þ 2ðN 2 4Þ cf
1 d ½/ ðxÞ þ ð1 cða þ ð1ÞN bÞ/N ðxÞ ¼ cð1 a bÞ 2 dx" N 1 # N N þ1 ð1 þ ð1Þ Þ ð1 þ ð1Þ Þ /0 ðxÞ 3b1 /1 ðxÞ ; m ¼ N 2ðN 2 1Þ 2ðN 2 4Þ
/m ðxÞ ¼ am ðmÞex=m
ð26Þ
Substituting Eq. (26) into Eq. (18), we get following systems of algebraic equations: a1 ¼ mð1 ðcs þ cf ÞÞa0
ð27Þ
a2 ¼ a0 2m½ð1 ða bÞ ðcs þ cf Þð1 a bÞb1 a1 ð28Þ m
amþ1 ¼ am1 2mð1 cs ða þ ð1Þ bÞam cs fð1 a bÞ " # m mþ1 ð1 þ ð1Þ Þ ð1 þ ð1Þ Þ a0 3b1 a1 2ðm2 1Þ 2ðm2 4Þ m
cf
ð1 þ ð1Þ Þ a0 ; 2ðm2 1Þ
36m6N
ð29Þ
to compute am with a0 = 1. If the equations given above are compatible, the determinant of coefficients must vanish. This is equivalent to use well-known closure relationship aN þ1 ¼ 0
ð20Þ
ðm ¼ 0; 1; . . . ; N Þ
ð30Þ
which enables us to compute eigenvalues. To obtain eigenvalue spectrum for isotropic scattering or for embedding approach, applying similar steps used above or just putting c = cs + cf in Eqs. (27)–(29), below system of algebraic equations could be obtained
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
a1 ¼ mð1 ðcs þ cf ÞÞa0 a2 ¼ a0 2m½ð1 cða bÞ cð1 a bÞb1 a1
ð31Þ ð32Þ
T N þ1 ðlÞ ¼ 2
ð33Þ
with the closure relationship aN þ1 ¼ 0
ð34Þ
which enables us to compute eigenvalues.
In this work, we used Mark and Marshak boundary conditions for the slab criticality. Mark boundary condition is merely used only for one illustrative case to demonstrate that Marshak boundary condition gives more accurate results than those of Mark for small N (Davison, 1958). Rest of the computations is based on Marshak boundary conditions. It could be inspected through Eqs. (27)–(29) that a‘ (m) = (1)‘a‘(m). Hence, for TN approximation with N odd, there are N + 1 roots ±m‘, corresponding to a solution for the mth Chebyshev moment of angular flux as /m ðxÞ ¼
m x=m‘
b‘ am ðm‘ Þ½ex=m‘ þ ð1Þ e
ð35Þ
‘¼0
When we use Eq. (35) in the TN approximation given by Eq. (11), angular flux becomes N X 1 wðx; lÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi ð2 dm0 ÞT m ðlÞ p 1 l2 m¼0
ðNX þ1Þ=2
ð36Þ
negative roots lk which correspond to the specified angles where incoming angular flux vanish on the boundary become
ð2k 1Þp lk ¼ cos ; ðN þ 1Þ=2 < k 6 N þ 1 ð40Þ 2ðN þ 1Þ or letting k ! k + (N + 1)/2 transformation
ð2k þ N Þp lk ¼ cos ; 1 6 k 6 ðN þ 1Þ=2 2ðN þ 1Þ
ð41Þ
Mkm ðaÞbm ¼ 0; k ¼ 1; . . . ; ðN þ 1Þ=2; m ¼ 1; . . . ; ðN þ 1Þ=2 ð42Þ where bm is the vector comprising elements bm and Mkm ðaÞ is the coefficient matrix with [(N + 1)/2]2 elements and whose components are given as Mkm ¼ 2
N X ‘¼0
1 ‘ a‘ ðmm Þðea=mm þ ð1Þ ea=mm ÞT ‘ ðlk Þ; d‘0 þ 1
k ¼ 1; . . . ; ðN þ 1Þ=2;
m ¼ 1; . . . ; ðN þ 1Þ=2
ð43Þ
Eq. (42) which is a linear system of equations in bm has a non-trivial solution only if det½Mkm ðaÞ ¼ 0 which is the criticality condition. If we use Marshak boundary condition which is based on zero incoming current on the vacuum boundary and could be approximated as Z 1 T 2k1 ðlÞwða; lÞ dl ¼ 0; k ¼ 1; 2; . . . ; ðN þ 1Þ=2
in the TN approximation using discrete eigenvalues given by Eq. (36), we obtain N Z 1 X ð2 dm0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi T 2k1 ðlÞT m ðlÞ dl 1 l2 m¼0 0 p
l>0
ðNX þ1Þ=2
ð45Þ
Using peculiarity property Tn(l) = (1)nTn(l) of Chebyshev polynomials, Eq. (45) could be rearranged as Z 1 N X ð2 dm0 Þ 2kþm1 pffiffiffiffiffiffiffiffiffiffiffiffiffi T 2k1 ðlÞT m ðlÞ dl ð1Þ 1 l2 0 p m¼0
ð38Þ
where lk are negative roots of TN+1(l) = 0. Since first type Chebyshev polynomials of order N + 1 could be represented as
m
b‘ am ðm‘ Þ½ea=m‘ þ ð1Þ ea=m‘ ¼ 0
‘¼0
ð37Þ
where a is the slab half thickness in units of mean free path. We could use the Mark boundary condition which implies zero incoming angular flux on the boundaries for specified values of l wða; lk Þ ¼ 0
ð39Þ
ð44Þ m
b‘ am ðm‘ Þ½ex=m‘ þ ð1Þ ex=m‘
where the coefficients b‘ should be determined from boundary conditions. Solution of transport equation with vacuum boundary condition must be symmetrical with respect to position as
wðx; lÞ ¼ wðx; lÞ
ð2k 1Þp l cos 2ðN þ 1Þ
0
‘¼0
wða; lÞ ¼ 0;
If we use boundary condition given by Eq. (38) in Eq. (36) we obtain following matrix equation for coefficients bm:
4. Criticality condition
ðNX þ1Þ=2
N þ1 Y k¼1
m
amþ1 ¼ am1 mf2ð1 cða þ ð1Þ bÞam ð1 þ ð1Þm Þ þc ð1 a bÞ a0 ðm2 1Þ ##) mþ1 ð1 þ ð1Þ Þ þ3b1 a1 ; 36m6N 2 ðm 4Þ
N
747
ðNX þ1Þ=2 ‘¼0
Since
m
b‘ am ðm‘ Þ½ea=m‘ þ ð1Þ ea=m‘ ¼ 0
ð46Þ
748
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
Z
1
comparison with separate kernels approach. Whenever such a discrepancy is encountered in the calculations it will be demonstrated that this anomaly can be alleviated using separate kernels approach. Table 1 presents calculated critical thicknesses using TN approximation in comparison with the ones obtained using PN in Yildiz’s work (1998) for different degrees of forward scattering and linearly anisotropic scattering. For the sake of comparability on the same basis, embedding approach or using the same kernel for the fission as the scattering kernel is employed initially in the TN method as in the work of Yildiz’s PN approximation. Results shown in Table 1 are for a collision parameter c = 1.2 and degree of backward scattering b = 0. Some important features of TN approximation can be drawn from Table 1. First, even though low order TN approximations are worse than the same order PN approximation in terms of convergence, higher order approximations are in almost complete agreement with the PN approximations. Inconsistencies encountered at low order approximations are primarily due to the different form of Fick’s Law or diffusion constant in two approximations which have already been explained in section three of this work. The other reason is the use of Mark boundary condition in TN approximations for this part, whereas Marshak boundary condition was used in PN calculations of Yildiz’s work. We deliberately used Mark condition instead of Marshak condition to see how the former works for different order approximations. As a result, it has been observed that Marshak condition produces faster convergent results than the Mark condition as seen from Table 1. Another important distinction appears for the greater values of backward scattering probability b, especially for the limiting value of cb ! 1 with b1 = 0.3 as seen from
ð2 dm0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi T n ðlÞT m ðlÞ dl 1 l2 0 p
1 sinððm þ nÞp=2Þ ðsinðm nÞp=2Þ þ ¼ pð1 þ dm0 Þ ðm þ nÞ ðm nÞ
I m;n ¼ ð1Þ
mþn
ð47Þ Eq. (45) takes the form N X
I m;2k1
m¼0
ðNX þ1Þ=2
b‘ am ðm‘ Þ½ea=m‘ þ ð1Þm ea=m‘ ¼ 0
ð48Þ
‘¼0
Eq. (48) could be written in the matrix form as following: Mkm ðaÞbm ¼ 0; k ¼ 1; 2; . . . ; ðN þ 1Þ=2; m ¼ 1; 2; . . . ; ðN þ 1Þ=2
ð49Þ Mkm ðaÞ
where bm is the vector comprising elements bm and is the coefficient matrix with [(N + 1)/2]2 elements and whose components are given as Mkm ¼
N X
‘
I ‘;2m1 a‘ ðmm Þðea=mm þ ð1Þ ea=mm Þ;
‘¼0
k ¼ 1; 2; . . . ; ðN þ 1Þ=2;
m ¼ 1; 2; . . . ; ðN þ 1Þ=2 ð50Þ
For non-trivial solution of Eq. (49), det½Mkm ðaÞ ¼ 0 that is the criticality condition. 5. Numerical results and conclusion We will present some critical thicknesses computed using TN approximation and comparison with available literature values. First of all, it should be noted that using an embedding approach for fission and scattering kernels would result completely different spectrums in nature in
Table 1 Critical thickness 2a computed by TN method using Mark boundary condition and comparison with PN method using Marshak boundary condition for c = 1.2, b = 0 and different degrees of forward and linearly anisotropic scattering a
b1
T3
P3 (Yildiz, 1998)
T9
P9 (Yildiz, 1998)
T15
P15 (Yildiz, 1998)
0.00
0.3 0.0 0.3
2.40930 2.65376 3.00758
2.37181 2.60400 2.93726
2.35384 2.58512 2.91670
2.35023 2.58076 2.91098
2.35042 2.58096 2.91124
2.34973 2.58009 2.91004
0.20
0.3 0.0 0.3
2.60514 2.85760 3.21840
2.54790 2.78353 3.11661
2.51752 2.75152 3.08210
2.51227 2.74530 3.07415
2.51252 2.74560 3.07450
2.51161 2.74442 3.07287
0.40
0.3 0.0 0.3
2.87114 3.13118 3.49570
2.77465 3.00915 3.33257
2.71520 2.94548 3.26320
2.70590 2.93484 3.25021
2.70614 2.93526 3.25078
2.70483 2.93354 3.24841
0.70
0.3 0.0 0.3
3.57882 3.84597 4.20043
3.28121 3.48333 3.73973
3.03698 3.20864 3.42392
2.97131 3.13710 3.34434
2.95452 3.12088 3.32920
2.94978 3.11533 3.32244
0.8333
0.3 0.0 0.3
4.26134 4.52760 4.86556a
3.60274 3.74994 3.92139
3.11326 3.20610 3.30976a
2.73464 2.79771 2.86635
2.72236 2.78482 2.89747a
2.60780 2.66235 2.72216
a
All eigenvalues of the spectrum are real in this case.
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
Table 1. This case differs from all previous ones given in Table 1 in that all eigenvalues of the spectrum of the one-speed neutron transport operator obtained using TN approximation are symmetric real numbers for this case. On the contrary, except for this limiting case, set of eigenvalues in previous calculations has consisted of symmetric real numbers and one purely imaginary conjugate pair. We will indicate later another kind of anomaly in which all eigenvalues are purely imaginary conjugate pairs. This unphysical set of spectrum eigenvalues encountered in the limiting case of cb ! 1 with b1 = 0.3 let us to inspect the use of embedded kernel approach or somehow equivalently using the same fission kernel as the scattering kernel instead of using an isotropic fission kernel which is the most proper choice beyond question. Of course, same considerations hold for the calculations based on spectrums consistent with the Case’s spectrum but obtained not using an isotropic fission kernel. Table 2 gives computed critical thicknesses for the limiting case cb ! 1 with b1 = 0.3 obtained using TN approximation for an isotropic fission kernel. The backward–forward-isotropic scattering kernel
749
is preserved in these calculations. The first two rows of Table 2 show the computed critical thicknesses using TN method with embedded kernel approach (c = cs + cf = 1.2) for Mark and Marshak boundary conditions, respectively, and it is not possible to say that convergency is achieved through the first 15 order approximations. In fact, monotonically decreasing values of critical thicknesses are obtained as the order of approximation is increased for embedded kernel approach in TN method for this limiting case in contrary to the PN values of Yildiz’s (1998) work given in the third row of Table 2. However, when we use a separate and isotropic kernel for the fission, it becomes possible to get rid of the unphysical spectrum of eigenvalues, all of which are real in this case, and convergency even in the low order approximations is achieved as seen from Table 2. Calculations are carried out from the purely absorbing medium (cs = 0, cf = c = 1.2) limit to purely scattering medium limit (cs = 0.99, cf = 0.21). As seen from Table 2, critical thickness increases as the scattering ratio cs increases, or conversely fission ratio cf decreases, this is the expected result. In addition to the better convergency
Table 2 Critical thickness 2a computed by TN method using an isotropic fission kernel when all eigenvalues are real and comparison with PN method for c = 1.2, b = 0, b1 = 0.3, a = 0.8333 c = 1.2 (TN)-Mark c = 1.2 (TN)-Marshak c = 1.2 (PN)-Marshak (Yildiz, 1998) cs = 0.0, cf = 1.2 (TN)-Marshak cs = 0.3, cf = 0.9 (TN)-Marshak cs = 0.6, cf = 0.6 (TN)-Marshak cs = 0.8, cf = 0.4 (TN)-Marshak cs = 0.9, cf = 0.3 (TN)-Marshak cs = 0.99, cf = 0.21 (TN)-Marshak
N=3
N=5
N=7
N=9
N = 11
N = 13
N = 15
4.86556 4.12903 3.92139 2.60184 2.80986 3.08679 3.33359 3.4854 3.64416
4.02143 3.48955 3.50097 2.58972 2.78165 3.0203 3.20866 3.30843 3.39534
3.58488 3.12908 3.15135 2.58190 2.76903 2.99449 3.17764 3.23209 3.27963
3.30976 2.91284 2.86635 2.58114 2.76771 2.99075 3.15657 3.21027 3.23766
3.11562 2.75363 2.72127 2.58021 2.76635 2.98821 3.14542 3.19896 3.21433
2.96868 2.63623 2.721126 2.57988 2.76589 2.98740 3.14116 3.19486 3.20393
2.85277 2.54105 2.72113 2.57984 2.76544 2.98665 3.13887 3.19229 3.19764
Table 3 Critical thickness 2a computed by TN method and comparison with PN method using Marshak boundary condition for c = 1.2, a = 0 and different degrees of backward and linearly anisotropic scattering b
b1
T3
P3 (Yildiz, 1998)
T9
P9 (Yildiz, 1998)
T15
P15 (Yildiz, 1998)
0.10
0.3 0.0 0.3
2.30192 2.49204 2.74762
2.30440 2.49189 2.74378
2.28054 2.46668 2.71646
2.28019 2.46631 2.71603
2.27924 2.46516 2.71454
2.27969 2.46568 2.71518
0.20
0.3 0.0 0.3
2.23742 2.38962 2.58324
2.23691 2.38662 2.57695
2.20910 2.36263 2.54640
2.20873 2.35729 2.54599
2.20776 2.35473 2.54461
2.20822 2.35668 2.54522
0.40
0.3 0.0 0.3
2.10352 2.19558 2.30281
2.09371 2.18334 2.28767
2.04942 2.13740 2.24111
2.04897 2.13755 2.24067
2.04781 2.13630 2.23929
2.04841 2.13693 2.23996
0.70
0.3 0.0 0.3
1.84737 1.87682 1.90794
1.79546 1.82199 1.84995
1.60211 1.62464 1.64837
1.59587 1.62114 1.64492
1.58794 1.61064 1.63456
1.59177 1.61447 1.63840
0.8333
0.3 0.0 0.3
1.67927 1.69053 1.70203
1.58106 1.58891 1.59781
1.06022 1.06114 1.06370
1.03458 1.03696 1.03936
0.84732 0.84862 0.84992
0.92482 0.92652 0.92823
750
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
of separate kernels approach, it should be addressed that critical thickness varies non-negligibly depending on the combination of the (cs, cf). In Table 3, critical thicknesses 2a computed by TN method and comparison with PN results of Yildiz (1998) work using Marshak boundary condition for c = 1.2, a = 0 and different degrees of backward and linearly anisotropic scattering are shown. Similar to the only forward scattering case with different degrees of linearly anisotropic scattering discussed previously, TN approximation produces results which are almost in complete agreement with the PN except for high backward probability ranges and at low order approximations. But, it could be observed from Table 3 that high-order approximations yield in very close values of critical thicknesses in both methods. Table 4 presents computed critical thickness 2a using embedded TN method (c = cs + cf) for various values of collision parameter c, backward and forward scattering probabilities and different degrees of linearly anisotropic scattering in comparison with literature values of Yildiz (1998), Sahni et al. (1992) and Dahl and Sjo¨strand (1979). A complete agreement with the literature values could be observed in Table 4, since the all the parameters fall within the range where any anomaly in the spectrum might be encountered or high backward or forward scattering probabilities would exist for which TN approximation may not converge or converge slowly as exemplified previously. However, it should be recalled that if a separate and isotropic kernel for fission were used all the calculated critical thicknesses had been different. To illustrate another anomaly which could be encountered in the solution of the one-speed neutron transport equation, we will give one more sample case in which all eigenvalues of the spectrum are complex conjugate pairs.
Such a case in the spectrum occurs in the limiting case when c(a + b) > 1 in the embedded kernel approach. Table 5 presents a case study in which critical thicknesses are computed using an isotropic fission kernel in TN approximation for c = 1.6, a = 0.42, and b = 0.42. First row in Table 1 shows the critical thicknesses obtained using embedded kernel approach for various order approximations in TN approximation. It is seen from the table that as the order of approximation increases, monotonically decreasing critical thicknesses are obtained which illustrates the non-convergent behavior of TN with embedded kernel approach for which all eigenvalues of the spectrum are of complex conjugate pairs. Whereas, as seen from the second and following rows of Table 5 in which a separate and isotropic fission kernel have been employed in TN approximation, convergency is achieved even in the low order approximations. It is also observed that in the limit of purely absorbing medium, that is (cs = 0, cf = c = 1.6), Chebyshev polynomial approximation is superior to spherical harmonics approximation. But, for the purely scattering limit convergency is weak and computations require higher order approximations. As a result of foregoing calculations we conclude that TN approximation applies equally well as PN approximation in most cases except for some limiting cases in which degree of forward scattering probability is high or sum of forward and backward scattering probabilities exceed certain limiting values. TN approximation gives convergent results for absorbing dominant cases even in the low order approximations, but it also converges in most instances requiring higher order approximations. More importantly, anomalies encountered in some limiting cases such as spectrums with all eigenvalues are symmetric real numbers or purely imaginary conjugate pairs could easily be overcome
Table 4 Comparison of computed critical thicknesses using TN with literature values for various degree of anisotropy and cross-section parameters a, b
b1
c
Sahni et al. (1992)
Dahl and Sjo¨strand (1979)
a = 0.3, b = 0.0
0.3
1.5857300 1.0808500 1.6183300 1.0936500 1.6557200 1.1116600
1.12550 5.00548 1.14623 5.01758 1.17456 5.03748
1.02786 5.00325 1.03207 5.00445 1.02405 5.00669
1.00861 5.00165 1.01001 5.00213 1.01196 5.00308
1.00549 5.00133 1.00634 5.00169 1.00754 5.00242
1.00720 5.00146 1.00841 5.00185 1.01009 5.00265
1.00 5.00 1.00 5.00 1.00 5.00
– – – – – –
1.4801000 1.0579700 1.4998100 1.0640700 1.5215200 1.0716700
1.10043 4.99629 1.11135 5.00061 1.12462 5.00582
1.01369 5.00285 1.01513 5.0024 1.01688 5.00277
1.00306 5.00133 1.00343 5.00118 1.03966 5.00124
1.00193 5.00078 1.00216 5.00095 1.00245 5.00095
1.00270 5.00061 1.00303 5.00104 1.00345 5.00106
1.00 5.00 1.00 5.00 1.00 5.00
– – – – – –
3.8303067 1.1631293 1.1084678 1.0071358 1.3162908 1.1921797 1.1309166 1.6726350
0.26464 3.01633 4.00560 19.9911 2.05025 3.02966 4.01747 1.07972
0.22351 3.00300 4.00282 20.0024 2.00513 3.00444 4.00418 1.00983
0.21134 3.00139 4.00131 20.0011 2.00217 3.00395 4.00187 1.00284
0.20835 3.00108 4.00103 20.0008 2.00166 3.00207 4.00145 1.00202
0.21010 3.00121 4.00113 20.00091 2.00190 3.00190 4.00113 1.00266
– – – – – – – –
0.20 3.00 4.00 20.00 2.00 3.00 4.00 1.00
0.0 0.3 a = 0.0, b = 0.3
0.3 0.0 0.3
a = 0.0, b = 0.0
0.0
0.3
T3
T7
T11
T13
P13-Yildiz (1998)
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
751
Table 5 Critical thickness 2a when all eigenvalues are complex for different combinations of scattering–fission ratios and for different degrees of linearly anisotropic scattering (a = 0.42, b = 0.42) b1
T1
T3
T7
T11
T13
T15
c = 1.6
0.3 0.0 0.3
0.77288 0.78536 0.79856
0.43772 0.44090 0.44415
0.25246 0.25283 0.25321
0.17810 0.17811 0.17829
0.15529 0.15534 0.15539
0.13762 0.13765 0.13768
cs = 0.0, cf = 1.6
0.3 0.0 0.3
1.57072 1.57072 1.57072
1.08778 1.08778 1.08778
1.03171 1.03171 1.03171
1.02608 1.02608 1.02608
1.02548 1.02548 1.02548
1.02503 1.02503 1.02503
cs = 0.3, cf = 1.3
0.3 0.0 0.3
1.56593 1.57072 1.57555
1.06520 1.06775 1.07033
0.98787 0.99451 0.99817
0.97807 0.98022 0.98239
0.97692 0.97908 0.98125
0.97620 0.97836 0.98053
cs = 0.6, cf = 1.0
0.3 0.0 0.3
1.56120 1.57072 1.58044
1.03721 1.04191 1.04668
0.92507 0.92865 0.93229
0.90931 0.90852 0.91208
0.90206 0.90556 0.90911
0.90033 0.90513 0.90861
cs = 0.9, cf = 0.7
0.3 0.0 0.3
1.55651 1.57072 1.58538
1.00168 1.00794 1.01433
0.83097 0.83474 0.83857
0.78241 0.78577 0.78918
0.77158 0.77484 0.77819
0.76436 0.76763 0.77347
cs = 0.99, cf = 0.61
0.3 0.0 0.3
1.55511 1.57072 1.58687
0.98916 0.99574 1.00247
0.79507 0.79823 0.80179
0.73063 0.73352 0.73647
0.71404 0.72438 0.72761
0.70216 0.70487 0.70762
using a separate and isotropic kernel for fission. This physically more reasonable approach in which scattering and fission ratios are handled as two different parameters could produce different critical sizes in comparison with embedded kernels approach not only in TN approximation but also in PN approximation. So, using an isotropic kernel for fission is inevitable for accurate calculations. Acknowledgement I would like to express my gratitude to Dr. Mehmet Tombakoglu of Nuclear Engineering Department of Hacettepe University for his valuable suggestions and help. References ¨ ztu¨rk, H., 2006a. TN approximation to Anli, F., Yasßa, F., Gu¨ngo¨r, S., O neutron transport equation and application to critical slab program. Journal of Quantitative Spectroscopy 101, 129. ¨ ztu¨rk, H., 2006b. TN approximation to Anli, F., Gu¨ngo¨r, S., Yasßa, F., O reflected slab and calculation of the critical half thicknesses. Journal of Quantitative Spectroscopy 101, 135.
Atalay, M.A., 1997. The reflected slab and sphere criticality problem with anisotropic scattering in one-speed neutron transport theory. Progress in Nuclear Energy 3, 229. Dahl, E.B., Sjo¨strand, N.G., 1979. Eigenvalue spectrum of multiplying slabs and spheres for monoenergetic neutrons with anisotropic scattering. Nuclear Science and Engineering 69, 114. Davison, B., 1958. Neutron Transport Theory. Oxford University Press, London. Sahni, D.C., Sjo¨strand, N.G., Garis, N.S., 1992. Criticality and time eigenvalues for one-speed neutrons in a slab with forward and backward scattering. Journal of Physics D: Applied Physics 25, 1381. Sahni, D.C., Kulkarni, M., Sjo¨strand, N.G., 2004. Criticality of spheres by PN method. Annals of Nuclear Energy 31, 991. Williams, M.M.R., 2002. The integral form of the neutron transport equation for backward–forward scattering in bare spheres. Annals of Nuclear Energy 29, 777. Yasßa, F., Anli, F., Gu¨ngo¨r, S., 2006. Eigenvalue spectrum with Chebyshev polynomial approximation of the transport equation in slab geometry. Journal of Quantitative Spectroscopy 97, 51. Yildiz, C., 1998. Variation of the critical slab thickness with the degree of strongly anisotropic scattering in one-speed neutron transport equation. Annals of Nuclear Energy 25, 529. Yildiz, C., Alcan, E., 1995. The effect of strong anisotropic scattering on the critical sphere problem in neutron transport theory using a synthetic kernel. Annals of Nuclear Energy 22, 671.
NUCLEAR ENERGY Annals of Nuclear Energy 34 (2007) 743–751 www.elsevier.com/locate/anucene
Solution of one-speed neutron transport equation for strongly anisotropic scattering by TN approximation: Slab criticality problem Ayhan Yilmazer
*
Department of Nuclear Engineering, Hacettepe University, Beytepe, Ankara 06800, Turkey Received 27 March 2006; received in revised form 17 March 2007; accepted 18 March 2007 Available online 9 May 2007
Abstract In this study, a recently proposed version of Chebyshev polynomial approximation which was used in spectrum and criticality calculations by one-speed neutron transport equation for slabs with isotropic scattering is further developed to slab criticality problems for strongly anisotropic scattering. Backward–forward-isotropic model is employed for the scattering kernel which is a combination of linearly anisotropic and strongly backward–forward kernels. Further to that, the common approaches of using the same functional form for scattering and fission kernels or embedding fission kernel into the scattering kernel even in strongly anisotropic scattering is questioned for TN approximation via taking an isotropic fission kernel in the transport equation. As a starting point, eigenvalue spectrum of one-speed neutron transport equation for a multiplying slab with different degrees of anisotropy in scattering and for different cross-section parameters is obtained using Chebyshev method. Later on, the spectra obtained for different degree of anisotropies and cross-section parameters are made use of in criticality problem of bare homogeneous slab with strongly anisotropic scattering. Calculated critical thicknesses by Chebysev method are almost in complete agreement with literature data except for some limiting cases. More importantly, it is observed that using a different kernel (isotropic) for fission rather than assuming it equal to the scattering kernel which is a more realistic physical approach yields in deviations in critical sizes in comparison with the values presented in literature. This separate kernel approach also eliminates the slow convergency and/or non-convergent behavior of high-order approximations arising from unphysical eigenspectrum calculations. Ó 2007 Elsevier Ltd. All rights reserved.
1. Introduction A large collection of works exists on the one-speed neutron transport equation with isotropic, linearly anisotropic, and strongly anisotropic scattering for spheres and slabs either bare or reflected. Among the most widely used methods in those works is spherical harmonics or PN method whose earlier applications could be found in the monograph of Davison (1958). The critical slab problem in one-speed neutron transport theory was studied by Yildiz (1998) using an admixture of linearly anisotropic and strongly peaked backward and/or forward scattering kernel, and it was concluded that spherical harmonics method *
Tel.: +90 312 2977300; fax: +90 312 2992122. E-mail address: [email protected]
0306-4549/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2007.03.010
is very suitable for the numerical solution of neutron transport equation with linearly anisotropic scattering. As seen from the data of critical thicknesses given in this study, the effect of existence of strongly peaked backward and forward scattering components together in the scattering kernel resulted in missing data when obtaining critical thicknesses for large backward and forward scattering probabilities and collision parameter. Similar behavior at higher order PN approximation was reported in the study of Sahni et al. (2004) in which the criticality problem of the one-speed transport equation with isotropic scattering for reflected spheres was studied using PN approximation. In this study, it was reported that the method diverges as one increases the order PN approximation, but, the method produces results that agree with those of the integral equation. Williams (2002) developed an integral form of the
744
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
transport equation for bare spheres incorporating fission plus and arbitrary proportion of isotropic, backward and forward scattering. In this study, resulting integral equation was solved numerically and obtained critical radii were compared with those obtained by Yildiz and Alcan (1995) who accounted for strong anisotropy using a synthetic kernel. William’s work should be particularly noted since it uses a separate and isotropic kernel for fission instead of taking its functional form equal to the scattering kernel or embedding it into the scattering kernel in the transport equation. We will discuss later that this physically more convenient choice could also be used in Chebyshev polynomial-TN approximation or in any proper polynomial approximation as it was used in the integral transport equation formulation of William’s study. Among other various methods employed in solving onespeed neutron transport equation is singular eigenfunction or Case’s normal-mode expansion technique which was successfully applied reflected slab and sphere criticality problems with linearly anisotropic scattering by Atalay (1997). Atalay, in his study, has shown that normal-mode expansion technique serves to evaluate the complete spectrum of a critical system. Recently, a new version of Chebyshev polynomial expansion for the angular flux has been proposed to obtain eigenspectrum of one-speed neutron transport operator with isotropic scattering in a slab (Yasßa et al., 2006). Then, the so-called TN method was applied to bare and reflected critical slab problems with isotropic scattering (Anli et al., 2006a,b). One of our main objectives in this study is to extend Chebyshev method in slab criticality problem with anisotropic scattering which is an admixture of linearly anisotropic and strong forward–backward parts (backward–forward-isotropic model). When doing this, to compare the critical thicknesses computed by TN method with literature values we will initially use the common approach of using same kernel for fission with anisotropic kernel. Later on, to observe the effect of fission kernel on spectrum or equivalently slab critical size, we will use an isotropic kernel for fission and the backward–forward-isotropic model for the scattering kernel. It was reported in studies (Anli et al., 2006a,b) that TN approximation in criticality of bare and reflected slab problems with isotropic scattering applies equally well as the PN approximation. One presumably can predict that TN approximation could equally be applicable in criticality of bare and reflected slab problems with anisotropic scattering as well as the PN approximation since both Legendre and Chebyshev polynomials belong to a general class of polynomials known as Jacobi polynomials. However, it is still worthwhile to adopt Chebyshev polynomials expansion for the angular part of the neutron angular flux in the solution of the transport equation because of its comparatively simple algebraic manipulation and it is appealing to investigate whether it could alleviate the slowly convergent or divergent cases encountered in PN approximation. Of course, a more rigorous study would be to employ the most general form of the Jacobi polynomials expansion for the
angular flux and to discuss qualitatively on purely mathematical basis for which Jacobi polynomial best convergent solution is obtainable for a particular problem. 2. Basic equations To obtain an accurate solution to transport equation when anisotropy is taken into account is more difficult and involved and effect of anisotropy was studied by many authors. The one-speed neutron transport equation could written as Xrwðr;XÞþRT wðr;XÞ Z Z ¼ RS fs ðXX0 Þwðr;X0 ÞdX0 þmRF ff ðXX0 Þwðr;X0 ÞdX0 4p
4p
ð1Þ where w(r, X) is the angular flux, RT is the total cross-section, RS is the scattering cross section, RF is the fission cross-section, m is the mean number of neutrons per fission, fs and ff stand for distribution functions (kernel) for scattered and fission neutrons, respectively. Since fission kernel is an isotropic event in the laboratory reference frame, it is reasonable to assume it in the form ff(X Æ X 0 ) = 1/4p. It is commonly assumed that scattering function fs(X 0 Æ X) which is representative of anisotropy in the interactions of neutrons with matter could be developed into a Legendre polynomial series (Davison, 1958) of the argument X 0 Æ X as follows: fs ðX0 XÞ ¼
N 1 X ð2n þ 1Þbn P n ðX0 XÞ 4p n¼0
ð2Þ
where X 0 and X stand for the direction of the neutron velocity before and after scattering, respectively, with the normalization condition, and b1 is the mean cosine of scattering angle. For the linearly anisotropic scattering, it is enough to take only the first two terms of series expansion. However, we will use backward–forward-isotropic model given as fs ðX0 XÞ ¼
1ab a ð1 þ 3b1 X0 XÞ þ dðX0 X 1Þ 4p 2p b ð3Þ þ dðX0 X þ 1Þ 2p
for scattering kernel where a and b are forward and backward scattering probabilities in a collision, respectively. Such a kind of anisotropic scattering has been used by many investigators. The first term on the right hand side of Eq. (3) represents linearly anisotropic scattering, second and third terms represent forward and backward scattering, respectively. The backward–forward-isotropic model that we have used in this study was proposed by Fermi for inclusion of anisotropy into the transport equation and has been used by many other investigators; for example, Yildiz (1998) used this kernel to calculate critical slab thickness for various degrees of anisotropy and cross section parameters in PN. Using scattering kernel given by
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
Eq. (3), we will first calculate eigenspectrum of one-speed neutron transport equation for a homogeneous slab with anisotropic scattering and use resulting spectrum in the slab criticality problems. For one-dimensional case Eq. (3) takes the following form: 1 fs ðll0 Þ ¼ ð1 a bÞð1 þ 3b1 ll0 Þ þ adðl l0 Þ 2 þ bdðl þ l0 Þ
ð4Þ
and the one-speed neutron transport equation given by Eq. (1) becomes l
o ½wðx;lÞ þ RT wðx;lÞ ox Z 1 Z 1 fs ðll0 Þwðx;l0 Þdl0 þmRF ff ðll0 Þwðr;l0 Þdl0 ¼ RS 1
ð5Þ
1
where ff(ll 0 ) is the one-dimensional fission kernel for a slab and is of the form ff(ll 0 ) = 1/2. If we scale the distance in units of mean free path and define collision and fission ratios as cs ¼
Rs Rt
and
cf ¼
mRf Rt
ð6Þ
the transport equation given by Eq. (5) takes the following form: Z 1 o l ½wðx; lÞ þ wðx; lÞ ¼ cs fs ðll0 Þwðx; l0 Þ dl0 ox 1 Z 1 þ cf ff ðll0 Þwðr; l0 Þ dl0 ð7Þ 1
If we use the conventional approaches of either using the same kernel for scattering and fission or embedding isotropic fission kernel into the scattering kernel, or the scattering event is isotropic Eq. (7) will become Z 1 o l ½wðx; lÞ þ wðx; lÞ ¼ c fs ðll0 Þwðx; l0 Þ dl0 ð8Þ ox 1 where total number secondary neutrons generated per collision c is related to cs and cf as c ¼ cs þ cf
ð9Þ
3. Chebyshev polynomial (TN) approximation As we have already touched on that Legendre and Chebyshev polynomials belong to a more general class of polynomials known as Jacobi or hypergeometric polynomials which are solutions to the Jacobi differential equation. Jacobi polynomials are represented by P ðg;jÞ ðlÞ and form n a complete set of orthogonal system in the interval [1, 1] with respect to the weighting function (1 l)g(1 + l)j. For g = j = 0, P nð0;0Þ ðlÞ reduces to Legendre polynomial, and for g = j = 1/2, P nð1=2;1=2Þ ðlÞ reduces to the Chebyshev polynomials of the first kind Tn(l). The orthogonality relation of Jacobi polynomials are of the following form:
Z
745
1 g
1
j
ð1 xÞ ð1 þ xÞ P mðg;jÞ ðlÞP nðg;gÞ ðlÞ dl
¼
2gþjþ1 Cðm þ g þ 1ÞCðm þ j þ 1Þ dmn m!Cðm þ g þ j þ 1Þ 2m þ g þ j þ 1
ð10Þ
Following previous information, it could be inferred that angular flux could be expanded in any Jacobi polynomial series since any function defined in the interval [1, 1] satisfying certain requirements could be expressed in series expansion of an orthogonal set of functions which is complete. Apart from some mathematical concerns such as convergency and possible singularities over the domain, it seems choice of Jacobi polynomials are completely arbitrary for the representation of the angular flux on purely mathematical basis. But, in addition to these mathematical annoyances, we will indicate some physical inconsistencies which could arise, etc. assigning Tn(l) for series representation of angular flux results in a difference form of diffusion coefficient or Fick’s Law obtained using Pn(l) expansion. Notwithstanding the mathematical and physical considerations, we could expand angular flux in terms of Chebyshev polynomials as proposed by Anli et al. (2006a) as 1 X 1 wðx; lÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi ð2 dn0 ÞT n ðlÞ/n ðxÞ ð11Þ p 1 l2 n¼0 where Tn(l) is nth order Chebyshev polynomial of the first kind. The mth Chebyshev moment of the angular flux can be defined as follows: Z 1 1 T m ðlÞ X pffiffiffiffiffiffiffiffiffiffiffiffiffi /m ðxÞ ¼ ð2 dn0 ÞT n ðlÞ/n ðxÞdl ð12Þ 1 l2 n¼0 1 p If we introduce the expansions given by Eq. (11) into Eq. (5) and use isotropic kernel for fission and the kernel given by Eq. (4) for scattering kernel transport equation takes the following form: 1 X ð2 dn0 ÞT n ðlÞ o pffiffiffiffiffiffiffiffiffiffiffiffiffi l /n ðxÞ þ ð1 cS ða þ ð1Þn bÞ/n ðxÞ ox p 1 l2 n¼0 cs cf ¼ ð1 a bÞ½/0 ðxÞ þ 3b1 l/1 ðxÞ þ /0 ðxÞ ð13Þ 2 2 Eq. (13) may be arranged if we use the recurrence relation for Chebyshev polynomials given below 1 lT n ðlÞ ¼ ðð1 þ dn;0 ÞT nþ1 ðlÞ þ T n1 ðlÞÞ 2 as " 1 X ð2 dn0 Þ½ð1 þ dn0 ÞT nþ1 ðlÞ þ T n1 ðlÞ o pffiffiffiffiffiffiffiffiffiffiffiffiffi ½/ ðxÞ ox n 2p 1 l2 n¼0 ð2 dn0 ÞT n ðlÞ n pffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 cs ða þ ð1Þ bÞ/n ðxÞ þ p 1 l2 cs cf ¼ ð1 a bÞ½/0 ðxÞ þ 3b1 l/1 ðxÞ þ /0 ðxÞ 2 2
ð14Þ
ð15Þ
After multiplying Eq. (15) by Tm(l) and integrating over l 2 (1, 1) and using the orthogonality relation of Chebyshev polynomials given below
746
Z
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
1
1
ð2 dn0 ÞT m ðlÞT n ðlÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi dl ¼ dnm p 1 l2
ð16Þ
we obtain 1d m ½ð1þdm0 Þ/mþ1 ðxÞþ/m1 ðxÞþð1cs ðaþð1Þ bÞ/m ðxÞ 2 dx ( " #) ð1þð1Þm Þ ð1þð1Þmþ1 Þ ¼cs ð1abÞ / ðxÞ3b1 /1 ðxÞ 2ðm2 1Þ 0 2ðm2 4Þ cf
ð1þð1Þm Þ / ðxÞ 2ðm2 1Þ 0
ð17Þ
The quantities /m(x) (m = 0, 1, 2 . . .) could be referred as the Chebyshev moments of the angular flux similar to spherical harmonics moments in the PN approximation. It is obvious that the first two moments /0(x) and /1(x) are identical with the scalar flux and the net current. A discrepancy merely appears in the definition or restatement of the Fick’s Law as will be shown later. This could be interpreted as if a new definition for diffusion coefficient were introduced in TN approximation. Although the exact solution of the differential equations system in Eq. (17) is impossible, an approximate solution can be found by assuming that /N+1(x), say, is negligibly small. The resulting equations are called the TN approximation: 1 d ½ð1 þ dm0 Þ/mþ1 ðxÞ þ /m1 ðxÞ 2 dx þ ð1 cs ða þ ð1Þm bÞ/m ðxÞ ¼ cs ð1 a bÞ " # ð1 þ ð1Þm Þ ð1 þ ð1Þmþ1 Þ / ðxÞ 3b1 /1 ðxÞ 2ðm2 1Þ 0 2ðm2 4Þ cf
ð1 þ ð1Þm Þ / ðxÞ; 2ðm2 1Þ 0
m ¼ 0; 1; 2; . . . ; N 1
m¼N
ð18Þ
ð19Þ
If we use the common approach of embedding fission kernel into the scattering kernel, or use same kernels for both events as in the isotropic scattering, c = cs + cf would be introduced as a lumped parameter for secondary neutrons, and the coupled set of equations given by Eqs. (18) and (19) take the form 1 d ð1 þ dm0 Þ/mþ1 ðxÞ þ /m1 ðxÞ 2 dx þ ð1 cða þ ð1Þm bÞ/m ðxÞ ¼ cð1 a bÞ " # ð1 þ ð1Þm Þ ð1 þ ð1Þmþ1 Þ / ðxÞ 3b1 /1 ðxÞ ; 2ðm2 1Þ 0 2ðm2 4Þ m ¼ 0; 1; 2; . . . ; N 1
It should be reminded that the approach just indicated before obtaining Eqs. (20) and (21) has also been used in PN approximation somehow in spectrum and criticality calculations of slab and sphere for anisotropic scattering problems. T1, namely taking N = 1, approximation results in two coupled differential equations for scalar flux /0(x) and net current /1(x) as following: d ½/ ðxÞ þ ð1 cÞ/0 ðxÞ ¼ 0 ð22Þ dx 1 1 d ½/ ðxÞ þ ð1 cfð1 a bÞb1 þ ða bÞgÞ/1 ðxÞ ¼ 0 2 dx 0 ð23Þ Eq. (23) could be rearranged for isotropic scattering as /1 ðxÞ ¼
1 d ½/ ðxÞ 2ð1 cb1 gÞ dx 0
ð24Þ
introducing another version of Fick’s Law stating that diffusion coefficient D is equal to 1 2ð1 cb1 Þ
ð25Þ
In order to obtain eigenvalue spectrum, we employ the well known procedure of seeking a solution of the form
N
ð1 þ ð1Þ Þ /0 ðxÞ; 2ðN 2 1Þ
ð21Þ
D¼
1 d N ½/ ðxÞ þ ð1 cs ða þ ð1Þ bÞ/N ðxÞ ¼ cs ð1 a bÞ 2 dx " N 1 # ð1 þ ð1ÞN Þ ð1 þ ð1ÞN þ1 Þ /0 ðxÞ 3b1 /1 ðxÞ 2ðN 2 1Þ 2ðN 2 4Þ cf
1 d ½/ ðxÞ þ ð1 cða þ ð1ÞN bÞ/N ðxÞ ¼ cð1 a bÞ 2 dx" N 1 # N N þ1 ð1 þ ð1Þ Þ ð1 þ ð1Þ Þ /0 ðxÞ 3b1 /1 ðxÞ ; m ¼ N 2ðN 2 1Þ 2ðN 2 4Þ
/m ðxÞ ¼ am ðmÞex=m
ð26Þ
Substituting Eq. (26) into Eq. (18), we get following systems of algebraic equations: a1 ¼ mð1 ðcs þ cf ÞÞa0
ð27Þ
a2 ¼ a0 2m½ð1 ða bÞ ðcs þ cf Þð1 a bÞb1 a1 ð28Þ m
amþ1 ¼ am1 2mð1 cs ða þ ð1Þ bÞam cs fð1 a bÞ " # m mþ1 ð1 þ ð1Þ Þ ð1 þ ð1Þ Þ a0 3b1 a1 2ðm2 1Þ 2ðm2 4Þ m
cf
ð1 þ ð1Þ Þ a0 ; 2ðm2 1Þ
36m6N
ð29Þ
to compute am with a0 = 1. If the equations given above are compatible, the determinant of coefficients must vanish. This is equivalent to use well-known closure relationship aN þ1 ¼ 0
ð20Þ
ðm ¼ 0; 1; . . . ; N Þ
ð30Þ
which enables us to compute eigenvalues. To obtain eigenvalue spectrum for isotropic scattering or for embedding approach, applying similar steps used above or just putting c = cs + cf in Eqs. (27)–(29), below system of algebraic equations could be obtained
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
a1 ¼ mð1 ðcs þ cf ÞÞa0 a2 ¼ a0 2m½ð1 cða bÞ cð1 a bÞb1 a1
ð31Þ ð32Þ
T N þ1 ðlÞ ¼ 2
ð33Þ
with the closure relationship aN þ1 ¼ 0
ð34Þ
which enables us to compute eigenvalues.
In this work, we used Mark and Marshak boundary conditions for the slab criticality. Mark boundary condition is merely used only for one illustrative case to demonstrate that Marshak boundary condition gives more accurate results than those of Mark for small N (Davison, 1958). Rest of the computations is based on Marshak boundary conditions. It could be inspected through Eqs. (27)–(29) that a‘ (m) = (1)‘a‘(m). Hence, for TN approximation with N odd, there are N + 1 roots ±m‘, corresponding to a solution for the mth Chebyshev moment of angular flux as /m ðxÞ ¼
m x=m‘
b‘ am ðm‘ Þ½ex=m‘ þ ð1Þ e
ð35Þ
‘¼0
When we use Eq. (35) in the TN approximation given by Eq. (11), angular flux becomes N X 1 wðx; lÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi ð2 dm0 ÞT m ðlÞ p 1 l2 m¼0
ðNX þ1Þ=2
ð36Þ
negative roots lk which correspond to the specified angles where incoming angular flux vanish on the boundary become
ð2k 1Þp lk ¼ cos ; ðN þ 1Þ=2 < k 6 N þ 1 ð40Þ 2ðN þ 1Þ or letting k ! k + (N + 1)/2 transformation
ð2k þ N Þp lk ¼ cos ; 1 6 k 6 ðN þ 1Þ=2 2ðN þ 1Þ
ð41Þ
Mkm ðaÞbm ¼ 0; k ¼ 1; . . . ; ðN þ 1Þ=2; m ¼ 1; . . . ; ðN þ 1Þ=2 ð42Þ where bm is the vector comprising elements bm and Mkm ðaÞ is the coefficient matrix with [(N + 1)/2]2 elements and whose components are given as Mkm ¼ 2
N X ‘¼0
1 ‘ a‘ ðmm Þðea=mm þ ð1Þ ea=mm ÞT ‘ ðlk Þ; d‘0 þ 1
k ¼ 1; . . . ; ðN þ 1Þ=2;
m ¼ 1; . . . ; ðN þ 1Þ=2
ð43Þ
Eq. (42) which is a linear system of equations in bm has a non-trivial solution only if det½Mkm ðaÞ ¼ 0 which is the criticality condition. If we use Marshak boundary condition which is based on zero incoming current on the vacuum boundary and could be approximated as Z 1 T 2k1 ðlÞwða; lÞ dl ¼ 0; k ¼ 1; 2; . . . ; ðN þ 1Þ=2
in the TN approximation using discrete eigenvalues given by Eq. (36), we obtain N Z 1 X ð2 dm0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi T 2k1 ðlÞT m ðlÞ dl 1 l2 m¼0 0 p
l>0
ðNX þ1Þ=2
ð45Þ
Using peculiarity property Tn(l) = (1)nTn(l) of Chebyshev polynomials, Eq. (45) could be rearranged as Z 1 N X ð2 dm0 Þ 2kþm1 pffiffiffiffiffiffiffiffiffiffiffiffiffi T 2k1 ðlÞT m ðlÞ dl ð1Þ 1 l2 0 p m¼0
ð38Þ
where lk are negative roots of TN+1(l) = 0. Since first type Chebyshev polynomials of order N + 1 could be represented as
m
b‘ am ðm‘ Þ½ea=m‘ þ ð1Þ ea=m‘ ¼ 0
‘¼0
ð37Þ
where a is the slab half thickness in units of mean free path. We could use the Mark boundary condition which implies zero incoming angular flux on the boundaries for specified values of l wða; lk Þ ¼ 0
ð39Þ
ð44Þ m
b‘ am ðm‘ Þ½ex=m‘ þ ð1Þ ex=m‘
where the coefficients b‘ should be determined from boundary conditions. Solution of transport equation with vacuum boundary condition must be symmetrical with respect to position as
wðx; lÞ ¼ wðx; lÞ
ð2k 1Þp l cos 2ðN þ 1Þ
0
‘¼0
wða; lÞ ¼ 0;
If we use boundary condition given by Eq. (38) in Eq. (36) we obtain following matrix equation for coefficients bm:
4. Criticality condition
ðNX þ1Þ=2
N þ1 Y k¼1
m
amþ1 ¼ am1 mf2ð1 cða þ ð1Þ bÞam ð1 þ ð1Þm Þ þc ð1 a bÞ a0 ðm2 1Þ ##) mþ1 ð1 þ ð1Þ Þ þ3b1 a1 ; 36m6N 2 ðm 4Þ
N
747
ðNX þ1Þ=2 ‘¼0
Since
m
b‘ am ðm‘ Þ½ea=m‘ þ ð1Þ ea=m‘ ¼ 0
ð46Þ
748
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
Z
1
comparison with separate kernels approach. Whenever such a discrepancy is encountered in the calculations it will be demonstrated that this anomaly can be alleviated using separate kernels approach. Table 1 presents calculated critical thicknesses using TN approximation in comparison with the ones obtained using PN in Yildiz’s work (1998) for different degrees of forward scattering and linearly anisotropic scattering. For the sake of comparability on the same basis, embedding approach or using the same kernel for the fission as the scattering kernel is employed initially in the TN method as in the work of Yildiz’s PN approximation. Results shown in Table 1 are for a collision parameter c = 1.2 and degree of backward scattering b = 0. Some important features of TN approximation can be drawn from Table 1. First, even though low order TN approximations are worse than the same order PN approximation in terms of convergence, higher order approximations are in almost complete agreement with the PN approximations. Inconsistencies encountered at low order approximations are primarily due to the different form of Fick’s Law or diffusion constant in two approximations which have already been explained in section three of this work. The other reason is the use of Mark boundary condition in TN approximations for this part, whereas Marshak boundary condition was used in PN calculations of Yildiz’s work. We deliberately used Mark condition instead of Marshak condition to see how the former works for different order approximations. As a result, it has been observed that Marshak condition produces faster convergent results than the Mark condition as seen from Table 1. Another important distinction appears for the greater values of backward scattering probability b, especially for the limiting value of cb ! 1 with b1 = 0.3 as seen from
ð2 dm0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi T n ðlÞT m ðlÞ dl 1 l2 0 p
1 sinððm þ nÞp=2Þ ðsinðm nÞp=2Þ þ ¼ pð1 þ dm0 Þ ðm þ nÞ ðm nÞ
I m;n ¼ ð1Þ
mþn
ð47Þ Eq. (45) takes the form N X
I m;2k1
m¼0
ðNX þ1Þ=2
b‘ am ðm‘ Þ½ea=m‘ þ ð1Þm ea=m‘ ¼ 0
ð48Þ
‘¼0
Eq. (48) could be written in the matrix form as following: Mkm ðaÞbm ¼ 0; k ¼ 1; 2; . . . ; ðN þ 1Þ=2; m ¼ 1; 2; . . . ; ðN þ 1Þ=2
ð49Þ Mkm ðaÞ
where bm is the vector comprising elements bm and is the coefficient matrix with [(N + 1)/2]2 elements and whose components are given as Mkm ¼
N X
‘
I ‘;2m1 a‘ ðmm Þðea=mm þ ð1Þ ea=mm Þ;
‘¼0
k ¼ 1; 2; . . . ; ðN þ 1Þ=2;
m ¼ 1; 2; . . . ; ðN þ 1Þ=2 ð50Þ
For non-trivial solution of Eq. (49), det½Mkm ðaÞ ¼ 0 that is the criticality condition. 5. Numerical results and conclusion We will present some critical thicknesses computed using TN approximation and comparison with available literature values. First of all, it should be noted that using an embedding approach for fission and scattering kernels would result completely different spectrums in nature in
Table 1 Critical thickness 2a computed by TN method using Mark boundary condition and comparison with PN method using Marshak boundary condition for c = 1.2, b = 0 and different degrees of forward and linearly anisotropic scattering a
b1
T3
P3 (Yildiz, 1998)
T9
P9 (Yildiz, 1998)
T15
P15 (Yildiz, 1998)
0.00
0.3 0.0 0.3
2.40930 2.65376 3.00758
2.37181 2.60400 2.93726
2.35384 2.58512 2.91670
2.35023 2.58076 2.91098
2.35042 2.58096 2.91124
2.34973 2.58009 2.91004
0.20
0.3 0.0 0.3
2.60514 2.85760 3.21840
2.54790 2.78353 3.11661
2.51752 2.75152 3.08210
2.51227 2.74530 3.07415
2.51252 2.74560 3.07450
2.51161 2.74442 3.07287
0.40
0.3 0.0 0.3
2.87114 3.13118 3.49570
2.77465 3.00915 3.33257
2.71520 2.94548 3.26320
2.70590 2.93484 3.25021
2.70614 2.93526 3.25078
2.70483 2.93354 3.24841
0.70
0.3 0.0 0.3
3.57882 3.84597 4.20043
3.28121 3.48333 3.73973
3.03698 3.20864 3.42392
2.97131 3.13710 3.34434
2.95452 3.12088 3.32920
2.94978 3.11533 3.32244
0.8333
0.3 0.0 0.3
4.26134 4.52760 4.86556a
3.60274 3.74994 3.92139
3.11326 3.20610 3.30976a
2.73464 2.79771 2.86635
2.72236 2.78482 2.89747a
2.60780 2.66235 2.72216
a
All eigenvalues of the spectrum are real in this case.
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
Table 1. This case differs from all previous ones given in Table 1 in that all eigenvalues of the spectrum of the one-speed neutron transport operator obtained using TN approximation are symmetric real numbers for this case. On the contrary, except for this limiting case, set of eigenvalues in previous calculations has consisted of symmetric real numbers and one purely imaginary conjugate pair. We will indicate later another kind of anomaly in which all eigenvalues are purely imaginary conjugate pairs. This unphysical set of spectrum eigenvalues encountered in the limiting case of cb ! 1 with b1 = 0.3 let us to inspect the use of embedded kernel approach or somehow equivalently using the same fission kernel as the scattering kernel instead of using an isotropic fission kernel which is the most proper choice beyond question. Of course, same considerations hold for the calculations based on spectrums consistent with the Case’s spectrum but obtained not using an isotropic fission kernel. Table 2 gives computed critical thicknesses for the limiting case cb ! 1 with b1 = 0.3 obtained using TN approximation for an isotropic fission kernel. The backward–forward-isotropic scattering kernel
749
is preserved in these calculations. The first two rows of Table 2 show the computed critical thicknesses using TN method with embedded kernel approach (c = cs + cf = 1.2) for Mark and Marshak boundary conditions, respectively, and it is not possible to say that convergency is achieved through the first 15 order approximations. In fact, monotonically decreasing values of critical thicknesses are obtained as the order of approximation is increased for embedded kernel approach in TN method for this limiting case in contrary to the PN values of Yildiz’s (1998) work given in the third row of Table 2. However, when we use a separate and isotropic kernel for the fission, it becomes possible to get rid of the unphysical spectrum of eigenvalues, all of which are real in this case, and convergency even in the low order approximations is achieved as seen from Table 2. Calculations are carried out from the purely absorbing medium (cs = 0, cf = c = 1.2) limit to purely scattering medium limit (cs = 0.99, cf = 0.21). As seen from Table 2, critical thickness increases as the scattering ratio cs increases, or conversely fission ratio cf decreases, this is the expected result. In addition to the better convergency
Table 2 Critical thickness 2a computed by TN method using an isotropic fission kernel when all eigenvalues are real and comparison with PN method for c = 1.2, b = 0, b1 = 0.3, a = 0.8333 c = 1.2 (TN)-Mark c = 1.2 (TN)-Marshak c = 1.2 (PN)-Marshak (Yildiz, 1998) cs = 0.0, cf = 1.2 (TN)-Marshak cs = 0.3, cf = 0.9 (TN)-Marshak cs = 0.6, cf = 0.6 (TN)-Marshak cs = 0.8, cf = 0.4 (TN)-Marshak cs = 0.9, cf = 0.3 (TN)-Marshak cs = 0.99, cf = 0.21 (TN)-Marshak
N=3
N=5
N=7
N=9
N = 11
N = 13
N = 15
4.86556 4.12903 3.92139 2.60184 2.80986 3.08679 3.33359 3.4854 3.64416
4.02143 3.48955 3.50097 2.58972 2.78165 3.0203 3.20866 3.30843 3.39534
3.58488 3.12908 3.15135 2.58190 2.76903 2.99449 3.17764 3.23209 3.27963
3.30976 2.91284 2.86635 2.58114 2.76771 2.99075 3.15657 3.21027 3.23766
3.11562 2.75363 2.72127 2.58021 2.76635 2.98821 3.14542 3.19896 3.21433
2.96868 2.63623 2.721126 2.57988 2.76589 2.98740 3.14116 3.19486 3.20393
2.85277 2.54105 2.72113 2.57984 2.76544 2.98665 3.13887 3.19229 3.19764
Table 3 Critical thickness 2a computed by TN method and comparison with PN method using Marshak boundary condition for c = 1.2, a = 0 and different degrees of backward and linearly anisotropic scattering b
b1
T3
P3 (Yildiz, 1998)
T9
P9 (Yildiz, 1998)
T15
P15 (Yildiz, 1998)
0.10
0.3 0.0 0.3
2.30192 2.49204 2.74762
2.30440 2.49189 2.74378
2.28054 2.46668 2.71646
2.28019 2.46631 2.71603
2.27924 2.46516 2.71454
2.27969 2.46568 2.71518
0.20
0.3 0.0 0.3
2.23742 2.38962 2.58324
2.23691 2.38662 2.57695
2.20910 2.36263 2.54640
2.20873 2.35729 2.54599
2.20776 2.35473 2.54461
2.20822 2.35668 2.54522
0.40
0.3 0.0 0.3
2.10352 2.19558 2.30281
2.09371 2.18334 2.28767
2.04942 2.13740 2.24111
2.04897 2.13755 2.24067
2.04781 2.13630 2.23929
2.04841 2.13693 2.23996
0.70
0.3 0.0 0.3
1.84737 1.87682 1.90794
1.79546 1.82199 1.84995
1.60211 1.62464 1.64837
1.59587 1.62114 1.64492
1.58794 1.61064 1.63456
1.59177 1.61447 1.63840
0.8333
0.3 0.0 0.3
1.67927 1.69053 1.70203
1.58106 1.58891 1.59781
1.06022 1.06114 1.06370
1.03458 1.03696 1.03936
0.84732 0.84862 0.84992
0.92482 0.92652 0.92823
750
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
of separate kernels approach, it should be addressed that critical thickness varies non-negligibly depending on the combination of the (cs, cf). In Table 3, critical thicknesses 2a computed by TN method and comparison with PN results of Yildiz (1998) work using Marshak boundary condition for c = 1.2, a = 0 and different degrees of backward and linearly anisotropic scattering are shown. Similar to the only forward scattering case with different degrees of linearly anisotropic scattering discussed previously, TN approximation produces results which are almost in complete agreement with the PN except for high backward probability ranges and at low order approximations. But, it could be observed from Table 3 that high-order approximations yield in very close values of critical thicknesses in both methods. Table 4 presents computed critical thickness 2a using embedded TN method (c = cs + cf) for various values of collision parameter c, backward and forward scattering probabilities and different degrees of linearly anisotropic scattering in comparison with literature values of Yildiz (1998), Sahni et al. (1992) and Dahl and Sjo¨strand (1979). A complete agreement with the literature values could be observed in Table 4, since the all the parameters fall within the range where any anomaly in the spectrum might be encountered or high backward or forward scattering probabilities would exist for which TN approximation may not converge or converge slowly as exemplified previously. However, it should be recalled that if a separate and isotropic kernel for fission were used all the calculated critical thicknesses had been different. To illustrate another anomaly which could be encountered in the solution of the one-speed neutron transport equation, we will give one more sample case in which all eigenvalues of the spectrum are complex conjugate pairs.
Such a case in the spectrum occurs in the limiting case when c(a + b) > 1 in the embedded kernel approach. Table 5 presents a case study in which critical thicknesses are computed using an isotropic fission kernel in TN approximation for c = 1.6, a = 0.42, and b = 0.42. First row in Table 1 shows the critical thicknesses obtained using embedded kernel approach for various order approximations in TN approximation. It is seen from the table that as the order of approximation increases, monotonically decreasing critical thicknesses are obtained which illustrates the non-convergent behavior of TN with embedded kernel approach for which all eigenvalues of the spectrum are of complex conjugate pairs. Whereas, as seen from the second and following rows of Table 5 in which a separate and isotropic fission kernel have been employed in TN approximation, convergency is achieved even in the low order approximations. It is also observed that in the limit of purely absorbing medium, that is (cs = 0, cf = c = 1.6), Chebyshev polynomial approximation is superior to spherical harmonics approximation. But, for the purely scattering limit convergency is weak and computations require higher order approximations. As a result of foregoing calculations we conclude that TN approximation applies equally well as PN approximation in most cases except for some limiting cases in which degree of forward scattering probability is high or sum of forward and backward scattering probabilities exceed certain limiting values. TN approximation gives convergent results for absorbing dominant cases even in the low order approximations, but it also converges in most instances requiring higher order approximations. More importantly, anomalies encountered in some limiting cases such as spectrums with all eigenvalues are symmetric real numbers or purely imaginary conjugate pairs could easily be overcome
Table 4 Comparison of computed critical thicknesses using TN with literature values for various degree of anisotropy and cross-section parameters a, b
b1
c
Sahni et al. (1992)
Dahl and Sjo¨strand (1979)
a = 0.3, b = 0.0
0.3
1.5857300 1.0808500 1.6183300 1.0936500 1.6557200 1.1116600
1.12550 5.00548 1.14623 5.01758 1.17456 5.03748
1.02786 5.00325 1.03207 5.00445 1.02405 5.00669
1.00861 5.00165 1.01001 5.00213 1.01196 5.00308
1.00549 5.00133 1.00634 5.00169 1.00754 5.00242
1.00720 5.00146 1.00841 5.00185 1.01009 5.00265
1.00 5.00 1.00 5.00 1.00 5.00
– – – – – –
1.4801000 1.0579700 1.4998100 1.0640700 1.5215200 1.0716700
1.10043 4.99629 1.11135 5.00061 1.12462 5.00582
1.01369 5.00285 1.01513 5.0024 1.01688 5.00277
1.00306 5.00133 1.00343 5.00118 1.03966 5.00124
1.00193 5.00078 1.00216 5.00095 1.00245 5.00095
1.00270 5.00061 1.00303 5.00104 1.00345 5.00106
1.00 5.00 1.00 5.00 1.00 5.00
– – – – – –
3.8303067 1.1631293 1.1084678 1.0071358 1.3162908 1.1921797 1.1309166 1.6726350
0.26464 3.01633 4.00560 19.9911 2.05025 3.02966 4.01747 1.07972
0.22351 3.00300 4.00282 20.0024 2.00513 3.00444 4.00418 1.00983
0.21134 3.00139 4.00131 20.0011 2.00217 3.00395 4.00187 1.00284
0.20835 3.00108 4.00103 20.0008 2.00166 3.00207 4.00145 1.00202
0.21010 3.00121 4.00113 20.00091 2.00190 3.00190 4.00113 1.00266
– – – – – – – –
0.20 3.00 4.00 20.00 2.00 3.00 4.00 1.00
0.0 0.3 a = 0.0, b = 0.3
0.3 0.0 0.3
a = 0.0, b = 0.0
0.0
0.3
T3
T7
T11
T13
P13-Yildiz (1998)
A. Yilmazer / Annals of Nuclear Energy 34 (2007) 743–751
751
Table 5 Critical thickness 2a when all eigenvalues are complex for different combinations of scattering–fission ratios and for different degrees of linearly anisotropic scattering (a = 0.42, b = 0.42) b1
T1
T3
T7
T11
T13
T15
c = 1.6
0.3 0.0 0.3
0.77288 0.78536 0.79856
0.43772 0.44090 0.44415
0.25246 0.25283 0.25321
0.17810 0.17811 0.17829
0.15529 0.15534 0.15539
0.13762 0.13765 0.13768
cs = 0.0, cf = 1.6
0.3 0.0 0.3
1.57072 1.57072 1.57072
1.08778 1.08778 1.08778
1.03171 1.03171 1.03171
1.02608 1.02608 1.02608
1.02548 1.02548 1.02548
1.02503 1.02503 1.02503
cs = 0.3, cf = 1.3
0.3 0.0 0.3
1.56593 1.57072 1.57555
1.06520 1.06775 1.07033
0.98787 0.99451 0.99817
0.97807 0.98022 0.98239
0.97692 0.97908 0.98125
0.97620 0.97836 0.98053
cs = 0.6, cf = 1.0
0.3 0.0 0.3
1.56120 1.57072 1.58044
1.03721 1.04191 1.04668
0.92507 0.92865 0.93229
0.90931 0.90852 0.91208
0.90206 0.90556 0.90911
0.90033 0.90513 0.90861
cs = 0.9, cf = 0.7
0.3 0.0 0.3
1.55651 1.57072 1.58538
1.00168 1.00794 1.01433
0.83097 0.83474 0.83857
0.78241 0.78577 0.78918
0.77158 0.77484 0.77819
0.76436 0.76763 0.77347
cs = 0.99, cf = 0.61
0.3 0.0 0.3
1.55511 1.57072 1.58687
0.98916 0.99574 1.00247
0.79507 0.79823 0.80179
0.73063 0.73352 0.73647
0.71404 0.72438 0.72761
0.70216 0.70487 0.70762
using a separate and isotropic kernel for fission. This physically more reasonable approach in which scattering and fission ratios are handled as two different parameters could produce different critical sizes in comparison with embedded kernels approach not only in TN approximation but also in PN approximation. So, using an isotropic kernel for fission is inevitable for accurate calculations. Acknowledgement I would like to express my gratitude to Dr. Mehmet Tombakoglu of Nuclear Engineering Department of Hacettepe University for his valuable suggestions and help. References ¨ ztu¨rk, H., 2006a. TN approximation to Anli, F., Yasßa, F., Gu¨ngo¨r, S., O neutron transport equation and application to critical slab program. Journal of Quantitative Spectroscopy 101, 129. ¨ ztu¨rk, H., 2006b. TN approximation to Anli, F., Gu¨ngo¨r, S., Yasßa, F., O reflected slab and calculation of the critical half thicknesses. Journal of Quantitative Spectroscopy 101, 135.
Atalay, M.A., 1997. The reflected slab and sphere criticality problem with anisotropic scattering in one-speed neutron transport theory. Progress in Nuclear Energy 3, 229. Dahl, E.B., Sjo¨strand, N.G., 1979. Eigenvalue spectrum of multiplying slabs and spheres for monoenergetic neutrons with anisotropic scattering. Nuclear Science and Engineering 69, 114. Davison, B., 1958. Neutron Transport Theory. Oxford University Press, London. Sahni, D.C., Sjo¨strand, N.G., Garis, N.S., 1992. Criticality and time eigenvalues for one-speed neutrons in a slab with forward and backward scattering. Journal of Physics D: Applied Physics 25, 1381. Sahni, D.C., Kulkarni, M., Sjo¨strand, N.G., 2004. Criticality of spheres by PN method. Annals of Nuclear Energy 31, 991. Williams, M.M.R., 2002. The integral form of the neutron transport equation for backward–forward scattering in bare spheres. Annals of Nuclear Energy 29, 777. Yasßa, F., Anli, F., Gu¨ngo¨r, S., 2006. Eigenvalue spectrum with Chebyshev polynomial approximation of the transport equation in slab geometry. Journal of Quantitative Spectroscopy 97, 51. Yildiz, C., 1998. Variation of the critical slab thickness with the degree of strongly anisotropic scattering in one-speed neutron transport equation. Annals of Nuclear Energy 25, 529. Yildiz, C., Alcan, E., 1995. The effect of strong anisotropic scattering on the critical sphere problem in neutron transport theory using a synthetic kernel. Annals of Nuclear Energy 22, 671.