- Email: [email protected]
Variational boundary conditions for the spherical harmonics approximation to the neutron transport equation
OF PHYSICS:
.iiX’bALS
27,
193-215
Variational Boundary Approximation
(1964)
Conditions for the Spherical Harmonics to the Neutron Transport Equation G.
General
Electt,ic
C’wnpccny,
C.
POMRANING
Vallecitos
9tomic
Labomtory,
Plrasanton,
C’alifornia
11) is shown that the monoenergetic neutron transport equation and the associated boundary conditions can be characterized by a Lagrangian. A proper choice of the trial function for this Lagrangian leads to the widely used spherical harmonics approximation as the Euler-Lagrange equations. Further, a set of boundary conditions for the spherical harmonic equations is the result of the logical application of the variational method. These variational boundary conditions appear to be significantly more accurate than the boundary conditions presently in general use. For example, the use of the variational boundary conditions at a free surface reduces the error (compared with the boundary conditions presently used) in the linear extrapolation distance for the Milne problem by several factors. In particular, the P-l (diffusion) approximation yields a value of 0.7071 (in units of mean free paths) and the P-3 approximation yields a value of 0.7118, both comparing quite favorably with t.he exact value of 0.7104. I. INTRODUCTION A commo~l method of solving the monoenergetic neutron transport equatiou is to expand the angular dependence of the directional flux in a truncated spherical harmonics series (1). This reduces the integro-differential equation to a finite number of coupled differential equatious for t,he spatially dependent expansion coefficients. While the boundary conditions for the transport problem are well set, these conditions cannot always be satisfied exactly by a truncated spherical harmonics expansion. In particular, at an interface between two media the proper transport boundary condition is that the directional flux be continuous iu all directions. It is not, in general, possible to satisfy this boundary condition in a truncated spherical hamonies expansion of the directional flux, and thus approximate boundary conditions may be needed. Further, at the outside houndaries of the system, which shall be referred to as outer surfaces, the general transport
boundary
incoming
hemisphere
condition
condition
exactly
in
is that
the
directional
flux
of solid angle. Again one cannot a truncated
spherical
harnlonics
mation is required. 193
be
satisfy
expansion
specified
in
the
this boundary and
an
approxi-
194
POMRANING
Two sets of boundary conditions at an outer surface for the spherical harmonics equations, due to Mark (2) and Marshak (3), are presently in general use. A detailed discussion of these boundary conditions can be found in ref. 1, and a brief discussion will be given in this paper at the appropriate point. At an interface between two media, the boundary conditions in general use are due to Davison (1) which he derives on physical grounds. The main point to be made here is that all of these boundary conditions are quite arbitrary, and a large number (also arbitrary) of other boundary conditions could equahy well be proposed. Previous authors (4, 5) have shown that the neutron transport equation can be properly characterized by a Lagrangian, and that the spherical harmonics approximation can be derived as the Euler-Lagrange equations of this Lagrangian by the proper choice of the trial function. In this paper we show that the variational method can also logically be used to derive a set of boundary couditions for the spherical harmonics equations. For the outer surface boundary conditions, this involves leaving the proper boundary terms in the Lagrangian. For the conditions at an interface, one cousiders the interface to be a limiting case of a medium whose properties are a continuous function of space. The a priori advantage of these bouudary conditions over those preseutly in use are that they do riot suffer from any arbitrariness. The only assumption required in their derivatiou is implicit in the equations themselves; i.e., the directional flux is assumed expandable in a truncated spherical harmonics series. We emphasize that the use of a Lagrangiau to describe neutron transport is not in any way au assumption, since reuderiug this Lagrangian statiouary is entirely equivalent to solving the neutron transport equation. Thus, if one prefers, the Lagrangian cau be taken as the fundamental equation of ueutrou transport theory. The a posteriori advant#age of the variatioual bouudary conditions is that they are 110 more complex t,han those presently in use, but appear to be significantly more accurate. For algebraic simplicity the formalism in this paper is restricted to slab geometry. Presumably the same formalism could be carried through in other geometries. The spherical harmonics reduce to Legendre polynomials in slab geometry, aud we shall refer to the flux expansiou carrying the first (N + 1) Legendre polynoniials (PO to PN) as the P-N approximation. Sections II aud III of the paper discuss briefly the general variational formulation of neutron transport theory and give the variational derivatiou of the P-N equations. Sectious IV and V discuss in some detail the forniulatiou of the variational bouudary conditions and compare the results with the boundary conditions currently in use. In particular, emphasis is giveu to the outer surface boundary conditions since the variational method yields results quite different from and
NEUTRON
significantly in use.
TRANSPORT
195
THEORY
more accurate than those of Mark and Marshak II.
The monoenergetic H(x, PM&
VARIATIONAL
transport
equation
/-d = m% co,
which are currently
FORMALISM
in slab geometry (-1
5 /.l 5 1;
can be written
(1)
a 5 x 5 b),
(1)
where
and
.fn= I: 4Jof(Po)P?z(Po).
(3)
The spatial region of interest has quite generally been assumed to be a s z 5 b. Here cp(z, +L) is the directional flux, z is the slab coordinate, P is the cosine of the angle between the velocity vector of the neutron and the x-axis, c is the mean number of secondaries per collision, Z is the macroscopic collision cross section, f( pO) dpO is the probability that a scattering event will scatter a neutron through an angle ~0 into C&Q,, S(z, p) is the external directjonal source, and P%(p) is the 71th Legendre polynomial. S, c, and fn have an arbitrary spatial dependence, although for simplicity in notation, we do not explicitly show this z-dependence. The general boundary conditions on Eq. (1) at the outer surfaces are specifications of the incoming flux, i.e., 4% /J) = ~+iP),
(0 < p 5 I),
(3:)
54, P) = A-(P),
(-1
(-5)
5 /.I < O),
where ,4+(p) and A-(P) are arbitrary, but known, functions. We now define the adjoint operator according t.o / dr+~(.r)H*(.r)+~*(~j
= [ d.rc~*(s)H(z)cp(lc)
+ (boundary
terms),
(6’)
where z represents the phase space variables z and p and the integration extends over the region of interest in r-space. Requiring that Eq. (Gj hold for any two functions &z) and (p*(z) yields
196
POMRANING
Thus the adjoint H*k,
equation is
P)cp*(.% PI = uz,
(-1
PI,
5 p 5 1;
a 5 z 5 b),
(8)
where ‘p*(z, p) is the adjoint directional flux, and T(z, p) is the adjoint directional source which, at this point, is arbitrary. The adjoint boundary conditions, also arbitrary at this point, are taken as cp*ca, IL) = B+(PL
(-1
cp*u4 Pu) = B-(P),
(0 -=c P 6 11,
where B+(P) and B-(p) are arbitrary, but specified, functions. evident from comparison of Eqs. (2) and (7) that
5 p < O),
source and boundary ux,
PI = St% -PI,
(11) properly.
That
is, if (13)
B+(P)
= A+(-d,
(13)
B-(P)
= A-(-EL),
(l-1)
then Eq. (11) is satisfied. We now postulate that a proper transport equations, (1) and (8)) Eqs. (4), (5), (91, and (101, is Fb, (P*I = s,’ c&zI: hw*(%
conditions
(10)
It is immediately
4% PI = cp*c% -PI if we choose the adjoint we set
19)
Lagrangian for the transport and adjoint and the associated boundary conditions,
Pcr)N(% cl)&,
P>
For Eq. (1.5) to be a correct Lagrangian, setting the first variation with respect to arbitrary and independent variations of q(x) and (p*(x) equal to zero must be equivalent to solving the transport and adjoint transport equations. It is
NEUTROK
TRANSPORT
197
THEORY
easily verified that setting the first variation of Eq. (15) with respect to (O*(X) equal to zero yields Eq. (1) as the Euler-Lagrange equation with Eqs. (4) and (5) as the subsidiary or bourldary conditions. Similarly, setting the first variation of Eq. ( 1.3) with respect to P(X) equal to zero and retaining the boundary terms from the integration by parts yields Eq. (8) as the Euler-Lagrange equation with Eqs. (9) and (10) as the boundary conditions. Thus we are assured that Eq. (15) is, in fact, a proper Lagrangian. We now assume that in any subsequent use of the Lagrangian, the set of t,rial functions for cp(.r) and P*(X) is restricted according to Eq. (11). This allows us t,o simplify the Lagrangian. In particular, the first variation of the Lagrangian can then be writ’ten, using Eqs. ( 11) through (14)) as
h-ow, since the variations in the interior of x-space are independent from those on the boundaries for any set of trial functions, setting 6F = 0 in Eq. (1Gj implies t,wo separate equations, i.e.,
ss b
a
1
dx
-1
d/.M(~, P)P(X, Co - X(x, EL)IBp(x, -p)
= 0,
(17)
Equations (17) and (18) are the starting expressionsfrom which we shall proceed. In particular, Eq. (17)) with the proper choice of the trial function, will yield the P-N equations and the boundary conditions on these equations at an interface between two media. Equation (18), with the same trial function, will yield the boundary conditions for the P-N equations at an outer surface. III.
DERIVATION
OF THE
P-iv EQUATIONS
We take as our trial function a truiicat,ed Legendre polynomial seriesin angle, i.e.,
198 Substituting [
POMRAKING
Eq. (19) into Eq. (18) yields
dz I: d/J [n(%
P)Y$
(“i”)
(Pm(Z)Prn(PcL) - S(% Pi] (20) ( -l)nPn(Ph%(Z)
.lgq) Now, a priori all of the &D~(x) are arbitrary implies (N + 1)Equations:
and independent.
= 0. Thus Eq. (20)
(0 5 n 5 N). We recognize Eq. (21) as the widely used P-N equations. That is, the classical derivation of the P-N equations is to expand the angular dependence of the directional flux in a truncated Legendre polynomial series and use orthogonality to obtain the differential equations for the spatially dependent. expansion coefficients. This is precisely the content of Eq. (21). Two points should be made here. First, it is not necessary to choose Legendre polynomials as the trial functions in order for the variational method to yield the P-N equations. Any set of polynomials, or even the simple integer powers of I*, would yield the the same result since the ingredients of the trial functions are the same in all cases, i.e., linear combinations of the integer powers of p. The Legendre polynomials were chosen so that the result, Eq. (21)) would be readily recognizable. Secondly, while a priori all of the &%(z) are independent, a posteriori this may or may not be the case depending upon the order of the expansion. As is well known (1) , the solution of Eq. (21), which contains (N + 1) unknowns, for cpn(z) contains (N + 1) constants if N is odd and N constants if N is even. Thus for N even, prior to applying the boundary conditions, one of the cpn(x) is linearly dependent upon the other N p,(x), and therefore a posteriori the Bv~(z) are not all independent. Let us explicitly write out the P-N equations as given by Eq. (21). Performing the indicated angular integrations, we have
+
&u(z) dz
Zcpl(X)
+ %h42)
=
Xl(X)
+
Wlcpl~~),
= A%(z) + =fNPN(Z),
NEVTRON
TRANSPORT
199
THEORY
where X,,(Z) is the ?Ith Legendre moment of the source, i.e., (33) IV.
BOUSI~lARY
CONDITIONS
AT
.4N
INTERFACE
If the medium properties have a discontinuity at some point within the system the usual method of solution of the P-N equations is to find the solution in the regions on either side of this discontinuity and join these solutions by prescribed boundary conditions at the interface (discontinuity). Since adding a region introduces another interface, the number of conditions that can (must) be satisfied at an int’erface is equal to the number of constant8s in the solution within a region. We know (1) that within a region the solution of the P-N equations contains (N + 1) constants if N is odd and N constants if N is even. For example, assuming a solution of the form (L&(z) = g,P
(24)
within a homogeneous region and setting the coefficient determinant of the g,,‘s equal to zero yields (N + 1) values for v if N is odd and N values for v if N is even. Thus we must construct (,N + 1) boundary conditions for N odd and N boundary conditions for N even. Kow, in deriving the P-N equations no assumptions were made concerning the spatial dependence of the mediunl properties, and thus the equations hold for a medium whose properties vary arbitrarily with Z. In particular, the P-N equations are valid at an interface and one should be able to determine the boundary conditions from the equations themselves. Because of the inherent differences in the even and odd order approximations (N even or odd ), we consider the two cases separately, taking the interface to be z = _Y(~ . ODD ORDER
We number the (N + 1) equations given by l?q. (22) from zero to N. Integrating the zeroth equation over x from (z c - t ) to (zO + E) and letting E approach zero, noting that in this limit only the derivative term contributes, yields the fact that cpl(x) is continuous at the interface, x = pg. A similar integration over the second equation, using the fact that cpl(zO) is continuous, yields the fact that (p3(x) is continuous at the interface. Continuing this procedure of integrating over the even numbered equations yields the fact that all the od’dmoments are continuolls at an interfare. The conditions on the even moments are found by integrating over the odd numbered equations. In particular, integration over the Nth ecluation yields the fact that (~~~~(2) is continuous at the interface. Integration over the (N - 2)th equation, using the fart that (P~-~(,Q)is continuous, yields the fact that (P~-~(x) is continuous at the interface. Continuing this pro-
200
POMRANING
cedure, we deduce that all the eoen moments are continuous at an interface. Thus we conclude that for an odd order expansion all moments are continuous at an interface. These are the boundary conditions universally used for odd order expansions, but the usual argument (1) to obtain these conditions is that the directional flux must be continuous in all directions at an interface, and thus all of the moments, cp,(x,,), must be continuous. EVEN
ORDER
As is well known (I), the argument that all moments must be continuous at an interface leads to one too many boundary conditions for an even order expansion. It is here that the procedure of obtaining the interface boundary conditions from the equations themselves is very useful and is, in fact, the only consistent treatment. As in the odd order analysis, subsequent integration of the zeroth, second, fourth, etc. equations leads to continuity of the odd moments at an interface. Integration of the first, third, fifth, etc. equations yields the fact that
cpzn(zo) + ( g$ >w&f&“),
(0 5 n I (N - 2)/2)
(25)
is continuous. By appropriately combining these continuity conditions, it is easily shown that Eq. (25) is equivalent to demanding continuity of
(0 5 n 5 (N - 2),'2).
(26)
Now, according to Sansone (B), a property of Legendre polynomials for n an integer is
Pz,(O) =(-1)”o1(s >.
(27)
Using Eq. (27) in Eq. (26) yields
(ozn(zo)- p2,(o)Po(~o),
(1 5 n 5 N/2)
(28)
is continuous. Noting that PZn+l(0) is zero for n an integer, we can demand continuity of the odd moments and continuity of Eq. (28) for the even moments by requiring con
- pn(o)cpo(~o),
(1 5 n 5 N)
(29)
be continuous. Equation (29) represents the N boundary conditions at an interface for an even order P-N expansion. We note that the scalar flux, cpo(z),is not continuous at an interface. As remarked by Davison (I), one should not be surprised to find a discontinuity in the scalar flux in an approximate solution of the transport equation.
NEUTRON
TRANSPORT
THEORY
201
At this point it is iuteresting to compare Eq. (29) with the even order interface boundary conditions proposed by other authors. Davison (1) has suggested the identical conditions as found in Eq. (29). His derivation, however, is conpletely different from the one giveu in this paper and is based on rather iuvolved physical arguments. The work by Rumyantsev (7), when specialized to slab geometry, also yields Eq. (29). His considerations were much the same as those in this paper, but were not in the context of the variational method. Finally, following the spirit of Marehak’s conditions at an outer boundary (discussed in ref. 1 and in t’he next section of this paper), the Ilarshak conditions at an interface are logically taken as continuity of s0
’ &JP*n+l(EL)pF(Xo, p),
t.0 5 11s (N - 2)/2),
(30)
(0 5 ‘n 5 (N - 2)/2).
(31)
and 0 s-1
dPPzn+lbFL)&O,
PFC),
That is, Marshak demands continuity of the odd Legendre half-moments of the directional flux. Setting n = 0 in Eqs. (30) and (31), we see Marshak’s conditions have the advantage of including the physical condition that the partial currents across an interface are continuous. The use of the expansion for the directional flux, Eq. (19)) in Eqs. (30) and (31) yields N conditions on the moments, p,(z), at the interface. Subsequent addition of Eqs. (30) and (31) and use of the orthogonality and parity properties of the Legendre polynomials yields the fact that rnn+1(~0),
(0 5 n 5 (N - 2)/2)
(32)
is continuous. Subtraction of Eq. (31) from IQ. (30) and use of the parity properties of the Legendre polynomials yields continuity of 1 N/2 dpP2,,+l(r)m~o (h + 1)P2m(P)(P2m(Z”), (0 5 n s (N - 2)/2). (33) s0
According to Bateman (8))
s0
l d.xP,(xP,(z)
= ; [(CT- Y)(U + lJ + l)]-’ .[A sin (UT/~) cos (VT/~) - .4-l sin (m/2) Cos(m/2)1,
where A4 = r(95 + v/2)r(l qj_; + u/a)r(l
+ u/2) + V/2) ’
(35)
202
POMRANING
and I’(X) is the gamma function,
defined as (36)
Setting v = 2n + 1 and c = 2~1 (both nl and n are integers) properties of the gamma function, we can simplify Eq. (34) to 1 s0
dxPz,+~(x)Pzm(x)
=
(-1y+m+l (2m)!(2n (71!l)b!)“2(2m+sn+1)(2,m - 2n -
Using Eq. (37) in Eq. (33), we demand continuity ji!
and using the
+ l)! l)(nk + n + 1)’
(37)
of
( -1)1)2(4m + 1)(2m)! [ (m!)222~(2m - 2n - l)(nz + n + 1)
1(0*&o),
(38)
(0 5 n 5 (N - 2>/2). Noting that Eq. (27) can be written P2n(0)
as =
(-1)“(2n)! yn(n!)2
’
Eq. (38) is identical to requiring continuity of g
1
(4m + l)PdO) (0 5 n 4 (N - 2)/2). [(2m - 2n - ~)(HL + n + 1 (P2m(Xo)7
(40)
Equations (32) and (40) are the N Marshak conditions at an interface. A comparison of these conditions with the variational conditions, Eq. (29)) shows that the two are not the same. For example, in the P-2 approximation the variational conditions, Eq. (29)) reduce to cpl(x0), (o~(20) + >~c,cJ,,( zo) being continuous, while the Marshak conditions, Eqs. (32) and (40)) reduce to (Ok, (p2(zo) + +&,(z,) being continuous. Since the variational boundary conditions (the same as those of Davison and Rumyantsev) are consistent with the P-N equations, whereas the Marshak conditions are not, it is suggested that the variational conditions, Eq. (29)) should be used. V. BOUlVDARY
CONDITIONS
AT
A?;
OUTER
SURFACE
We now consider the boundary conditions at an outer surface. It is here that the variational method yields entirely new results which are no more complex to apply than the boundary conditions of I\Iark and Xarshak, currently in use, but appear to be significantly more accurate. Before proceeding, however, let us briefly discuss the Mark and Marshak boundary conditions. To be specific, we consider the left hand boundary of a system, taken as x = a. The right hand boundary is handled similarly. We know from the number of constants in the
NEUTROS
TRANSPORT
203
THEORY
solution of the P-N equations that we require (N + 1)/2 outer boundary conditions at z = a for N odd and N/2 conditions for N even. At a free surface (no neutrons entering the system) the rigorous transport boundary condition is da, ~1 = 0,
(0 5 /.Lf
1).
(41)
It is obvious that Eq. (41) cannot be satisfied exactly by a truncated Legendre polynomial expansion of the directional flux. As an approximation, Mark (2) has suggestedsetting (43)
da, h) = 0,
where the /*i are the R,/2 positive roots of Pn(pi) = 0, and R = N + 1 for N odd and R = N for N even. Using the expansion for t#hedirectional flux, Eq. (19)) in Eq. (42) yields R/2 linear relationships, which is the required number of boundary conditions, on the moments at the free surface. Davison (1) has proved that Eq. (42) is equivalent to replacing the vacuum to the left of the free surface with a pure absorber and treating the resulting interface according to his boundary conditions discussedin the last section. Marshak (3) has suggested
(0 -I n 5 (R - 2)/2),
(4.3)
to approximate Eq. (41). Setting 71 = 0 in Ii:q. (43) shows that Rlarshak’s boundary conditions include the physical condition that the current into the system is zero. As with the l\lark conditions, the llse of Eq. (19) in Eq. (43)
The use of the directional flux expansion, Eq. (19), in E:q. (44) allows one to carry out all of the angular integrations, since 4+(p) is presumed known. This leads to R/Y inhomogeneouslinear relationships between the moments at z = a. Based on the above discussion,it is evident that the Marshak homidary condi-
204
POMRANING
tions at au outer surface enjoy much more usage than do the Mark conditions. Accordingly, the emphasis here will be on the comparison of the accuracy of the variational boundary conditions with that of the Marshak conditions. AIark’s conditions will be mentioned only in an incidental mauner. Before proceeding with the variational analysis for a general N, we consider the special case of N = 1. This is a particularly important case since it corresponds to the widely used diffusion approximation. Further, the consideration of this special case will serve to point the way toward the more general analysis. We consider, within the framework of diffusion theory, the quite general problem of a free surface at the left hand face, x = 0, and a ( -p)“ entrant flux distribution at the right hand face, x = zO. Thus we set A+(P)
= 0,
(45)
A-(/l)
= S(-/.4Ln,
(46)
where S is a constant indicating the magnitude of the entrant flux and n is an arbitrary nonnegative integer. The first variation of the boundary terms of the functional, Eq. (18), for this problem is 1’ &wP(O, co&40, -PI 0
+ s: &P[S( -PI”
- (0(X”, P)l&h”,
-pj
= 0.
(47)
Using the directional flux expansion, Eq. (19), with N = 1 in Eq. (47) and performing all angular integrations yields
+[’p (0)+~m16J(O) [& (0)+;m16540) (2”) - ;Jcz”, +[ -s + +;$s 1&o(b(Z”) 2(n
2)
+[2(-y3) +;cp (20) - ; J(s,i] GJ(Z") =0,
where we have defined cpO(z)= p(x), the scalar flux, and cpl(z) = J(x), the current, so as to use the usual diffusion theory notation. Xow, in the use of the variational method one usually assumesthat all the variations on the boundaries of the system are independent and arbitrary. Thus Eq. (48) implies J(O) ~ (o(O)
zz
--
1
2’
J(O) __ = -- 4 40) 9:
(50)
SEVTRON
TRANSPORT
205
THEORY
(52) where we have eliminated use of the relationship
S in favor of a more physical
quantity,
Jin(xo),
by
Prom the discussion earlier in this section, we desire one boundary condition at each outer surface. Equations (49) through (52) represent two boundary conditions at each surface. Further, it is evident that these boundary conditions are contradictory. The reason that this trouble arises is apparently a peculiarity in the application of the variational method to the transport equation. In particular for this problem once one has derived the P-l equations as the Euler-Lagrange equations, the four variations of the trial functions on the boundaries of the system are no longer independent since the P-l solution contains only two constants. Thus we must impose a restriction on the trial functions at each boundary to reflect this fact. At x = 0 we demand, quite generally, that q(O) and J(0) be linearly related according to .1(O) = --rep(O). (54) At x = z. , the general linear form can be taken as cicp(X”) - pd(x,)
(55)
= Ji,(ZO).
r, 01, and p are known from physical considerations to be positive. (49) through (52) also suggest that these quantities are positive. (54) and (55) in Eq. (48) to eliminate all the current terms, and coefficients of (o(O)&(O), &(x0), and cp(zO)d~(xO) equal to zero yields tions for CY,0, and r : ? 6 - ?fsI” ?$j - K(y (a/p)’ 1 1 -~ + 2)P - a(n + 2j 4(n The solution of Eqs. (-56) through
8
= 0,
(56)
= 0,
(57)
= 0.
(58)
(58) is
r = d2/3, a = 2/2 p/3,
@= (n+3)
Equations Using Eqs. setting the three equa-
(59) (60)
32/2+4 01+3)+@j(,?%+2)
1 -
(f-31)
206 For isotropic ing to
POMRANI
scattering,
h-c;
E’-1 theory relates the current
to the scalar flux accord-
J(Z) = - L 4&J 3c T’ Thus the linear extrapolation
distance at .Z = Ois,fromKqs.
d = [,b,!$g’ where the transport
(54), (59),and
= 0.7071 At,,
mean free path, Xtr , is defined for isotropic At, = l/Z.
(62), ( 63 )
scattering
as (,ti4)
We note that the linear extrapolation distance is independent of both z. , the size of the system, and c, the mean number of secondaries per collision. This is due to the crudeness of diffusion theory. Comparison of Eq. (63) with the P-l Marshak result of 0.6667 At, and the exact Milne value (for a pure scatterer) of 0.7104 At, (9) shows that the variational result is significantly more accurate than the ,\Iarshak result. The P-l R’lark value for the linear extrapolation distance is o..i774 x1, . We now turn our attention to the boundary condition at x = x0 , given by l+:qs. (M), (60), and (61). We compare this with the Xlarshak result, which is
dd -
4
J(xo) - __ = JLn(Z”), 2
(65)
i.e., Eq. (35) with cy = 4; and p = 3.6. We note that the variational values for (Y and p contain a dependence on n, the index of the entrant angular distribution, whereas the Marshak values do not. Since any exact solution of the transport equation will depend upon the entrant angular distribution, one might espect the variational boundary condition to be superior to the Marshak condition. Table I lists some numerical results for the parameters 01and /3 as a function of 11. To obtain an explicit indication of the accuracy of the variational boundary condition, a diffusion theory solution for P, the capture probability of a semiinfinite slab with a pn entrant directional flux, was obtained. Using Eq. (55) as the boundary condition, the solution was found to be p = &3(1
“e;‘:(
1 - c) .
For c = 0.99, use of the Marshak values for ty and p in E:q. (66) yields a value of P equal to 0.2070 for all n, while the use of the variational values for (Y and fi, Eqs. (60) and (61), yields 0.20.58 for n = 0 and 0.2183 for 7~ = 1. Comparison
n (angular
index) .-.
U The
a
.___
.~.~~~
B
0 (isotropic 1 2 3
flus)
0.2500 0.2357 0.23i8 0 ‘I”20 Iti .
0.5303 0.5000 0.4834 0.1730
0) (normal
heam j
0.2011
0.42ti7
Marshak
boundary
condition
uses oL = 0.2500
and p = ().5()(K) for
all 12.
of these result#s with the exact I-alues ( IO), 020.54 for n = 0 and 021cici for 71 = 1, indeed indicates that the variational boundary condition is superior. For smaller values of c, both the nlarshak and variational boundary conditions yield rather poor quantitative results for P, due t’o the inadequacy of diffusion theory for anything but weak absorbers. Qualitatively, however, the variational boundary condition always ( for all c) shows the proper trend, i.e., a greates value of I’ for a larger 1~ ( more normal heanl), while the ILIarshak conditions always yield a result independent, of n. We now consider the extension of the treatment of the outer surface boundary conditions to higher order expansions. WC first reduce EC]. (18), the first variation of the boundary terms, to the equivalent form in terms of the moments for an Nth order expansion. We have seen in t,he P-l case that the treatment of the outer hollndaries separates itself into two independent treatments, one for each boundary. This is also true in the general P-N expansion and thus, for simplicity, we consider the left hand boundary only. The right hand boundary is handled similarly. We assunw that the entrant’ directional flux is represented by a Legendre polynomial series awording to (.67) ITsing Eqs. l.19) and i(i7) in the left hand boundary yields
temls (z = a) of Ii;q. (18)
208
POMRANING
In order to evaluate integral :
the angular integrals
in Eq. ( 68) we consider t’he following
(69) Use of the recurrence
relationship
(n + 1) P,,+1(s) allows us to write
for the Legendre polynomials,
i.e.,
= !2n + 1) .zP,(.Ec) - ?11>,‘-1(2),
(70)
Eq. (69) as
Equation (71) is the sum of two integrals, each of the type given by Eq. (34). Using Eq. (34) in the special case that v and g are integers, we can reduce Eq. (71) to
APn Bz,,, U,, + Ann+1 Bzm+l Vrm (,nL +
+ ,z [(Zw
1
1)
+ 1)(21x + 3)
1
(72)
[A, B,+I + Am+1 Bml,
where
u,, =
(--I) (‘+“~+l’(2m)!(2n)!(2n” (.!,,,!)za(*n+3m+l)(2t)1 - 271- 1)(2w
+ n + 2m’ + m) ~ - 2n + l)(?i + 112+ l)(n + m) ’
(73)
( -l)‘“fm+“(2n + 1)!(2m + l)! TV,,, = (7~!n2!)z2(Zn+Zm+1)(2v~- 2n - 1)(2w - 211.+ l)(n + 112+ l)(% + m + 2) ’ (74) Proper identification of -4,, and B,, in Eq. (72 ) with cpn(a), a~,( a), and A,+ in 1Cq. (68) completes the reduction of the boundary terms to the equivalent form in terms of the moments. Kow, as in the P-1 case, setting the coefficients of 6pn( a) equal to zero yields too many boundary conditions. As before, the reason is, due to the structure of the P-N equations, that the 6~,,( a) are not all independent., and thus we need to restrict the trial functions at the boundaries. Because of the inherent differences in the even and odd order treatments, we consider each separately. We discuss first the odd order treatment, as it is the most straightforward.
Since an odd order I’-N expansion requires (N + I)/2 boundary conditions at each outer surface, we write (N + 1)/2 general linear relationships between
XIX-TROX
the nlonlents
TRANSPORT
209
THEORY
at z = a, i.e., (1 5 II2 5 (N + 1)/Q).
The usual Irarshak
conditions
(7.5)
have this form where
2n+1 AI,,, = 3 -
l s0
&.lr??,( p ) I’, ( /J‘1)
(76)
1
N,
= s0
(77‘)
&&&.&‘1+(p~.
We shall use the variational method to determine the constants in the linear relationships. Appropriately conlbining (linearly) the (N + 1 J/2 equations given by Eq. (75) we can in general write (N-1) 2 a+-m(a) + c Rnm(on(a) + Tm = 0, (1 5 m 5 (N + 1)/2), (78) ,i=o where R,, and T, are composite constants. Equation (78) represents the final forln of the (N + 1)/2 restrictions we place on the nloments at the outer boundary. Substitution of Eq. (78) into the boundary terms and setting to zero the coefficients of the renlaining variations will yield (N + 1) (N + 3)/4 equations for the sanle number of unknowns, R,, and T, . EVEN
ORDER
Since an even order I’-N expansion requires N/2 boundary conditions at each outer surface, we write N/2 general linear relationships between the rnonlents at, ,7 = a, i.e., 2 ~~fnm~n(a)= Nm , ?I=”
(1 5 m 5 N/2).
(79)
Sow, nunlbering the P-N equations given by Eq. (22) from zero to N, we find front the even nunlbered equations,
where &,, and er,, depend upon the various fz,, , 8, and c. Thus for an even order expansion the equations themselves place an additional restriction on the MOrnrnts. CTsingEq. (80) in Eq. (79) yields
where ML, and Nk represent cotnpositc constants. We see that Eq. (SO), derived from the P-N equations themselves, has allowed us to reduce the sum in 1Xq. (79) by one term. Appropriately combining (linearly) the N/2 equatiotts represented by ISq. (81) me can in getteral write (N-2) :2 (1 4 )lz 5 ND), (82) (PN-m(u) + c n,&,‘(a) + ?‘“, = 0, ,i=0
where R,,, aud 1’, are composite constants. 1Squations ( 80) and (82) represent the final forttt of the N/2 restrictions we place ou the ntotnettts at the outer boundary. Substitution of Kqs. (80 j and (82) into the boundary terms and setting to zero the coefficients of the remaining variations will yield (N)(N + a)/-1 equations for the saute number of unknowns, R,,, attd T, . As an example of the use of these outer bouttdary considerations, we consider a free surface (no entrattt neutrons) at the left hand boundary, taken to be z = 0, of a sourceless system with isotropic scat’tering in the P-2 and P-3 approximations. We first consider the I’-2 approximation which, for this problem, Las the follomitig equations : (y
(83j
+ Z(1 - c)cpo(,x) = 0,
(81) 2 d,,(x) sdi The appropriate
boundary
term, Eq. (68) xit’h N = 2 and A’(p) 1 TO
9 +
$1
+
3
1
,cp:! ->
cpz(O) = ,“$(I
0, is
&l
(86)
where the argutnettt of t,he tnotttettts, z = 0, has been ontitted P’rotu Eqs. (83) and (85 ) we deduce
as the specialization oiily one equation
(8.5)
+ Z$@(z) = 0.
for simplicity.
- c) cpo(O)
of Eq. (80) for this problem.
For N = 2, Eq. (82) yields
VI (0) + Go (0) + T = 0, where we have suppressed
the subscripts
(87)
on R and T for sitnplicit,y.
(88) Using F:qs.
NEUTRON
TRANSPORT
211
THEORY
(87‘) and (88) toelinlinatecpl(0),6~l(0),cp2(0),and6~2(0) inEq. (86) andsetting the resultiug coefficients of 8p0(0) and cp0(0)6(0~(0) equal to zero yields two equations for R and T, i.e., ,j&T - B/isRT
- ,?:i(l - c)T = 0,
;g [I + (1 - c) + (1 - c)‘]
(89)
= “;&.
(90)
EIquation (8:)) yields 7’ = 0. This could have been anticipated since with 110 entrant flux one wollld expect no inhomogeneous term in ~:q. (82). Equatiou i 90 ) yields [l
+
(1
-
C)
+
(1
-
d’Il’.“~
(91)
Physical considerations (,the current is negative at 2 = 0) required taking positive root of Eq. (90 ) to ohtain Eq. (911. Xow, for isotropic scattering, P-2 equations, ( 83) through (85), yield
4m(O) 1
&j = $-II i + i ii - C) o’z. Thus the linear extrapolation (92),
d =
1
d&awe
at 2 = 0 is, fro111 Eqs. (88))
1 [
dp”(O) -I =
[ mclz
0.7071il + fs(1 - C)) (I + (1 - cj + (,I - c)s)liz
the the
(92) (91)) and
(,93 )
1
xi,.
We note that the 1’.2 expression for d depends upon the composition of the lnedium, whereas the I’-1 result did not. This inlprovement in detail is due to the addit,ional complexity of the P-2 trial function (p’ term) over that of the P-l trial function. However, we still have not folmd a dependence on the size of t,he system, indicating that, at least p3 detail (1’.3 expansion) is required for this. We uow consider t’he 1’.3 variat,ional boundary conditions at a free surface. The appropriate boundary term, IQ. (68) with IV = 3 aud =1+(p) = 0, is
where the argument, of the monlcuts, x = 0, has been omitt,ed for For N =z 3, RI. (78‘) yields two equations $n(O) + &%(O~ + &l(O) do)
+ PdO)
silllplicity.
= 0,
(95)
+ JJp’lO‘) = 0,
(96)
222
POMRANlNC;
where we have set Rol = A, RI1 = R, I& = C, and R1?= D for simplicity in notation. Further, we have set Y1 = Ii2 = 0 since the problem is that of a free surface. [If T, and T, are retained in Eqs. (9.5) and (S(j), as they could be, the analysis would yield T1 = T2 = 0.) Using Eys. (93) and (96) in Eq. (94) to eliminate (02, y3 , 8~~ , aud &Q and setting to zero the coefficients of ~&Q , ~,,s~~ , (prS(p,, and +&Q yields four nonlinear equations for the four unknowns, A through
D: Js’ - >1/fy2 + ‘g& fi
-
,5+Jj
-
>i(f
+
25’32c’D
+
3@C
+
L&4
-
- rA7izsL4”
= 0,
“iAD
- 1'7fz8AB = 0, -$i + >f&'+ :,&A - y&D + zp&CD + BiAD - '?JBC - 13T{2sAB = 0, - %is + pi6B + 25&D' - 1"&8B2 = 0. By appropriately manipulating Eqs. (97) through ( 100 j, it can be shown they are equivalent
(97)
(98) (99) (100) that
to 16 - 40C + 100C2 - 1478’ = 0,
(101)
2OCD - 70AC + 49A - 16D = 0,
(102)
3BC - 3AD = 0, -72 + ;iBB + 100D2 - 147B' = 0.
(103)
I-2C+
(104)
Equations (101 j through ( 104) are in a convenient form for a numerical solution. If one assumes a value for C, Eq. (101) can be solved for A, then Eq. (102) cau be solved for D, and finally Eq. (103) can be solved for B. Equation (104) is then used as a cousistency check. The results of this numerical solution are:
A = -0.532591, c = +0.744949,
B = -0.868925, D = +1..522003.
(105)
A second solution of Eqs. (101) to (104) exists (corresponding to a positive value of A), but is rejected on physical grounds. For comparison, we list the Alarshak values, which are: A = - /es,/ 8,
B= -1,
(7 = +g,
D = +3$
(106)
We note that the variational values, Eq. (105), follow the same trend as do the Marshak values, Eq. (106)) but are c.onsistently smaller (in absolute value). It, is not possible from the boundary conditions alone, Eqs. (95) and (96), to
NEUTRON
TRANSPORT
213
THEORY
compute the linear extrapolation distance. One must solve a specific problem using these boundary conditions and, from this solution, compute the linear extrapolation dist’ance. Thus the l’-3 approximation gives both a composition and size dependence to the linear extrapolation distance. So as to compare the accuracy of the Rlarshak and variational P-3 free surface boundary conditions, we consider the RIilne problem, i.e., a pure isotropic scatterer occupying the halfspace x > 0 with a source at infinity. The P-3 equations for this problem are 107)
The general solution of Eqs. (107) through at infinity is
(110) which has the proper behavior
cpo(z) = CUB.2+ /3 + ye?,
(111)
q(z)
(112)
= -;$Y,
p2(x) = +~pz,
(113)
where v = d:<:i/3 and o(, 6, and y are constants to be determined from the boundary conditions. rsing l~~qs. (85) aud (96) as the boundary conditions at z = 0 yields (115)
B= 3vCD
- 7RC’ ldC(SD
l-lC(Dil - RC)
- RC)
y.
(116)
Sow, the linear extrapolation distance, d, is defined as the reciprocal logarithmic derivative of the asymptotic(in Ohis case linear) portion scalar flux, cpU, at the free surface. Thus we have 1 3vCID - 7BC - lX(DA - BC) x cl = P ; Xt, = yjyYc :W
-
'74
1tr.
of the of the (117)
214
POMRANING
TABLE THE
LINEAR
EXTRAPOLATION
DISTAXCE THE
P-l P-2 P-3 P-4 P-5 P-6 P-7 P-8
&kLNE
II (IN
ITNITS
OK >~G.N
FREE
PATHS)
~011
kkoBLEM'
Mark
Marshal
0.5774
0.6667 o.mx7 0.7051 0. i02 0.7082
0.7i-16
n.ci940 0.7297 0. i039 0.7198 0.7069 0.7159
Variational 0.7071 0.7071 0.7118
0 The exact result is 0.7104. Using the Marshak values for 9 through D, Eq. ( 106)) in Eq. ( 117) yields the familiar result d = 0.7051 Xt, . This value, when compared to the exact WienerHopf result of 0.7104 Xt, , is seen to be 0.75% in error. Using the variational values for A through D, Eq. (105), in Eq. (117) yields rl = 0.7118 Xt, , which is 0.20% in error. Thus the use of the variational boundary conditions has reduced the error by almost a factor of four. For completeness we mention the P-3 Mark result for d which is 0.6940 Xt, , 3.31% in error. Table II summarizes the linear extrapolation distance results found in the literature for the Xlne problem as computed by various approximations. It is evident that the P-3 variational result is more accurate than higher order P-N results using either the 1\Iark or lIarshak boundary conditions. VI.
CONCLUSIONS
It has been shown that nlonoenergetic transport theory can properly be characterized by a Lagrangian, and that the use of a truncated spherical harmonics series (Legendre polynomials in slab geometry) as a trial function leads to the widely used P-N approximation. Further, the required boundary conditions for the P-N differential equations are the result of the logical application of the variational method. In particular, at an interface between two media one obtains the boundary conditions by considering an interface as the limiting case of a medium with continuously varying (in space) properties. Using this procedure one obtains rigorously the same interface boundary conditions currently in use, but which were previously derived using rather involved (and somewhat arbitrary) physical considerations. At an outer surface, the boundary conditions are derived by retaining the proper boundary terms in the functional. This leads to boundary conditions which are significantly more accurate, but no more complex, than those currently ill use. In short, the variational calculus has been
SETTROX
TRANSPORT
used t’o give a consistent and logical derivation of all aspects of the tions, and has resulted ill a more accurate P-N approximation. RE(~EIVE:I):
October
31.’
THEORY
P-N
eyua-
14. 196.1 REFERENCES
J.B.SYKES, “Seutron Transport Theory,” pp. 116-156. The Clarendon Press, Oxford, 1957. i. J. C. MARK, The Spherical Harmonic Method, I, CRT-340. Atomic Energy of Canada, Ltd., Ontario, 1947. 3. R. MAISHAS, Phgs. Rev. 71, 443 (1947). 4. G. C. PUMILANING AND bl. CLARK. JK., Sucl. Sci. E.ng.16, 147 (1963). 5. 1.1. S. SELENGTT, Hanford Quart. Rept. HW-59126, pp. 899124. Richland, Wash., 1959. 6. G. SANSONE, “Orthogonal Polynomials,” p. 172. Int’erscience, New Tork, 1959. Y. G. Y. Rl-MI-ANTSEV, J. .VUC~. Ener.gy, Ports A/U, 16, 111 (1902). 8. H. BATEMAN, “Higher Transcendental Functions,” Vol. II, p. 173. Compiled by the staff of the Bnteman Manuscript Project, A. Frdelyi, Director. McGraw-Hill, New Tork, 1953. 9. K.&I. CASE, F. DE HOFFMAS, ASU G. PLACZEB, Introductiontothe Theory of Neutron Diffusion, p. 136. Los Alamos Scientific Laboratory, Los illamos, Sew Mexico, 1953. 10. D. P. SELENGUT, private communication, July, 1963. 1. B.~)AVISON.~XD
.iiX’bALS
27,
193-215
Variational Boundary Approximation
(1964)
Conditions for the Spherical Harmonics to the Neutron Transport Equation G.
General
Electt,ic
C’wnpccny,
C.
POMRANING
Vallecitos
9tomic
Labomtory,
Plrasanton,
C’alifornia
11) is shown that the monoenergetic neutron transport equation and the associated boundary conditions can be characterized by a Lagrangian. A proper choice of the trial function for this Lagrangian leads to the widely used spherical harmonics approximation as the Euler-Lagrange equations. Further, a set of boundary conditions for the spherical harmonic equations is the result of the logical application of the variational method. These variational boundary conditions appear to be significantly more accurate than the boundary conditions presently in general use. For example, the use of the variational boundary conditions at a free surface reduces the error (compared with the boundary conditions presently used) in the linear extrapolation distance for the Milne problem by several factors. In particular, the P-l (diffusion) approximation yields a value of 0.7071 (in units of mean free paths) and the P-3 approximation yields a value of 0.7118, both comparing quite favorably with t.he exact value of 0.7104. I. INTRODUCTION A commo~l method of solving the monoenergetic neutron transport equatiou is to expand the angular dependence of the directional flux in a truncated spherical harmonics series (1). This reduces the integro-differential equation to a finite number of coupled differential equatious for t,he spatially dependent expansion coefficients. While the boundary conditions for the transport problem are well set, these conditions cannot always be satisfied exactly by a truncated spherical harmonics expansion. In particular, at an interface between two media the proper transport boundary condition is that the directional flux be continuous iu all directions. It is not, in general, possible to satisfy this boundary condition in a truncated spherical hamonies expansion of the directional flux, and thus approximate boundary conditions may be needed. Further, at the outside houndaries of the system, which shall be referred to as outer surfaces, the general transport
boundary
incoming
hemisphere
condition
condition
exactly
in
is that
the
directional
flux
of solid angle. Again one cannot a truncated
spherical
harnlonics
mation is required. 193
be
satisfy
expansion
specified
in
the
this boundary and
an
approxi-
194
POMRANING
Two sets of boundary conditions at an outer surface for the spherical harmonics equations, due to Mark (2) and Marshak (3), are presently in general use. A detailed discussion of these boundary conditions can be found in ref. 1, and a brief discussion will be given in this paper at the appropriate point. At an interface between two media, the boundary conditions in general use are due to Davison (1) which he derives on physical grounds. The main point to be made here is that all of these boundary conditions are quite arbitrary, and a large number (also arbitrary) of other boundary conditions could equahy well be proposed. Previous authors (4, 5) have shown that the neutron transport equation can be properly characterized by a Lagrangian, and that the spherical harmonics approximation can be derived as the Euler-Lagrange equations of this Lagrangian by the proper choice of the trial function. In this paper we show that the variational method can also logically be used to derive a set of boundary couditions for the spherical harmonics equations. For the outer surface boundary conditions, this involves leaving the proper boundary terms in the Lagrangian. For the conditions at an interface, one cousiders the interface to be a limiting case of a medium whose properties are a continuous function of space. The a priori advantage of these bouudary conditions over those preseutly in use are that they do riot suffer from any arbitrariness. The only assumption required in their derivatiou is implicit in the equations themselves; i.e., the directional flux is assumed expandable in a truncated spherical harmonics series. We emphasize that the use of a Lagrangiau to describe neutron transport is not in any way au assumption, since reuderiug this Lagrangian statiouary is entirely equivalent to solving the neutron transport equation. Thus, if one prefers, the Lagrangian cau be taken as the fundamental equation of ueutrou transport theory. The a posteriori advant#age of the variatioual bouudary conditions is that they are 110 more complex t,han those presently in use, but appear to be significantly more accurate. For algebraic simplicity the formalism in this paper is restricted to slab geometry. Presumably the same formalism could be carried through in other geometries. The spherical harmonics reduce to Legendre polynomials in slab geometry, aud we shall refer to the flux expansiou carrying the first (N + 1) Legendre polynoniials (PO to PN) as the P-N approximation. Sections II aud III of the paper discuss briefly the general variational formulation of neutron transport theory and give the variational derivatiou of the P-N equations. Sectious IV and V discuss in some detail the forniulatiou of the variational bouudary conditions and compare the results with the boundary conditions currently in use. In particular, emphasis is giveu to the outer surface boundary conditions since the variational method yields results quite different from and
NEUTRON
significantly in use.
TRANSPORT
195
THEORY
more accurate than those of Mark and Marshak II.
The monoenergetic H(x, PM&
VARIATIONAL
transport
equation
/-d = m% co,
which are currently
FORMALISM
in slab geometry (-1
5 /.l 5 1;
can be written
(1)
a 5 x 5 b),
(1)
where
and
.fn= I: 4Jof(Po)P?z(Po).
(3)
The spatial region of interest has quite generally been assumed to be a s z 5 b. Here cp(z, +L) is the directional flux, z is the slab coordinate, P is the cosine of the angle between the velocity vector of the neutron and the x-axis, c is the mean number of secondaries per collision, Z is the macroscopic collision cross section, f( pO) dpO is the probability that a scattering event will scatter a neutron through an angle ~0 into C&Q,, S(z, p) is the external directjonal source, and P%(p) is the 71th Legendre polynomial. S, c, and fn have an arbitrary spatial dependence, although for simplicity in notation, we do not explicitly show this z-dependence. The general boundary conditions on Eq. (1) at the outer surfaces are specifications of the incoming flux, i.e., 4% /J) = ~+iP),
(0 < p 5 I),
(3:)
54, P) = A-(P),
(-1
(-5)
5 /.I < O),
where ,4+(p) and A-(P) are arbitrary, but known, functions. We now define the adjoint operator according t.o / dr+~(.r)H*(.r)+~*(~j
= [ d.rc~*(s)H(z)cp(lc)
+ (boundary
terms),
(6’)
where z represents the phase space variables z and p and the integration extends over the region of interest in r-space. Requiring that Eq. (Gj hold for any two functions &z) and (p*(z) yields
196
POMRANING
Thus the adjoint H*k,
equation is
P)cp*(.% PI = uz,
(-1
PI,
5 p 5 1;
a 5 z 5 b),
(8)
where ‘p*(z, p) is the adjoint directional flux, and T(z, p) is the adjoint directional source which, at this point, is arbitrary. The adjoint boundary conditions, also arbitrary at this point, are taken as cp*ca, IL) = B+(PL
(-1
cp*u4 Pu) = B-(P),
(0 -=c P 6 11,
where B+(P) and B-(p) are arbitrary, but specified, functions. evident from comparison of Eqs. (2) and (7) that
5 p < O),
source and boundary ux,
PI = St% -PI,
(11) properly.
That
is, if (13)
B+(P)
= A+(-d,
(13)
B-(P)
= A-(-EL),
(l-1)
then Eq. (11) is satisfied. We now postulate that a proper transport equations, (1) and (8)) Eqs. (4), (5), (91, and (101, is Fb, (P*I = s,’ c&zI: hw*(%
conditions
(10)
It is immediately
4% PI = cp*c% -PI if we choose the adjoint we set
19)
Lagrangian for the transport and adjoint and the associated boundary conditions,
Pcr)N(% cl)&,
P>
For Eq. (1.5) to be a correct Lagrangian, setting the first variation with respect to arbitrary and independent variations of q(x) and (p*(x) equal to zero must be equivalent to solving the transport and adjoint transport equations. It is
NEUTROK
TRANSPORT
197
THEORY
easily verified that setting the first variation of Eq. (15) with respect to (O*(X) equal to zero yields Eq. (1) as the Euler-Lagrange equation with Eqs. (4) and (5) as the subsidiary or bourldary conditions. Similarly, setting the first variation of Eq. ( 1.3) with respect to P(X) equal to zero and retaining the boundary terms from the integration by parts yields Eq. (8) as the Euler-Lagrange equation with Eqs. (9) and (10) as the boundary conditions. Thus we are assured that Eq. (15) is, in fact, a proper Lagrangian. We now assume that in any subsequent use of the Lagrangian, the set of t,rial functions for cp(.r) and P*(X) is restricted according to Eq. (11). This allows us t,o simplify the Lagrangian. In particular, the first variation of the Lagrangian can then be writ’ten, using Eqs. ( 11) through (14)) as
h-ow, since the variations in the interior of x-space are independent from those on the boundaries for any set of trial functions, setting 6F = 0 in Eq. (1Gj implies t,wo separate equations, i.e.,
ss b
a
1
dx
-1
d/.M(~, P)P(X, Co - X(x, EL)IBp(x, -p)
= 0,
(17)
Equations (17) and (18) are the starting expressionsfrom which we shall proceed. In particular, Eq. (17)) with the proper choice of the trial function, will yield the P-N equations and the boundary conditions on these equations at an interface between two media. Equation (18), with the same trial function, will yield the boundary conditions for the P-N equations at an outer surface. III.
DERIVATION
OF THE
P-iv EQUATIONS
We take as our trial function a truiicat,ed Legendre polynomial seriesin angle, i.e.,
198 Substituting [
POMRAKING
Eq. (19) into Eq. (18) yields
dz I: d/J [n(%
P)Y$
(“i”)
(Pm(Z)Prn(PcL) - S(% Pi] (20) ( -l)nPn(Ph%(Z)
.lgq) Now, a priori all of the &D~(x) are arbitrary implies (N + 1)Equations:
and independent.
= 0. Thus Eq. (20)
(0 5 n 5 N). We recognize Eq. (21) as the widely used P-N equations. That is, the classical derivation of the P-N equations is to expand the angular dependence of the directional flux in a truncated Legendre polynomial series and use orthogonality to obtain the differential equations for the spatially dependent. expansion coefficients. This is precisely the content of Eq. (21). Two points should be made here. First, it is not necessary to choose Legendre polynomials as the trial functions in order for the variational method to yield the P-N equations. Any set of polynomials, or even the simple integer powers of I*, would yield the the same result since the ingredients of the trial functions are the same in all cases, i.e., linear combinations of the integer powers of p. The Legendre polynomials were chosen so that the result, Eq. (21)) would be readily recognizable. Secondly, while a priori all of the &%(z) are independent, a posteriori this may or may not be the case depending upon the order of the expansion. As is well known (1) , the solution of Eq. (21), which contains (N + 1) unknowns, for cpn(z) contains (N + 1) constants if N is odd and N constants if N is even. Thus for N even, prior to applying the boundary conditions, one of the cpn(x) is linearly dependent upon the other N p,(x), and therefore a posteriori the Bv~(z) are not all independent. Let us explicitly write out the P-N equations as given by Eq. (21). Performing the indicated angular integrations, we have
+
&u(z) dz
Zcpl(X)
+ %h42)
=
Xl(X)
+
Wlcpl~~),
= A%(z) + =fNPN(Z),
NEVTRON
TRANSPORT
199
THEORY
where X,,(Z) is the ?Ith Legendre moment of the source, i.e., (33) IV.
BOUSI~lARY
CONDITIONS
AT
.4N
INTERFACE
If the medium properties have a discontinuity at some point within the system the usual method of solution of the P-N equations is to find the solution in the regions on either side of this discontinuity and join these solutions by prescribed boundary conditions at the interface (discontinuity). Since adding a region introduces another interface, the number of conditions that can (must) be satisfied at an int’erface is equal to the number of constant8s in the solution within a region. We know (1) that within a region the solution of the P-N equations contains (N + 1) constants if N is odd and N constants if N is even. For example, assuming a solution of the form (L&(z) = g,P
(24)
within a homogeneous region and setting the coefficient determinant of the g,,‘s equal to zero yields (N + 1) values for v if N is odd and N values for v if N is even. Thus we must construct (,N + 1) boundary conditions for N odd and N boundary conditions for N even. Kow, in deriving the P-N equations no assumptions were made concerning the spatial dependence of the mediunl properties, and thus the equations hold for a medium whose properties vary arbitrarily with Z. In particular, the P-N equations are valid at an interface and one should be able to determine the boundary conditions from the equations themselves. Because of the inherent differences in the even and odd order approximations (N even or odd ), we consider the two cases separately, taking the interface to be z = _Y(~ . ODD ORDER
We number the (N + 1) equations given by l?q. (22) from zero to N. Integrating the zeroth equation over x from (z c - t ) to (zO + E) and letting E approach zero, noting that in this limit only the derivative term contributes, yields the fact that cpl(x) is continuous at the interface, x = pg. A similar integration over the second equation, using the fact that cpl(zO) is continuous, yields the fact that (p3(x) is continuous at the interface. Continuing this procedure of integrating over the even numbered equations yields the fact that all the od’dmoments are continuolls at an interfare. The conditions on the even moments are found by integrating over the odd numbered equations. In particular, integration over the Nth ecluation yields the fact that (~~~~(2) is continuous at the interface. Integration over the (N - 2)th equation, using the fart that (P~-~(,Q)is continuous, yields the fact that (P~-~(x) is continuous at the interface. Continuing this pro-
200
POMRANING
cedure, we deduce that all the eoen moments are continuous at an interface. Thus we conclude that for an odd order expansion all moments are continuous at an interface. These are the boundary conditions universally used for odd order expansions, but the usual argument (1) to obtain these conditions is that the directional flux must be continuous in all directions at an interface, and thus all of the moments, cp,(x,,), must be continuous. EVEN
ORDER
As is well known (I), the argument that all moments must be continuous at an interface leads to one too many boundary conditions for an even order expansion. It is here that the procedure of obtaining the interface boundary conditions from the equations themselves is very useful and is, in fact, the only consistent treatment. As in the odd order analysis, subsequent integration of the zeroth, second, fourth, etc. equations leads to continuity of the odd moments at an interface. Integration of the first, third, fifth, etc. equations yields the fact that
cpzn(zo) + ( g$ >w&f&“),
(0 5 n I (N - 2)/2)
(25)
is continuous. By appropriately combining these continuity conditions, it is easily shown that Eq. (25) is equivalent to demanding continuity of
(0 5 n 5 (N - 2),'2).
(26)
Now, according to Sansone (B), a property of Legendre polynomials for n an integer is
Pz,(O) =(-1)”o1(s >.
(27)
Using Eq. (27) in Eq. (26) yields
(ozn(zo)- p2,(o)Po(~o),
(1 5 n 5 N/2)
(28)
is continuous. Noting that PZn+l(0) is zero for n an integer, we can demand continuity of the odd moments and continuity of Eq. (28) for the even moments by requiring con
- pn(o)cpo(~o),
(1 5 n 5 N)
(29)
be continuous. Equation (29) represents the N boundary conditions at an interface for an even order P-N expansion. We note that the scalar flux, cpo(z),is not continuous at an interface. As remarked by Davison (I), one should not be surprised to find a discontinuity in the scalar flux in an approximate solution of the transport equation.
NEUTRON
TRANSPORT
THEORY
201
At this point it is iuteresting to compare Eq. (29) with the even order interface boundary conditions proposed by other authors. Davison (1) has suggested the identical conditions as found in Eq. (29). His derivation, however, is conpletely different from the one giveu in this paper and is based on rather iuvolved physical arguments. The work by Rumyantsev (7), when specialized to slab geometry, also yields Eq. (29). His considerations were much the same as those in this paper, but were not in the context of the variational method. Finally, following the spirit of Marehak’s conditions at an outer boundary (discussed in ref. 1 and in t’he next section of this paper), the Ilarshak conditions at an interface are logically taken as continuity of s0
’ &JP*n+l(EL)pF(Xo, p),
t.0 5 11s (N - 2)/2),
(30)
(0 5 ‘n 5 (N - 2)/2).
(31)
and 0 s-1
dPPzn+lbFL)&O,
PFC),
That is, Marshak demands continuity of the odd Legendre half-moments of the directional flux. Setting n = 0 in Eqs. (30) and (31), we see Marshak’s conditions have the advantage of including the physical condition that the partial currents across an interface are continuous. The use of the expansion for the directional flux, Eq. (19)) in Eqs. (30) and (31) yields N conditions on the moments, p,(z), at the interface. Subsequent addition of Eqs. (30) and (31) and use of the orthogonality and parity properties of the Legendre polynomials yields the fact that rnn+1(~0),
(0 5 n 5 (N - 2)/2)
(32)
is continuous. Subtraction of Eq. (31) from IQ. (30) and use of the parity properties of the Legendre polynomials yields continuity of 1 N/2 dpP2,,+l(r)m~o (h + 1)P2m(P)(P2m(Z”), (0 5 n s (N - 2)/2). (33) s0
According to Bateman (8))
s0
l d.xP,(xP,(z)
= ; [(CT- Y)(U + lJ + l)]-’ .[A sin (UT/~) cos (VT/~) - .4-l sin (m/2) Cos(m/2)1,
where A4 = r(95 + v/2)r(l qj_; + u/a)r(l
+ u/2) + V/2) ’
(35)
202
POMRANING
and I’(X) is the gamma function,
defined as (36)
Setting v = 2n + 1 and c = 2~1 (both nl and n are integers) properties of the gamma function, we can simplify Eq. (34) to 1 s0
dxPz,+~(x)Pzm(x)
=
(-1y+m+l (2m)!(2n (71!l)b!)“2(2m+sn+1)(2,m - 2n -
Using Eq. (37) in Eq. (33), we demand continuity ji!
and using the
+ l)! l)(nk + n + 1)’
(37)
of
( -1)1)2(4m + 1)(2m)! [ (m!)222~(2m - 2n - l)(nz + n + 1)
1(0*&o),
(38)
(0 5 n 5 (N - 2>/2). Noting that Eq. (27) can be written P2n(0)
as =
(-1)“(2n)! yn(n!)2
’
Eq. (38) is identical to requiring continuity of g
1
(4m + l)PdO) (0 5 n 4 (N - 2)/2). [(2m - 2n - ~)(HL + n + 1 (P2m(Xo)7
(40)
Equations (32) and (40) are the N Marshak conditions at an interface. A comparison of these conditions with the variational conditions, Eq. (29)) shows that the two are not the same. For example, in the P-2 approximation the variational conditions, Eq. (29)) reduce to cpl(x0), (o~(20) + >~c,cJ,,( zo) being continuous, while the Marshak conditions, Eqs. (32) and (40)) reduce to (Ok, (p2(zo) + +&,(z,) being continuous. Since the variational boundary conditions (the same as those of Davison and Rumyantsev) are consistent with the P-N equations, whereas the Marshak conditions are not, it is suggested that the variational conditions, Eq. (29)) should be used. V. BOUlVDARY
CONDITIONS
AT
A?;
OUTER
SURFACE
We now consider the boundary conditions at an outer surface. It is here that the variational method yields entirely new results which are no more complex to apply than the boundary conditions of I\Iark and Xarshak, currently in use, but appear to be significantly more accurate. Before proceeding, however, let us briefly discuss the Mark and Marshak boundary conditions. To be specific, we consider the left hand boundary of a system, taken as x = a. The right hand boundary is handled similarly. We know from the number of constants in the
NEUTROS
TRANSPORT
203
THEORY
solution of the P-N equations that we require (N + 1)/2 outer boundary conditions at z = a for N odd and N/2 conditions for N even. At a free surface (no neutrons entering the system) the rigorous transport boundary condition is da, ~1 = 0,
(0 5 /.Lf
1).
(41)
It is obvious that Eq. (41) cannot be satisfied exactly by a truncated Legendre polynomial expansion of the directional flux. As an approximation, Mark (2) has suggestedsetting (43)
da, h) = 0,
where the /*i are the R,/2 positive roots of Pn(pi) = 0, and R = N + 1 for N odd and R = N for N even. Using the expansion for t#hedirectional flux, Eq. (19)) in Eq. (42) yields R/2 linear relationships, which is the required number of boundary conditions, on the moments at the free surface. Davison (1) has proved that Eq. (42) is equivalent to replacing the vacuum to the left of the free surface with a pure absorber and treating the resulting interface according to his boundary conditions discussedin the last section. Marshak (3) has suggested
(0 -I n 5 (R - 2)/2),
(4.3)
to approximate Eq. (41). Setting 71 = 0 in Ii:q. (43) shows that Rlarshak’s boundary conditions include the physical condition that the current into the system is zero. As with the l\lark conditions, the llse of Eq. (19) in Eq. (43)
The use of the directional flux expansion, Eq. (19), in E:q. (44) allows one to carry out all of the angular integrations, since 4+(p) is presumed known. This leads to R/Y inhomogeneouslinear relationships between the moments at z = a. Based on the above discussion,it is evident that the Marshak homidary condi-
204
POMRANING
tions at au outer surface enjoy much more usage than do the Mark conditions. Accordingly, the emphasis here will be on the comparison of the accuracy of the variational boundary conditions with that of the Marshak conditions. AIark’s conditions will be mentioned only in an incidental mauner. Before proceeding with the variational analysis for a general N, we consider the special case of N = 1. This is a particularly important case since it corresponds to the widely used diffusion approximation. Further, the consideration of this special case will serve to point the way toward the more general analysis. We consider, within the framework of diffusion theory, the quite general problem of a free surface at the left hand face, x = 0, and a ( -p)“ entrant flux distribution at the right hand face, x = zO. Thus we set A+(P)
= 0,
(45)
A-(/l)
= S(-/.4Ln,
(46)
where S is a constant indicating the magnitude of the entrant flux and n is an arbitrary nonnegative integer. The first variation of the boundary terms of the functional, Eq. (18), for this problem is 1’ &wP(O, co&40, -PI 0
+ s: &P[S( -PI”
- (0(X”, P)l&h”,
-pj
= 0.
(47)
Using the directional flux expansion, Eq. (19), with N = 1 in Eq. (47) and performing all angular integrations yields
+[’p (0)+~m16J(O) [& (0)+;m16540) (2”) - ;Jcz”, +[ -s + +;$s 1&o(b(Z”) 2(n
2)
+[2(-y3) +;cp (20) - ; J(s,i] GJ(Z") =0,
where we have defined cpO(z)= p(x), the scalar flux, and cpl(z) = J(x), the current, so as to use the usual diffusion theory notation. Xow, in the use of the variational method one usually assumesthat all the variations on the boundaries of the system are independent and arbitrary. Thus Eq. (48) implies J(O) ~ (o(O)
zz
--
1
2’
J(O) __ = -- 4 40) 9:
(50)
SEVTRON
TRANSPORT
205
THEORY
(52) where we have eliminated use of the relationship
S in favor of a more physical
quantity,
Jin(xo),
by
Prom the discussion earlier in this section, we desire one boundary condition at each outer surface. Equations (49) through (52) represent two boundary conditions at each surface. Further, it is evident that these boundary conditions are contradictory. The reason that this trouble arises is apparently a peculiarity in the application of the variational method to the transport equation. In particular for this problem once one has derived the P-l equations as the Euler-Lagrange equations, the four variations of the trial functions on the boundaries of the system are no longer independent since the P-l solution contains only two constants. Thus we must impose a restriction on the trial functions at each boundary to reflect this fact. At x = 0 we demand, quite generally, that q(O) and J(0) be linearly related according to .1(O) = --rep(O). (54) At x = z. , the general linear form can be taken as cicp(X”) - pd(x,)
(55)
= Ji,(ZO).
r, 01, and p are known from physical considerations to be positive. (49) through (52) also suggest that these quantities are positive. (54) and (55) in Eq. (48) to eliminate all the current terms, and coefficients of (o(O)&(O), &(x0), and cp(zO)d~(xO) equal to zero yields tions for CY,0, and r : ? 6 - ?fsI” ?$j - K(y (a/p)’ 1 1 -~ + 2)P - a(n + 2j 4(n The solution of Eqs. (-56) through
8
= 0,
(56)
= 0,
(57)
= 0.
(58)
(58) is
r = d2/3, a = 2/2 p/3,
@= (n+3)
Equations Using Eqs. setting the three equa-
(59) (60)
32/2+4 01+3)+@j(,?%+2)
1 -
(f-31)
206 For isotropic ing to
POMRANI
scattering,
h-c;
E’-1 theory relates the current
to the scalar flux accord-
J(Z) = - L 4&J 3c T’ Thus the linear extrapolation
distance at .Z = Ois,fromKqs.
d = [,b,!$g’ where the transport
(54), (59),and
= 0.7071 At,,
mean free path, Xtr , is defined for isotropic At, = l/Z.
(62), ( 63 )
scattering
as (,ti4)
We note that the linear extrapolation distance is independent of both z. , the size of the system, and c, the mean number of secondaries per collision. This is due to the crudeness of diffusion theory. Comparison of Eq. (63) with the P-l Marshak result of 0.6667 At, and the exact Milne value (for a pure scatterer) of 0.7104 At, (9) shows that the variational result is significantly more accurate than the ,\Iarshak result. The P-l R’lark value for the linear extrapolation distance is o..i774 x1, . We now turn our attention to the boundary condition at x = x0 , given by l+:qs. (M), (60), and (61). We compare this with the Xlarshak result, which is
dd -
4
J(xo) - __ = JLn(Z”), 2
(65)
i.e., Eq. (35) with cy = 4; and p = 3.6. We note that the variational values for (Y and p contain a dependence on n, the index of the entrant angular distribution, whereas the Marshak values do not. Since any exact solution of the transport equation will depend upon the entrant angular distribution, one might espect the variational boundary condition to be superior to the Marshak condition. Table I lists some numerical results for the parameters 01and /3 as a function of 11. To obtain an explicit indication of the accuracy of the variational boundary condition, a diffusion theory solution for P, the capture probability of a semiinfinite slab with a pn entrant directional flux, was obtained. Using Eq. (55) as the boundary condition, the solution was found to be p = &3(1
“e;‘:(
1 - c) .
For c = 0.99, use of the Marshak values for ty and p in E:q. (66) yields a value of P equal to 0.2070 for all n, while the use of the variational values for (Y and fi, Eqs. (60) and (61), yields 0.20.58 for n = 0 and 0.2183 for 7~ = 1. Comparison
n (angular
index) .-.
U The
a
.___
.~.~~~
B
0 (isotropic 1 2 3
flus)
0.2500 0.2357 0.23i8 0 ‘I”20 Iti .
0.5303 0.5000 0.4834 0.1730
0) (normal
heam j
0.2011
0.42ti7
Marshak
boundary
condition
uses oL = 0.2500
and p = ().5()(K) for
all 12.
of these result#s with the exact I-alues ( IO), 020.54 for n = 0 and 021cici for 71 = 1, indeed indicates that the variational boundary condition is superior. For smaller values of c, both the nlarshak and variational boundary conditions yield rather poor quantitative results for P, due t’o the inadequacy of diffusion theory for anything but weak absorbers. Qualitatively, however, the variational boundary condition always ( for all c) shows the proper trend, i.e., a greates value of I’ for a larger 1~ ( more normal heanl), while the ILIarshak conditions always yield a result independent, of n. We now consider the extension of the treatment of the outer surface boundary conditions to higher order expansions. WC first reduce EC]. (18), the first variation of the boundary terms, to the equivalent form in terms of the moments for an Nth order expansion. We have seen in t,he P-l case that the treatment of the outer hollndaries separates itself into two independent treatments, one for each boundary. This is also true in the general P-N expansion and thus, for simplicity, we consider the left hand boundary only. The right hand boundary is handled similarly. We assunw that the entrant’ directional flux is represented by a Legendre polynomial series awording to (.67) ITsing Eqs. l.19) and i(i7) in the left hand boundary yields
temls (z = a) of Ii;q. (18)
208
POMRANING
In order to evaluate integral :
the angular integrals
in Eq. ( 68) we consider t’he following
(69) Use of the recurrence
relationship
(n + 1) P,,+1(s) allows us to write
for the Legendre polynomials,
i.e.,
= !2n + 1) .zP,(.Ec) - ?11>,‘-1(2),
(70)
Eq. (69) as
Equation (71) is the sum of two integrals, each of the type given by Eq. (34). Using Eq. (34) in the special case that v and g are integers, we can reduce Eq. (71) to
APn Bz,,, U,, + Ann+1 Bzm+l Vrm (,nL +
+ ,z [(Zw
1
1)
+ 1)(21x + 3)
1
(72)
[A, B,+I + Am+1 Bml,
where
u,, =
(--I) (‘+“~+l’(2m)!(2n)!(2n” (.!,,,!)za(*n+3m+l)(2t)1 - 271- 1)(2w
+ n + 2m’ + m) ~ - 2n + l)(?i + 112+ l)(n + m) ’
(73)
( -l)‘“fm+“(2n + 1)!(2m + l)! TV,,, = (7~!n2!)z2(Zn+Zm+1)(2v~- 2n - 1)(2w - 211.+ l)(n + 112+ l)(% + m + 2) ’ (74) Proper identification of -4,, and B,, in Eq. (72 ) with cpn(a), a~,( a), and A,+ in 1Cq. (68) completes the reduction of the boundary terms to the equivalent form in terms of the moments. Kow, as in the P-1 case, setting the coefficients of 6pn( a) equal to zero yields too many boundary conditions. As before, the reason is, due to the structure of the P-N equations, that the 6~,,( a) are not all independent., and thus we need to restrict the trial functions at the boundaries. Because of the inherent differences in the even and odd order treatments, we consider each separately. We discuss first the odd order treatment, as it is the most straightforward.
Since an odd order I’-N expansion requires (N + I)/2 boundary conditions at each outer surface, we write (N + 1)/2 general linear relationships between
XIX-TROX
the nlonlents
TRANSPORT
209
THEORY
at z = a, i.e., (1 5 II2 5 (N + 1)/Q).
The usual Irarshak
conditions
(7.5)
have this form where
2n+1 AI,,, = 3 -
l s0
&.lr??,( p ) I’, ( /J‘1)
(76)
1
N,
= s0
(77‘)
&&&.&‘1+(p~.
We shall use the variational method to determine the constants in the linear relationships. Appropriately conlbining (linearly) the (N + 1 J/2 equations given by Eq. (75) we can in general write (N-1) 2 a+-m(a) + c Rnm(on(a) + Tm = 0, (1 5 m 5 (N + 1)/2), (78) ,i=o where R,, and T, are composite constants. Equation (78) represents the final forln of the (N + 1)/2 restrictions we place on the nloments at the outer boundary. Substitution of Eq. (78) into the boundary terms and setting to zero the coefficients of the renlaining variations will yield (N + 1) (N + 3)/4 equations for the sanle number of unknowns, R,, and T, . EVEN
ORDER
Since an even order I’-N expansion requires N/2 boundary conditions at each outer surface, we write N/2 general linear relationships between the rnonlents at, ,7 = a, i.e., 2 ~~fnm~n(a)= Nm , ?I=”
(1 5 m 5 N/2).
(79)
Sow, nunlbering the P-N equations given by Eq. (22) from zero to N, we find front the even nunlbered equations,
where &,, and er,, depend upon the various fz,, , 8, and c. Thus for an even order expansion the equations themselves place an additional restriction on the MOrnrnts. CTsingEq. (80) in Eq. (79) yields
where ML, and Nk represent cotnpositc constants. We see that Eq. (SO), derived from the P-N equations themselves, has allowed us to reduce the sum in 1Xq. (79) by one term. Appropriately combining (linearly) the N/2 equatiotts represented by ISq. (81) me can in getteral write (N-2) :2 (1 4 )lz 5 ND), (82) (PN-m(u) + c n,&,‘(a) + ?‘“, = 0, ,i=0
where R,,, aud 1’, are composite constants. 1Squations ( 80) and (82) represent the final forttt of the N/2 restrictions we place ou the ntotnettts at the outer boundary. Substitution of Kqs. (80 j and (82) into the boundary terms and setting to zero the coefficients of the remaining variations will yield (N)(N + a)/-1 equations for the saute number of unknowns, R,,, attd T, . As an example of the use of these outer bouttdary considerations, we consider a free surface (no entrattt neutrons) at the left hand boundary, taken to be z = 0, of a sourceless system with isotropic scat’tering in the P-2 and P-3 approximations. We first consider the I’-2 approximation which, for this problem, Las the follomitig equations : (y
(83j
+ Z(1 - c)cpo(,x) = 0,
(81) 2 d,,(x) sdi The appropriate
boundary
term, Eq. (68) xit’h N = 2 and A’(p) 1 TO
9 +
$1
+
3
1
,cp:! ->
cpz(O) = ,“$(I
0, is
&l
(86)
where the argutnettt of t,he tnotttettts, z = 0, has been ontitted P’rotu Eqs. (83) and (85 ) we deduce
as the specialization oiily one equation
(8.5)
+ Z$@(z) = 0.
for simplicity.
- c) cpo(O)
of Eq. (80) for this problem.
For N = 2, Eq. (82) yields
VI (0) + Go (0) + T = 0, where we have suppressed
the subscripts
(87)
on R and T for sitnplicit,y.
(88) Using F:qs.
NEUTRON
TRANSPORT
211
THEORY
(87‘) and (88) toelinlinatecpl(0),6~l(0),cp2(0),and6~2(0) inEq. (86) andsetting the resultiug coefficients of 8p0(0) and cp0(0)6(0~(0) equal to zero yields two equations for R and T, i.e., ,j&T - B/isRT
- ,?:i(l - c)T = 0,
;g [I + (1 - c) + (1 - c)‘]
(89)
= “;&.
(90)
EIquation (8:)) yields 7’ = 0. This could have been anticipated since with 110 entrant flux one wollld expect no inhomogeneous term in ~:q. (82). Equatiou i 90 ) yields [l
+
(1
-
C)
+
(1
-
d’Il’.“~
(91)
Physical considerations (,the current is negative at 2 = 0) required taking positive root of Eq. (90 ) to ohtain Eq. (911. Xow, for isotropic scattering, P-2 equations, ( 83) through (85), yield
4m(O) 1
&j = $-II i + i ii - C) o’z. Thus the linear extrapolation (92),
d =
1
d&awe
at 2 = 0 is, fro111 Eqs. (88))
1 [
dp”(O) -I =
[ mclz
0.7071il + fs(1 - C)) (I + (1 - cj + (,I - c)s)liz
the the
(92) (91)) and
(,93 )
1
xi,.
We note that the 1’.2 expression for d depends upon the composition of the lnedium, whereas the I’-1 result did not. This inlprovement in detail is due to the addit,ional complexity of the P-2 trial function (p’ term) over that of the P-l trial function. However, we still have not folmd a dependence on the size of t,he system, indicating that, at least p3 detail (1’.3 expansion) is required for this. We uow consider t’he 1’.3 variat,ional boundary conditions at a free surface. The appropriate boundary term, IQ. (68) with IV = 3 aud =1+(p) = 0, is
where the argument, of the monlcuts, x = 0, has been omitt,ed for For N =z 3, RI. (78‘) yields two equations $n(O) + &%(O~ + &l(O) do)
+ PdO)
silllplicity.
= 0,
(95)
+ JJp’lO‘) = 0,
(96)
222
POMRANlNC;
where we have set Rol = A, RI1 = R, I& = C, and R1?= D for simplicity in notation. Further, we have set Y1 = Ii2 = 0 since the problem is that of a free surface. [If T, and T, are retained in Eqs. (9.5) and (S(j), as they could be, the analysis would yield T1 = T2 = 0.) Using Eys. (93) and (96) in Eq. (94) to eliminate (02, y3 , 8~~ , aud &Q and setting to zero the coefficients of ~&Q , ~,,s~~ , (prS(p,, and +&Q yields four nonlinear equations for the four unknowns, A through
D: Js’ - >1/fy2 + ‘g& fi
-
,5+Jj
-
>i(f
+
25’32c’D
+
3@C
+
L&4
-
- rA7izsL4”
= 0,
“iAD
- 1'7fz8AB = 0, -$i + >f&'+ :,&A - y&D + zp&CD + BiAD - '?JBC - 13T{2sAB = 0, - %is + pi6B + 25&D' - 1"&8B2 = 0. By appropriately manipulating Eqs. (97) through ( 100 j, it can be shown they are equivalent
(97)
(98) (99) (100) that
to 16 - 40C + 100C2 - 1478’ = 0,
(101)
2OCD - 70AC + 49A - 16D = 0,
(102)
3BC - 3AD = 0, -72 + ;iBB + 100D2 - 147B' = 0.
(103)
I-2C+
(104)
Equations (101 j through ( 104) are in a convenient form for a numerical solution. If one assumes a value for C, Eq. (101) can be solved for A, then Eq. (102) cau be solved for D, and finally Eq. (103) can be solved for B. Equation (104) is then used as a cousistency check. The results of this numerical solution are:
A = -0.532591, c = +0.744949,
B = -0.868925, D = +1..522003.
(105)
A second solution of Eqs. (101) to (104) exists (corresponding to a positive value of A), but is rejected on physical grounds. For comparison, we list the Alarshak values, which are: A = - /es,/ 8,
B= -1,
(7 = +g,
D = +3$
(106)
We note that the variational values, Eq. (105), follow the same trend as do the Marshak values, Eq. (106)) but are c.onsistently smaller (in absolute value). It, is not possible from the boundary conditions alone, Eqs. (95) and (96), to
NEUTRON
TRANSPORT
213
THEORY
compute the linear extrapolation distance. One must solve a specific problem using these boundary conditions and, from this solution, compute the linear extrapolation dist’ance. Thus the l’-3 approximation gives both a composition and size dependence to the linear extrapolation distance. So as to compare the accuracy of the Rlarshak and variational P-3 free surface boundary conditions, we consider the RIilne problem, i.e., a pure isotropic scatterer occupying the halfspace x > 0 with a source at infinity. The P-3 equations for this problem are 107)
The general solution of Eqs. (107) through at infinity is
(110) which has the proper behavior
cpo(z) = CUB.2+ /3 + ye?,
(111)
q(z)
(112)
= -;$Y,
p2(x) = +~pz,
(113)
where v = d:<:i/3 and o(, 6, and y are constants to be determined from the boundary conditions. rsing l~~qs. (85) aud (96) as the boundary conditions at z = 0 yields (115)
B= 3vCD
- 7RC’ ldC(SD
l-lC(Dil - RC)
- RC)
y.
(116)
Sow, the linear extrapolation distance, d, is defined as the reciprocal logarithmic derivative of the asymptotic(in Ohis case linear) portion scalar flux, cpU, at the free surface. Thus we have 1 3vCID - 7BC - lX(DA - BC) x cl = P ; Xt, = yjyYc :W
-
'74
1tr.
of the of the (117)
214
POMRANING
TABLE THE
LINEAR
EXTRAPOLATION
DISTAXCE THE
P-l P-2 P-3 P-4 P-5 P-6 P-7 P-8
&kLNE
II (IN
ITNITS
OK >~G.N
FREE
PATHS)
~011
kkoBLEM'
Mark
Marshal
0.5774
0.6667 o.mx7 0.7051 0. i02 0.7082
0.7i-16
n.ci940 0.7297 0. i039 0.7198 0.7069 0.7159
Variational 0.7071 0.7071 0.7118
0 The exact result is 0.7104. Using the Marshak values for 9 through D, Eq. ( 106)) in Eq. ( 117) yields the familiar result d = 0.7051 Xt, . This value, when compared to the exact WienerHopf result of 0.7104 Xt, , is seen to be 0.75% in error. Using the variational values for A through D, Eq. (105), in Eq. (117) yields rl = 0.7118 Xt, , which is 0.20% in error. Thus the use of the variational boundary conditions has reduced the error by almost a factor of four. For completeness we mention the P-3 Mark result for d which is 0.6940 Xt, , 3.31% in error. Table II summarizes the linear extrapolation distance results found in the literature for the Xlne problem as computed by various approximations. It is evident that the P-3 variational result is more accurate than higher order P-N results using either the 1\Iark or lIarshak boundary conditions. VI.
CONCLUSIONS
It has been shown that nlonoenergetic transport theory can properly be characterized by a Lagrangian, and that the use of a truncated spherical harmonics series (Legendre polynomials in slab geometry) as a trial function leads to the widely used P-N approximation. Further, the required boundary conditions for the P-N differential equations are the result of the logical application of the variational method. In particular, at an interface between two media one obtains the boundary conditions by considering an interface as the limiting case of a medium with continuously varying (in space) properties. Using this procedure one obtains rigorously the same interface boundary conditions currently in use, but which were previously derived using rather involved (and somewhat arbitrary) physical considerations. At an outer surface, the boundary conditions are derived by retaining the proper boundary terms in the functional. This leads to boundary conditions which are significantly more accurate, but no more complex, than those currently ill use. In short, the variational calculus has been
SETTROX
TRANSPORT
used t’o give a consistent and logical derivation of all aspects of the tions, and has resulted ill a more accurate P-N approximation. RE(~EIVE:I):
October
31.’
THEORY
P-N
eyua-
14. 196.1 REFERENCES
J.B.SYKES, “Seutron Transport Theory,” pp. 116-156. The Clarendon Press, Oxford, 1957. i. J. C. MARK, The Spherical Harmonic Method, I, CRT-340. Atomic Energy of Canada, Ltd., Ontario, 1947. 3. R. MAISHAS, Phgs. Rev. 71, 443 (1947). 4. G. C. PUMILANING AND bl. CLARK. JK., Sucl. Sci. E.ng.16, 147 (1963). 5. 1.1. S. SELENGTT, Hanford Quart. Rept. HW-59126, pp. 899124. Richland, Wash., 1959. 6. G. SANSONE, “Orthogonal Polynomials,” p. 172. Int’erscience, New Tork, 1959. Y. G. Y. Rl-MI-ANTSEV, J. .VUC~. Ener.gy, Ports A/U, 16, 111 (1902). 8. H. BATEMAN, “Higher Transcendental Functions,” Vol. II, p. 173. Compiled by the staff of the Bnteman Manuscript Project, A. Frdelyi, Director. McGraw-Hill, New Tork, 1953. 9. K.&I. CASE, F. DE HOFFMAS, ASU G. PLACZEB, Introductiontothe Theory of Neutron Diffusion, p. 136. Los Alamos Scientific Laboratory, Los illamos, Sew Mexico, 1953. 10. D. P. SELENGUT, private communication, July, 1963. 1. B.~)AVISON.~XD