- Email: [email protected]
Variation of the critical slab thickness with the degree of strongly anisotropic scattering in one-speed neutron transport theory
Pergamon
Ann. Nucl. Energy, Vol. 25, No. 8, pp. 529 540, 1998 © 1998 ElsevierScience Ltd. All rights reserved PiI: S0306-4549(97)00114-X Printed in Great Britain 0306-4549/98 $19.00+ 0.00
VARIATION OF THE CRITICAL SLAB THICKNESS WITH THE DEGREE OF STRONGLY ANISOTROPIC SCATTERING IN ONE-SPEED NEUTRON TRANSPORT THEORY C.Yfldiz Department of Physics, lstanbul Technical University, 80626, Maslak, Istanbul-TURKEY
(Received 30 September 1997)
Abstraet--The critical slab problem is studied in one-speed neutron transport theory using a linearly anisotropic kernel which combines forward and backward scattering. It is shown that, the recently observed non-monotonic variation of the thickness also exists in this strongly anisotropic case. In addition, the influence of the linear anisotropy on the critical thickness is analysed in detail. Numerical analysis for the critical thickness are performed using the spherical harmonics method and results are tabulated for selected illustrative cases as a function of different degrees of anisotropic scattering. Finally, some results are discussed and compared with those already obtained by other methods, the agreement is satisactory. The spherical harmonic method gives generally accurate results in one dimensional geometry, and it is very suitable for the numerical solution of the neutron transport equation with linearly anisotropic scattering. © 1998 Elsevier Science Ltd
I. INTRODUCTION For many years the neutron transport problems have been subject to several theoretical and numerical studies. Several workers (see, e.g. Carlvik, 1968; Kaper et al.,1974; Dahl and Sjostrand, 1979;Sahni and SjOstrand, 1990; Garis,1991 ) have studied the criticality problems of one-speed neutrons in homogeneous slabs and spheres using various methods which are applicable for exact analysis. On the other hand, higher order c eigenvalues ( e.g., the mean number of neutrons emitted per collision), have been calculated by Kemer et al (1967), and Kschwendt (1971) for critical slabs, and Dabl and SjOstrand (1979) for both critical slabs and spheres. Gaffs and Sjrstrand (1989) have given an approximate expression for the criticality eigenvalues Cn of order n for the case when cn is large. In all of these works, the results have been obtained with great accuracy in terms of the criticality factor c as a function of the dimensions of the systems, or vice versa. Two of the more popular methods of solving onedimensional critical problems are the Wiener - Hopf techniques (Davison, 1958; Williams, 1971) and the singular eigenfunction method (Case and Zweifel,1967). However in neutron transport theory explicit solution can be found only for very few problems. A number of general methods have been developed for solving the transport equation, and in fact, these methods have been applied to various problems. For some further investigation see Sahni and SjOstrand (1990). It is usually more difficult to obtain very accurate solution with anisotropic scattering. The neutron distribution problem can be solved by inserting into transport equation an appropriate scattering function, which represents the neutron interaction probability. In many of the above mentioned works, the scattering is generally assumed to be isotropic. However, many authors were interested in the effects of anisotropic scattering. (see Sahni
529
530
C. Ylldlz
and SjOsWand,1990). The problem of anisotropy in the interactions of neutrons with matter and of its effects on the critical size is one of the most classic problems of neutron transport theory. It is commonly assumed that the scattering function f(f~'.fl) can be developed into a Legendre polynomial series of the argumentll'.fl, i.e. (Davison, 1958) 1
N
f(n'.n)=~-~(2n~ + l) b~Pdn'.n).
(1)
Here ~ ' and f~ denote the direction of the neutron velocity before and after a scattering collision, respectively, with the normalisation condition bo=l. Further we note that bj is the average cosine of the scattering angle. If the anisotropy is assumed to be linear, only the first two terms are used in equation (1). The first calculation with linearly anisotropic scattering has been studied by Dahl and SjOstrand (1979) who applied Carlvik's method to the integral equation with vacuum boundary condition and have given eigenvalues up to the sixth order. Very little work has been done in the case of quadratic anisotropy, i.e. when only the first three terms are used in the series development of equation (1). This case has been studied by Lathrop and Leonard (1965). The most detailed calculations with linearly and quadratically anisotropic scattering have been described by Sj6strand (1976, 1978, 1980) and more fundamental studies were made by Protopopescu and SjOstrand (1981a, 1981b). As well known, the inclusion of anisotropic scattering effect into neutron transport equation leads to great mathematical complexity. For example, for linearly anisotropic scattering, it is in general difficult to obtain an integral equation for neutron flux alone. In addition, analytical studies of amsotropic scattering effects are usually limited for accurate critical calculations. However a different form of anisotropic scattering can be introduced as the following form
fin'.n)
1-a-p a 4a: ~ 8 ( n ' . n
-1)+
8(n'.n+l)
(2)
where a and fl are parameters which vary from 0 to 1 and give the fractions of neutrons which emerge from a collision in the forward and backward directions, respectively. Such a scattering model (backward-forwardisotropic model) was originally proposed by Fermi to describe neutrons (see e.g. Williams, 1966,1978). This model has been developed to a high degree of sophistication by InOnii (1973, 1976) and has been exploited in a practical manner by Siewert and Williams (1977) for solving critical slab problems, and by Williams (1985) for a multigroup formalism. It was pointed out by lnOnOand Tezcan (1978) and more accurately shown by Pomraning (1978) that increasing the degree of forward anisotropy a , the critical thickness first increases, and reaches a maximum then decreases. This type of anisotropic scattering has also been used to study the critical slab problem by Spiga and Vestrucci (1981). Similar works with the scattering function given in equation (2) were done by many workers (Garcia and Siewert, 1981; Tezcan, 1979, 1981; Tezcan and Yddiz, 1986; Yddiz 1990, 1994; Sahni and SjOstrand, 1991; Sharma et al., 1994; Ganapol and Komreich, 1996). Strongly forward scattering combined with ordinary linearly anisotropic scattering was studied by Sahni and SjOstrand (1991). In a later work, the full problem including also strongly backward scattering was treated in slab geometry by Sahni et a1.(1992). These authors have presented numerical results for the variation of the c-eigenvalues as a function of anisotropic scattering. Finally some exact solutions have been given by SjOstrand (1992). The purpose of the present work is to extend the above investigation of the critical slab problems to include linearly anisotropic scattering together with both forward and backward scattering, and to search the effect of the linearly anisotropic scattering on the critical size. Usually, the results have been expressed in terms of the criticality factor c as a function of the dimensions of the system, or vice versa. In this paper, we fix c as a real quantity and then look for the critical size of the system as a function of the anisotropy parameters. The scattering function is assumed to be of the form f(~'.f~)=~(l+3bj
f ~ ' . £ ' 1 ) + ~ ( £ ' 1 ' .f~ - 1 ) + ~ 8 ( [ 1 ' . [ 1 + 1 ) .
(3)
Here 0 < a ,fl< 1, a + f l < 1. One usually assumes that the scattering function f ( f g . ~ ) is non negative. This means that the coefficient bl in equation (3) must fulfil the condition [bl [ < 1/3. Then, the c-eigenvalues are real for negative values of b/and always real for large systems, about 2mfp thick.. On other hand, many c-eigenvalues
Variation of the critical slab thickness
531
are complex (see e.g. Dalai and SjOstr~d, 1979,1989; Sj6stmnd, 1989). More fundamental studies of the bel~viour of the criticality eigenvalues of one-speed transport operator with linearly anisotropic scattering were carried out by Sabni et al. (1995, 1997). For the criticality problems, generally, analytical techniques based on the integro-differential form of the transport equation are used. These methods yield highly accurate results, but they can only be employed for low order anisotropic scattering. For reactor calculations several terms may be considered in the series development of equation (1). The great complication arises when the higher order anisotropy is considered. In addition, when strongly backward scattering is included, it is neeessa~j to employ a tnmsformation (In6n~ 1973) to convert the equation into canonical form for solving exactly by analytical procedure. As a result of this difficulty, when there is no possibility of using analytical method, one can try more accurate methods for the solution. The complexity of the equations usually forces one to implement numerical methods. Thus such approaches have been made through a direct numerical attack on the integro-differential form of the neutron Iransport equation by using the eigenfunetion techniques. Among these, like the Sn method and the spherical harmonic method (also known PN method) have thus been developed for the design of reactors. For higher-order calculations, the SN methods becomes more popular (Bell and Glasstone, 1972; SjOstrand, 1980; Makai and Sj6strand, 1996). Nevertheless, the PN method is the most direct as well as one of the most powerful appruaehes to numerical transport calculations. One converts the integrodifferential form of the transport equation into a linear system of algebraic equations. These equations are most amenable to solution by a digital computer. Usually only a few orders in the Legendre expansion series given in equation (1) are manageable for exact analysis. For this reason, in the present work, the computations were performed by using the Pr~ method for the critical thickness. In addition, the present method is certainly more accurate in one-dimensional geometries for sufficient large spatial regions. We report some papers on the results of a variety of one-speed calculations in the PN method ( Davison, 1958; Nesthat et al., 1977; Lee and Dias, 1984; Lee et al., 1985; Garcia and Siewert, 1996; Khouaja et al., 1997).
H. THE BASIC EQUATIONS For the investigation of critical slabs, assuming the slab material to be homogeneous, we call 2a the slab thickness and choose x-coordinate axis perpendicular to the slab surface. The origin of the x-axis is chosen such that the slab extends from x=-a to x - a . The system is surrounded by vacuum, and there are no external sources. Then, with conventional notation, the transport equation for neutrons of one-speed can be written as (see e.g. Davison, 1958; Williams, 1971; or Duderstadt and Martin, 1979)
/~-~-x
+ E, Ud(x,/a) = cE, .[ S f(,u',,u)tP(X,la')dla'dqg 1
(4)
o
Here W(x,/a) is the neutron flux density depending on x and on the cosine of the angle between the positive x-axis and the neutron velocity vector. E t is the total cross section and c is the number of secondary neutrons per collision. We assume that E t and c do not vary in space or direction. The multiplicity c is related to the material cross sections by equation cE, = o E : + E,
(5)
where E : is the fission cross section, E, the scattering cross section and v the number of neutrons per fission. It should be noted that in equation (4) we have assumed that the same scattering function is applicable to both fission and the scattering process. The fission source term can be taken to be isotropic in the laboratory system of o o o r ~ s . However, these two problems are not very different. We measure distance in units of total mean free path. The dimensions characterising changes in reactor core composition are usually comparable to a neutron mf-p. If we scale, that is, define a dimensionless space variable such that x ~ xE,, we can write W(xg,,/a) ~ Ud(x,/a). Now, when the scattering function given in equation (3) is inserted into equation (4) we obtain
532
C. Ylldlz ~ I ] ( x .X //~
t2
r+1
subject to the following free space boundary and symmetry conditions W(a, # ) = 0
(7a)
q'(x,/~) = ~e(-x,/~), /~ >0.
(7b)
For an infinite symmetric slab the angular flux can be expanded in terms of the Legendre polynomials 2n+l
• (x, u) = ~ ( - T - ) ~ ( x ) e . (/~)
(8)
where N is odd number andkt the cosine of the angle between~ and the major axis. When the expansion of equation (8) is inserted into equation (6), and using orthogenality properties and recursion relations of the Legendre polynomial, it is found that, after some rearrangement of terms, the PN equations have the following form d~. , d~.+ l n - - - - ~ + (n + 1 ) - - ~ + (2n + 1)dO (x) = (2n + l){c(1 - a - fl)~0a.o + (1 - ct - fl)b,~,8.,
+ a c ~ . +,Bc(-l)"~,},
0 < n _
(9)
One may employ the well-known procedure of seeking a solution of the homogeneous equation (9) in the form (Davison, 1958)
• ,,(x) = G,,(o)exp(-x/u)
(10)
where the functions Gn( o ) are some constants. Each of the ~ . ( x ) defined in equation (10) will satisfy equation (9) provided that the characteristic PNequations
(n+ 1)G,,(v)+nG,,_,(o) - ( 2 n + 1){1-c(1- a - fl)(6.o +b,6,1)- 4 a +(-')'fl]}G.(u)=0
(11)
are satisfied. Equation (11) has a homogeneous matrix form [M(v)IG=0
(12)
where the coefficient m a t r i x is (N+1)x(N+I), and G=[Go,G~,...G,]r. The unknown constant vector G is determined from the normalisation G_~(o) = 0 and Go(u ) = I. The general solutions of equation (6) in the present approximation can be expressed as N+I
q).(x) = £ G . C o , ) [ p j F . ( - x / o , ) + q j F . ( x / u 2 ) ] j=l
(13)
Variation of the critical slab thickness
533
where Fn(x/vj)-=exp(-x/vi)and Pi and qi are arbitrary constants to be determined from the boundm3t eoaditiom of the problem. The essential idea of the P~ method is that l~q(v)=0, i.e., the permissible eigeavalue vj is the j-th positive zero of GN+l(v). Therefore, this should be identical with the detexmin~tal equation (12). Non-trivial solution will be obtained only when deqM(v)l=0.There are L=(N+I)/2 roots. All the roots vj ( for c>l ) are either purely real or purely imaginary. The spatial variation of the angular flux can be evaluated from the moment expansion in terms of Legendre polynomials given in equation (8). The numerical determination of these roots is done using the Newton-Raphson iterative technique. The computation time is extremely short. All calculations reported here were performed in about 2 see. on a pentium machine running at 120 MHz.
m . Tm~ SLAB c R r r I c A L r r Y C O N D m O N
The exact solution of the transport equation (4) must satisfy certain boundary conditions at the surface boundary between two media and the free- surface condition on the external surface. However, it is not possible to satisfy this requirement exactly with PN approximation for finite N. The difficulty arises from the fact that the exact boundary conditions are imposed over half the angular range - 1
fP~ (-/.t)W(a,-p)d/z= 0,
i=1,3,5,...N
(14)
0
The form of the criticality condition can be deduced from equation (14) for a finite slab of half-thickness a. Thus, we obtain the critical equation in the matrix form as follows [ M] (a) ]Pi =[ 0 ]
(15)
where Pi is a vector with elements Pi i= 1,2,...(N+ 1)/2 and the [(N + 1) / 2]: elements of matrix Mtj (a) are given by (N-l)/2
M~(a) = E tp(n,i)G2,(uj)cosh(a / uj) n=0
+ G2,_,(oj)sirda(a / oj) ,
(16)
where (4n + 1X-1)"÷J(2i - l)!(2n)!
q~(n,i) = 22"2"-J(2n- 2i+ l)(i +n)[n!(i- 1)!]2 "
(17)
For the matrix system in equation (15) to have a unique solution, the determinant of M~(a) must vanish, which is the criticality condition. To obtain the one-group neutron flux, an arbitrary normalisation is necessarily required to evaluate the eigenfunction O.(x). After determining the values of G~ the complete angular solution PN can be obtained.
534
C. Ylldiz
IV. N U M E R I C A L R E S U L T S AND DISCUSSION In order to examine the variation o f the critical thickness as a function o f the anisotropy parameters, we have c a m e d out the calculation using vacuum boundary condition from the criticality condition given in equation (15). The permissible root vj, for a given different values o f c is calculated from the determinantal equation(12). Our numerical results are presented in Tables 1-4 for different combinations o f the forward and backward scattering parameters. We give the critical size for various values o f bl ranging from -0.3 to +0.3 for chosen values o f 0=1.2 and 0=2. In some situations, one can also fix the critical size o f the system and then look for the quantity c which is the number o f secondaries per collision. (see e.g.Sahni and SjOstrand, 1990; Sahni et al., 1992). Table 1. The critical thickness 2a for 0=1.2 with degree o f forward scattering a and linearly anisotropic scattering bl as calculated P s approximation.
a
bl
P3
P5
P7
P9
Pll
P13
P15
0.00
-0.3 0.0 0.3
2.37181 2.60400 2.93726
2.35366 2.58453 2.91562
2,35106 2,58177 2,91237
2.35023 2.58076 2.91098
2.34975 2.58014 2.91014
2.34967 2.58002 2.90995
2.34973 2.58009 2.91004
0.20
--0.3 0.0 0.3
2.54790 2.78353 3.11661
2.51866 2.75196 3.08163
2,51363 2.74684 3,07611
2.51227 274530 3.07415
2.51157 2.74444 3.07299
2.51150 2.74431 3.07275
2.51161 2.74442 3.07287
0.40
--0.3 0.0 0.3
2.77465 3.00915 3.33257
2.72186 2.95105 3.26716
2.70928 2.93828 3,25398
2.70590 293484 3.25021
2.70455 2.93333 3,24836
2.70457 2.93328 3.24815
2.70483 2.93354 3.24841
0.65
--0.3 0.0 0.3
3.17711 3.39006 3.66598
3.03995 3.23431 3.48433
2.97765 3A6655 3.40948
2.95103 3.13867 3.38004
2.94021 3.12761 3.36869
2.93950 3.12720 3.36863
2.94052 3.12854 3.37035
0.70
--0.3 0.0 0.3
3.28121 3.48333 3.73973
3.10924 3.28762 3.51124
3.01728 3.18634 3.39769
2.97131 3.13710 3.34434
2.95115 3.11595 3.32199
2.94897 3.11405 3.32051
2.94978 3.11533 3.32244
0.8333
--0.3 0.0 0.3
3.60274 3.74994 3.92139
3.26866 3.37765 3.50097
2.97875 3~06067 3.15135
2.73464 2.79771 2.86635
2.60780 2.66236 2.72127
2.60780 2.66235 2.72126
2.60780 2.66235 2.72126
In Table 5 the results o f the critical thickness are given for a+fl ranging from 0 towards the limit a + f f =l/c as a function o f the different combination o f various anisotropy parameters. In order to see the performance and the accuracy o f the present computational method, the critical thickness for some selected value o f c are given in Table 6. We therefore take as reference values, the values given in the papers by Sahni et al. (1992) and Dahl and Sj6strand (1989). For the scattering function under investigation here the numerical and theoretical criticality calculation have been performed by many workers (see Sahni and Sj6strand, 1990) using various methods. R can be seen that our results for critical thickness versus a,fl and bt, follow the same general trend as those o f In6nO and Tezcan (1981), Pommning (1978), Sahni and Sj6strand (1991) and Sahni et al. (1992). The agreement was very good. It has been confirmed that the earlier observed non-monotonic variation o f the
Variation of the critical slab thickness
535
critical thickness exists in this case as well. The smprising trend previously reported is that the critical thickness for a fixed c first increases with forward anisotmpy parameters a(fl= 0) and then decreases. However, the variation of the critical thickness with a and linear anisotropy parameters bl depend upon how the value c has been chosen (Sahni and SjOstrand, 1991). We see from Table 1-4, as mentioned above, the variation will be non monotonic only
Table 2. The critical thickness 2a for c=2 with degree of forward scattering a and linearly anisotropic scattering b~ as calculated Ps approximation. t/
bl
P3
Ps
P7
P9
PH
P13
PLs
0.00
-0.3 0.0 0.3
0.64215 0.68410 0.73809
0.61261 0,65052 0,69890
0.59847 0.63506 0.68172
0.59211 0.62835 0.67457
0.58946 0.62562 0.67175
0.58925 0.62547 0.67168
0.58945 0.62575 0.67205
0.10
-- 0.3 0.0 0.3
0.65624 0.69667 0.74795
0,62185 0.65752 0.70225
0.60346 0.63727 0.67954
0.59426 0.62742 0.66887
0.59023 0.62319 0.66440
0.58979 0.62281 0.66410
0.58996 0.62307 0.66449
0.20
-- 0.3 0.0 0.3
0.67119 0.70971 0.75777
0.63091 0.66391 0.70443
0.60675 0.63716 0.67429
0.59316 0.62244 0.65814
0.58689 0.61576 0.65093
0.58615 0.61503 0.65024
0.58622 0.61519 0.65054
0.30
-- 0.3 0.0 0.3
0.68699 0.72313 076735
0.63955 0.66937 0.70510
0.60746 0.63380 0.66505
0.58695 0.61145 0.64039
0.57704 0.60076 0.62873
0.57599 0.59967 0.62761
0.57594 0.59969 0.62772
0.49
-- 0.3 0.0 0.3
0.71885 0.74870 0.78356
0.65319 0.67543 0.70068
0.59690 0.61382 0.63261
0.55012 0.56335 0.57781
0.52586 0.53742 0.54996
0.52571 0.53726 0.54979
0.52571 0.53726 0.54979
0.4999
-- 0.3 0.0 0.3
0.72055 0.74999 0.78428
065373 0.67553 0.70019
0.59575 0.61213 0.63027
0.54693 0.55954 0.57327
0.52156 0.53247 0.54425
0.52156 0.53247 054425
0.52156 0.53247 0.54425
if the critical thickness is larger than 0.67 for isotropic scattering ( a , fl = 0). This condition is equivalent to the requirement about c<2. For a certain degree of linearly anisotropic scattering hi, the critical size 2a is a function of two quantities, viz. the quantity c and the degree of forward scattering a . This is seen very clearly from the results in Table 1 and 2. As an example, for c=l.2 the critical thickness first increases with increasing forward scattering a for all values of bl ranging from -0.3 to +0.3 and reaches a maximum at about a = 0.71 and the thickness 2a=2.95. In that case, that is, if one wants to have a maximum at the value bf-0.3 and a =0.65 it is necessary to change the range of hi (from b1=-0.3 to b~=+0.3, in this case ) or one must reduce the quantity c. From the above discussion we see that for a given value of c=2, in Table 2, the maximum of thickness 2a (2a=_=_0.58996) turns out to be a =0.1 for b1=-0.3. While for bT0.0 and b1=+0.3, the critical thickness reduces monotonically with increasing a . That is, the behaviour of the critical thickness depending on bl is still non monotonic for c=2 and b~=-0.3. From this result, the variation of the critical thickness is monotonic for all values ofbt only when the quantity c is larger than about 1.9 or 2. In the case of backward scattering f l ; e 0 ( a = 0 ) the critical sizes for c=l.2 and c=2. decrease monotonically with increasing aniso~ropy. It was pointed out by Sahni et al. (1992), if the slab is thicker the presence of fl ;e 0 shows that the behaviour can be non-monotonic for some scattering. Similar conclusion are also confirmed from the numerical results of In6nO and Tezcan (1978) and Salmi and Sj6sUand (1991) and also Sahni ct al.(1992). In6m3 and Tezean (1978) obtained the critical slab thickness for c=1.2 using a variational method. The size parameter 2q corresponding to the maximum 2a works out be 3.11 with about a = 0 . 6 5 . We tried to arrive at these values from the curves given in their papers. Sahni and Sj6strand (1991) obtained a useful criterion for the monotonic or non-monotonic variation of c with a at fixed values of system dimensions. They found that the variation of the criticality factor is non-monotonic for all values of bl only when the thickness is larger than about
536
C. Ylldlz
0.7 mfp. This condition is equivalent to the requirement c<1.9. Consequently, we conclude that their results are in perfect agreement with our values in all considered cases. Another interesting ease is that the critical thickness seems to go to zero in the limit o f a =l/c. Note that when a becomes equal to l/c, the neutron dislribution is entirely concentrated in the normal plan x=0 (InrnO and Tezean, 1978). As can be seen from Tables 3 and 4, for c=2, when a or ,6 approaches the value l/c, the critical slab thickness is attained in the computational accuracy to the value o f 2a---0.5 and 0.3 for forward and backward scattering ease respectively, independently from the value o f b j . In the limit o f a,fl=l/c the critical thickness falls sharply to zero (In0na and Tezcan, 1978; Sahni et al., 1992). It is not possible to arrive at this sharp value in this limit using the present method. This limit value is computed with a precision o f nine significant figures for tt and ft. For example given a value o f fl = 0.4999 .... we obtain 2a---0.328. The other surprising result is the behaviour o f the critical thickness with anisotropy parameters a + f l . We see from Table 5 that the variation o f the critical thickness with increasing a + fl may sometimes be nonmonotonic. There is maximum at about a = fl=0.01 for c=1.2 and bi = -0.3 while for bt = -0.3 and 0.0, the variation is still monotonic. This depends on the values o f a ,fl and b I and also on the value c. However, for such systems below a certain size the variation is monotonic. It can be seen that the change occurs at about 2a-~ 1. One might think that there is certain value o f the thickness 2a for which the variation does not depend on the anisotropy parameters. For the thick slab, the variation o f the thickness for c<1.5 is still non-monotonic while for about c>1.6 the variation is monotonic. Such a behaviour o f the critical thickness is practically independent o f a + fl and bt-
Table 3. The critical slab thickness 2a for c=1.2 with various degree o f backward scattering fl and linearly anisotropic scattering bi as calculated Pr~ approximation. fl
bl
P3
P5
P7
P9
Pll
P13
P15
0.10
- 0.3 0.0 0.3
230440 249189 2.74378
2.28399 2.47034 2.72060
2.28106 2.46732 2.71730
2.28019 2.46631 2.71603
2.27970 2.46571 2.71526
2.27962 2.46560 2.71509
2.27969 2.46568 2.71518
0.20
-0.3 0.0 0.3
2.23691 2.38662 2.57695
2.21316 2.36185 2.55085
2.20968 2.35834 2.54722
2.20873 2.35729 2.54599
2.20822 2.35670 2.54527
2.20815 2.35660 2.54513
2.20822 2.35668 2.54522
0.40
-0.3 0.0 0.3
2.09371 2.18334 2.28767
205670 2.14521 2.24829
2.05042 2.13902 2.24218
2.04897 2.13755 2.24067
2.04834 2.13687 2.23993
204830 2.13682 2.23985
2.04841 2.13693 2.23996
0.70
-0.3 0.0 0.3
1.79546 1.82199 1.84995
1.66590 1.68907 1.71346
1.61681 1.63937 166315
1.59857 1.62114 1.64492
1.59171 1.61435 1.63820
1.59120 1.61387 1.63777
1.59177 1.61447 1.63840
0.80
-0.3 0.0 0.3
1.64053 1.65289 1.66556
1.41613 1.42436 1.43275
1.28466 1.29110 1.29765
1.20064 1.20616 1.21178
1.14702 1.15204 1.15715
113489 1.13983 1.14485
1.13110 1.13603 1.14105
0.8333
-0.3 0.0 0.3
1.58016 1.58891 1.59781
1.31841 1.32340 1.32845
1.15208 1.15538 1.15871
1.03458 1.03696 1.03936
0 95199 0.95384 0.95569
0.93253 0.93427 0.93601
0.92482 0.92652 0.92823
This was already indicated in earlier paper by Sahni et al. (1992) .They have given the value o f c<1.9 as the condition for the appearance o f the non-monotonic behaviour o f the criticality factor c. One example is shown in the paper o f Salmi." et al. (1992) that the variation o f c with a = ,8 is monotonic for a slab o f thickness 5 mfp, for various value o f bl ranging from 0 to +0.3. That is, the c value decreases with increasing ct+fl. But for b~ = 43.3 there is maximum in the curve at about at +fl=0.47. This depends on the values o f the a , p and bl and also on the
Variation of the critical slab thickness
537
thickness of the slab. For I mfp the variation of the criticality factor is always monotonic with increasing a + fl and for all values of bj. This should be eoml~tred to the value c=1.6 and the thickness about 2a<. 1.1 from our calculations using the P~ method. Thus the agreement between our values and with those of In6nO and Tezean (1978) and Sahni et al. (1992)is very good. Our results follow the same general trend as those of them
Table 4. The critical thickness 2a for c=2 with various degree of backward scattering fl and linearly anisotropic scattering bl as calculated Pr~ approximation.
P7
P9
PII
O10
--0.3 0.0 0.3
0.61856 0.65127 0.69127
0.58338 0.61202 0.64671
0.56559 0.59280 062568
0.55706 0.58382 061617
0.55332 0.57996 0.61218
0.55294 0.57963 0.61193
0.55314 0.57992 0.61231
0.20
--0.3 0.0 0.3
0.59453 0.61951 0.64874
0.55238 0.57327 0.59744
0.52939 0.54866 0.57087
0.51727 0.53592 0.55743
0.51151 0.52996 0.55123
0.51078 0.52926 0.55059
0.51094 0.52949 0.55091
0.30
--0.3 0.0 0.3
0.56974 0.58827 0.60914
0.51901 0.53350 0.54960
0.48875 0.50140 0.51539
0.47073 0.48253 0.49556
0.46115 0.47259 0.48520
0.45970 0.47112 0.48373
0.45964 0.47112 0.48380
0.40
--0.3 0.0 0.3
0.54391 0.55712 0.57153
0.48299 0.49237 0.50243
0.44319 0.45064 0.45858
0.41625 0.42264 0.42941
0.40000 0.40583 0.41199
0.39726 0.40301 0.40910
0.39671 0.40248 0.40858
0.49
--0.3 0.0 0.3
0.51965 0.52898 0.53889
0.44857 0.45444 0.46059
0.39904 0.40315 0.40740
0.36239 0.36547 0.36863
0.33843 0.34093 0.34349
0.33440 0.33681 0.33928
0.33329 0.33569 0.33815
0.4999
-- 0.3 0.0 0.3
0.51692 0.52587 0.53536
0.44470 0.45025 0.45605
0.39410 0.39792 0.40187
0.35641 0.35922 0.36210
0.33164 0.33389 0.33618
0.32750 0.32966 0.33187
0.32635 0.32849 0.33069
fl
bl
P3
P5
Pl3
PI5
We note that also from results in Table 5, for a given values of aandfl, the maximum points shift towards increasing a + fl for the decreasing value of c. Such behaviour of the maximum point may also be explained by the change in the magnitude of the value of c. The decrease of c increases the size, then leakage is favoured by the anisotropy to arrive at the critical values. For c=1.2, when the values a + f l approach the limit value of about 2a = 1.2; 1.7 and 2.0 for fl = a ; 0.15; 0.30 respectively. This happens independently of the value of bl. Finally, we see from Table 6 that the present method gives very accurate results and is well suited for the analysis of such systems. Numerical results of Dalai and SjOstrand (1979) and Sahni et al. (1992) practically coincide with our results to 4 or 5 significant figures especially for the thickness larger than 1.
V. CONCLUSION In this paper we have considered the slab criticality problem of the linear transport equation with forward, backward and linearly amsotropic scattering in a homogeneous slab. It is pointed out by many workers (see e.g. Sabni and SjOstrand,1990 ) that for fixed value of the factor c, the slab thickness varies non-monotonically with increasing the forward scattering under certain conditions. We observe that the addition of the linearly amsotropic scattering increases the critical slab thickness. That is, the leakage is less favoured by isotropic scattering than by the linear anisotropy. One characteristic feature of the linear anisotropie scattering is that the critical dimension of a
538
C. Ylldlz
system for a fixed c- value increases with the degree of anisotropy (see e.g. Sahni and SjOstrand,1990 ). In these thin slabs, the angular distribution of neutrons is very strongly peaked in the direction along the slab. However, for such systems below a certain size, especially for the thinner slab, (about 2a--0.67 in this work ) the variation of the critical thickness is monotonic with increasing a and for all values of b~. That is criticality factor c becomes only the controlling factor for the size when the thickness is larger than about c_=_1.9. In the ease of backward scattering (i.e. a = 0), the thickness decreases uniformly with increasing fl independently the value of bl. A special case of the scattering function in equation (1) has been studied in detail by Sahni e~ al. (1992). This surprising variation of the critical slab thickness may be understood by considering the change in the angular distribution (InOnO and Tezcan, 1978; and Pomraning, 1978; Sabni et al. 1992 ).
Table 5. The critical slab thickness 2a as a function of some selected anisotropy parameters a + fl and linearly anisotropie scattering bt for three value ofc as calculated Pl3 approximation. c=1.2 c~
c=1.5
c=1.6
bl fl=a
fl=0.15
fl=0.30
fl=a
fl=0.15 fl=0.30
fl=a
fl=0.15 fl 0.30 =
0.00
- 0.3 0.0 0.3
2.34967 2.58002 2.90995
2.24421 2.41067 2.62769
2.13252 2.24885 2.38984
1.12371 1.21194 1.33109
1.05134 1.11212 1.18771
0.96393 1.00265 1.04771
0.95516 1.02598 1.12028
0.88728 0.93511 0.99388
0.80304 0.83232 0.86599
0.01
- 0.3 0.0 0.3
2.34996 2.57565 2.89640
2.24984 2.41583 2.63195
2.13594 2.25137 2.39104
1.12113 1.20705 1.32225
1.05192 1.11220 1.18704
0.96218 1.00023 1.04440
0.95226 1.02109 111206
0.88723 0.93458 0.99265
0.80065 0.82930 0.86217
0.05
- 0.3 0.0 0.3
2.34938 2.55616 2.84102
2.27251 2.43642 2.64868
2.14887 2.26044 2.39450
1.10847 1.18500 1.28463
1.05356 1.11170 1.18337
0.95339 0.98856 1.02903
093825 0.99894 1.07683
0.88628 0.93157 0.98667
0.78921 0.81519 0.84471
0.33
- 0.3 0.0 0.3
2.05201 2.10806 2.17005
2.41891 2.55539 2.72314
2.13176 2.19949 227558
0.66691 0.67168 0.67658
0.99992 1.03356 1.07199
0,73009 0.73770 0.74562
0.42
- 0.3 0.0 0.3
1.21760 1.22168 1.22580
2.43809 2.55800 2.70183
1.98612 2.02982 2.07734
0.68
- 0.3 0.0 0.3
0.92990 0.95144 0.97506
0.81156 0.83417 0.85945 0.74085 0.75387 0.76779
1.71137 1.72439 1.73776
Another important conclusion of this calculation is that the PN method provides accurate results with very short computing time. All calculations reported here were performed in about 3 rain. on a pontium machine runnin~ at 120 MI-Iz. Inspection o f these numerical results in tables shows that the results are correct at least in the first 3 and 4 digits. The calculated critical thicknesses are rounded off at the 5th decimal point. The accuracy depends on the order o f PN approximation. The convergence of the method is very fast. Four stable decimal figures are always achieved with N _<11,13. The P15 approximation is relatively accurate than the P3 approximation. But there are no great improvements o f the P15 approximation over the Pi3. Therefore some numerical results are given in the P13 approximation. It is not necessary to go to further approximation. In addition, the accuracy oll the critical thickness depends also on the accuracy on the fixed values o f c. We have compared some o f our results with those obtained by means o f various methods. These are in very good agreement with those o f Dahl and Sj6strand (1979) and Sahin et al. (1992) up to the 5 or 6 digit, as shown in Table 6.
V a r i a t i o n o f the critical slab t h i c k n e s s
539
Consequently, we conclude that the present method is suited and also rapid for this type o f problems with less computational effort. For large systems the lower order approximations yield somewhat better results. In this ease the method works reasonably well. It is applicable to large system, and also to very small assemblies. In addition it can also be extended to more general problems with higher order anisotropy and to multiregion energy dependent problems.
Table 6. The some critical slab thickness obtained in various PN approximation with different value o f c, a , fl and degree o f linearly anisotropic scattering bl comparison with results o f Sahni et al. (1992)andDahlandSj0sWand(1979).
a,fl
b]
-0.30 ~Z= 0.3 0.00 fl= 0.0 0.30
-0.30 IOr= 0.0 0.00 i fl=~ 0.3 0.30
P3
P9
PIl
Pt3
8ahni
Dahl
1.5857300 1.0808500 1.6183300 1.0936500 1.6557200 1.1116600
1.10876 5.02327 1.12337 5.02910 1.14309 5.03908
1.01337 5.00225 1.01532 5.00291 1.01798 5.00417
1.00775 5.00160 1.00895 5.00204 1.01062 5.00294
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.08251 5.01155 1.08968 5.01365 1.09834 5.01595
1.00521 5.00106 1.00575 5.00157 1.00642 5.00170
1.00274 5.00069 1.00305 5.00113 1.00346 5.00117
1.00270 5.00061 1.00303 5.00104 1.00345 5.00106
1.00 5.00 1.00 5.00 100 5.00
-----
-----
--
--
0,00
3.8303067 0.25661 0.21472 0.21050 0.21010 1.1631293 3.02170 3.00193 3.00133 3.00121 1.1084678 4.01675 4.00182 4.00126 4.00113 1.0071358 20.01077 20.00152 20.00104 20.00091
-----
0.20 3.00 4.00 20.00
0.30
1.3162908 1.1921797 1.1309166 1.6726350
-----
2.00 3.00 4.00 1.00
5 = 0.0
/~=00
c
2.04010 3.02825 4.02268 1.06249
2.00297 3.00271 4.00257 1.00422
200206 3.00187 4.00178 1.00260
2.00190 3.00168 4.00159 1.00266
REFERENCES Bell G.I. and Glasstone S. (1972) Nuclear Reactor Theory. Van Nostrand-Reinhold,New York, NY. Carlvik I. (1968) Nucl. Sci. Engng., 31,295. Case K. M. and Zweifel P.F. (1967) Linear Transport Theory. Addison-Wesley, Reading, Mass. Dahl E.B. and SjOstrand N.G. (1979) Nuel. Sci. Engng., 69, 114. Dahl E.B. and Sj6strand N.G. (1989) Ann. Nucl.Energy 16, 527. Davison B. (1958) Neutron Transport Theory. Oxford University Press, London. Duderstadt J.J. and Martin W.R. (1979) Transport Theory. Wiley, New York. Ganapol B.D. and Komreieh D.E. (1996) Ann. Nucl. Energy 23,301. Gareia R . D . M and Siewert C.E. (1981) J. Phys. D: Appl. Phys.,14, L65. Gareia R D . M . and Siewert C.E. (1996) Ann. Nuel. Energy 23, 321. Garis N.S. (1991)Nuel. Sei. Engng., 107, 343. Gaffs N.S. and SjSslrand N,G. (1989) Kernteehnik 54, 183. InOnO E. (1973)Transp. Theory Statist. Phys., 3, 107. In6n• E. (1976) Phys. Fluids 19, 1332. In~nt~ E. and Tezean C. (1978) Transp. Theory Statist. Phys., 7, 51.
540
C. Yddlz
Kaper H.G., Lindeman A.J. and LeafG. K. (1974) Nuel. Sci. Engng., 54, 94. Kerner I.O., Kiesewetter H. and von Weber S. (1967) Kernenergie 60, 299. Khouaja H., Edwards D. R. and Tsoulfaridis N. (1997) Ann. blue. Energy 24, 518. Kschwendt H. (1971) Nucl. Sci. Engng., 44,423. Lathrop K.D. and Leonard A. (1965) Nucl. Sci. Engng., 22,115. Lee C. E. and Dias M. P. (1984) Ann. Nucl. Energy 11,515. Lee C. E., Pan W.C.P. and Dias M. P. (1985) Ann.Nucl. Energy 12, 613. Makai M. and SjOstrandN.G. (1996) Ann. Nucl. Energy 23, 397. Neshat K., Siewert C. E. and Ishiguro Y. (1977) Nu¢l. Sci. Engng., 62, 330. Pomraning G.C. (1978) Transp. Theory Statist. Phys.,7, 161. Protopopescu V. and Sj6strand N.G. (1981a) Progr. Nucl. Energy 7, 47. Protopopeseu V. and SjOstrand N.G. (1981b) On the solution of the Dispersion Equation for Monoenergetic Neutron Transport with Quadratically Anisotropie Scattering. CTH-RF-36, Chalmers Uni., Sweden. Sahni D.C. and Sj6strand N.G. (1990) Progr. Nuci. Energy 23, 241. Sahni D.C. and Sj6strand N.G. (1991) Transp. Theory Statist. Phys., 20, 339. Sahni D.C., Dahl B. and SjOstrandN.G. (1997) Ann. Nucl. Energy 24,135. Sahni D.C., Garis N.S. and Sj6strand N.G. (1995) Transp. Theory Statist. Phys., 24, 629. Sahni D.C., Sj6strand N.G. and Gaffs N.S. (1992) J. Phys. D: Appl. Phys., 25, 1381. Sharma A., Paranjape S.D., Kumar A. and Dwivedi S.R. (1994) Ann. Nucl. Energy 21,267. Siewert C.E and Williams M.M.R (1977)J. Phys. D: Appl. Phys., 10, 2031. Sj6strand N.G. (1976)J. Nucl.Sci.Tech., 13,81. SjOstrand NG. (1978) Atomkemenergie 31, 16. Sj6strand N.G. (1980) Ann. Nucl. Energy 7, 435. SjOstrand N.G. (1989) Approximate High-Order Eigenvalues in Two-Medium One-Speed Neutron Transport. CTH-RF-67, Chalmers Uni., Sweden. SjOstrand N.G. (1992) Exact Solution of Some One-Speed Neutron Transport Problems. CTH-RF-94., Chalmers Uni. Sweden.. Spiga G. and Vestrucci P. (1981) Ann. Nucl Energy 8, 165. Tezcan C. (1979) Transp. Theory Statist. Phys., 8, 187. Tezean C. (1981) Transp. Theory Statist. Phys., 10, 105. Tezcan C. and Yddlz C. (1986) Ann.Nucl. Energy 13, 345. Williams M.M.R. (1966) The Slowing Down and Thermalization of Neutrons. North - Holland. Amsterdam. Williams M.M.R. (1971) Mathematical Methods in Particle Transport Theory. Butterworths, London. Williams M.M.R. (1978)Math. Proc. Camb. Phil. Sot., 84, 549. Williams M.M.R. (1985) Ann. Nucl Energy 12, 167. Yfldlz C. (1990) IL Nuovo Cimento 105B, 1215. Ydchz C. (1995) Ann. Nucl. Energy 22, 671.
Ann. Nucl. Energy, Vol. 25, No. 8, pp. 529 540, 1998 © 1998 ElsevierScience Ltd. All rights reserved PiI: S0306-4549(97)00114-X Printed in Great Britain 0306-4549/98 $19.00+ 0.00
VARIATION OF THE CRITICAL SLAB THICKNESS WITH THE DEGREE OF STRONGLY ANISOTROPIC SCATTERING IN ONE-SPEED NEUTRON TRANSPORT THEORY C.Yfldiz Department of Physics, lstanbul Technical University, 80626, Maslak, Istanbul-TURKEY
(Received 30 September 1997)
Abstraet--The critical slab problem is studied in one-speed neutron transport theory using a linearly anisotropic kernel which combines forward and backward scattering. It is shown that, the recently observed non-monotonic variation of the thickness also exists in this strongly anisotropic case. In addition, the influence of the linear anisotropy on the critical thickness is analysed in detail. Numerical analysis for the critical thickness are performed using the spherical harmonics method and results are tabulated for selected illustrative cases as a function of different degrees of anisotropic scattering. Finally, some results are discussed and compared with those already obtained by other methods, the agreement is satisactory. The spherical harmonic method gives generally accurate results in one dimensional geometry, and it is very suitable for the numerical solution of the neutron transport equation with linearly anisotropic scattering. © 1998 Elsevier Science Ltd
I. INTRODUCTION For many years the neutron transport problems have been subject to several theoretical and numerical studies. Several workers (see, e.g. Carlvik, 1968; Kaper et al.,1974; Dahl and Sjostrand, 1979;Sahni and SjOstrand, 1990; Garis,1991 ) have studied the criticality problems of one-speed neutrons in homogeneous slabs and spheres using various methods which are applicable for exact analysis. On the other hand, higher order c eigenvalues ( e.g., the mean number of neutrons emitted per collision), have been calculated by Kemer et al (1967), and Kschwendt (1971) for critical slabs, and Dabl and SjOstrand (1979) for both critical slabs and spheres. Gaffs and Sjrstrand (1989) have given an approximate expression for the criticality eigenvalues Cn of order n for the case when cn is large. In all of these works, the results have been obtained with great accuracy in terms of the criticality factor c as a function of the dimensions of the systems, or vice versa. Two of the more popular methods of solving onedimensional critical problems are the Wiener - Hopf techniques (Davison, 1958; Williams, 1971) and the singular eigenfunction method (Case and Zweifel,1967). However in neutron transport theory explicit solution can be found only for very few problems. A number of general methods have been developed for solving the transport equation, and in fact, these methods have been applied to various problems. For some further investigation see Sahni and SjOstrand (1990). It is usually more difficult to obtain very accurate solution with anisotropic scattering. The neutron distribution problem can be solved by inserting into transport equation an appropriate scattering function, which represents the neutron interaction probability. In many of the above mentioned works, the scattering is generally assumed to be isotropic. However, many authors were interested in the effects of anisotropic scattering. (see Sahni
529
530
C. Ylldlz
and SjOsWand,1990). The problem of anisotropy in the interactions of neutrons with matter and of its effects on the critical size is one of the most classic problems of neutron transport theory. It is commonly assumed that the scattering function f(f~'.fl) can be developed into a Legendre polynomial series of the argumentll'.fl, i.e. (Davison, 1958) 1
N
f(n'.n)=~-~(2n~ + l) b~Pdn'.n).
(1)
Here ~ ' and f~ denote the direction of the neutron velocity before and after a scattering collision, respectively, with the normalisation condition bo=l. Further we note that bj is the average cosine of the scattering angle. If the anisotropy is assumed to be linear, only the first two terms are used in equation (1). The first calculation with linearly anisotropic scattering has been studied by Dahl and SjOstrand (1979) who applied Carlvik's method to the integral equation with vacuum boundary condition and have given eigenvalues up to the sixth order. Very little work has been done in the case of quadratic anisotropy, i.e. when only the first three terms are used in the series development of equation (1). This case has been studied by Lathrop and Leonard (1965). The most detailed calculations with linearly and quadratically anisotropic scattering have been described by Sj6strand (1976, 1978, 1980) and more fundamental studies were made by Protopopescu and SjOstrand (1981a, 1981b). As well known, the inclusion of anisotropic scattering effect into neutron transport equation leads to great mathematical complexity. For example, for linearly anisotropic scattering, it is in general difficult to obtain an integral equation for neutron flux alone. In addition, analytical studies of amsotropic scattering effects are usually limited for accurate critical calculations. However a different form of anisotropic scattering can be introduced as the following form
fin'.n)
1-a-p a 4a: ~ 8 ( n ' . n
-1)+
8(n'.n+l)
(2)
where a and fl are parameters which vary from 0 to 1 and give the fractions of neutrons which emerge from a collision in the forward and backward directions, respectively. Such a scattering model (backward-forwardisotropic model) was originally proposed by Fermi to describe neutrons (see e.g. Williams, 1966,1978). This model has been developed to a high degree of sophistication by InOnii (1973, 1976) and has been exploited in a practical manner by Siewert and Williams (1977) for solving critical slab problems, and by Williams (1985) for a multigroup formalism. It was pointed out by lnOnOand Tezcan (1978) and more accurately shown by Pomraning (1978) that increasing the degree of forward anisotropy a , the critical thickness first increases, and reaches a maximum then decreases. This type of anisotropic scattering has also been used to study the critical slab problem by Spiga and Vestrucci (1981). Similar works with the scattering function given in equation (2) were done by many workers (Garcia and Siewert, 1981; Tezcan, 1979, 1981; Tezcan and Yddiz, 1986; Yddiz 1990, 1994; Sahni and SjOstrand, 1991; Sharma et al., 1994; Ganapol and Komreich, 1996). Strongly forward scattering combined with ordinary linearly anisotropic scattering was studied by Sahni and SjOstrand (1991). In a later work, the full problem including also strongly backward scattering was treated in slab geometry by Sahni et a1.(1992). These authors have presented numerical results for the variation of the c-eigenvalues as a function of anisotropic scattering. Finally some exact solutions have been given by SjOstrand (1992). The purpose of the present work is to extend the above investigation of the critical slab problems to include linearly anisotropic scattering together with both forward and backward scattering, and to search the effect of the linearly anisotropic scattering on the critical size. Usually, the results have been expressed in terms of the criticality factor c as a function of the dimensions of the system, or vice versa. In this paper, we fix c as a real quantity and then look for the critical size of the system as a function of the anisotropy parameters. The scattering function is assumed to be of the form f(~'.f~)=~(l+3bj
f ~ ' . £ ' 1 ) + ~ ( £ ' 1 ' .f~ - 1 ) + ~ 8 ( [ 1 ' . [ 1 + 1 ) .
(3)
Here 0 < a ,fl< 1, a + f l < 1. One usually assumes that the scattering function f ( f g . ~ ) is non negative. This means that the coefficient bl in equation (3) must fulfil the condition [bl [ < 1/3. Then, the c-eigenvalues are real for negative values of b/and always real for large systems, about 2mfp thick.. On other hand, many c-eigenvalues
Variation of the critical slab thickness
531
are complex (see e.g. Dalai and SjOstr~d, 1979,1989; Sj6stmnd, 1989). More fundamental studies of the bel~viour of the criticality eigenvalues of one-speed transport operator with linearly anisotropic scattering were carried out by Sabni et al. (1995, 1997). For the criticality problems, generally, analytical techniques based on the integro-differential form of the transport equation are used. These methods yield highly accurate results, but they can only be employed for low order anisotropic scattering. For reactor calculations several terms may be considered in the series development of equation (1). The great complication arises when the higher order anisotropy is considered. In addition, when strongly backward scattering is included, it is neeessa~j to employ a tnmsformation (In6n~ 1973) to convert the equation into canonical form for solving exactly by analytical procedure. As a result of this difficulty, when there is no possibility of using analytical method, one can try more accurate methods for the solution. The complexity of the equations usually forces one to implement numerical methods. Thus such approaches have been made through a direct numerical attack on the integro-differential form of the neutron Iransport equation by using the eigenfunetion techniques. Among these, like the Sn method and the spherical harmonic method (also known PN method) have thus been developed for the design of reactors. For higher-order calculations, the SN methods becomes more popular (Bell and Glasstone, 1972; SjOstrand, 1980; Makai and Sj6strand, 1996). Nevertheless, the PN method is the most direct as well as one of the most powerful appruaehes to numerical transport calculations. One converts the integrodifferential form of the transport equation into a linear system of algebraic equations. These equations are most amenable to solution by a digital computer. Usually only a few orders in the Legendre expansion series given in equation (1) are manageable for exact analysis. For this reason, in the present work, the computations were performed by using the Pr~ method for the critical thickness. In addition, the present method is certainly more accurate in one-dimensional geometries for sufficient large spatial regions. We report some papers on the results of a variety of one-speed calculations in the PN method ( Davison, 1958; Nesthat et al., 1977; Lee and Dias, 1984; Lee et al., 1985; Garcia and Siewert, 1996; Khouaja et al., 1997).
H. THE BASIC EQUATIONS For the investigation of critical slabs, assuming the slab material to be homogeneous, we call 2a the slab thickness and choose x-coordinate axis perpendicular to the slab surface. The origin of the x-axis is chosen such that the slab extends from x=-a to x - a . The system is surrounded by vacuum, and there are no external sources. Then, with conventional notation, the transport equation for neutrons of one-speed can be written as (see e.g. Davison, 1958; Williams, 1971; or Duderstadt and Martin, 1979)
/~-~-x
+ E, Ud(x,/a) = cE, .[ S f(,u',,u)tP(X,la')dla'dqg 1
(4)
o
Here W(x,/a) is the neutron flux density depending on x and on the cosine of the angle between the positive x-axis and the neutron velocity vector. E t is the total cross section and c is the number of secondary neutrons per collision. We assume that E t and c do not vary in space or direction. The multiplicity c is related to the material cross sections by equation cE, = o E : + E,
(5)
where E : is the fission cross section, E, the scattering cross section and v the number of neutrons per fission. It should be noted that in equation (4) we have assumed that the same scattering function is applicable to both fission and the scattering process. The fission source term can be taken to be isotropic in the laboratory system of o o o r ~ s . However, these two problems are not very different. We measure distance in units of total mean free path. The dimensions characterising changes in reactor core composition are usually comparable to a neutron mf-p. If we scale, that is, define a dimensionless space variable such that x ~ xE,, we can write W(xg,,/a) ~ Ud(x,/a). Now, when the scattering function given in equation (3) is inserted into equation (4) we obtain
532
C. Ylldlz ~ I ] ( x .X //~
t2
r+1
subject to the following free space boundary and symmetry conditions W(a, # ) = 0
(7a)
q'(x,/~) = ~e(-x,/~), /~ >0.
(7b)
For an infinite symmetric slab the angular flux can be expanded in terms of the Legendre polynomials 2n+l
• (x, u) = ~ ( - T - ) ~ ( x ) e . (/~)
(8)
where N is odd number andkt the cosine of the angle between~ and the major axis. When the expansion of equation (8) is inserted into equation (6), and using orthogenality properties and recursion relations of the Legendre polynomial, it is found that, after some rearrangement of terms, the PN equations have the following form d~. , d~.+ l n - - - - ~ + (n + 1 ) - - ~ + (2n + 1)dO (x) = (2n + l){c(1 - a - fl)~0a.o + (1 - ct - fl)b,~,8.,
+ a c ~ . +,Bc(-l)"~,},
0 < n _
(9)
One may employ the well-known procedure of seeking a solution of the homogeneous equation (9) in the form (Davison, 1958)
• ,,(x) = G,,(o)exp(-x/u)
(10)
where the functions Gn( o ) are some constants. Each of the ~ . ( x ) defined in equation (10) will satisfy equation (9) provided that the characteristic PNequations
(n+ 1)G,,(v)+nG,,_,(o) - ( 2 n + 1){1-c(1- a - fl)(6.o +b,6,1)- 4 a +(-')'fl]}G.(u)=0
(11)
are satisfied. Equation (11) has a homogeneous matrix form [M(v)IG=0
(12)
where the coefficient m a t r i x is (N+1)x(N+I), and G=[Go,G~,...G,]r. The unknown constant vector G is determined from the normalisation G_~(o) = 0 and Go(u ) = I. The general solutions of equation (6) in the present approximation can be expressed as N+I
q).(x) = £ G . C o , ) [ p j F . ( - x / o , ) + q j F . ( x / u 2 ) ] j=l
(13)
Variation of the critical slab thickness
533
where Fn(x/vj)-=exp(-x/vi)and Pi and qi are arbitrary constants to be determined from the boundm3t eoaditiom of the problem. The essential idea of the P~ method is that l~q(v)=0, i.e., the permissible eigeavalue vj is the j-th positive zero of GN+l(v). Therefore, this should be identical with the detexmin~tal equation (12). Non-trivial solution will be obtained only when deqM(v)l=0.There are L=(N+I)/2 roots. All the roots vj ( for c>l ) are either purely real or purely imaginary. The spatial variation of the angular flux can be evaluated from the moment expansion in terms of Legendre polynomials given in equation (8). The numerical determination of these roots is done using the Newton-Raphson iterative technique. The computation time is extremely short. All calculations reported here were performed in about 2 see. on a pentium machine running at 120 MHz.
m . Tm~ SLAB c R r r I c A L r r Y C O N D m O N
The exact solution of the transport equation (4) must satisfy certain boundary conditions at the surface boundary between two media and the free- surface condition on the external surface. However, it is not possible to satisfy this requirement exactly with PN approximation for finite N. The difficulty arises from the fact that the exact boundary conditions are imposed over half the angular range - 1
fP~ (-/.t)W(a,-p)d/z= 0,
i=1,3,5,...N
(14)
0
The form of the criticality condition can be deduced from equation (14) for a finite slab of half-thickness a. Thus, we obtain the critical equation in the matrix form as follows [ M] (a) ]Pi =[ 0 ]
(15)
where Pi is a vector with elements Pi i= 1,2,...(N+ 1)/2 and the [(N + 1) / 2]: elements of matrix Mtj (a) are given by (N-l)/2
M~(a) = E tp(n,i)G2,(uj)cosh(a / uj) n=0
+ G2,_,(oj)sirda(a / oj) ,
(16)
where (4n + 1X-1)"÷J(2i - l)!(2n)!
q~(n,i) = 22"2"-J(2n- 2i+ l)(i +n)[n!(i- 1)!]2 "
(17)
For the matrix system in equation (15) to have a unique solution, the determinant of M~(a) must vanish, which is the criticality condition. To obtain the one-group neutron flux, an arbitrary normalisation is necessarily required to evaluate the eigenfunction O.(x). After determining the values of G~ the complete angular solution PN can be obtained.
534
C. Ylldiz
IV. N U M E R I C A L R E S U L T S AND DISCUSSION In order to examine the variation o f the critical thickness as a function o f the anisotropy parameters, we have c a m e d out the calculation using vacuum boundary condition from the criticality condition given in equation (15). The permissible root vj, for a given different values o f c is calculated from the determinantal equation(12). Our numerical results are presented in Tables 1-4 for different combinations o f the forward and backward scattering parameters. We give the critical size for various values o f bl ranging from -0.3 to +0.3 for chosen values o f 0=1.2 and 0=2. In some situations, one can also fix the critical size o f the system and then look for the quantity c which is the number o f secondaries per collision. (see e.g.Sahni and SjOstrand, 1990; Sahni et al., 1992). Table 1. The critical thickness 2a for 0=1.2 with degree o f forward scattering a and linearly anisotropic scattering bl as calculated P s approximation.
a
bl
P3
P5
P7
P9
Pll
P13
P15
0.00
-0.3 0.0 0.3
2.37181 2.60400 2.93726
2.35366 2.58453 2.91562
2,35106 2,58177 2,91237
2.35023 2.58076 2.91098
2.34975 2.58014 2.91014
2.34967 2.58002 2.90995
2.34973 2.58009 2.91004
0.20
--0.3 0.0 0.3
2.54790 2.78353 3.11661
2.51866 2.75196 3.08163
2,51363 2.74684 3,07611
2.51227 274530 3.07415
2.51157 2.74444 3.07299
2.51150 2.74431 3.07275
2.51161 2.74442 3.07287
0.40
--0.3 0.0 0.3
2.77465 3.00915 3.33257
2.72186 2.95105 3.26716
2.70928 2.93828 3,25398
2.70590 293484 3.25021
2.70455 2.93333 3,24836
2.70457 2.93328 3.24815
2.70483 2.93354 3.24841
0.65
--0.3 0.0 0.3
3.17711 3.39006 3.66598
3.03995 3.23431 3.48433
2.97765 3A6655 3.40948
2.95103 3.13867 3.38004
2.94021 3.12761 3.36869
2.93950 3.12720 3.36863
2.94052 3.12854 3.37035
0.70
--0.3 0.0 0.3
3.28121 3.48333 3.73973
3.10924 3.28762 3.51124
3.01728 3.18634 3.39769
2.97131 3.13710 3.34434
2.95115 3.11595 3.32199
2.94897 3.11405 3.32051
2.94978 3.11533 3.32244
0.8333
--0.3 0.0 0.3
3.60274 3.74994 3.92139
3.26866 3.37765 3.50097
2.97875 3~06067 3.15135
2.73464 2.79771 2.86635
2.60780 2.66236 2.72127
2.60780 2.66235 2.72126
2.60780 2.66235 2.72126
In Table 5 the results o f the critical thickness are given for a+fl ranging from 0 towards the limit a + f f =l/c as a function o f the different combination o f various anisotropy parameters. In order to see the performance and the accuracy o f the present computational method, the critical thickness for some selected value o f c are given in Table 6. We therefore take as reference values, the values given in the papers by Sahni et al. (1992) and Dahl and Sj6strand (1989). For the scattering function under investigation here the numerical and theoretical criticality calculation have been performed by many workers (see Sahni and Sj6strand, 1990) using various methods. R can be seen that our results for critical thickness versus a,fl and bt, follow the same general trend as those o f In6nO and Tezcan (1981), Pommning (1978), Sahni and Sj6strand (1991) and Sahni et al. (1992). The agreement was very good. It has been confirmed that the earlier observed non-monotonic variation o f the
Variation of the critical slab thickness
535
critical thickness exists in this case as well. The smprising trend previously reported is that the critical thickness for a fixed c first increases with forward anisotmpy parameters a(fl= 0) and then decreases. However, the variation of the critical thickness with a and linear anisotropy parameters bl depend upon how the value c has been chosen (Sahni and SjOstrand, 1991). We see from Table 1-4, as mentioned above, the variation will be non monotonic only
Table 2. The critical thickness 2a for c=2 with degree of forward scattering a and linearly anisotropic scattering b~ as calculated Ps approximation. t/
bl
P3
Ps
P7
P9
PH
P13
PLs
0.00
-0.3 0.0 0.3
0.64215 0.68410 0.73809
0.61261 0,65052 0,69890
0.59847 0.63506 0.68172
0.59211 0.62835 0.67457
0.58946 0.62562 0.67175
0.58925 0.62547 0.67168
0.58945 0.62575 0.67205
0.10
-- 0.3 0.0 0.3
0.65624 0.69667 0.74795
0,62185 0.65752 0.70225
0.60346 0.63727 0.67954
0.59426 0.62742 0.66887
0.59023 0.62319 0.66440
0.58979 0.62281 0.66410
0.58996 0.62307 0.66449
0.20
-- 0.3 0.0 0.3
0.67119 0.70971 0.75777
0.63091 0.66391 0.70443
0.60675 0.63716 0.67429
0.59316 0.62244 0.65814
0.58689 0.61576 0.65093
0.58615 0.61503 0.65024
0.58622 0.61519 0.65054
0.30
-- 0.3 0.0 0.3
0.68699 0.72313 076735
0.63955 0.66937 0.70510
0.60746 0.63380 0.66505
0.58695 0.61145 0.64039
0.57704 0.60076 0.62873
0.57599 0.59967 0.62761
0.57594 0.59969 0.62772
0.49
-- 0.3 0.0 0.3
0.71885 0.74870 0.78356
0.65319 0.67543 0.70068
0.59690 0.61382 0.63261
0.55012 0.56335 0.57781
0.52586 0.53742 0.54996
0.52571 0.53726 0.54979
0.52571 0.53726 0.54979
0.4999
-- 0.3 0.0 0.3
0.72055 0.74999 0.78428
065373 0.67553 0.70019
0.59575 0.61213 0.63027
0.54693 0.55954 0.57327
0.52156 0.53247 0.54425
0.52156 0.53247 054425
0.52156 0.53247 0.54425
if the critical thickness is larger than 0.67 for isotropic scattering ( a , fl = 0). This condition is equivalent to the requirement about c<2. For a certain degree of linearly anisotropic scattering hi, the critical size 2a is a function of two quantities, viz. the quantity c and the degree of forward scattering a . This is seen very clearly from the results in Table 1 and 2. As an example, for c=l.2 the critical thickness first increases with increasing forward scattering a for all values of bl ranging from -0.3 to +0.3 and reaches a maximum at about a = 0.71 and the thickness 2a=2.95. In that case, that is, if one wants to have a maximum at the value bf-0.3 and a =0.65 it is necessary to change the range of hi (from b1=-0.3 to b~=+0.3, in this case ) or one must reduce the quantity c. From the above discussion we see that for a given value of c=2, in Table 2, the maximum of thickness 2a (2a=_=_0.58996) turns out to be a =0.1 for b1=-0.3. While for bT0.0 and b1=+0.3, the critical thickness reduces monotonically with increasing a . That is, the behaviour of the critical thickness depending on bl is still non monotonic for c=2 and b~=-0.3. From this result, the variation of the critical thickness is monotonic for all values ofbt only when the quantity c is larger than about 1.9 or 2. In the case of backward scattering f l ; e 0 ( a = 0 ) the critical sizes for c=l.2 and c=2. decrease monotonically with increasing aniso~ropy. It was pointed out by Sahni et al. (1992), if the slab is thicker the presence of fl ;e 0 shows that the behaviour can be non-monotonic for some scattering. Similar conclusion are also confirmed from the numerical results of In6nO and Tezcan (1978) and Salmi and Sj6sUand (1991) and also Sahni ct al.(1992). In6m3 and Tezean (1978) obtained the critical slab thickness for c=1.2 using a variational method. The size parameter 2q corresponding to the maximum 2a works out be 3.11 with about a = 0 . 6 5 . We tried to arrive at these values from the curves given in their papers. Sahni and Sj6strand (1991) obtained a useful criterion for the monotonic or non-monotonic variation of c with a at fixed values of system dimensions. They found that the variation of the criticality factor is non-monotonic for all values of bl only when the thickness is larger than about
536
C. Ylldlz
0.7 mfp. This condition is equivalent to the requirement c<1.9. Consequently, we conclude that their results are in perfect agreement with our values in all considered cases. Another interesting ease is that the critical thickness seems to go to zero in the limit o f a =l/c. Note that when a becomes equal to l/c, the neutron dislribution is entirely concentrated in the normal plan x=0 (InrnO and Tezean, 1978). As can be seen from Tables 3 and 4, for c=2, when a or ,6 approaches the value l/c, the critical slab thickness is attained in the computational accuracy to the value o f 2a---0.5 and 0.3 for forward and backward scattering ease respectively, independently from the value o f b j . In the limit o f a,fl=l/c the critical thickness falls sharply to zero (In0na and Tezcan, 1978; Sahni et al., 1992). It is not possible to arrive at this sharp value in this limit using the present method. This limit value is computed with a precision o f nine significant figures for tt and ft. For example given a value o f fl = 0.4999 .... we obtain 2a---0.328. The other surprising result is the behaviour o f the critical thickness with anisotropy parameters a + f l . We see from Table 5 that the variation o f the critical thickness with increasing a + fl may sometimes be nonmonotonic. There is maximum at about a = fl=0.01 for c=1.2 and bi = -0.3 while for bt = -0.3 and 0.0, the variation is still monotonic. This depends on the values o f a ,fl and b I and also on the value c. However, for such systems below a certain size the variation is monotonic. It can be seen that the change occurs at about 2a-~ 1. One might think that there is certain value o f the thickness 2a for which the variation does not depend on the anisotropy parameters. For the thick slab, the variation o f the thickness for c<1.5 is still non-monotonic while for about c>1.6 the variation is monotonic. Such a behaviour o f the critical thickness is practically independent o f a + fl and bt-
Table 3. The critical slab thickness 2a for c=1.2 with various degree o f backward scattering fl and linearly anisotropic scattering bi as calculated Pr~ approximation. fl
bl
P3
P5
P7
P9
Pll
P13
P15
0.10
- 0.3 0.0 0.3
230440 249189 2.74378
2.28399 2.47034 2.72060
2.28106 2.46732 2.71730
2.28019 2.46631 2.71603
2.27970 2.46571 2.71526
2.27962 2.46560 2.71509
2.27969 2.46568 2.71518
0.20
-0.3 0.0 0.3
2.23691 2.38662 2.57695
2.21316 2.36185 2.55085
2.20968 2.35834 2.54722
2.20873 2.35729 2.54599
2.20822 2.35670 2.54527
2.20815 2.35660 2.54513
2.20822 2.35668 2.54522
0.40
-0.3 0.0 0.3
2.09371 2.18334 2.28767
205670 2.14521 2.24829
2.05042 2.13902 2.24218
2.04897 2.13755 2.24067
2.04834 2.13687 2.23993
204830 2.13682 2.23985
2.04841 2.13693 2.23996
0.70
-0.3 0.0 0.3
1.79546 1.82199 1.84995
1.66590 1.68907 1.71346
1.61681 1.63937 166315
1.59857 1.62114 1.64492
1.59171 1.61435 1.63820
1.59120 1.61387 1.63777
1.59177 1.61447 1.63840
0.80
-0.3 0.0 0.3
1.64053 1.65289 1.66556
1.41613 1.42436 1.43275
1.28466 1.29110 1.29765
1.20064 1.20616 1.21178
1.14702 1.15204 1.15715
113489 1.13983 1.14485
1.13110 1.13603 1.14105
0.8333
-0.3 0.0 0.3
1.58016 1.58891 1.59781
1.31841 1.32340 1.32845
1.15208 1.15538 1.15871
1.03458 1.03696 1.03936
0 95199 0.95384 0.95569
0.93253 0.93427 0.93601
0.92482 0.92652 0.92823
This was already indicated in earlier paper by Sahni et al. (1992) .They have given the value o f c<1.9 as the condition for the appearance o f the non-monotonic behaviour o f the criticality factor c. One example is shown in the paper o f Salmi." et al. (1992) that the variation o f c with a = ,8 is monotonic for a slab o f thickness 5 mfp, for various value o f bl ranging from 0 to +0.3. That is, the c value decreases with increasing ct+fl. But for b~ = 43.3 there is maximum in the curve at about at +fl=0.47. This depends on the values o f the a , p and bl and also on the
Variation of the critical slab thickness
537
thickness of the slab. For I mfp the variation of the criticality factor is always monotonic with increasing a + fl and for all values of bj. This should be eoml~tred to the value c=1.6 and the thickness about 2a<. 1.1 from our calculations using the P~ method. Thus the agreement between our values and with those of In6nO and Tezean (1978) and Sahni et al. (1992)is very good. Our results follow the same general trend as those of them
Table 4. The critical thickness 2a for c=2 with various degree of backward scattering fl and linearly anisotropic scattering bl as calculated Pr~ approximation.
P7
P9
PII
O10
--0.3 0.0 0.3
0.61856 0.65127 0.69127
0.58338 0.61202 0.64671
0.56559 0.59280 062568
0.55706 0.58382 061617
0.55332 0.57996 0.61218
0.55294 0.57963 0.61193
0.55314 0.57992 0.61231
0.20
--0.3 0.0 0.3
0.59453 0.61951 0.64874
0.55238 0.57327 0.59744
0.52939 0.54866 0.57087
0.51727 0.53592 0.55743
0.51151 0.52996 0.55123
0.51078 0.52926 0.55059
0.51094 0.52949 0.55091
0.30
--0.3 0.0 0.3
0.56974 0.58827 0.60914
0.51901 0.53350 0.54960
0.48875 0.50140 0.51539
0.47073 0.48253 0.49556
0.46115 0.47259 0.48520
0.45970 0.47112 0.48373
0.45964 0.47112 0.48380
0.40
--0.3 0.0 0.3
0.54391 0.55712 0.57153
0.48299 0.49237 0.50243
0.44319 0.45064 0.45858
0.41625 0.42264 0.42941
0.40000 0.40583 0.41199
0.39726 0.40301 0.40910
0.39671 0.40248 0.40858
0.49
--0.3 0.0 0.3
0.51965 0.52898 0.53889
0.44857 0.45444 0.46059
0.39904 0.40315 0.40740
0.36239 0.36547 0.36863
0.33843 0.34093 0.34349
0.33440 0.33681 0.33928
0.33329 0.33569 0.33815
0.4999
-- 0.3 0.0 0.3
0.51692 0.52587 0.53536
0.44470 0.45025 0.45605
0.39410 0.39792 0.40187
0.35641 0.35922 0.36210
0.33164 0.33389 0.33618
0.32750 0.32966 0.33187
0.32635 0.32849 0.33069
fl
bl
P3
P5
Pl3
PI5
We note that also from results in Table 5, for a given values of aandfl, the maximum points shift towards increasing a + fl for the decreasing value of c. Such behaviour of the maximum point may also be explained by the change in the magnitude of the value of c. The decrease of c increases the size, then leakage is favoured by the anisotropy to arrive at the critical values. For c=1.2, when the values a + f l approach the limit value of about 2a = 1.2; 1.7 and 2.0 for fl = a ; 0.15; 0.30 respectively. This happens independently of the value of bl. Finally, we see from Table 6 that the present method gives very accurate results and is well suited for the analysis of such systems. Numerical results of Dalai and SjOstrand (1979) and Sahni et al. (1992) practically coincide with our results to 4 or 5 significant figures especially for the thickness larger than 1.
V. CONCLUSION In this paper we have considered the slab criticality problem of the linear transport equation with forward, backward and linearly amsotropic scattering in a homogeneous slab. It is pointed out by many workers (see e.g. Sabni and SjOstrand,1990 ) that for fixed value of the factor c, the slab thickness varies non-monotonically with increasing the forward scattering under certain conditions. We observe that the addition of the linearly amsotropic scattering increases the critical slab thickness. That is, the leakage is less favoured by isotropic scattering than by the linear anisotropy. One characteristic feature of the linear anisotropie scattering is that the critical dimension of a
538
C. Ylldlz
system for a fixed c- value increases with the degree of anisotropy (see e.g. Sahni and SjOstrand,1990 ). In these thin slabs, the angular distribution of neutrons is very strongly peaked in the direction along the slab. However, for such systems below a certain size, especially for the thinner slab, (about 2a--0.67 in this work ) the variation of the critical thickness is monotonic with increasing a and for all values of b~. That is criticality factor c becomes only the controlling factor for the size when the thickness is larger than about c_=_1.9. In the ease of backward scattering (i.e. a = 0), the thickness decreases uniformly with increasing fl independently the value of bl. A special case of the scattering function in equation (1) has been studied in detail by Sahni e~ al. (1992). This surprising variation of the critical slab thickness may be understood by considering the change in the angular distribution (InOnO and Tezcan, 1978; and Pomraning, 1978; Sabni et al. 1992 ).
Table 5. The critical slab thickness 2a as a function of some selected anisotropy parameters a + fl and linearly anisotropie scattering bt for three value ofc as calculated Pl3 approximation. c=1.2 c~
c=1.5
c=1.6
bl fl=a
fl=0.15
fl=0.30
fl=a
fl=0.15 fl=0.30
fl=a
fl=0.15 fl 0.30 =
0.00
- 0.3 0.0 0.3
2.34967 2.58002 2.90995
2.24421 2.41067 2.62769
2.13252 2.24885 2.38984
1.12371 1.21194 1.33109
1.05134 1.11212 1.18771
0.96393 1.00265 1.04771
0.95516 1.02598 1.12028
0.88728 0.93511 0.99388
0.80304 0.83232 0.86599
0.01
- 0.3 0.0 0.3
2.34996 2.57565 2.89640
2.24984 2.41583 2.63195
2.13594 2.25137 2.39104
1.12113 1.20705 1.32225
1.05192 1.11220 1.18704
0.96218 1.00023 1.04440
0.95226 1.02109 111206
0.88723 0.93458 0.99265
0.80065 0.82930 0.86217
0.05
- 0.3 0.0 0.3
2.34938 2.55616 2.84102
2.27251 2.43642 2.64868
2.14887 2.26044 2.39450
1.10847 1.18500 1.28463
1.05356 1.11170 1.18337
0.95339 0.98856 1.02903
093825 0.99894 1.07683
0.88628 0.93157 0.98667
0.78921 0.81519 0.84471
0.33
- 0.3 0.0 0.3
2.05201 2.10806 2.17005
2.41891 2.55539 2.72314
2.13176 2.19949 227558
0.66691 0.67168 0.67658
0.99992 1.03356 1.07199
0,73009 0.73770 0.74562
0.42
- 0.3 0.0 0.3
1.21760 1.22168 1.22580
2.43809 2.55800 2.70183
1.98612 2.02982 2.07734
0.68
- 0.3 0.0 0.3
0.92990 0.95144 0.97506
0.81156 0.83417 0.85945 0.74085 0.75387 0.76779
1.71137 1.72439 1.73776
Another important conclusion of this calculation is that the PN method provides accurate results with very short computing time. All calculations reported here were performed in about 3 rain. on a pontium machine runnin~ at 120 MI-Iz. Inspection o f these numerical results in tables shows that the results are correct at least in the first 3 and 4 digits. The calculated critical thicknesses are rounded off at the 5th decimal point. The accuracy depends on the order o f PN approximation. The convergence of the method is very fast. Four stable decimal figures are always achieved with N _<11,13. The P15 approximation is relatively accurate than the P3 approximation. But there are no great improvements o f the P15 approximation over the Pi3. Therefore some numerical results are given in the P13 approximation. It is not necessary to go to further approximation. In addition, the accuracy oll the critical thickness depends also on the accuracy on the fixed values o f c. We have compared some o f our results with those obtained by means o f various methods. These are in very good agreement with those o f Dahl and Sj6strand (1979) and Sahin et al. (1992) up to the 5 or 6 digit, as shown in Table 6.
V a r i a t i o n o f the critical slab t h i c k n e s s
539
Consequently, we conclude that the present method is suited and also rapid for this type o f problems with less computational effort. For large systems the lower order approximations yield somewhat better results. In this ease the method works reasonably well. It is applicable to large system, and also to very small assemblies. In addition it can also be extended to more general problems with higher order anisotropy and to multiregion energy dependent problems.
Table 6. The some critical slab thickness obtained in various PN approximation with different value o f c, a , fl and degree o f linearly anisotropic scattering bl comparison with results o f Sahni et al. (1992)andDahlandSj0sWand(1979).
a,fl
b]
-0.30 ~Z= 0.3 0.00 fl= 0.0 0.30
-0.30 IOr= 0.0 0.00 i fl=~ 0.3 0.30
P3
P9
PIl
Pt3
8ahni
Dahl
1.5857300 1.0808500 1.6183300 1.0936500 1.6557200 1.1116600
1.10876 5.02327 1.12337 5.02910 1.14309 5.03908
1.01337 5.00225 1.01532 5.00291 1.01798 5.00417
1.00775 5.00160 1.00895 5.00204 1.01062 5.00294
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.08251 5.01155 1.08968 5.01365 1.09834 5.01595
1.00521 5.00106 1.00575 5.00157 1.00642 5.00170
1.00274 5.00069 1.00305 5.00113 1.00346 5.00117
1.00270 5.00061 1.00303 5.00104 1.00345 5.00106
1.00 5.00 1.00 5.00 100 5.00
-----
-----
--
--
0,00
3.8303067 0.25661 0.21472 0.21050 0.21010 1.1631293 3.02170 3.00193 3.00133 3.00121 1.1084678 4.01675 4.00182 4.00126 4.00113 1.0071358 20.01077 20.00152 20.00104 20.00091
-----
0.20 3.00 4.00 20.00
0.30
1.3162908 1.1921797 1.1309166 1.6726350
-----
2.00 3.00 4.00 1.00
5 = 0.0
/~=00
c
2.04010 3.02825 4.02268 1.06249
2.00297 3.00271 4.00257 1.00422
200206 3.00187 4.00178 1.00260
2.00190 3.00168 4.00159 1.00266
REFERENCES Bell G.I. and Glasstone S. (1972) Nuclear Reactor Theory. Van Nostrand-Reinhold,New York, NY. Carlvik I. (1968) Nucl. Sci. Engng., 31,295. Case K. M. and Zweifel P.F. (1967) Linear Transport Theory. Addison-Wesley, Reading, Mass. Dahl E.B. and SjOstrand N.G. (1979) Nuel. Sci. Engng., 69, 114. Dahl E.B. and Sj6strand N.G. (1989) Ann. Nucl.Energy 16, 527. Davison B. (1958) Neutron Transport Theory. Oxford University Press, London. Duderstadt J.J. and Martin W.R. (1979) Transport Theory. Wiley, New York. Ganapol B.D. and Komreieh D.E. (1996) Ann. Nucl. Energy 23,301. Gareia R . D . M and Siewert C.E. (1981) J. Phys. D: Appl. Phys.,14, L65. Gareia R D . M . and Siewert C.E. (1996) Ann. Nuel. Energy 23, 321. Garis N.S. (1991)Nuel. Sei. Engng., 107, 343. Gaffs N.S. and SjSslrand N,G. (1989) Kernteehnik 54, 183. InOnO E. (1973)Transp. Theory Statist. Phys., 3, 107. In6n• E. (1976) Phys. Fluids 19, 1332. In~nt~ E. and Tezean C. (1978) Transp. Theory Statist. Phys., 7, 51.
540
C. Yddlz
Kaper H.G., Lindeman A.J. and LeafG. K. (1974) Nuel. Sci. Engng., 54, 94. Kerner I.O., Kiesewetter H. and von Weber S. (1967) Kernenergie 60, 299. Khouaja H., Edwards D. R. and Tsoulfaridis N. (1997) Ann. blue. Energy 24, 518. Kschwendt H. (1971) Nucl. Sci. Engng., 44,423. Lathrop K.D. and Leonard A. (1965) Nucl. Sci. Engng., 22,115. Lee C. E. and Dias M. P. (1984) Ann. Nucl. Energy 11,515. Lee C. E., Pan W.C.P. and Dias M. P. (1985) Ann.Nucl. Energy 12, 613. Makai M. and SjOstrandN.G. (1996) Ann. Nucl. Energy 23, 397. Neshat K., Siewert C. E. and Ishiguro Y. (1977) Nu¢l. Sci. Engng., 62, 330. Pomraning G.C. (1978) Transp. Theory Statist. Phys.,7, 161. Protopopescu V. and Sj6strand N.G. (1981a) Progr. Nucl. Energy 7, 47. Protopopeseu V. and SjOstrand N.G. (1981b) On the solution of the Dispersion Equation for Monoenergetic Neutron Transport with Quadratically Anisotropie Scattering. CTH-RF-36, Chalmers Uni., Sweden. Sahni D.C. and Sj6strand N.G. (1990) Progr. Nuci. Energy 23, 241. Sahni D.C. and Sj6strand N.G. (1991) Transp. Theory Statist. Phys., 20, 339. Sahni D.C., Dahl B. and SjOstrandN.G. (1997) Ann. Nucl. Energy 24,135. Sahni D.C., Garis N.S. and Sj6strand N.G. (1995) Transp. Theory Statist. Phys., 24, 629. Sahni D.C., Sj6strand N.G. and Gaffs N.S. (1992) J. Phys. D: Appl. Phys., 25, 1381. Sharma A., Paranjape S.D., Kumar A. and Dwivedi S.R. (1994) Ann. Nucl. Energy 21,267. Siewert C.E and Williams M.M.R (1977)J. Phys. D: Appl. Phys., 10, 2031. Sj6strand N.G. (1976)J. Nucl.Sci.Tech., 13,81. SjOstrand NG. (1978) Atomkemenergie 31, 16. Sj6strand N.G. (1980) Ann. Nucl. Energy 7, 435. SjOstrand N.G. (1989) Approximate High-Order Eigenvalues in Two-Medium One-Speed Neutron Transport. CTH-RF-67, Chalmers Uni., Sweden. SjOstrand N.G. (1992) Exact Solution of Some One-Speed Neutron Transport Problems. CTH-RF-94., Chalmers Uni. Sweden.. Spiga G. and Vestrucci P. (1981) Ann. Nucl Energy 8, 165. Tezcan C. (1979) Transp. Theory Statist. Phys., 8, 187. Tezean C. (1981) Transp. Theory Statist. Phys., 10, 105. Tezcan C. and Yddlz C. (1986) Ann.Nucl. Energy 13, 345. Williams M.M.R. (1966) The Slowing Down and Thermalization of Neutrons. North - Holland. Amsterdam. Williams M.M.R. (1971) Mathematical Methods in Particle Transport Theory. Butterworths, London. Williams M.M.R. (1978)Math. Proc. Camb. Phil. Sot., 84, 549. Williams M.M.R. (1985) Ann. Nucl Energy 12, 167. Yfldlz C. (1990) IL Nuovo Cimento 105B, 1215. Ydchz C. (1995) Ann. Nucl. Energy 22, 671.