- Email: [email protected]
TN method for the critical thickness of one-speed neutrons in a slab with forward and backward scattering
ARTICLE IN PRESS
Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216 www.elsevier.com/locate/jqsrt
TN method for the critical thickness of one-speed neutrons in a slab with forward and backward scattering Hakan O¨ztu¨rka,, Fikret Anlib, Su¨leyman Gu¨ngo¨ra a
C - ukurova University, Faculty of Science and Letters, Department of Physics, Adana 01330, Turkey b KSU., Faculty of Science and Letters, Department of Physics, K.Maras- 46100, Turkey Received 5 May 2006; received in revised form 1 December 2006; accepted 1 December 2006
Abstract The critical slab problem in the case of combination of forward and backward scattering with usual isotropic scattering is studied for one-speed neutrons in a uniform finite slab by using TN method based on Chebyshev polynomial approximation and Marshak boundary conditions. It is shown that TN method gives accurate results in one-dimensional geometry and is very efficient both in derivation of equations and rapid convergence. Numerical results obtained by TN method are compared against the PN method in tabular form, which agreed quite well. r 2006 Elsevier Ltd. All rights reserved. Keywords: Critical slab; Chebyshev polynomials; Neutron transport; Forward and backward scattering
1. Introduction Neutron transport theory is concerned with the migration of neutrons through bulk media. The transfer of radiation in stellar atmospheres and radiative equilibrium in them, the penetration of X-rays and g-rays through scattering media, cosmic ray showers, and other such problems are of interest to the physicists. These processes differ only in that the collision laws are different, and even these are very similar for some of the processes mentioned. This is particularly true of neutron transport theory and radiative transfer in stellar atmospheres. The equations of the constant cross-section approximation of neutron transport are identical to those for the corresponding case in radiative transfer theory. Many methods recently developed in neutron transport theory could be used in radiative transfer problem [1]. The determination of the critical size is usually the most important problem of neutron transport theory. Neutron population of the system is directly proportional to the dimension of the system, x in one dimension. There should be a value of x, let a, such that for xoa, the change in the neutron population is negative, and in the absence of independent sources as in this work, the neutron density will decay exponentially. For x4a, however, the change in the neutron population is positive, and the neutron density will increase exponentially.
Corresponding author. Tel.: +90 322 3386084; fax: +90 322 338 6070.
E-mail address: [email protected] (H. O¨ztu¨rk). 0022-4073/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2006.12.002
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
212
If x ¼ a, the system at which nuclear reactions are self-sustaining is then said to be critical and a is called the critical size of the dimension of the system [1]. The Chebyshev polynomial expansion to one-dimensional neutron transport in a slab was first used by Aspelund [4], Conkie [3], and Yabushita [5]. They explained the advantages of using Chebyshev polynomials instead of Legendre polynomials in neutron transport theory. Additional discussions of this subject can also be found in [14,17,18]. Since the neutrons are migrated anisotropically in a reactor, isotropic representation of them is naturally not sufficient. The forward and backward scattering of neutrons in a system can be characterized as anisotropic scattering, because the probabilities of scattering of neutrons to all directions are not equal to each other as in isotropic scattering. However, in the preceding papers, we have considered the isotropic scattering for computation of the critical slab thicknesses by using TN method [14,17,18]. As a first step to approach real system we investigate in this study, the effect of forward and backward scattering other than isotropic scattering on the critical thickness of the slab for a fixed criticality factor c ¼ 2.0 (number of secondary neutrons per collision). Thus, this work is an extension of those studies on the subject previously carried out by the same authors [17,18]. In the current study, Marshak-type vacuum boundary conditions are used. The results obtained by TN expansion are compared against those obtained by PN expansion and other works. Note that an alternative to PN method, known as the Case method, is extensively used in the literature for criticality calculations. Here, we do not present the details of PN method as they are available in [11] or elsewhere [6–10]. 2. Theory With conventional notation, the transport equation for neutrons of one speed in a homogeneous critical slab with the scattering kernel given in [12,13] and with no sources is Z qcðx; mÞ sS ð1 a bÞ 1 þ sT cðx; mÞ ¼ m cðx; m0 Þ dm0 þ sS ½a cðx; mÞ þ b cðx; mÞ, (1) qx 2 1 subject to following free space boundary and symmetry conditions: cða; mÞ ¼ 0,
(2a)
cðx; mÞ ¼ cðx; mÞ; m40.
(2b)
where a and b are the parameters varying over the range of 0pa, bp1, a+b p1 and give the fractions of neutrons which emerge from a collision in forward and backward directions, respectively. c(x,m) is the angular flux of the neutrons at position x traveling in direction m, cosine of the angle between the neutron velocity vector and the positive x-axis. sT and sS are the macroscopic total and scattering differential cross sections, respectively; c is the number of secondary neutrons per collision, csT ¼ sS. It is assumed that the slab material is homogeneous, 2a is the slab thickness, and the origin of the x-axis is chosen such that the slab extends from x ¼ a to x ¼ a. Since the relation between neutron transport and radiative transfer has been given briefly in the introduction: Eq. (1) is indeed identical to the radiative transfer equation. That is, neutron flux in Eq. (1) corresponds to the intensity (radiance) of the radiation field, where coordinate x can be interpreted as the optical variable, m as the cosine of the polar angle that specify the direction of propagation of the radiation in the medium, and c as the albedo for single scattering and the scattering kernel to the phase function. 2.1. Chebyshev polynomial expansion method (TN method) We consider the angular flux in terms of Chebyshev polynomials given in [14,17,18]. Inserting it into Eq. (1) by using recurrence and orthogonality relations of Chebyshev polynomials of the first kind [15], and multiplying the resulting equation with Tn(m) Chebyshev polynomial, and finally integrating the resulting equation over dm in the interval [1,1], one can obtain the equations of moments: df1 ðxÞ þ sT f0 ðxÞ ¼ sS f0 ðxÞ, dx
(3a)
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
213
df2 ðxÞ df0 ðxÞ þ þ 2½sT sS ða bÞ f1 ðxÞ ¼ 0, dx dx
(3b)
and in general, dfnþ1 ðxÞ dfn1 ðxÞ ½1 þ ð1Þn þ þ 2½sT sS ða þ ð1Þn bÞ fn ðxÞ ¼ sS ð1 a bÞ f0 ðxÞ; dx dx n2 1
nX2.
(3c)
The well-known solutions to Eqs. (3) are in the form of [1] fn ðxÞ ¼ G n ðvÞ expðsT x=vÞ,
(4)
and by substituting this into Eqs. (3), it is possible to obtain analytic expressions for all Gn(n)’s: G 1 ðvÞ ¼ vð1 cÞ,
(5a)
G 2 ðvÞ ¼ 2n2 ð1 cÞ ½1 cða bÞ 1,
(5b)
G nþ1 ðvÞ þ 2n½1 cða þ ð1Þn bÞ G n ðnÞ þ G n1 ðvÞ ¼
½1 þ ð1Þn vc ð1 a bÞ; nX2, n2 1
(5c)
where the normalizations G1(n) ¼ 0 and G0(n) ¼ 1 are used. As seen from Eqs. (5), the functions Gn(n) depend on n and one can obtain the discrete n eigenvalues by setting GN+1(n) ¼ 0 for various c, a and b values. Since fN+1(n) ¼ 0 in the TN approximation as in the PN approximation, we should have GN+1(nj) ¼ 0 in our method. With so computed discrete eigenvalues, the general solutions of Eq. (1) in the present approximation can be expressed as fn ðxÞ ¼
ðNþ1Þ=2 X
bk Gn ðvk Þ expðsT x=vk Þ þ ð1Þn expðsT x=vk Þ ,
(6)
k¼1
where bk’s are constants and can be determined from the boundary conditions, and the property of Gn(n) ¼ (1)n Gn(n) is used. 3. Critical slab Finding the critical size of a system means to solve the eigenvalue problem. A reactor cannot be critical if fewer neutrons are emitted than absorbed following each collision (co1). Then, it can be said that for a slab of finite thickness with c41, a reactor may be subcritical or critical, in which case there will be physical solutions to the transport equation. Thus, in this section we carry out a good estimate of critical thickness by requiring that the asymptotic flux go to zero at the extrapolated boundary [2]. For the criticality of the slab, we use Marshak-type boundary condition, and by following the same procedure as in [17,18], we obtain the critical equation: ( ) ðNþ1Þ=2 N X X n sT a=vj sT a=vj I k G 0 ðvj Þbj coshðsT a=vj Þ þ bj G n ðvj Þ e þ ð1Þ e I n;k ¼ 0, (7) j¼1
n¼1
k ¼ 1; 3; 5; . . . ; N, where I k and I n;k are given in [17,18]. For example for T1 approximation, we obtain an analytic solution for the critical half thickness by setting N ¼ 1 in Eq. (7): pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 2 2½1 cða bÞ 1 1 pffiffiffiffiffiffiffiffiffiffiffi a ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tan h . (8) p 1c sT 2ð1 cÞ½1 cða bÞ
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
214
4. Numerical results We carry out the calculations using vacuum boundary condition from the criticality condition given in Eq. (7) in order to determine the critical thickness as a function of the anisotropy parameters. All numerical results reported here are performed by Maple software in about 1 min on a Celeron 1100 MHz machine. For all calculations we take the total differential cross section to its normalized value sT ¼ 1 cm1. The permissible eigenvalues, nj, are calculated from Eqs. (5). Numerical results of our calculations are presented in Tables 1–3 for different values of anisotropy parameters a and b, and for c ¼ 2.0. Meanwhile, one can also fix the critical size of the system and seek the quantity of c. We also give the results of the critical thickness for different combinations of a+b in Table 3. In order to see the performance and the accuracy of the present computational method, the critical thickness for some selected value c is given in Table 3. In addition, in Table 3, we use the criticality factors c calculated by Sahni et al. [9] and Dahl and Sjo¨strand [16] for slabs of three thicknesses which are assumed as benchmark values, 1, 4
Table 1 The critical thickness 2a for c ¼ 2.0 with the degree of forward scattering a a
T1
P1
T3
P3
T5
P5
T7
P7
T9
P9
T11
P11
T13
P13
T15
P15
0.00 0.10 0.20 0.30 0.40 0.49
1.0366 1.0719 1.1118 1.1574 1.2104 1.2664
0.8241 0.8508 0.8808 0.9149 0.9541 0.9950
0.6962 0.7112 0.7274 0.7448 0.7635 0.7815
0.6841 0.6967 0.7097 0.7231 0.7367 0.7487
0.6556 0.6627 0.6689 0.6736 0.6758 0.6745
0.6477 0.6533 0.6575 0.6595 0.6583 0.6531
0.6361 0.6384 0.6382 0.6342 0.6248 0.6099
0.6342 0.6357 0.6344 0.6288 0.6170 0.5991
0.6297 0.6293 0.6252 0.6155 0.5976 0.5718
0.6284 0.6274 0.6224 0.6115 0.5916 0.5634
0.6261 0.6239 0.6170 0.6029 0.5780 0.5431
0.6256 0.6232 0.6158 0.6008 0.5745 0.5374
0.6246 0.6214 0.6127 0.5954 0.5651 0.5223
0.6243 0.6209 0.6118 0.5938 0.5622 0.5175
0.6237 0.6199 0.6098 0.5901 0.5551 0.5052
0.6236 0.6196 0.6094 0.5891 0.5531 0.5016
Table 2 The critical thickness 2a for c ¼ 2.0 with the degree of backward scattering b b
T1
P1
T3
P3
T5
P5
T7
P7
T9
P9
T11
P11
T13
P13
T15
P15
0.10 0.20 0.30 0.40 0.49
1.0050 0.9765 0.9506 0.9268 0.9071
0.8001 0.7784 0.7585 0.7403 0.7251
0.6661 0.6374 0.6097 0.5825 0.5584
0.6513 0.6195 0.5883 0.5571 0.5290
0.6178 0.5800 0.5414 0.5016 0.4649
0.6088 0.5693 0.5286 0.4861 0.4469
0.5943 0.5509 0.5046 0.4550 0.4087
0.5917 0.5471 0.4991 0.4473 0.3985
0.5855 0.5384 0.4862 0.4283 0.3732
0.5838 0.5359 0.4825 0.4226 0.3655
0.5806 0.5310 0.4746 0.4095 0.3465
0.5800 0.5300 0.4726 0.4058 0.3410
0.5784 0.5274 0.4678 0.3966 0.3263
0.5780 0.5267 0.4664 0.3938 0.3217
0.5771 0.5251 0.4633 0.3870 0.3098
0.5769 0.5248 0.4625 0.3850 0.3061
Table 3 Some critical slab thicknesses 2a with different values of c, a and b compared with the results of Sahni et al. [12] and Dahl and Sjo¨strand [16] a, b
c
T3
P3
T7
P7
T11
P11
T13
P13
Sahni
Dahl
a ¼ 0.2 b¼0 a¼0 b ¼ 0.2 a¼0 b¼0 a+b ¼ 0.1
1.6226400 1.0876700 1.5391200 1.0681200 1.1084678 1.0775700 1.5977300 1.0771300 1.4664300 1.0737600
1.1083 5.0080 1.0904 4.9997 4.0056 4.9997 1.0791 5.0008 1.2417 5.0180
1.0904 5.0217 1.0730 5.0133 4.0168 5.0145 1.0650 5.0153 1.2034 5.0265
1.0191 5.0035 1.0115 5.0022 4.0028 5.0028 1.0106 5.0027 1.0600 5.0031
1.0165 5.0036 1.0097 5.0022 4.0029 5.0029 1.0091 5.0027 1.0520 5.0031
1.0053 5.0016 1.0028 5.0009 4.0013 5.0014 1.0028 5.0012 1.0186 5.0013
1.0047 5.0015 1.0025 5.0009 4.0013 5.0013 1.0025 5.0011 1.0164 5.0012
1.0034 5.0012 1.0018 5.0007 4.0010 5.0011 1.0019 5.0009 1.0112 5.0009
1.0030 5.0011 1.0016 5.0006 4.0009 5.0010 1.0017 5.0008 1.0098 5.0008
1.00 5.00 1.00 5.00 — 5.00 1.00 5.00 1.00 5.00
— — — — 4.00 — — — — —
a+b ¼ 0.5
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
215
and 5 mfp (mean free path). As one can see from that table, the variation of the criticality factor c with a ¼ b (a+b) is still monotonic. That is, the c value decreases with increasing a+b. Furthermore, the critical size increases with the decrease of c, then leakage is supported by forward and backward scattering of neutrons to arrive at the critical values. The same manner is valid for the other values of a and b in Table 3. As a result, the critical thicknesses 2a calculated from our TN method are in good agreement with the reference values. 5. Conclusions In this work, we considered the critical slab problem for one-speed neutrons with forward and backward scattering in a homogeneous slab. We discussed the TN method and compared it against the well established PN approximation results. Results presented in Tables 1 and 2 show that for fixed value of the criticality factor c, the monotonic behavior of the critical thickness is observed in two cases. First, the criticality factor satisfies the condition c42, and second when all neutrons tend to scatter in backward direction along the slab. In contrast, it is known that slab thickness varies non-monotonically with increasing forward scattering. In these cases, the angular distribution of neutrons is strongly peaked in the direction along the slab. However, for the thinner slab (2a ¼ 0.62 in this work) the variation of the critical thickness is monotonic with increasing a. That means the criticality factor c, becomes only the controlling factor for the size when the thickness is larger. For the reasons mentioned above, in this work, the monotonic behavior of the critical thickness is observed both for forward and backward scattering. However, the presence of b 6¼ 0 may cause non-monotonic behavior for some scattering when the slab is thicker. This condition is equivalent to the requirement about co2. This surprising variation of the critical slab thickness may be understood by considering the change in the angular distribution in scattering kernel. In Tables 1 and 2, when a or b approaches to the value 1/c, the critical slab thickness is attained to the value of 2a ¼ 0.50 and 0.30 for forward and backward scattering, respectively. This is a sign that in the limit of a, b ¼ 1/c, the critical thickness falls sharply to zero, i.e. the distribution of neutrons in the system is completely dense in the normal plane x ¼ 0. In this work, we did not arrive at this sharp value in T15 approximation limit. But it is realized from Tables 1 and 2 that if the approximation order is increased N415, then the critical size would approach to zero in the limit of a, b ¼ 1/c. The calculated critical thicknesses are given in the tables with four decimal points. Although the T15 approximation is relatively accurate than the T13 approximation, there are no great improvements between them. Therefore, we find enough to give some numerical results with T13 approximation in Table 3. We also compare some of our results with those obtained by means of PN and other methods in Table 3. It is seen that our results are in very good agreement with those of Dahl and Sjo¨strand [16] and Sahni et al. [9] up to three or four digits. We conclude that the present TN method converges rapidly. It is relatively simple in derivation, and requires less computational effort for this type of transport problems. The lower order approximations of TN method yield better results for large systems. Furthermore, it can also be extended to more general problems with anisotropic scattering and to multiregion energy and time-dependent problems. In the near future, we plan to apply this method to time dependent problems in curvilinear geometries with linearly anisotropic scattering and different boundary conditions. References [1] [2] [3] [4] [5] [6]
Davison B. Neutron transport theory. New York: Oxford University Press; 1958. p. 1–142. Bell GI, Glasstone S. Nuclear reactor theory. New York: VNR Company; 1972. p. 64–126. Conkie WR. Polynomial approximations in neutron transport theory. Nucl Sci Eng 1959;6:260–6. Aspelund O. On a new method for solving the (Boltzmann) equation in neutron transport theory. PICG 1958;16:530–4. Yabushita S. Tschebyscheff polynomials approximation method of the neutron transport equation. J Math Phys 1961;2:543–9. Lee CE, Dias MP. Analytical solutions to the moment transport equations-I; one-group one-region slab and sphere criticality. Ann Nucl Energy 1984;11:515–30. [7] Aranson R. Critical problems for bare and reflected slabs and spheres. Nucl Sci Eng 1984;86:150–6. [8] Tezcan C, Sever R. The critical slab with the backward–forward-isotropic scattering model. Ann Nucl Energy 1985;12:573–6.
ARTICLE IN PRESS 216
H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
[9] Sahni DC, Sjo¨strand NG, Garis NS. Criticality and time eigenvalues for one-speed neutrons in a slab with forward and backward scattering. J Phys D: Appl Phys 1992;25:1381–9. [10] Atalay MA. The critical slab problem for reflecting boundary conditions in one-speed neutron transport theory. Ann Nucl Energy 1996;23:183–93. [11] Yildiz C. Variation of the critical slab thickness with the degree of strongly anisotropic scattering in one-speed neutron transport theory. Ann Nucl Energy 1998;25:529–40. [12] Ino¨nu¨ E. A theorem on anisotropic scattering. Transp Theory Stat Phys 1973;3:137–46. [13] Siewert CE, Williams MMR. The effect of anisotropic scattering on the critical slab problem in neutron transport theory using a synthetic kernel. J Phys D: Appl Phys 1977;10:2031–40. [14] Yasa F, Anli F, Gu¨ngo¨r S. Eigenvalue spectrum with Chebyshev polynomial approximation of the transport equation in slab geometry. JQSRT 2006;97:51–7. [15] Bell WW. Special functions for scientists and engineers. London: D. Van Nostrand Company; 1968. p. 187–96. [16] Dahl EB, Sjo¨strand NG. Eigenvalue spectrum of multiplying slabs and spheres for monoenergetic neutrons with anisotropic scattering. Nucl Sci Eng 1979;69:114–25. [17] Anli F, Yasa F, Gu¨ngo¨r S, O¨ztu¨rk H. TN approximation to neutron transport equation and application to critical slab problem. JQSRT 2006;101:129–34. [18] Anli F, Yasa F, Gu¨ngo¨r S, O¨ztu¨rk H. TN approximation to reflected slab and computation of the critical half thickness. JQSRT 2006;101:135–40.
Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216 www.elsevier.com/locate/jqsrt
TN method for the critical thickness of one-speed neutrons in a slab with forward and backward scattering Hakan O¨ztu¨rka,, Fikret Anlib, Su¨leyman Gu¨ngo¨ra a
C - ukurova University, Faculty of Science and Letters, Department of Physics, Adana 01330, Turkey b KSU., Faculty of Science and Letters, Department of Physics, K.Maras- 46100, Turkey Received 5 May 2006; received in revised form 1 December 2006; accepted 1 December 2006
Abstract The critical slab problem in the case of combination of forward and backward scattering with usual isotropic scattering is studied for one-speed neutrons in a uniform finite slab by using TN method based on Chebyshev polynomial approximation and Marshak boundary conditions. It is shown that TN method gives accurate results in one-dimensional geometry and is very efficient both in derivation of equations and rapid convergence. Numerical results obtained by TN method are compared against the PN method in tabular form, which agreed quite well. r 2006 Elsevier Ltd. All rights reserved. Keywords: Critical slab; Chebyshev polynomials; Neutron transport; Forward and backward scattering
1. Introduction Neutron transport theory is concerned with the migration of neutrons through bulk media. The transfer of radiation in stellar atmospheres and radiative equilibrium in them, the penetration of X-rays and g-rays through scattering media, cosmic ray showers, and other such problems are of interest to the physicists. These processes differ only in that the collision laws are different, and even these are very similar for some of the processes mentioned. This is particularly true of neutron transport theory and radiative transfer in stellar atmospheres. The equations of the constant cross-section approximation of neutron transport are identical to those for the corresponding case in radiative transfer theory. Many methods recently developed in neutron transport theory could be used in radiative transfer problem [1]. The determination of the critical size is usually the most important problem of neutron transport theory. Neutron population of the system is directly proportional to the dimension of the system, x in one dimension. There should be a value of x, let a, such that for xoa, the change in the neutron population is negative, and in the absence of independent sources as in this work, the neutron density will decay exponentially. For x4a, however, the change in the neutron population is positive, and the neutron density will increase exponentially.
Corresponding author. Tel.: +90 322 3386084; fax: +90 322 338 6070.
E-mail address: [email protected] (H. O¨ztu¨rk). 0022-4073/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2006.12.002
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
212
If x ¼ a, the system at which nuclear reactions are self-sustaining is then said to be critical and a is called the critical size of the dimension of the system [1]. The Chebyshev polynomial expansion to one-dimensional neutron transport in a slab was first used by Aspelund [4], Conkie [3], and Yabushita [5]. They explained the advantages of using Chebyshev polynomials instead of Legendre polynomials in neutron transport theory. Additional discussions of this subject can also be found in [14,17,18]. Since the neutrons are migrated anisotropically in a reactor, isotropic representation of them is naturally not sufficient. The forward and backward scattering of neutrons in a system can be characterized as anisotropic scattering, because the probabilities of scattering of neutrons to all directions are not equal to each other as in isotropic scattering. However, in the preceding papers, we have considered the isotropic scattering for computation of the critical slab thicknesses by using TN method [14,17,18]. As a first step to approach real system we investigate in this study, the effect of forward and backward scattering other than isotropic scattering on the critical thickness of the slab for a fixed criticality factor c ¼ 2.0 (number of secondary neutrons per collision). Thus, this work is an extension of those studies on the subject previously carried out by the same authors [17,18]. In the current study, Marshak-type vacuum boundary conditions are used. The results obtained by TN expansion are compared against those obtained by PN expansion and other works. Note that an alternative to PN method, known as the Case method, is extensively used in the literature for criticality calculations. Here, we do not present the details of PN method as they are available in [11] or elsewhere [6–10]. 2. Theory With conventional notation, the transport equation for neutrons of one speed in a homogeneous critical slab with the scattering kernel given in [12,13] and with no sources is Z qcðx; mÞ sS ð1 a bÞ 1 þ sT cðx; mÞ ¼ m cðx; m0 Þ dm0 þ sS ½a cðx; mÞ þ b cðx; mÞ, (1) qx 2 1 subject to following free space boundary and symmetry conditions: cða; mÞ ¼ 0,
(2a)
cðx; mÞ ¼ cðx; mÞ; m40.
(2b)
where a and b are the parameters varying over the range of 0pa, bp1, a+b p1 and give the fractions of neutrons which emerge from a collision in forward and backward directions, respectively. c(x,m) is the angular flux of the neutrons at position x traveling in direction m, cosine of the angle between the neutron velocity vector and the positive x-axis. sT and sS are the macroscopic total and scattering differential cross sections, respectively; c is the number of secondary neutrons per collision, csT ¼ sS. It is assumed that the slab material is homogeneous, 2a is the slab thickness, and the origin of the x-axis is chosen such that the slab extends from x ¼ a to x ¼ a. Since the relation between neutron transport and radiative transfer has been given briefly in the introduction: Eq. (1) is indeed identical to the radiative transfer equation. That is, neutron flux in Eq. (1) corresponds to the intensity (radiance) of the radiation field, where coordinate x can be interpreted as the optical variable, m as the cosine of the polar angle that specify the direction of propagation of the radiation in the medium, and c as the albedo for single scattering and the scattering kernel to the phase function. 2.1. Chebyshev polynomial expansion method (TN method) We consider the angular flux in terms of Chebyshev polynomials given in [14,17,18]. Inserting it into Eq. (1) by using recurrence and orthogonality relations of Chebyshev polynomials of the first kind [15], and multiplying the resulting equation with Tn(m) Chebyshev polynomial, and finally integrating the resulting equation over dm in the interval [1,1], one can obtain the equations of moments: df1 ðxÞ þ sT f0 ðxÞ ¼ sS f0 ðxÞ, dx
(3a)
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
213
df2 ðxÞ df0 ðxÞ þ þ 2½sT sS ða bÞ f1 ðxÞ ¼ 0, dx dx
(3b)
and in general, dfnþ1 ðxÞ dfn1 ðxÞ ½1 þ ð1Þn þ þ 2½sT sS ða þ ð1Þn bÞ fn ðxÞ ¼ sS ð1 a bÞ f0 ðxÞ; dx dx n2 1
nX2.
(3c)
The well-known solutions to Eqs. (3) are in the form of [1] fn ðxÞ ¼ G n ðvÞ expðsT x=vÞ,
(4)
and by substituting this into Eqs. (3), it is possible to obtain analytic expressions for all Gn(n)’s: G 1 ðvÞ ¼ vð1 cÞ,
(5a)
G 2 ðvÞ ¼ 2n2 ð1 cÞ ½1 cða bÞ 1,
(5b)
G nþ1 ðvÞ þ 2n½1 cða þ ð1Þn bÞ G n ðnÞ þ G n1 ðvÞ ¼
½1 þ ð1Þn vc ð1 a bÞ; nX2, n2 1
(5c)
where the normalizations G1(n) ¼ 0 and G0(n) ¼ 1 are used. As seen from Eqs. (5), the functions Gn(n) depend on n and one can obtain the discrete n eigenvalues by setting GN+1(n) ¼ 0 for various c, a and b values. Since fN+1(n) ¼ 0 in the TN approximation as in the PN approximation, we should have GN+1(nj) ¼ 0 in our method. With so computed discrete eigenvalues, the general solutions of Eq. (1) in the present approximation can be expressed as fn ðxÞ ¼
ðNþ1Þ=2 X
bk Gn ðvk Þ expðsT x=vk Þ þ ð1Þn expðsT x=vk Þ ,
(6)
k¼1
where bk’s are constants and can be determined from the boundary conditions, and the property of Gn(n) ¼ (1)n Gn(n) is used. 3. Critical slab Finding the critical size of a system means to solve the eigenvalue problem. A reactor cannot be critical if fewer neutrons are emitted than absorbed following each collision (co1). Then, it can be said that for a slab of finite thickness with c41, a reactor may be subcritical or critical, in which case there will be physical solutions to the transport equation. Thus, in this section we carry out a good estimate of critical thickness by requiring that the asymptotic flux go to zero at the extrapolated boundary [2]. For the criticality of the slab, we use Marshak-type boundary condition, and by following the same procedure as in [17,18], we obtain the critical equation: ( ) ðNþ1Þ=2 N X X n sT a=vj sT a=vj I k G 0 ðvj Þbj coshðsT a=vj Þ þ bj G n ðvj Þ e þ ð1Þ e I n;k ¼ 0, (7) j¼1
n¼1
k ¼ 1; 3; 5; . . . ; N, where I k and I n;k are given in [17,18]. For example for T1 approximation, we obtain an analytic solution for the critical half thickness by setting N ¼ 1 in Eq. (7): pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 2 2½1 cða bÞ 1 1 pffiffiffiffiffiffiffiffiffiffiffi a ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tan h . (8) p 1c sT 2ð1 cÞ½1 cða bÞ
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
214
4. Numerical results We carry out the calculations using vacuum boundary condition from the criticality condition given in Eq. (7) in order to determine the critical thickness as a function of the anisotropy parameters. All numerical results reported here are performed by Maple software in about 1 min on a Celeron 1100 MHz machine. For all calculations we take the total differential cross section to its normalized value sT ¼ 1 cm1. The permissible eigenvalues, nj, are calculated from Eqs. (5). Numerical results of our calculations are presented in Tables 1–3 for different values of anisotropy parameters a and b, and for c ¼ 2.0. Meanwhile, one can also fix the critical size of the system and seek the quantity of c. We also give the results of the critical thickness for different combinations of a+b in Table 3. In order to see the performance and the accuracy of the present computational method, the critical thickness for some selected value c is given in Table 3. In addition, in Table 3, we use the criticality factors c calculated by Sahni et al. [9] and Dahl and Sjo¨strand [16] for slabs of three thicknesses which are assumed as benchmark values, 1, 4
Table 1 The critical thickness 2a for c ¼ 2.0 with the degree of forward scattering a a
T1
P1
T3
P3
T5
P5
T7
P7
T9
P9
T11
P11
T13
P13
T15
P15
0.00 0.10 0.20 0.30 0.40 0.49
1.0366 1.0719 1.1118 1.1574 1.2104 1.2664
0.8241 0.8508 0.8808 0.9149 0.9541 0.9950
0.6962 0.7112 0.7274 0.7448 0.7635 0.7815
0.6841 0.6967 0.7097 0.7231 0.7367 0.7487
0.6556 0.6627 0.6689 0.6736 0.6758 0.6745
0.6477 0.6533 0.6575 0.6595 0.6583 0.6531
0.6361 0.6384 0.6382 0.6342 0.6248 0.6099
0.6342 0.6357 0.6344 0.6288 0.6170 0.5991
0.6297 0.6293 0.6252 0.6155 0.5976 0.5718
0.6284 0.6274 0.6224 0.6115 0.5916 0.5634
0.6261 0.6239 0.6170 0.6029 0.5780 0.5431
0.6256 0.6232 0.6158 0.6008 0.5745 0.5374
0.6246 0.6214 0.6127 0.5954 0.5651 0.5223
0.6243 0.6209 0.6118 0.5938 0.5622 0.5175
0.6237 0.6199 0.6098 0.5901 0.5551 0.5052
0.6236 0.6196 0.6094 0.5891 0.5531 0.5016
Table 2 The critical thickness 2a for c ¼ 2.0 with the degree of backward scattering b b
T1
P1
T3
P3
T5
P5
T7
P7
T9
P9
T11
P11
T13
P13
T15
P15
0.10 0.20 0.30 0.40 0.49
1.0050 0.9765 0.9506 0.9268 0.9071
0.8001 0.7784 0.7585 0.7403 0.7251
0.6661 0.6374 0.6097 0.5825 0.5584
0.6513 0.6195 0.5883 0.5571 0.5290
0.6178 0.5800 0.5414 0.5016 0.4649
0.6088 0.5693 0.5286 0.4861 0.4469
0.5943 0.5509 0.5046 0.4550 0.4087
0.5917 0.5471 0.4991 0.4473 0.3985
0.5855 0.5384 0.4862 0.4283 0.3732
0.5838 0.5359 0.4825 0.4226 0.3655
0.5806 0.5310 0.4746 0.4095 0.3465
0.5800 0.5300 0.4726 0.4058 0.3410
0.5784 0.5274 0.4678 0.3966 0.3263
0.5780 0.5267 0.4664 0.3938 0.3217
0.5771 0.5251 0.4633 0.3870 0.3098
0.5769 0.5248 0.4625 0.3850 0.3061
Table 3 Some critical slab thicknesses 2a with different values of c, a and b compared with the results of Sahni et al. [12] and Dahl and Sjo¨strand [16] a, b
c
T3
P3
T7
P7
T11
P11
T13
P13
Sahni
Dahl
a ¼ 0.2 b¼0 a¼0 b ¼ 0.2 a¼0 b¼0 a+b ¼ 0.1
1.6226400 1.0876700 1.5391200 1.0681200 1.1084678 1.0775700 1.5977300 1.0771300 1.4664300 1.0737600
1.1083 5.0080 1.0904 4.9997 4.0056 4.9997 1.0791 5.0008 1.2417 5.0180
1.0904 5.0217 1.0730 5.0133 4.0168 5.0145 1.0650 5.0153 1.2034 5.0265
1.0191 5.0035 1.0115 5.0022 4.0028 5.0028 1.0106 5.0027 1.0600 5.0031
1.0165 5.0036 1.0097 5.0022 4.0029 5.0029 1.0091 5.0027 1.0520 5.0031
1.0053 5.0016 1.0028 5.0009 4.0013 5.0014 1.0028 5.0012 1.0186 5.0013
1.0047 5.0015 1.0025 5.0009 4.0013 5.0013 1.0025 5.0011 1.0164 5.0012
1.0034 5.0012 1.0018 5.0007 4.0010 5.0011 1.0019 5.0009 1.0112 5.0009
1.0030 5.0011 1.0016 5.0006 4.0009 5.0010 1.0017 5.0008 1.0098 5.0008
1.00 5.00 1.00 5.00 — 5.00 1.00 5.00 1.00 5.00
— — — — 4.00 — — — — —
a+b ¼ 0.5
ARTICLE IN PRESS H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
215
and 5 mfp (mean free path). As one can see from that table, the variation of the criticality factor c with a ¼ b (a+b) is still monotonic. That is, the c value decreases with increasing a+b. Furthermore, the critical size increases with the decrease of c, then leakage is supported by forward and backward scattering of neutrons to arrive at the critical values. The same manner is valid for the other values of a and b in Table 3. As a result, the critical thicknesses 2a calculated from our TN method are in good agreement with the reference values. 5. Conclusions In this work, we considered the critical slab problem for one-speed neutrons with forward and backward scattering in a homogeneous slab. We discussed the TN method and compared it against the well established PN approximation results. Results presented in Tables 1 and 2 show that for fixed value of the criticality factor c, the monotonic behavior of the critical thickness is observed in two cases. First, the criticality factor satisfies the condition c42, and second when all neutrons tend to scatter in backward direction along the slab. In contrast, it is known that slab thickness varies non-monotonically with increasing forward scattering. In these cases, the angular distribution of neutrons is strongly peaked in the direction along the slab. However, for the thinner slab (2a ¼ 0.62 in this work) the variation of the critical thickness is monotonic with increasing a. That means the criticality factor c, becomes only the controlling factor for the size when the thickness is larger. For the reasons mentioned above, in this work, the monotonic behavior of the critical thickness is observed both for forward and backward scattering. However, the presence of b 6¼ 0 may cause non-monotonic behavior for some scattering when the slab is thicker. This condition is equivalent to the requirement about co2. This surprising variation of the critical slab thickness may be understood by considering the change in the angular distribution in scattering kernel. In Tables 1 and 2, when a or b approaches to the value 1/c, the critical slab thickness is attained to the value of 2a ¼ 0.50 and 0.30 for forward and backward scattering, respectively. This is a sign that in the limit of a, b ¼ 1/c, the critical thickness falls sharply to zero, i.e. the distribution of neutrons in the system is completely dense in the normal plane x ¼ 0. In this work, we did not arrive at this sharp value in T15 approximation limit. But it is realized from Tables 1 and 2 that if the approximation order is increased N415, then the critical size would approach to zero in the limit of a, b ¼ 1/c. The calculated critical thicknesses are given in the tables with four decimal points. Although the T15 approximation is relatively accurate than the T13 approximation, there are no great improvements between them. Therefore, we find enough to give some numerical results with T13 approximation in Table 3. We also compare some of our results with those obtained by means of PN and other methods in Table 3. It is seen that our results are in very good agreement with those of Dahl and Sjo¨strand [16] and Sahni et al. [9] up to three or four digits. We conclude that the present TN method converges rapidly. It is relatively simple in derivation, and requires less computational effort for this type of transport problems. The lower order approximations of TN method yield better results for large systems. Furthermore, it can also be extended to more general problems with anisotropic scattering and to multiregion energy and time-dependent problems. In the near future, we plan to apply this method to time dependent problems in curvilinear geometries with linearly anisotropic scattering and different boundary conditions. References [1] [2] [3] [4] [5] [6]
Davison B. Neutron transport theory. New York: Oxford University Press; 1958. p. 1–142. Bell GI, Glasstone S. Nuclear reactor theory. New York: VNR Company; 1972. p. 64–126. Conkie WR. Polynomial approximations in neutron transport theory. Nucl Sci Eng 1959;6:260–6. Aspelund O. On a new method for solving the (Boltzmann) equation in neutron transport theory. PICG 1958;16:530–4. Yabushita S. Tschebyscheff polynomials approximation method of the neutron transport equation. J Math Phys 1961;2:543–9. Lee CE, Dias MP. Analytical solutions to the moment transport equations-I; one-group one-region slab and sphere criticality. Ann Nucl Energy 1984;11:515–30. [7] Aranson R. Critical problems for bare and reflected slabs and spheres. Nucl Sci Eng 1984;86:150–6. [8] Tezcan C, Sever R. The critical slab with the backward–forward-isotropic scattering model. Ann Nucl Energy 1985;12:573–6.
ARTICLE IN PRESS 216
H. O¨ztu¨rk et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 105 (2007) 211–216
[9] Sahni DC, Sjo¨strand NG, Garis NS. Criticality and time eigenvalues for one-speed neutrons in a slab with forward and backward scattering. J Phys D: Appl Phys 1992;25:1381–9. [10] Atalay MA. The critical slab problem for reflecting boundary conditions in one-speed neutron transport theory. Ann Nucl Energy 1996;23:183–93. [11] Yildiz C. Variation of the critical slab thickness with the degree of strongly anisotropic scattering in one-speed neutron transport theory. Ann Nucl Energy 1998;25:529–40. [12] Ino¨nu¨ E. A theorem on anisotropic scattering. Transp Theory Stat Phys 1973;3:137–46. [13] Siewert CE, Williams MMR. The effect of anisotropic scattering on the critical slab problem in neutron transport theory using a synthetic kernel. J Phys D: Appl Phys 1977;10:2031–40. [14] Yasa F, Anli F, Gu¨ngo¨r S. Eigenvalue spectrum with Chebyshev polynomial approximation of the transport equation in slab geometry. JQSRT 2006;97:51–7. [15] Bell WW. Special functions for scientists and engineers. London: D. Van Nostrand Company; 1968. p. 187–96. [16] Dahl EB, Sjo¨strand NG. Eigenvalue spectrum of multiplying slabs and spheres for monoenergetic neutrons with anisotropic scattering. Nucl Sci Eng 1979;69:114–25. [17] Anli F, Yasa F, Gu¨ngo¨r S, O¨ztu¨rk H. TN approximation to neutron transport equation and application to critical slab problem. JQSRT 2006;101:129–34. [18] Anli F, Yasa F, Gu¨ngo¨r S, O¨ztu¨rk H. TN approximation to reflected slab and computation of the critical half thickness. JQSRT 2006;101:135–40.