- Email: [email protected]
The collision probability code—minos
Journal of Nuclear Eneqy, 1967, Vol. 21. pp. 767 ‘o 785. Pqamon
Prcsa Ltd. Printed in Northem Ma&
THE COLLISION PROBABILITY CODE-MINOS * R. G. HARPW Rolls-Royce Limited, Derby (First received 1 August 1966 and injkalform
14 April 1967)
Abstract-A method is developed to solve the energy-dependenttransport equation with anisotropic scattering, by the method of collision probabilities, in carte&n pmetry. The program MPM)G has been written to produce the resulting collision probability matices in a typical lattice cell. These are then inverted using other programs. The running time of the scheme is comparable with multigroup two-dimensional diffusion theory. Comparisons of the method with other theoretical techniques are reported. 1. INTRODUCTION THE most frequently used technique for the prediction of the neutron flux distribution in reactor systems is the method of diffusion theory. Frequently this method is used in situations where the conditions for the applicability of the method are violated, e.g. close to material discontinuities and throughout heavily absorbing regions. Recourse to correlated data allows predictions to be made of specific parameters at the expense of generality in application. Alternatively, solutions Qf the transport equation in the integro-differential form lead to very timeconsuming programs, especially if more than one-dimension is considered. However, another approximation, applied to the integral form of the transport equation, is receiving considerable attention recently-this is known as the method of collision probabilities. In this approach the total reaction rate in any region of the system is expressed in terms of the sources, either direct or scattering, in all other regions of the system multiplied by the probability that the next collision will be made in the region under Now, in general, the collision probabilities depend upon the consideration. flux shape in the source region, but under certain circumstances it has been demonstrated (see for example STUART and WOODRUFF, 1957 or WEXLER,1963) that the collision probabilities are insensitive to the actual flux shape used. STACE(1962) shows that the criterion for the applicability of this assumption varies with the ratio of the scattering to total cross sections in the region being considered. His conclusions are discussed below. If, for all regions considered, the above criteria are satisfied, then we can choose a flux shape which is spatially independent within each region. Consequently the collision probabilities can be determined to sufficient accuracy ab initio, without recourse to any iterative approach. This step leads to significant saving of computer time when compared with finite difference methods. Finally, if some alternative technique is available for the treatment of large extents of moderator-where the use of collision probability techniques would lead to an unrealistically large number of regions-then the method provides a rapid and accurate prediction of the neutron spectrum, upon which considerable reliance can be placed from a design point of view. The collision probability approach is developed below for the case of anisotropic 767
.
768
R G. HARPER
neutron scattering-this is an extension of the more usual treatment which assumes isotropic sources and fluxes (see however the work of HYSLOP(1963) which considers anisotropic fluxes with isotropic neutron scattering and the paper by TAEMIASHI (1966) considering anisotropic flux with anisotropic scattering in a cylindrical system). Also CARL~IK(1966) has suggested that the extension of collision probability methods to take account of anisotropic neutron scattering events can be effected using Gaussian quadrature methods. Results obtained using the code are compared with a onedimensional anisotropic transport code TET-Pl and also some comparisons with a two-dimensional isotropic transport code 2D X-Y are reported. 2. DERIVATION
OF THE COLLISION
PROBABILITY
EQUATIONS
The collision probability approach consists essentially of describing neutron transport in terms of the probability that neutrons born, or scattered, in one region of the cell should have their next collision in some other region. The theory of the method as usually developed makes, as fundamental, the assumption that all emission sources are isotropic and that neutron scattering is isotropic in the Laboratory System. Any treatment of anisotropic scattering will result in anisotropic fluxes and consequently anisotropic sources and will therefore require that explicit information relating to the angular dependence of the neutron flux be available. These two requirements are not compatible. However using plausible physical assumptions, it is possible to obtain solutions of the equations within the framework of collision probability techniques and also to take account to some extent of the angular dependence of the neutron flux and neutron scattering events. 2.1 The collisionprobabilityform of the integral Boltzmannequation withanisotropic scattering Consider the point S at position r relative to some origin 0. Then if a direction P is chosen through this point S, the position of any point P on the line in the negative direction of S2can be described as (r - R&2),where R is the distance from S to P and Sz is a unit vector. The number of collisions in dR at r, for neutrons moving in direction S2along the line, is given by: E)+(r, E, Sz) dE dR = [ dR’PI-k!k
&(r,
JR
‘S(r - R’S& E, a) dE + *, d&2’&, cl(r - R’S2, E’S2’--+ En)+@ - R’QE’, Sz’) dE’ dE s I
X
$1 _
(1)
dS&’ F(r - R’Q, JI) sn s E’ x y(r - R’B, E’) c,(r - R’B, E’)#(r - R’P, E’, a’) dE’dE _
Here the element of solid angle is taken to be unity. Where the term in brackets represents the source of neutrons moving in direction Szat the point (r - RQ), and where: P,-%+
= &(r, E) dR exp (--LB’&(S)
dS)
769
The collision probability code-rmos
represents the probability that neutrons born at r - R’S2 and moving in direction n should suffer their next collisions in dR at r. Ais an eigenvalue, the inverse of which is the criticality factor of Ken. S(r, E, $2) is a term containing all sources external to the system. Consider the above equation applied to a multi-region system of N regions, with constant cross sections in each region, then we may write for the collision rate at r in direction S2:
&(r, E)&, E, Q) dR dE =
s
5 ,dR’P,_&2 ~
i=l
Rt
‘S(r - R’O, E, a) X
+ n dL2’1, x,,(E’Q’ + ES&j0 - R’Q, E’, Q’) dE’
s
+Iz
dE
s
sn
d&Y F(E)v,(E’) z,,(E’)#r - R’S& E’, a’) dE’ s 1’
Multiplying both sides by dA, we may replace the integrals over the line elements, dR, by volume integrals, dr. Performing the volume integrals over dr’, i.e. along the line through r in direction SZ,and also the integrals over the energy group intervals, we obtain: &,(+A&
a) dr = jr K(r, n)%,(r, WL,,(~)
dr
+ 47A 5 &,K(r, a) 2 vi,. &,k
Q)+&r, Q2,WP”(Q) dr i-w,0
d-1
i-1
{g =
1,2, . . . , G). (2)
Where Vi(r, SQ = dA . Ri
(34
huh Q, W =
V,(r,WUr, W (a) dr = i-+r.n
W
dES(r’, E, rR)P (E) r_R’Ll -)r
.
(W
770
R. G. HARPER
Equation (2) above has expressed the collision rate in dr at point r for neutrons in group g moving in direction P, in terms of the angularly dependent flux and sources along a line passing through r direction 51. In the next section the group-averaged parameters above will be obtained, subject to two simplifying assumptions. These are: (i) That both the flux and source can be considered to be spatially constant in any region. (ii) That the collision probability sections.
can be obtained using group-averaged
cross
2.2 Production of group-averaged parameters In order to include explicit information on the angular dependence of the scattering events, the flux and scattering cross section will be expanded in Legendre polynomials. The coefficients of the polynomials in the flux expansion will be the required solution. Consider the expansion of the differential scattering cross section into surface harmonics : a@ --t E’, e) = $ $: (2z+ l)uz((E + E’)P,(cos e)
(4)
where 19is the angle through which the neutron is scattered. Let Q be the unit vector in the direction of the incident neutron Q’ be the unit vector in the direction of the emergent neutron. Then using the addition theorem for Legendre polynomials we obtain: P,(COSe) = ~,(s2. s2’) = P,(COSe)P,(cos en) + 2 i (I - m)! m=~(2 + m)!
x
P,~os
B)P,~cos
e’)
cos
[m(p, - y’>]. (5)
Here axes have been chosen such that cos 0 is the projection of &Zonto the Z-axis and 8 is the angle between a fixed line in the plane 2 = 0 and the projection of Sz onto this plane. Consider also the expansion of the flux into surface harmonics: a,M, #(r,E,Q
E)p,(cos 0) +$1[a,“x,~9a(r, W,Q(cos 0) cos (49)
= Ta
(6) [
+ b,QXPa(r,JW,a(cos 0) sin (c~$lI
The co-ordinate system is the same as in (5) above. Here we have:
a,x,(r, E) = 5,(r, E) =
#r,EW,QdSk
(7)
Given equations (5) and (6), together with the subsidiary relations (7), we can now construct the group-averaged parameters in the equations (3a)-(3e) above.
The collisionprobability code-mms Consider for example equation (3~):
dE &(E’sZ’ --t ELl)+(r,l, E’, Sk’)
co@’ - JX,(r,, E%(Q) + cl@’-+ E)&(r,,E’)h(Ql [ + 1~3’- E)Pl’ (Q)[Ell(ri, E’) ~0s v + rl?(r,,E’) sin~11 where only first-degree terms have been retained. Finally, performing the integration over a;ri, gives: V,(r, Q)Ur,
Q, a’) CS,,*,, (r, Q, S2’) = V,(r,S2)/zfldE’c+1dE
&(E’ -+ EM,
1
fi2, E’)P&os 0) + &(E’ -+ E)E,,(r, a, E’)P,‘(cos 0)
[ + c14E’ -+ E>P~(cos @[&t(r, Q, E’) cos q + qlt(r, S2, E’) sin ~1
where the l&r, P, E’) etc. define mean components of the flux along the line through r in direction P given by: l&r, S2E) 3 s VikW
& hh
E’)/Ur, a).
Let us now make the two assumptions below: (i) that the collision probability in an energy group g, say, can be obtained using effective cross sections in that group, and
(ii) that the scalar and angular components of the neutron flux and source are constant over each region.
The first assumption follows from using monoenergetic collision probabifity theory within each neutron energy group. The second assumption has been extensively examined by several authors-STUART and WOODRUFF,WEXLER,and STAR in the
772
R a. HARPER
isotropic case. STACE,in particular, concludes, that independent of the ratio of scattering to total cross section, the flat fhrx assumption should be satisfactory for regions where the mean chord is less than about two mean free paths. He also derives a second condition, necessary for the solution to be independent of the flux shape in the region, namely: (&,P&,[$,(r)l
Q1
where P&,(r)] is the self-collision probability in regionj determined from the actual flux shape in the region. Hence, in moderating materials, this is a much more stringent criterion and is true only if P,,< 1, i.e. if the regions are very small compared with a mean free path. Making use of the above assumptions, and integrating equations (2) over dr, for a flxed line direction P we obtain:
I ‘=1,2,...,N
g = 1,2,. ..,G Here we have expressed the mean angular flux over a region j in group g in terms of the mean angular flux in all other regions and groups. The analysis so far has been general, in that no assumption has been made concerning the form of the scattering cross section or the weighting spectrum to be used for producing group-averaged parameters. In a typical reactor physics problem, it is well known that the flux in a homogeneous material region large in terms of mean free paths is characteristic of the material of the region and is, therefore, separable in space and energy. Consequently a method of performing consistent group averaging of equations suggests itself. Let us suppose that the characteristic scalar spectrum of the material of region i has been determined, and let the spectrum be represented by #CO(E). Then we can average the zero-moment cross section &&!?’+ E) in equation (8) relative to this scalar flux, i.e. the zero-moment of the expansion of the flux &,@) has the same energy dependence as this spectrum. To perform the averages of the first moment of the cross section, we note the relation between the zero and Crst moments of the neutron flux, given by WB~VBERO and WIGNER(1958)p. 233 for the PIspherical harmonics approximation in the case of anisotropic scattering:
The collisionprobability code-rams
773
Bearing in mind the fact that &Z, B) is separable in space and energy, we have that the energy dependence of the first moment of the flux expansion &(Z, E) is the same as that of MZ J9lM9 - %lbm Hence the energy dependence of &IQ, &$!Z) and q?(E) is as above. These two conclusions allow the group averaging to be performed knowing the equilibrium scalar spectrum characteristic of the region under consideration. We thus obtain a set of coupled simultaneous equations; in the next section we examine a method of solution of these equations. 3. SOLUTION OF THE GROUP-AVERAGED COLLISION PROBABILITY EQUATIONS Consider the group-averaged collision probability equations (8) above:
If we now expand the L.H.S. of the above equation in terms of surface harmonics, and then multiply both sides by successive Lengendre polynomials, and integrate over dfi bearing in mind the definition of group-averaged parameters as laid out in Section 2.2 we obtain the equations below:
Plcos 8
Several of the integrals involved in the determination of the terms P&2) in the above equations can be performed analytically before calculation of the numerical values is attempted. We shall assume the system considered is effectively infinite in axial extent. Furthermore, we shall write Pm,, the probability that a neutron born in region m and moving in direction Q has its next collision in region n, as the difference between the probability of a neutron born in region m, reaching n and leaving region n, i.e. Here
where a(Q(=JZ(s)&) is the optical path from ds to the surface of region n in the plane z = 0 and in the direction of the projection of P into this plane. Consider the application of these relations to the first of the above equations.
The collisionprobability code+mms
Where we have performed the integral over the polar angle, using the relation:
s UP
K,,(x) =
e-“/sine(sin B)n-l.d&
0
It should be noted that:
s U
e-a~sinesine(c0s ep+l
de = 0.
0
Similar equations for the other two orders in the expansion are given below. PI (cos e)
pli tcOs
e)
cos
v
PI1 (~0s e) sin 9
775
R. G. I-Ihpm
776
where
A=
It should be noted that, due to the assumption of invariance in the z-dimension, any modes in the expansion of the flux which are not symmetric in 0 about the plane z = 0 should be excluded from the solution. However, we shall retain these terms for the present and return to the point later. By an exactly equivalent treatment, higher order expansions could be considered, leading to more equations and considerable increase in complexity. Consider now the integrals over area in each of the above equations. Assume a mesh of equally spaced lines has been placed across the cell in direction q. i dr,=dx,dZ Let : where dxi is!an elemental increase in length in the direction 9 in region i where dl is the mesh spacing for lines placed across the cell. / Now we have also that: &,a+&)
=~&nW
dt.
Using this relation, it may be shown that:
Consequently we obtain finally the equations 4~ &,oA,&,o = 25
5 -
1
i-w-1~~~~
(9)
(12)
777
The collision probability cod-mm
where B = [4rraF,,y,,~~,~r~~~*o~ + SW,+ ~r,~~,~ot,fI[K~~(ortS3+ + KG,cm v where
+ %f,sin Q,+ CW-&~~ ~0sv + ~4,~sin ~~1PG&J1* K&d * = &&iinfon)- L(~iontgoJ- JL(aiinjJ + JGnbi,.t~ii.).
The L.H.S.‘s of the above equations are open to a simple interpretation. Consider the expression for the total flux and the currents in each of the coordinate directions respectively. We have:
jrio. = Jvir = -
&(sZ) cos 0 dS2 s &(S2) sin 8 sin y dn
= T rj&,
s &,,(S2) sin 8 cos cpda = $[;,,.
jar, = s
Consequently, since we are considering an axially infinite system, we would expect the current in the z-direction to be zero. Since the equation for j, contains the coefficient of the mode which is antisymmetric with respect to the plane z = 0, this is in agreement with the conclusion drawn earlier that these modes should vanish from the solution and is also reflected in equation (10) above, the solution of which is indeterminate relative to the solutions of the other equations, since the equation contains no eigenvalue. We shall therefore disregard this equation and solve the remaining set of simultaneous equations. Solutions to the resulting set of equations can be obtained using the method programmed by CLAYTON(1964a) in the code PIP-l, provided the equations are recast into a form in which the resulting matrix possesses Property A (see for example Fox, 1962; YOUNG, 1954). 4. PROGRAMMING
The
ORGANIZATION
program MrNoshas been written to obtain the coefficients for the above equations.
The code applies to a rectangular cell and albedo boundary conditions are applied at each cell boundary. The restrictions in the use of the code are discussed fully below. 4.1 The matrixform of the equations The equations developed in Section 3 are first of all recast into a form suitable for an iterative solution of the type used by the PIP code. We may write:
(Q:&&
+ $:')P$
1
+:',' =i;gi, + (Qi:!-,+,6"! + $,"'>P,!;: + (Q&@
+ S,‘;‘)P,;;
(13)
178
R. G. HARPER
where
Equations (13~(15) are therefore of the form:
i.e. where I is a unit matrix of order NG. Here the matrix
possesses ‘Property A’, hence the method of successive over-relaxation employed by the PIP-l code is applicable to these equations. Here 4 is the vector ((#, j = 1, N), (#, j = 1, N), (c@, j = 1, N), g = 1, G)
The collision probability code-mos
779
and S is the vector ((Soil, i = 1, N), (St,, i = 1, N),
L 1 p(1)
p(2)
p(3)
I
0
0
pf2) a Pa) p?’ L7
pm, I7
p(4) 0
p(K) I7
where the component matrices are of order N and are defined above. Q is a G x G array of diagonal sub-matrices, each of order 3N. The sub-matrix corresponding to row g’ and column g represents the isotropic and anisotropic scattering from group g’ to group g. The terms, along the diagonal, are given by:
<(Q$,, j = 1, N);
= 1, WI.
The emission vector is defined by the above relations together with equations (14). For consideration of isotropic problems, the above equations are simplified. Variables 4t2) and #s) are not considered and the equations relating to these elements are discarded-consequently only matrices of the type P(l) and Q(l) are retained, thus reducing the order of the matrices to be inverted by a factor of three. The equations are now in a suitable form for solution by the program PIP-l. 4.2 Treatment of the geometry problem For solution of this problem, the rectangular lattice cell is sub-divided into smaller rectangular regions in each of which the mean group fluxes are to be determined. The program checks the input information to ensure that there is no overlapping of spatial regions, and that no gaps have been left by the geometry specification. The next stage in the evaluation of the double integrals in each of the P& (i = 1, . . . . 6) is to place a parallel set of lines of constant spacing I and at an angle q, say, across the cell. Along each line, the distances cut off by each rectangular region traversed are determined and stored. The boundary regions cut by each line are also stored. This mesh of lines is then rotated through a small angle AT, and once again the distances cut off by each region traversed are stored. This process is repeated until the mesh of lines has been rotated through 180”, and all geometrical information has been recorded on magnetic tape. Both the mesh spacing and the angular increment are input to the program; and, in order to ensure that adequate information is obtained on small geometrical regions, the facility exists to specify certain ‘fine’ regions. The mesh spacing for all lines crossing these %ne’ regions is reduced to one-tenth of the basic line mesh. All this information is now available as input to the routine for determining collision probability matrices, and is independent of the cross sections in any region. 4.3 Evahution of the region-to-region collision probabilities The rectangular sub-regions of the cell are considered to be numbered sequentially up to N, and associated with these regions is a vector of total cross sections for each group, The cell boundary is itself divided into smaller boundary regions. 3
780
R. G. HARPER
Each line of each angle in turn is considered and each pair of rectangles cut by the line are inspected. The optical paths between each pair of rectangular interfaces are evaluated and then for each pair of rectangles K and J, where K < J, the contribution to the isotropic collision probability P,, in group g is given by:
where Al is the basic line mesh. Agl is the angular direction mesh. If either K or J is specified as a fine region, then the value of AI used in the contribution to the integral is one-tenth the basic mesh spacing. Similar expressions hold for PCz)to PCs). Summing these contributions over all lines at all angles gives an estimate of the collision probabilities. Tabular values of K, (n = 2, . . . , 5) are stored in the machine-these tables have been derived from rational approximations due to CLAYTON (1964b). For a general Consepoint X, say, not at a tabular point, the derivative of K,,(x) is -K&&(x). quently the above tables are sufficient to allow the first two terms of the Taylor expansion for Ki,(n = 3,4 or 5) to be determined, which enable approximations at non-tabular points to be found. Three ranges of the tabulations are considered, the range of argument points being 0(0.01)4, 4(0*1)14and 14(0.5)30. At arguments greater than 30, the value of the function is set to zero. The relative errors in each of the ranges are less than 1~2x 10m6,1.2 x 109 and 0.03 respectively. In parallel with the calculation of region-to-region collision probabilities, the region to boundary probabilities are also determined. These terms represent those neutrons which escape from the cell on tist flight and are defined in exactly the same fashion as the region-to-region probabilities. Consequently there results for each region and in each group six component probabilities each associated with separate modes of the anisotropic source terms. The self-collision terms, i.e. the elements Pit,, in the sub-matrices, are treated separately. These elements are calculated by summing the probability of neutrons, born in a particular region, suffering collision in all other regions, including escape from the cell and subtracting the resulting probability from the normalization value. This normalization value is not necessarily unity for the component probabilities, since the terms being evaluated contain the angularly dependent sources. Should the cell under consideration be a member of a lattice of similar cells, it is now necessary to modify the collision probabilities to allow for this. These modilkations and the method of allowing a macroscopic variation in flux over the cell are discussed in the next section.
4.4 Boundary modzfications Consider the cell under consideration to be a member of a lattice of similar cells, then neutrons leaving this cell by crossing the boundary, will enter a similar cell. Associated with each neutron path will be the probabilities of collisions in the regions crossed until another boundary is encountered, whereupon the process is repeated.
781
The collision probability code-hmos
This behaviour can be expressed as follows: P*ij = Pi* + Pia s P,j $ P$baPpb f P,* $ Pib * P’bl~? Pbj + * * * = Pi, + Pib * P,j{l =
pdj+
Pib
(1
+ pbb + p2bb + ’. ‘1
- pbj
-
Pbb) *
In matrix notation, this equation generalizes to: Lp,*l = Lp,l +
[prbI[l
-
Pbbl-l[Pbrl
where, given several boundary regions round the perimeter of the cell, then: [Prb] Represents the matrix of probabilities of neutrons born in any region in the
cell crossing any boundary region. [P,,] Represents the probabilities of neutrons entering the cell by any boundary region and leaving by any other. [Prb] Represents the probabilities of neutrons entering the cell by any boundary region having collisions in any region in the cell. If the boundary sub-regions are chosen such that the flux along any region is reasonably constant, then the above process will give a reasonable approximation to the collision probabilities in directly reflected, inhomogeneous cells provided that the angular distribution of the incoming neutrons is the same as that of the neutrons leaving the cell. In terms of the component probabilities discussed above for the region-to-region matrices, this means that the boundary-to-boundary and boundaryto-region matrices should be generated for each of the modal angular distributions of the sources. However, in the comparisons discussed below only isotropic distribution of the re-entrant current has been considered, and the boundary-to-boundary and boundary-to-rectangle matrices appropriate to this distribution have been used to modify all component matrices of the region-to-region probabilities. The elements of these boundary term matrices have similar forms to those of the region-to-region matrices. Typically, the elements of the matrix associated with the second mode are given by relations of the form:
dl sin (gl - 7r)
or
dl cos (9 -
7r) ’
tk+a
for x and y boundaries respectively. An alternative, and rather more exact, treatment of the boundary has been developed. In this approach neutrons are tracked until they reach a boundary. Reflection takes place at this boundary and the tracking procedure is continued until a further boundary is reached. This procedure is repeated up to a maximum of 7 reflexes. It is assumed that further first-flight probabilities are negligibly small. The
182
R. G. HARPER
small residual is added to the self-collision probability of the source region to preserve neutron balance. Results using this technique are also discussed below and compared with the first approach. To allow for non-symmetrical cells, the facility exists to allow neutrons to escape from a boundary region on one side of the cell and re-enter the cell at a boundary region located at a similar position on some other side. One further facility is incorporated in the treatment of boundary conditions-this applies to the situation arising when the cell under consideration is placed in material of finite extent. Consequently a macroscopic flux shape is imposed on the cell, and reflective boundary conditions at the cell boundaries are no longer appropriate. In this situation, the ratio of the number of neutrons leaving a particular boundary to the number entering the cell across the same boundary can be found from the macroscopic flux. Hence a knowledge of this ratio allows the probability of entering the cell to be evaluated. Given the macroscopic flux as &(x, y), we have, on a boundary,
g=-= JIN JoUT where K = x or y as appropriate. If now the boundary to rectangle probabilities are multiplied by the appropriate a-values determined from, say, a full-core diffusion theory calculation, then the solution of the cell calculation will be appropriate to the true physical situation. Sufficient information is now available to enable the complete modified collision probability matrices to be generated. The size of problem which can be treated by the present program on the IBM 7094 is limited by the storage facilities required for the PIP-l code. In order to ensure rapid solution of the matrix equations, all work is performed in the core, tapes being used only for input and output of information. Consequently, up to 40 spatial regions may be considered in the isotropic code, and up to 13 in the anisotropic version. Up to 16 boundary regions are allowed. The maximum number of groups must be less than 40, but is related to the number of regions used in the calculation. Fully dynamic storage is available, and the total storage available, 25,000 locations, must be greater than the required problem storage, given by: KG(K+ 6)+ N&G + 4) + K(G, + 4) + G, where K the number of regions in the isotropic case or three times the number of regions in the anisotropic case G the number of groups N, the number of materials G, the lowest energy group with non-zero fixed source G, the lowest energy group with fission source. This expression is identical to that quoted in the PIP-1 report by snwrm computer.
CLAYTON for the
The collisionprobabilitycode-rmm
783
The running time for a typical isotropic problem of 40 regions 6 groups with 14 boundary regions is about 5-6 min. For an equivalent anisotropic problem of 13 regions and 6 groups, the running time is only marginally greater. 5. COMPARISON OF THE PRESENT METHOD WITH OTHER TECHNIQUES
Two simple situations are discussed and the method is compared in one case with other two-dimensional solutions of the Boltzmann equation and in the second with a one-dimensional solution of the Boltxmann equation capable of treating anisotropic scattering. Case 1
The first problem examined is a two-region cell consisting of two similarly orientated rectangles with a common centre. The inner rectangle is a fuel region with dimension 3.048 x 0.103 cm, and the outer region is moderating material with outside dimension 34086 x 0.33036 cm. Cross sections for both materials in seven groups are supplied, and the resulting seven-group eigenvalue problem is solved with reflecting boundaries along the outside of the cell. The isotropic option of the MINOS calculation has been compared with an isotropic collision probability code PIJ (BEARDWOOD, 1965) and with a two-dimensional integration of the differential Boltxmann equation in the S,, approximation using the Los Alamos code 2D X-Y (BENGSTON et al., 1961). Results are given below. These indicate good agreement between both collision probability methods, and show a tendency for the 2D X-Y results to approach the collision probability solution as the order of S, approximation is increased. Case 2
The second comparison was made on a slab cell composed of two material regions ; a fuel region of thickness O-4cm and a moderator region of thickness O-5cm were considered. The moderator was assumed to be capable of anisotropic scattering. Free boundaries were assumed at each side of the cell. Calculations were performed assuming a uniform source density in each region into the top group of a four-thermalgroup scheme. MINOScalculations were compared with results obtained from a one-dimensional integration of the differential Boltzmann equation in the double P, approximation using the code TET-Pl (DAWSON, 1961). This code is capable of treating anisotropic scattering. Results are presented below in general form for calculations using these codes in both the isotropic mode and in the anisotropic mode up to the Pl term in the scattering cross section. Figure 1 shows predictions of the flux distribution in the fourth neutron group using various options of the isotropic MINOS code. The results should all be presented in histogram form, but to avoid confusion this has been done only for the TET-PI case, other points indicate histogram levels. The two calculations using isotropic boundary conditions, at the top and bottom boundaries, were performed on cells with heights 2-Ocm and 8-Ocm, these being the small and large cells respectively. This comparison indicates the effect of the somewhat
784
R. G. HARP= TABLE1.--coMpARIsoN
Method Collision probability
OF RRACITWTY AND DISADVANTAGE FACIQR
Measure of complexity 60 Lines 32 Angles
Code MINOS
(Reflected boundaries)
100 Lines 32 Angles 100 Lines 16 Angles &
CdiSiOIl
(rsotrop~zLkIaries)
pgb&;z probability
PIJ
Discrete ordinates Discrete ordinates
2D-XY 2D-XY
S‘
km
Fuel absorptions Moderator absorptions
1.1667
1.6822
l-1675*
1*6905+
1.1698
1.6941
1.1775
1.7218
1.1830
1.7463
+ These results apply to a 9-0~11problem. 9.01
D
a Fro.
l.-Xomparison
of
MINOS
with TET-Pl for slab calculation-isotropic
case.
+
6.0
z-tf
7.0
?f
-+-n-
ZiYJ ,&,
0
0.1
FIG. 2.-Comparison
,
0.2
0.3
,
04
-TET o ET
PI PI
isotropic atisotropic
;gJig;~;s 0.5
0.6
0.7
0.6
0.9
cm of MINED with TET-PI for slab calculation-anisotropic
case.
The collision probability code-mm
785
incorrect boundary conditions, the discrepancy between MINOSand TET-Pl being reduced when the horizontal boundaries are moved further apart. The correct reflecting boundary conditions produce a flux in good agreement with TET-Pl . Figure 2 iUustrates the same comparison for the anisotropic situation. For comparison the TET-Pl result for the isotropic case is included. The same general remarks apply to this comparison as previously, except that it should be noted that the discrepancy between the small and large cell calculations with isotropic boundary conditions is greater, indicating that the boundary conditions are less appropriate in this case. As in the isotropic case, the MINOScalculation with reflecting boundaries is in reasonable agreement with the TET-PI calculation. 6. CONCLUSIONS
A method has been described which solves the Boltxmam equation in twodimensional (x-u) geometry, allowing anisotropic scattering, by the use of the firstflight-collision probability approximation. Comparisons of the method with other two-dimensional solutions of the Boltxmann equation and with a one-dimensional solution capable of treating anisotropic scattering have been presented. Results in all cases indicate close agreement between the present method and other techniques. The computing time required for running this program is comparable with dii% sion theory methods and considerably less than equivalent transport theory solutions. Acknowledgment-The author is indebted to Mr. M. J. TERRY AEE Winfrith for permission to quote the results of his work with the PIJ and 2D X-Y codes. REFERENCES BEARDWOODJ. E. CLAYTONA. J. and PULLI. C. (1965)Pmcee&gs
ofthe Conference on the Applicution of Computing Methoak to Reactor Problems, Report ANL 7050. BENGSTON J. et al. (1961) Report AGN TM-392. CARLW~I. (1966) Private communication. CLAYTON A. J. (1964a) Report AEEW-R 326, &WTON A. J. (1964b) J. nucl. Energy (Parts A/B Reactor Sci. Technol.) l&82. DAWSON C. W. (1961) Report AML 120. Fox L. (1962) Numerical Solution of Ordinary and Partial D@vential Equations, Pergamon Press, Oxford. HYS~P J. (1963) 1. nucf. Energy (Parts A/B Reactor Sci. Tech&.) 17,237. MACROBERT T. M. (1947) Spherical Harmonics, Methueu. STACER. H. W. (1962) Report W/AT 909. STUARTG. W. and WOODRUFP R. W. (1957) Nucl. Sci. Engng 3,339. TAKAHAS~H. (1966) Nucl. Sci. Engng 24,60. WEINBERG A. M. and WIGNBRE. P. (1958) 27tePhysical Z%eoryof Neutron Chain Reactors, University of Chicago Press. WWLER G. (1963) Report TNPG/PDl. 571. YO~G D. M. (1954) Trans. Am. math. Sot. 76,92.
Prcsa Ltd. Printed in Northem Ma&
THE COLLISION PROBABILITY CODE-MINOS * R. G. HARPW Rolls-Royce Limited, Derby (First received 1 August 1966 and injkalform
14 April 1967)
Abstract-A method is developed to solve the energy-dependenttransport equation with anisotropic scattering, by the method of collision probabilities, in carte&n pmetry. The program MPM)G has been written to produce the resulting collision probability matices in a typical lattice cell. These are then inverted using other programs. The running time of the scheme is comparable with multigroup two-dimensional diffusion theory. Comparisons of the method with other theoretical techniques are reported. 1. INTRODUCTION THE most frequently used technique for the prediction of the neutron flux distribution in reactor systems is the method of diffusion theory. Frequently this method is used in situations where the conditions for the applicability of the method are violated, e.g. close to material discontinuities and throughout heavily absorbing regions. Recourse to correlated data allows predictions to be made of specific parameters at the expense of generality in application. Alternatively, solutions Qf the transport equation in the integro-differential form lead to very timeconsuming programs, especially if more than one-dimension is considered. However, another approximation, applied to the integral form of the transport equation, is receiving considerable attention recently-this is known as the method of collision probabilities. In this approach the total reaction rate in any region of the system is expressed in terms of the sources, either direct or scattering, in all other regions of the system multiplied by the probability that the next collision will be made in the region under Now, in general, the collision probabilities depend upon the consideration. flux shape in the source region, but under certain circumstances it has been demonstrated (see for example STUART and WOODRUFF, 1957 or WEXLER,1963) that the collision probabilities are insensitive to the actual flux shape used. STACE(1962) shows that the criterion for the applicability of this assumption varies with the ratio of the scattering to total cross sections in the region being considered. His conclusions are discussed below. If, for all regions considered, the above criteria are satisfied, then we can choose a flux shape which is spatially independent within each region. Consequently the collision probabilities can be determined to sufficient accuracy ab initio, without recourse to any iterative approach. This step leads to significant saving of computer time when compared with finite difference methods. Finally, if some alternative technique is available for the treatment of large extents of moderator-where the use of collision probability techniques would lead to an unrealistically large number of regions-then the method provides a rapid and accurate prediction of the neutron spectrum, upon which considerable reliance can be placed from a design point of view. The collision probability approach is developed below for the case of anisotropic 767
.
768
R G. HARPER
neutron scattering-this is an extension of the more usual treatment which assumes isotropic sources and fluxes (see however the work of HYSLOP(1963) which considers anisotropic fluxes with isotropic neutron scattering and the paper by TAEMIASHI (1966) considering anisotropic flux with anisotropic scattering in a cylindrical system). Also CARL~IK(1966) has suggested that the extension of collision probability methods to take account of anisotropic neutron scattering events can be effected using Gaussian quadrature methods. Results obtained using the code are compared with a onedimensional anisotropic transport code TET-Pl and also some comparisons with a two-dimensional isotropic transport code 2D X-Y are reported. 2. DERIVATION
OF THE COLLISION
PROBABILITY
EQUATIONS
The collision probability approach consists essentially of describing neutron transport in terms of the probability that neutrons born, or scattered, in one region of the cell should have their next collision in some other region. The theory of the method as usually developed makes, as fundamental, the assumption that all emission sources are isotropic and that neutron scattering is isotropic in the Laboratory System. Any treatment of anisotropic scattering will result in anisotropic fluxes and consequently anisotropic sources and will therefore require that explicit information relating to the angular dependence of the neutron flux be available. These two requirements are not compatible. However using plausible physical assumptions, it is possible to obtain solutions of the equations within the framework of collision probability techniques and also to take account to some extent of the angular dependence of the neutron flux and neutron scattering events. 2.1 The collisionprobabilityform of the integral Boltzmannequation withanisotropic scattering Consider the point S at position r relative to some origin 0. Then if a direction P is chosen through this point S, the position of any point P on the line in the negative direction of S2can be described as (r - R&2),where R is the distance from S to P and Sz is a unit vector. The number of collisions in dR at r, for neutrons moving in direction S2along the line, is given by: E)+(r, E, Sz) dE dR = [ dR’PI-k!k
&(r,
JR
‘S(r - R’S& E, a) dE + *, d&2’&, cl(r - R’S2, E’S2’--+ En)+@ - R’QE’, Sz’) dE’ dE s I
X
$1 _
(1)
dS&’ F(r - R’Q, JI) sn s E’ x y(r - R’B, E’) c,(r - R’B, E’)#(r - R’P, E’, a’) dE’dE _
Here the element of solid angle is taken to be unity. Where the term in brackets represents the source of neutrons moving in direction Szat the point (r - RQ), and where: P,-%+
= &(r, E) dR exp (--LB’&(S)
dS)
769
The collision probability code-rmos
represents the probability that neutrons born at r - R’S2 and moving in direction n should suffer their next collisions in dR at r. Ais an eigenvalue, the inverse of which is the criticality factor of Ken. S(r, E, $2) is a term containing all sources external to the system. Consider the above equation applied to a multi-region system of N regions, with constant cross sections in each region, then we may write for the collision rate at r in direction S2:
&(r, E)&, E, Q) dR dE =
s
5 ,dR’P,_&2 ~
i=l
Rt
‘S(r - R’O, E, a) X
+ n dL2’1, x,,(E’Q’ + ES&j0 - R’Q, E’, Q’) dE’
s
+Iz
dE
s
sn
d&Y F(E)v,(E’) z,,(E’)#r - R’S& E’, a’) dE’ s 1’
Multiplying both sides by dA, we may replace the integrals over the line elements, dR, by volume integrals, dr. Performing the volume integrals over dr’, i.e. along the line through r in direction SZ,and also the integrals over the energy group intervals, we obtain: &,(+A&
a) dr = jr K(r, n)%,(r, WL,,(~)
dr
+ 47A 5 &,K(r, a) 2 vi,. &,k
Q)+&r, Q2,WP”(Q) dr i-w,0
d-1
i-1
{g =
1,2, . . . , G). (2)
Where Vi(r, SQ = dA . Ri
(34
huh Q, W =
V,(r,WUr, W (a) dr = i-+r.n
W
dES(r’, E, rR)P (E) r_R’Ll -)r
.
(W
770
R. G. HARPER
Equation (2) above has expressed the collision rate in dr at point r for neutrons in group g moving in direction P, in terms of the angularly dependent flux and sources along a line passing through r direction 51. In the next section the group-averaged parameters above will be obtained, subject to two simplifying assumptions. These are: (i) That both the flux and source can be considered to be spatially constant in any region. (ii) That the collision probability sections.
can be obtained using group-averaged
cross
2.2 Production of group-averaged parameters In order to include explicit information on the angular dependence of the scattering events, the flux and scattering cross section will be expanded in Legendre polynomials. The coefficients of the polynomials in the flux expansion will be the required solution. Consider the expansion of the differential scattering cross section into surface harmonics : a@ --t E’, e) = $ $: (2z+ l)uz((E + E’)P,(cos e)
(4)
where 19is the angle through which the neutron is scattered. Let Q be the unit vector in the direction of the incident neutron Q’ be the unit vector in the direction of the emergent neutron. Then using the addition theorem for Legendre polynomials we obtain: P,(COSe) = ~,(s2. s2’) = P,(COSe)P,(cos en) + 2 i (I - m)! m=~(2 + m)!
x
P,~os
B)P,~cos
e’)
cos
[m(p, - y’>]. (5)
Here axes have been chosen such that cos 0 is the projection of &Zonto the Z-axis and 8 is the angle between a fixed line in the plane 2 = 0 and the projection of Sz onto this plane. Consider also the expansion of the flux into surface harmonics: a,M, #(r,E,Q
E)p,(cos 0) +$1[a,“x,~9a(r, W,Q(cos 0) cos (49)
= Ta
(6) [
+ b,QXPa(r,JW,a(cos 0) sin (c~$lI
The co-ordinate system is the same as in (5) above. Here we have:
a,x,(r, E) = 5,(r, E) =
#r,EW,QdSk
(7)
Given equations (5) and (6), together with the subsidiary relations (7), we can now construct the group-averaged parameters in the equations (3a)-(3e) above.
The collisionprobability code-mms Consider for example equation (3~):
dE &(E’sZ’ --t ELl)+(r,l, E’, Sk’)
co@’ - JX,(r,, E%(Q) + cl@’-+ E)&(r,,E’)h(Ql [ + 1~3’- E)Pl’ (Q)[Ell(ri, E’) ~0s v + rl?(r,,E’) sin~11 where only first-degree terms have been retained. Finally, performing the integration over a;ri, gives: V,(r, Q)Ur,
Q, a’) CS,,*,, (r, Q, S2’) = V,(r,S2)/zfldE’c+1dE
&(E’ -+ EM,
1
fi2, E’)P&os 0) + &(E’ -+ E)E,,(r, a, E’)P,‘(cos 0)
[ + c14E’ -+ E>P~(cos @[&t(r, Q, E’) cos q + qlt(r, S2, E’) sin ~1
where the l&r, P, E’) etc. define mean components of the flux along the line through r in direction P given by: l&r, S2E) 3 s VikW
& hh
E’)/Ur, a).
Let us now make the two assumptions below: (i) that the collision probability in an energy group g, say, can be obtained using effective cross sections in that group, and
(ii) that the scalar and angular components of the neutron flux and source are constant over each region.
The first assumption follows from using monoenergetic collision probabifity theory within each neutron energy group. The second assumption has been extensively examined by several authors-STUART and WOODRUFF,WEXLER,and STAR in the
772
R a. HARPER
isotropic case. STACE,in particular, concludes, that independent of the ratio of scattering to total cross section, the flat fhrx assumption should be satisfactory for regions where the mean chord is less than about two mean free paths. He also derives a second condition, necessary for the solution to be independent of the flux shape in the region, namely: (&,P&,[$,(r)l
Q1
where P&,(r)] is the self-collision probability in regionj determined from the actual flux shape in the region. Hence, in moderating materials, this is a much more stringent criterion and is true only if P,,< 1, i.e. if the regions are very small compared with a mean free path. Making use of the above assumptions, and integrating equations (2) over dr, for a flxed line direction P we obtain:
I ‘=1,2,...,N
g = 1,2,. ..,G Here we have expressed the mean angular flux over a region j in group g in terms of the mean angular flux in all other regions and groups. The analysis so far has been general, in that no assumption has been made concerning the form of the scattering cross section or the weighting spectrum to be used for producing group-averaged parameters. In a typical reactor physics problem, it is well known that the flux in a homogeneous material region large in terms of mean free paths is characteristic of the material of the region and is, therefore, separable in space and energy. Consequently a method of performing consistent group averaging of equations suggests itself. Let us suppose that the characteristic scalar spectrum of the material of region i has been determined, and let the spectrum be represented by #CO(E). Then we can average the zero-moment cross section &&!?’+ E) in equation (8) relative to this scalar flux, i.e. the zero-moment of the expansion of the flux &,@) has the same energy dependence as this spectrum. To perform the averages of the first moment of the cross section, we note the relation between the zero and Crst moments of the neutron flux, given by WB~VBERO and WIGNER(1958)p. 233 for the PIspherical harmonics approximation in the case of anisotropic scattering:
The collisionprobability code-rams
773
Bearing in mind the fact that &Z, B) is separable in space and energy, we have that the energy dependence of the first moment of the flux expansion &(Z, E) is the same as that of MZ J9lM9 - %lbm Hence the energy dependence of &IQ, &$!Z) and q?(E) is as above. These two conclusions allow the group averaging to be performed knowing the equilibrium scalar spectrum characteristic of the region under consideration. We thus obtain a set of coupled simultaneous equations; in the next section we examine a method of solution of these equations. 3. SOLUTION OF THE GROUP-AVERAGED COLLISION PROBABILITY EQUATIONS Consider the group-averaged collision probability equations (8) above:
If we now expand the L.H.S. of the above equation in terms of surface harmonics, and then multiply both sides by successive Lengendre polynomials, and integrate over dfi bearing in mind the definition of group-averaged parameters as laid out in Section 2.2 we obtain the equations below:
Plcos 8
Several of the integrals involved in the determination of the terms P&2) in the above equations can be performed analytically before calculation of the numerical values is attempted. We shall assume the system considered is effectively infinite in axial extent. Furthermore, we shall write Pm,, the probability that a neutron born in region m and moving in direction Q has its next collision in region n, as the difference between the probability of a neutron born in region m, reaching n and leaving region n, i.e. Here
where a(Q(=JZ(s)&) is the optical path from ds to the surface of region n in the plane z = 0 and in the direction of the projection of P into this plane. Consider the application of these relations to the first of the above equations.
The collisionprobability code+mms
Where we have performed the integral over the polar angle, using the relation:
s UP
K,,(x) =
e-“/sine(sin B)n-l.d&
0
It should be noted that:
s U
e-a~sinesine(c0s ep+l
de = 0.
0
Similar equations for the other two orders in the expansion are given below. PI (cos e)
pli tcOs
e)
cos
v
PI1 (~0s e) sin 9
775
R. G. I-Ihpm
776
where
A=
It should be noted that, due to the assumption of invariance in the z-dimension, any modes in the expansion of the flux which are not symmetric in 0 about the plane z = 0 should be excluded from the solution. However, we shall retain these terms for the present and return to the point later. By an exactly equivalent treatment, higher order expansions could be considered, leading to more equations and considerable increase in complexity. Consider now the integrals over area in each of the above equations. Assume a mesh of equally spaced lines has been placed across the cell in direction q. i dr,=dx,dZ Let : where dxi is!an elemental increase in length in the direction 9 in region i where dl is the mesh spacing for lines placed across the cell. / Now we have also that: &,a+&)
=~&nW
dt.
Using this relation, it may be shown that:
Consequently we obtain finally the equations 4~ &,oA,&,o = 25
5 -
1
i-w-1~~~~
(9)
(12)
777
The collision probability cod-mm
where B = [4rraF,,y,,~~,~r~~~*o~ + SW,+ ~r,~~,~ot,fI[K~~(ortS3+ + KG,cm v where
+ %f,sin Q,+ CW-&~~ ~0sv + ~4,~sin ~~1PG&J1* K&d * = &&iinfon)- L(~iontgoJ- JL(aiinjJ + JGnbi,.t~ii.).
The L.H.S.‘s of the above equations are open to a simple interpretation. Consider the expression for the total flux and the currents in each of the coordinate directions respectively. We have:
jrio. = Jvir = -
&(sZ) cos 0 dS2 s &(S2) sin 8 sin y dn
= T rj&,
s &,,(S2) sin 8 cos cpda = $[;,,.
jar, = s
Consequently, since we are considering an axially infinite system, we would expect the current in the z-direction to be zero. Since the equation for j, contains the coefficient of the mode which is antisymmetric with respect to the plane z = 0, this is in agreement with the conclusion drawn earlier that these modes should vanish from the solution and is also reflected in equation (10) above, the solution of which is indeterminate relative to the solutions of the other equations, since the equation contains no eigenvalue. We shall therefore disregard this equation and solve the remaining set of simultaneous equations. Solutions to the resulting set of equations can be obtained using the method programmed by CLAYTON(1964a) in the code PIP-l, provided the equations are recast into a form in which the resulting matrix possesses Property A (see for example Fox, 1962; YOUNG, 1954). 4. PROGRAMMING
The
ORGANIZATION
program MrNoshas been written to obtain the coefficients for the above equations.
The code applies to a rectangular cell and albedo boundary conditions are applied at each cell boundary. The restrictions in the use of the code are discussed fully below. 4.1 The matrixform of the equations The equations developed in Section 3 are first of all recast into a form suitable for an iterative solution of the type used by the PIP code. We may write:
(Q:&&
+ $:')P$
1
+:',' =i;gi, + (Qi:!-,+,6"! + $,"'>P,!;: + (Q&@
+ S,‘;‘)P,;;
(13)
178
R. G. HARPER
where
Equations (13~(15) are therefore of the form:
i.e. where I is a unit matrix of order NG. Here the matrix
possesses ‘Property A’, hence the method of successive over-relaxation employed by the PIP-l code is applicable to these equations. Here 4 is the vector ((#, j = 1, N), (#, j = 1, N), (c@, j = 1, N), g = 1, G)
The collision probability code-mos
779
and S is the vector ((Soil, i = 1, N), (St,, i = 1, N),
L 1 p(1)
p(2)
p(3)
I
0
0
pf2) a Pa) p?’ L7
pm, I7
p(4) 0
p(K) I7
where the component matrices are of order N and are defined above. Q is a G x G array of diagonal sub-matrices, each of order 3N. The sub-matrix corresponding to row g’ and column g represents the isotropic and anisotropic scattering from group g’ to group g. The terms, along the diagonal, are given by:
<(Q$,, j = 1, N);
= 1, WI.
The emission vector is defined by the above relations together with equations (14). For consideration of isotropic problems, the above equations are simplified. Variables 4t2) and #s) are not considered and the equations relating to these elements are discarded-consequently only matrices of the type P(l) and Q(l) are retained, thus reducing the order of the matrices to be inverted by a factor of three. The equations are now in a suitable form for solution by the program PIP-l. 4.2 Treatment of the geometry problem For solution of this problem, the rectangular lattice cell is sub-divided into smaller rectangular regions in each of which the mean group fluxes are to be determined. The program checks the input information to ensure that there is no overlapping of spatial regions, and that no gaps have been left by the geometry specification. The next stage in the evaluation of the double integrals in each of the P& (i = 1, . . . . 6) is to place a parallel set of lines of constant spacing I and at an angle q, say, across the cell. Along each line, the distances cut off by each rectangular region traversed are determined and stored. The boundary regions cut by each line are also stored. This mesh of lines is then rotated through a small angle AT, and once again the distances cut off by each region traversed are stored. This process is repeated until the mesh of lines has been rotated through 180”, and all geometrical information has been recorded on magnetic tape. Both the mesh spacing and the angular increment are input to the program; and, in order to ensure that adequate information is obtained on small geometrical regions, the facility exists to specify certain ‘fine’ regions. The mesh spacing for all lines crossing these %ne’ regions is reduced to one-tenth of the basic line mesh. All this information is now available as input to the routine for determining collision probability matrices, and is independent of the cross sections in any region. 4.3 Evahution of the region-to-region collision probabilities The rectangular sub-regions of the cell are considered to be numbered sequentially up to N, and associated with these regions is a vector of total cross sections for each group, The cell boundary is itself divided into smaller boundary regions. 3
780
R. G. HARPER
Each line of each angle in turn is considered and each pair of rectangles cut by the line are inspected. The optical paths between each pair of rectangular interfaces are evaluated and then for each pair of rectangles K and J, where K < J, the contribution to the isotropic collision probability P,, in group g is given by:
where Al is the basic line mesh. Agl is the angular direction mesh. If either K or J is specified as a fine region, then the value of AI used in the contribution to the integral is one-tenth the basic mesh spacing. Similar expressions hold for PCz)to PCs). Summing these contributions over all lines at all angles gives an estimate of the collision probabilities. Tabular values of K, (n = 2, . . . , 5) are stored in the machine-these tables have been derived from rational approximations due to CLAYTON (1964b). For a general Consepoint X, say, not at a tabular point, the derivative of K,,(x) is -K&&(x). quently the above tables are sufficient to allow the first two terms of the Taylor expansion for Ki,(n = 3,4 or 5) to be determined, which enable approximations at non-tabular points to be found. Three ranges of the tabulations are considered, the range of argument points being 0(0.01)4, 4(0*1)14and 14(0.5)30. At arguments greater than 30, the value of the function is set to zero. The relative errors in each of the ranges are less than 1~2x 10m6,1.2 x 109 and 0.03 respectively. In parallel with the calculation of region-to-region collision probabilities, the region to boundary probabilities are also determined. These terms represent those neutrons which escape from the cell on tist flight and are defined in exactly the same fashion as the region-to-region probabilities. Consequently there results for each region and in each group six component probabilities each associated with separate modes of the anisotropic source terms. The self-collision terms, i.e. the elements Pit,, in the sub-matrices, are treated separately. These elements are calculated by summing the probability of neutrons, born in a particular region, suffering collision in all other regions, including escape from the cell and subtracting the resulting probability from the normalization value. This normalization value is not necessarily unity for the component probabilities, since the terms being evaluated contain the angularly dependent sources. Should the cell under consideration be a member of a lattice of similar cells, it is now necessary to modify the collision probabilities to allow for this. These modilkations and the method of allowing a macroscopic variation in flux over the cell are discussed in the next section.
4.4 Boundary modzfications Consider the cell under consideration to be a member of a lattice of similar cells, then neutrons leaving this cell by crossing the boundary, will enter a similar cell. Associated with each neutron path will be the probabilities of collisions in the regions crossed until another boundary is encountered, whereupon the process is repeated.
781
The collision probability code-hmos
This behaviour can be expressed as follows: P*ij = Pi* + Pia s P,j $ P$baPpb f P,* $ Pib * P’bl~? Pbj + * * * = Pi, + Pib * P,j{l =
pdj+
Pib
(1
+ pbb + p2bb + ’. ‘1
- pbj
-
Pbb) *
In matrix notation, this equation generalizes to: Lp,*l = Lp,l +
[prbI[l
-
Pbbl-l[Pbrl
where, given several boundary regions round the perimeter of the cell, then: [Prb] Represents the matrix of probabilities of neutrons born in any region in the
cell crossing any boundary region. [P,,] Represents the probabilities of neutrons entering the cell by any boundary region and leaving by any other. [Prb] Represents the probabilities of neutrons entering the cell by any boundary region having collisions in any region in the cell. If the boundary sub-regions are chosen such that the flux along any region is reasonably constant, then the above process will give a reasonable approximation to the collision probabilities in directly reflected, inhomogeneous cells provided that the angular distribution of the incoming neutrons is the same as that of the neutrons leaving the cell. In terms of the component probabilities discussed above for the region-to-region matrices, this means that the boundary-to-boundary and boundaryto-region matrices should be generated for each of the modal angular distributions of the sources. However, in the comparisons discussed below only isotropic distribution of the re-entrant current has been considered, and the boundary-to-boundary and boundary-to-rectangle matrices appropriate to this distribution have been used to modify all component matrices of the region-to-region probabilities. The elements of these boundary term matrices have similar forms to those of the region-to-region matrices. Typically, the elements of the matrix associated with the second mode are given by relations of the form:
dl sin (gl - 7r)
or
dl cos (9 -
7r) ’
tk+a
for x and y boundaries respectively. An alternative, and rather more exact, treatment of the boundary has been developed. In this approach neutrons are tracked until they reach a boundary. Reflection takes place at this boundary and the tracking procedure is continued until a further boundary is reached. This procedure is repeated up to a maximum of 7 reflexes. It is assumed that further first-flight probabilities are negligibly small. The
182
R. G. HARPER
small residual is added to the self-collision probability of the source region to preserve neutron balance. Results using this technique are also discussed below and compared with the first approach. To allow for non-symmetrical cells, the facility exists to allow neutrons to escape from a boundary region on one side of the cell and re-enter the cell at a boundary region located at a similar position on some other side. One further facility is incorporated in the treatment of boundary conditions-this applies to the situation arising when the cell under consideration is placed in material of finite extent. Consequently a macroscopic flux shape is imposed on the cell, and reflective boundary conditions at the cell boundaries are no longer appropriate. In this situation, the ratio of the number of neutrons leaving a particular boundary to the number entering the cell across the same boundary can be found from the macroscopic flux. Hence a knowledge of this ratio allows the probability of entering the cell to be evaluated. Given the macroscopic flux as &(x, y), we have, on a boundary,
g=-= JIN JoUT where K = x or y as appropriate. If now the boundary to rectangle probabilities are multiplied by the appropriate a-values determined from, say, a full-core diffusion theory calculation, then the solution of the cell calculation will be appropriate to the true physical situation. Sufficient information is now available to enable the complete modified collision probability matrices to be generated. The size of problem which can be treated by the present program on the IBM 7094 is limited by the storage facilities required for the PIP-l code. In order to ensure rapid solution of the matrix equations, all work is performed in the core, tapes being used only for input and output of information. Consequently, up to 40 spatial regions may be considered in the isotropic code, and up to 13 in the anisotropic version. Up to 16 boundary regions are allowed. The maximum number of groups must be less than 40, but is related to the number of regions used in the calculation. Fully dynamic storage is available, and the total storage available, 25,000 locations, must be greater than the required problem storage, given by: KG(K+ 6)+ N&G + 4) + K(G, + 4) + G, where K the number of regions in the isotropic case or three times the number of regions in the anisotropic case G the number of groups N, the number of materials G, the lowest energy group with non-zero fixed source G, the lowest energy group with fission source. This expression is identical to that quoted in the PIP-1 report by snwrm computer.
CLAYTON for the
The collisionprobabilitycode-rmm
783
The running time for a typical isotropic problem of 40 regions 6 groups with 14 boundary regions is about 5-6 min. For an equivalent anisotropic problem of 13 regions and 6 groups, the running time is only marginally greater. 5. COMPARISON OF THE PRESENT METHOD WITH OTHER TECHNIQUES
Two simple situations are discussed and the method is compared in one case with other two-dimensional solutions of the Boltzmann equation and in the second with a one-dimensional solution of the Boltxmann equation capable of treating anisotropic scattering. Case 1
The first problem examined is a two-region cell consisting of two similarly orientated rectangles with a common centre. The inner rectangle is a fuel region with dimension 3.048 x 0.103 cm, and the outer region is moderating material with outside dimension 34086 x 0.33036 cm. Cross sections for both materials in seven groups are supplied, and the resulting seven-group eigenvalue problem is solved with reflecting boundaries along the outside of the cell. The isotropic option of the MINOS calculation has been compared with an isotropic collision probability code PIJ (BEARDWOOD, 1965) and with a two-dimensional integration of the differential Boltxmann equation in the S,, approximation using the Los Alamos code 2D X-Y (BENGSTON et al., 1961). Results are given below. These indicate good agreement between both collision probability methods, and show a tendency for the 2D X-Y results to approach the collision probability solution as the order of S, approximation is increased. Case 2
The second comparison was made on a slab cell composed of two material regions ; a fuel region of thickness O-4cm and a moderator region of thickness O-5cm were considered. The moderator was assumed to be capable of anisotropic scattering. Free boundaries were assumed at each side of the cell. Calculations were performed assuming a uniform source density in each region into the top group of a four-thermalgroup scheme. MINOScalculations were compared with results obtained from a one-dimensional integration of the differential Boltzmann equation in the double P, approximation using the code TET-Pl (DAWSON, 1961). This code is capable of treating anisotropic scattering. Results are presented below in general form for calculations using these codes in both the isotropic mode and in the anisotropic mode up to the Pl term in the scattering cross section. Figure 1 shows predictions of the flux distribution in the fourth neutron group using various options of the isotropic MINOS code. The results should all be presented in histogram form, but to avoid confusion this has been done only for the TET-PI case, other points indicate histogram levels. The two calculations using isotropic boundary conditions, at the top and bottom boundaries, were performed on cells with heights 2-Ocm and 8-Ocm, these being the small and large cells respectively. This comparison indicates the effect of the somewhat
784
R. G. HARP= TABLE1.--coMpARIsoN
Method Collision probability
OF RRACITWTY AND DISADVANTAGE FACIQR
Measure of complexity 60 Lines 32 Angles
Code MINOS
(Reflected boundaries)
100 Lines 32 Angles 100 Lines 16 Angles &
CdiSiOIl
(rsotrop~zLkIaries)
pgb&;z probability
PIJ
Discrete ordinates Discrete ordinates
2D-XY 2D-XY
S‘
km
Fuel absorptions Moderator absorptions
1.1667
1.6822
l-1675*
1*6905+
1.1698
1.6941
1.1775
1.7218
1.1830
1.7463
+ These results apply to a 9-0~11problem. 9.01
D
a Fro.
l.-Xomparison
of
MINOS
with TET-Pl for slab calculation-isotropic
case.
+
6.0
z-tf
7.0
?f
-+-n-
ZiYJ ,&,
0
0.1
FIG. 2.-Comparison
,
0.2
0.3
,
04
-TET o ET
PI PI
isotropic atisotropic
;gJig;~;s 0.5
0.6
0.7
0.6
0.9
cm of MINED with TET-PI for slab calculation-anisotropic
case.
The collision probability code-mm
785
incorrect boundary conditions, the discrepancy between MINOSand TET-Pl being reduced when the horizontal boundaries are moved further apart. The correct reflecting boundary conditions produce a flux in good agreement with TET-Pl . Figure 2 iUustrates the same comparison for the anisotropic situation. For comparison the TET-Pl result for the isotropic case is included. The same general remarks apply to this comparison as previously, except that it should be noted that the discrepancy between the small and large cell calculations with isotropic boundary conditions is greater, indicating that the boundary conditions are less appropriate in this case. As in the isotropic case, the MINOScalculation with reflecting boundaries is in reasonable agreement with the TET-PI calculation. 6. CONCLUSIONS
A method has been described which solves the Boltxmam equation in twodimensional (x-u) geometry, allowing anisotropic scattering, by the use of the firstflight-collision probability approximation. Comparisons of the method with other two-dimensional solutions of the Boltxmann equation and with a one-dimensional solution capable of treating anisotropic scattering have been presented. Results in all cases indicate close agreement between the present method and other techniques. The computing time required for running this program is comparable with dii% sion theory methods and considerably less than equivalent transport theory solutions. Acknowledgment-The author is indebted to Mr. M. J. TERRY AEE Winfrith for permission to quote the results of his work with the PIJ and 2D X-Y codes. REFERENCES BEARDWOODJ. E. CLAYTONA. J. and PULLI. C. (1965)Pmcee&gs
ofthe Conference on the Applicution of Computing Methoak to Reactor Problems, Report ANL 7050. BENGSTON J. et al. (1961) Report AGN TM-392. CARLW~I. (1966) Private communication. CLAYTON A. J. (1964a) Report AEEW-R 326, &WTON A. J. (1964b) J. nucl. Energy (Parts A/B Reactor Sci. Technol.) l&82. DAWSON C. W. (1961) Report AML 120. Fox L. (1962) Numerical Solution of Ordinary and Partial D@vential Equations, Pergamon Press, Oxford. HYS~P J. (1963) 1. nucf. Energy (Parts A/B Reactor Sci. Tech&.) 17,237. MACROBERT T. M. (1947) Spherical Harmonics, Methueu. STACER. H. W. (1962) Report W/AT 909. STUARTG. W. and WOODRUFP R. W. (1957) Nucl. Sci. Engng 3,339. TAKAHAS~H. (1966) Nucl. Sci. Engng 24,60. WEINBERG A. M. and WIGNBRE. P. (1958) 27tePhysical Z%eoryof Neutron Chain Reactors, University of Chicago Press. WWLER G. (1963) Report TNPG/PDl. 571. YO~G D. M. (1954) Trans. Am. math. Sot. 76,92.