- Email: [email protected]
Finite element methods for neutron transport based on maximum and minimum principles for discontinuous trial functions
Ann. nut/. Enerq~, Vol. 19, No. 1&i 2, pp. 565-592, 1992 Printed in Greal Britain. All rights reserved
0306-4549/92 SS.OO+O.OO Copyright Q I992 Pergmon Press Ltd
FINITEELEMENTMETHODSFORNEUTRON TRANSPORTBASEDONMAXIMUMANDMINIMUM PRINCIPLESFORDISCONTINUOUSTRIAL FUNCTIONS R.T. ACKROY D Energy Systems Section, Mechanical Engineering Imperial College, Exhibition Rd, London SW7 2BX
Department
Abstract Variational principles are given for the first order Boltzmann equation for neutron transport and its equivalent second order forms for even and odd-parity transport. For these extremum principles there are no boundary conditions to be imposed, and across interfaces the trial functions can be discontinuous. Numerical solutions are presented for the even-parity principle. The scale of interface discontinuities is controlled by a penalty parameter A. As A is increased for a given basis the degree of satisfaction of the transport equation by the optimised trial function decreases to that for a classical trial function, whereas the’interface discontinuities decrease. These trends lend themselves to a simple interpretation in terms of vectors in a Hilbert space. In applications the spatial dependence of a trial function is represented by finite elements, and a spherical harmonic expansion gives the directional dependence of the trial function. The principles are applied in the method of composite solutions in which the order of the spherical harmonic expansion can be varied according to the physical needs of the regions of a system. Large values of X are used in the method of composite solutions. On the other hand specific values of X when used with the simplest finite element lead to the well-known central difference equations for neutron transport.
1
Introduction
Finite element methods for neutron transport have achieved a geometrical flexibility comparable to that of Monte Carlo methods. The progress of the finite element method has been reviewed in of Williams and Goddard (1981), Goddard and Williams (1986), and Goddard (1991). Exposition 565
566
R . T . ACKROYD
variational treatments of finite element methods for neutron transport have been given by Lewis (1981) and Ackroyd et al (1987). The former is based on the even-parity transport equation, whereas the latter treats both even and odd-parity equations and the first order of Boltzmann transport equation. In the sequel current treatments using continuous variational solutions are extended so that discontinuous variational solutions can be used to increase the efficiency of computations. Variational treatments based on both the even and odd-parity transport equations can be used to provide efficient solutions and to establish rigorous numerical benchmarks for test problems. The maximum principles for even and odd-parity transport give independent variational approximations for the fluxes. The average of these fluxes for relatively coarse approximations can be as good as or better than refined calculation for the even-parity flux. An example of Adam's work is shown in section 7.2. The complementary principles for even and odd-parity transport can be used to establish benchmarks for both local and global characteristics of a system, because the principles are of the maximum and minimum kinds. Ackroyd and Splawksi (1982) showed that close upper and lower bounds for a local characteristic, such as the capture rate in a small region, can be obtained from complementary principles. Tight bounds on a global characteristic, viz the disadvantage factor, have been obtained by Ackroyd and Nanneh (1988). Bounds for a local and a global characteristic are illustrated in section 4.3. In almost all applications of the finite element methods for neutron transport the specification of the directional dependence of the angular flux is the same for all regions. For many applications strong directional dependence of the angular flux is confined to a fraction of the system, but its effect is significant. Then a conventional treatment, whether using spherical harmonic or discrete ordinate directional representations, is wasteful computationally. The method of composite solutions is a technique of using in a variational solution directional specifications that can be chosen independently for each region according to the physical need of the region. Thus the method of composite solutions can lead to a significant reduction in the number of coefficients to be evaluated in a variational solution. Ackroyd and Wilson (1988) implemented the method with the K +- ( ¢ + , d - ) maximum principle. The even-parity trial function and the odd-parity trial function in this principle can be discontinuous across interfaces. It is this property that permits the directional specification of a trial function in the method of composite solutions to be chosen independently for each region. The efficiency of the method is improved significantly by basing it on the K +° (¢+) principle, Ackroyd (1986), Ackroyd, Mirza, Pourzand, (1992), which uses a single trial function ¢+. In the K +° (~b+) functional there is a surface integral which is a least-squares measure of the scale of the discontinuities in at interfaces. This integral involves a weighting function and a penalty parameter A. The measure of the discontinuities for the variational solution yielded by the K +° (¢+) principle varies as 1/As. In the method of composite solutions A is taken to be large relative to a specific value of A. For this specific value of A, and a natural choice for the weighting function, the variational solution given by the K +° (¢+) principle for an infinite slab of an iostropic scattering medium is the central difference equation for neutron transport. The finite elements used are of uniform pitch and the spatial dependence of the trial function over an element is a constant. This result can be extended to a medium with an arbitrary scattering law on using the K~ (¢+) maximum principle employing discontinuous trial functions for the even-parity angular flux. The concept of complementary principles can be extended to discontinuous trial functions by introducing the K~- (¢-) maximum principle using discontinuous trial functions 4 for the odd-parity angular flux.
Finite elements for neutron transport
567
A simple means of establishing the K~ ($+) and K~" ($- )principles is to establish them, as in section 4, as special cases of the KA (¢) maximum principle for the first order Boltzmann equation. The derivation is given in section 3 for discontinuous trial functions $ for the angular flux. The I~'~ ($+) and K~- ( $ - ) principles have properties of a more general kind than the classical principles for even and odd-parity transport. This difference can be illustrated geometrically with the aid of an appropriate Hilbert space, as shown for the K + (¢+) principle in section 5. The ways in which the penalty parameter affects the variational solution, with respect to the scale of discontinuities at interfaces and the degree or satisfaction of the transport equation, are reviewed in section 6. The numerical solutions of section 7.1 illustrates the effect on flux profiles of the penalty parameter for a lattice cell. The basic theory in the sequel is for neutron transport in a typical energy group of a system with impressed sources. For simplicity the boundary condition is confined to that for a surface source. The theory can be extended readily to cover bare surface, perfect reflector and albedo boundary conditions. This exposition in the sequel follows that of Ackroyd and Goddard (1991) and Ackroyd (1992).
2 2.1
A review Boltzmann
of some
transport
first order
principles
equation
The steady state angular flux ¢o (_r,~_) at r for the direction ~ is given for each energy group by an equation of the kind n_. V¢o (_~,n_) + ~¢o (_~,n__)= S(r,n_) (1) for an absorbing and scattering system of volume V with a distributed source S. The operator is given in terms of the total cross section a(_r) and the differential scattering cross section
., (~,n_ n_n_')by
r.u (_~,n_) = ~ (_~) - f~, ~s (_~,n__.n__')u (_~,_n') dn'
(2)
On the boundary 6V of V there is a source T (_r,2) for which
¢o (_r,fl__)= T (_r,~ ) ,
n.n < 0
(3)
with L~tile outward normal to 6 V at _r. Across an interface of physically distinct regions $o (_r_r,2 ) is continuous for all those directions ~ which cross the interface at _r, but it is in general discontinuous for directions lying in the tangent plane at _r to the interface.
568 2.2
R.T. ACKROYD Parity
second
order
equations
The Boltzmann equation (1) is transformed into the pair of simultaneous equations
Q. v ~ : + c¢o+ ~ . v4+ + G - ' 40
= s+ } S-
(4)
by making the Vladimirov substitutions
¢0~ (r, ~) = [¢o (r, ~) ± ¢o (r, - ~ ) ] / 2 s + (r, n_) = [ s (r, fi) + s (r, -a_)]/2 f
(5)
The odd-parity angular flux 4,o is eliminated from the parity pair (4) to yield the second order even-parity equation - _~. v [Ga_. V¢~+] + C¢+~ = s + - a_. v G S (e) The odd-parity equation
- Q_. v [c-1_~. r e 0 ] +
G-~:
= S - - ~_. V C - ~ S +
(7)
is obtained in a similar way. The operators C and G -1 are defined by
c . (r,_~) = ~ (r). (r, ~) - / 4 . ~+ (r,_a- a_') ~ (r, a__')da'
(s)
C-~,, (,-,_~) = ,, (,.),, (,-,_~) - ~ , ,,7 (,-,~ ._~') ,, (,-, n') de'
(9)
and
Expressions are available for the operators C -x, G and ~-1 in terms of a, and a. The operators ~, C and G, and their inverses, are self-adjoint. All of these operators are positive definite. Useful identities are
F,u
= Cu + + G-lu -
~-lu
= C-Xu + + Gu-
]
i
The boundary conditions for the parity angular fluxes on 6V are
(10)
Finite elements for neutron transport
569
5Z /5'
Subregions I~,l~ etc of system volume V Subregional surfaces (i) Interior to the system , , Si for ~ , (ii) A part of the surface enclosing the system Boundary surface 6V
Sj for I~. S~ for Vi
Figure h Subregions and specified surfaces of a system
~t + G [S- - ~. V4o+] = T (r,_~) } 4I+C-*[S+--~.V4~ -] =T(_r,~)
ft.n<0
(11)
and
,/,;'- - G [S- - ~. VC, t] 4,; - C -~ [S+-- a_. v¢,;']
= T (r,-_a) } -T (r, -~_)
a"" > 0
(12)
with n the outward normal to 6V. Across an interface of physically distinct regions ~ (_r,~_) are continuous for all those directions ~_ which cross the interface at _r.
570
2.3
R.T.
ACKROYD
Classical principles
Associated with the second order even and odd-parity transport equations are the classical maximum principles K + (¢+) and K - ( ¢ - ) for ~b+ a trial function for ¢+, and~b- a trial function for ¢~-. [In the sequel ~b+ (r, fl) and ¢ - (_r,fl_) are even and odd functions of _fl respectively] The g + principle is classical because (i) the Euler-Lagrange equation for the K + (¢+) functional
V 4,r
V 4~-
+
4/r
/" in"
+ (_',m
d dS
6v _q.n_<0
_/I,._.._,(::.°,,
.,,
6 V 4T
is the even-parity equation (6); (ii) ~+ has to be continuous across interfaces in the manner (c) below. A trial function ¢+ (r_r,~_) is admissible in the i f + principle if :- (a) It is a continuous function of L' and ~_ almost everywhere within the subregions ~ into which V is partitioned by the specified surfaces of Fig.(1). These surfaces comprise the boundaries of physicM regions, the edges of neutron shadows (Ackroyd 1978), and the artificial boundaries introduced for computational convenience, eg. the faces of finite elements. (b) Within each subregion 2 " V¢ + is almost everywhere a continuous function of r and l). (c) On each interface ¢+ (r, 12) is continuous at _r for all directions ~_ which cross the interface. Note that the trial function is not required to satisfy on 6V the boundary condition for a surface source. The functional K - (~b-) for the K - maximum principle is
V 4~r
-
dV V 4~r
+
4f f
6v _n.n
-
4//,~_.nl(¢-)2dfmS
(14)
6V 4~t
The conditions for ¢- to be admissible in the K - principle are of the same type as the conditions (a)-(c) above for ~b+ to be admissible in the K + principle. For finite element work a trial function ¢+ or ¢ - is expressed as a linear combination of basis functions. Each basis function is the product of a function of position _r and a function of direction l). The functions of position are the shape functions for the finite elements into which the volume
Finite elements for neutron transport
571
V has been partitioned. Either spherical harmonics or discrete ordinates are the usual means for the specification of the functions of direction. Spherical harmonics have the advantage of being the analytical companions of the Legendre functions used to specify the differential scattering cross section. Even spherical harmonics are used for ¢+ and odd spherical harmonics for ¢ - , and either requires half the number of spherical harmonics to specify the angular flux ~b= ¢+ + ¢ - . In the maximum principles g +(¢+) K - (¢-)
_
(15)
equality holds if, and only if, the trial function is identical with the corresponding parity component of ¢o. If ¢+ is the even-parity trial function that maximizes K + (¢+) for trial functions with a given basis then the corresponding odd-parity angular flux approximation is usually taken to be
(16)
:
by analogy with the exact result (4). Similarly for a maximizing trial function ¢ - the corresponding even-parity angular flux approximation is (17)
=c-'
The principles K + and K - are complementary in the following sense. It can be shown that
[[{s+c-is + +s-Gs-}d v
:
V 4~
=
2a
(say)
(18)
On defining 2/3 = K + (¢+) - K - (¢~-)
(19)
there follows for the functionals K + (¢+) - a and a - K - ( ¢ - ) the complementary principle g + (¢+) - ~ < ~ < ~ - K - ( ¢ - )
(20)
with the equalities holding if, and only if, ¢+ = ¢+ and ¢ - = ¢~-. This result is shown to be of some value in section 4.3 for the benchmarking of transport solutions.
2.4
A neo-classical
principle
Although the Euler-Lagrange method is restricted to differential equations of even order, Axelsson and Barker (1984), a maximum principle for the first order Boltzmann equation can be obtained, Ackroyd (1983 a, b), using a generalized least-squares procedure. This method is illustrated in section 3.3 by deriving a maximum principle for the first order Boltzmann equation for which discontinuous trial functions are admissible. Here the neo-classical principle f~" for the first order transport equation is stated.
572
R . T . ACKROYD
Special choices for the continuous trial function 6 yield the classical K + and K - principles. The continuity conditions for 6 are of the same kinds as those for 6 + in the K + (6 +) and K - (6-) principles. T h e / ~ (6) functional is
V 4=
-- //{(~-.~'V6)r-1~--'V6-6r6}d a d V v 4~t
6 v fi-n_<0
-
4 / / I n , a l ( 6 ) 2 dadS
(21)
6V 4~
The maximum principle is
(6) _< k (6,) = 2~
(22)
with equality holding if, and only if, 6 = 60. Substituting the special choices of trial functions
(a)
6
(b)
6 -- 6 + + 6 -
=
6 + + 6; (23)
in the/~ (6) principle leads to the K + (6 +) and K - (6-) principles respectively. This simple means of deriving complementary principles is illustrated in section 4 by deriving maximum principles for even and odd-parity transport that accept discontinuous trial functions.
2.5
P r i n c i p l e s for d i s c o n t i n u o u s
trial functions
K (6) principle A generalization of t h e / ( ( 6 ) principle is the K (6) principle
g ( 6 ) < 2~
(24)
Admissible trial functions 6 can be discontinuous across interfaces, and they do not have to satisfy the boundary condition. Equality holds if, and only if, 6 = 6o. The R"(6) functional is given by
USi 4~¢
u ( . % n S j ) 4~-
Finite elements for neutron transport
573
Jumps at r in the trial functions are ~:k ( ~ , ~ ) _ ~- (rt,fi)
~
" ~
Interface S../'IS.
Figure 2: Jumps in trial functions at r on the interface S~ FI 8j of elements I~ and where the definition (21) for/{ (~b) is extended to discontinuous trial functions. The USi is the set of all element faces which lie within the volume V, and n~ is the outward normal to the element surface Si. The set of all element interfaces is U(S~ NSj) and ~ (~.~,i~) - ~ (r~,fl) is the jump at r on the interface St N Sj in the trial function in going from the element with volume I~ to the element with volume Vj in Fig(2). The weight ~7 (£, ~) is non-negative and is usually chosen to be positive for ft. nl ~ 0, where n= is the normal to Sl. The penalty parameter A is positive. Rules for selecting ~17 and A are obtained in the same w~t.vas the rules derived in section 6 for the maximum principles for even and odd-parity transport which agimit discontinuous trial functions. K +- (~+, ~ - ) principle The K +- ( ~ + , ~ - ) principle can be used to derive all of the other principles mentioned in this section. In the principle ~+ and ~ - are trial functions for ~b+ sad ~" respectively. They can be discontinuous at interfaces, and they need not satisfy the boundary condition. For the principle
K+-
<
equality holds if, and only if, ~b+ = ~+ and ~b- - ~ .
The functional K +- (~+, ~-) is
= K÷(,+)+K-(,-)-./ f t.t$ i 4 w
u(s~ns~) 4.
with l ~ (_r,fl) non-negativeweights,usuallychosen to be positiveforQ . m # O. Putting
=W-=W
(26)
574
R . T . ACKROYD
=
+
(28)
gives the I'~ (~) principle. Using the substitutions (28) and letting A ~ oo gives t h e / ~ (~) principle. The substition = ~bo (29) transforms the principle (27) into the even-parity principle, Ackroyd (1986), r ~+° (~+) < K+ (~o+)
(30)
with equality holding if~b+ = ~b+. The way in which this principle operates can be illustrated geometrically by means of a Hilbert space, in very much the same as that give in section 5 for the related even-parity maximum principle K + (~b+). Composite solutions Ackroyd and Wilson (1988) used the K + - (~b+, ~ - ) principle for the method of composite solutions, with A large relative to a specific value for this penalty parameter. For the specific value of A , and a uniform mesh of the primitive elements of section 6.2, the optimised trial functions satisfy the central difference equations for neutron transport in an infinite slab of an isotropic scattering medium. T h e drawback of the K + - (~b+,~b- ) principle in the use of two trial functions is removed for the K +° (~b+) and K ~ - (~b') principles. The K + (~b+) principle is derived in section 4.1 It is superior to the K +° (~+) principle, because with the same kind of trial function for an infinite slab, and the same specific value of A , it gives the central difference transport equation for an arbitrary scattering medium. The K + (~+)principle is illustrated in section T with numerical solutions for a lattice cell.
3
3.1
Maximum tion Error
p r i n c i p l e for t h e first o r d e r t r a n s p o r t
equa-
measures
In constructing a maximum principle some measures are required of the errors made by a trial function ~ (_r, ~_) for the angular flux #0 (__r,~_). A measure of the error made by ¢ in not satisfying the transport everywhere within the volume V, and for all directions, is defined to be
V 4r
Since ~-1 is a positive definite operator, the volume error e~ vanishes if, and only if, ~_- V~b + E~b - S = 0
(32)
Finite elements for neutron transport
575
ie.ff is a solution of the transport equation. The failure to satisfy the boundary condition given is measured by
~,~-2/ / In.fil[~(r.,fi)-T(,',fi)12df~dS
(33)
6V fl.n_<0 where tile normal n is outward to the boundary surface 6V. To define the error incurred by the failure of ~bto satisfy the interface continuity condition, consider in Fig 2 a typical interface Si f3Sj of two spatial finite elements V~ and Vj. For a point r on S, N Sj let ~ (r~, f/) denote the value of the trial function as given by its expression for the element ~ . Similarly ~b(_rj, fl_) is the value of the trial function at r as given by the trial function's expression for I,~. The error incurred at the set U (Si N Sj) of all interfaces is specified as
¢/=A
f
{fl •n~{_~..~ [~ (r~,f~) - ~b(_r,,fl)] } ×~_1 [~(~,~)_ ~ (_.~.~)] }
f
dfldS
(34)
u(S,nS~) 4.
where ni is the outward normal at r to Si. The penalty parameter A serves to bias the interface error relative to the other errors. The cumulative error is v~ (~)
=
S])dV
/(ft. v ~ + r.~ _ s, ~-1 L~. v ~ + ~ V
6v ~-n_<0
f f
{"-"-[*("'"-)-*('-""-)1} ×r.-, {._... [~ (.., a_) _ ¢ (_._.._)] }
o(s,.s,)'For brevity the notation
f f uvdfldV = f(u, v)dV is used in the V 41r
sequel. The functionl Vx (~b) has
V
the property that Vx(if) >_ 0
(36)
because )"~--1 is a positive definite operator.
3.2
A least-squares
identity
Expansion of the least-squares error functional Vx (O~) gives the identity Kx(~b)+Vx(~b)
=
/(S,E-1S, dV + 2 / V
=
/
[f~'n[T2dfldS
6V fl.n_<0
2a
(37)
The functional Kx (d) can be expressed as
K~ (~) =
[ [2(~,s) + 2(ft.v~, r.-'s)]dv - P (~,~) V
6V fl.n<0
USi 4~r
576
R.T. ACKROYD
where the bi-linear functional F (0, ~b) is defined as ~' (0, ~b) --
/ { ( ~ . V0, ~-]-1[~. V~) "dr (0, Y]~b)} dV v
+ f/[~,nlO~bdf~dS=
F(~b,0)
(39)
6V 4x
The expression (38) is obtained by using the divergence theorem to convert the term 2f(ff, f~ • V V~)dV=f f div [Q_~b2] d~dV in the expansion of the volume component of V into surface integrals v 4~ for USi, the set of all element faces inside V, and the boundary 6V.
3.3
The
K~ (~b) m a x i m u m
principle
It follows from the definition (35) that Vx (~b,) vanishes to give
Kx (¢o) = 2~ The condition
(a)
K~ (¢0)
implies (b)
(40) = 2~ ] /
~
(41)
= ~0
as follows. Putting
x (r, n__)= ~ (_~,9_) - ~0 (r, 9_)
(42)
the condition (41a) and the identity (37) give
0= Vx(if)= f [(_Q.vx, x~-lfl__,vx)+ (x, ~cx)] dV+ f f I,,nlxZdftdS V
(43)
6V 4x
on expanding the volume integrand of Vx (~b), noting that ~-1 is self-adjoint, and applying the divergence theorem to the term f f [2. ~7(X~) dfldV. The expression (43) implies the vanishing V ,ix of X, ie.~b = ~b,, because E and its inverse are positive definite operators. The conditions (40) and (41) and the property (36) establish the maximum principle Kx (~b) _< 2~
(44)
with equality holding if, and only if, ~b= ~b. If ~b(r, _~) is taken to be continuous across each interface, for all directions D crossing that interface, then Kx (5) =/~" (~b) (45)
Finite elementsfor neutrontransport
577
where
k(¢)
=
f [2(¢,s)+2(~.v,,r='s)]d r -
P(¢,¢)
V
SV N-n_
and
k (~b)_<2a
(47)
is the maximum principle of Ackroyd (1983 a,b). For this principle admissible trial function have to satisfy the interface continuity condition.
3.4
M a x i m u m principle in t e r m s of the parity c o m p o n e n t s of the trial function
On expressing q~as the sum of the parity components ¢+ and ¢- the functionals Vx (¢) and Kx (~b) become, with the aid of the properties (10) for I[i-1 and E ,
:
+(_~. v¢+
f{
+ G-l,-
- S-.G
[ft. V, + + a - ' , - - s - l ) }
dV
V
6V fl.n_.< 0
( {ll. ~ [~+ (,~,a)- ~+ (~i,a)] } +
A
f
/
)
x G { a ' n ~ [ ¢ + ( r - / ' a ) - ¢ + (_r,,f~)]}
+ {a-~ [¢-(,~,a)-,- (,_~,a)] }
dadS
(48)
×c-1 {~. ~ [,- (~,a)_,- (,~,~_)1}
o(s,osj),,
and
+,-) =
K+(,+)+K-(,-)-2f U S i 4~"
_ ~ f
f
o(~,°~,)'-
×a{,.o_-.~ [,+(r__i, (r~,_,)-, (~_,,._)1} _i_ { fI. n_./[,;b_ a) _ ~;b_+(r_./, a) ] } ×c-' 1~,_.~_,[,- (~..,_)- ,- (~_,a)] }
dadS
(49)
578
4
4.1
R.T.
Complementary port Maximum
principles
ACKROYD
for even and odd-parity
principle for even-parity
trans-
transport
For the special choice of trial function
¢ = ¢+ + ¢:
(50)
the functions (48) and (49), taken with the parity equations (4) and the boundary conditions (11) (12), give
V~ (¢+ +¢2)
=
r + (~,+-¢+,,~,+-¢t)
+
A
f f
G{~'n~[dz+(r~'~)-ek+(£J'~)]}× { a ~ [~+ (~,~) - ~+ (_~,~_)] } dadS
(51)
u ( S , n S i ) 4~r
and Kx(¢++¢:)
-
K-(¢:)
-
[J
=K +(¢+)
( ~ {a ,~ [~+(~,a)- ~+ (_~,a)] (-~,-~)]} }
[J
×al._.n~ [~+ (.~,._)-~+ +2~. ~ . [¢+ (a,~) - ¢+ .(,~,~)] ¢: (,,~) . .
u ( S , n S ~ ) 4~r
)
d~dS (52)
Here F + (u +, v+) is the bi-linear functional
r+ (u+,v+) = [ {(~. w+,G~, w+) + (~+,Cv+)} av v
+ f f I"-.nl,,+v+d.dS
(53)
6V 4~
with the property that F + (u+,u +) is positive definite. The bilinear function F - ( u - , v - ) is
F-(u-,,~-) = [ { l ~ v u - , c - ~ v,,-)+lu-,C-~v-)}dv v
8 V 41r
with the property that F - (u-, u - ) is positive definite. In obtaining the last result from the expression (49) the integral over US/, the set of all element faces within V, has been written as an integral over the set of all interfaces U(Si I-I
St)
On making use of the results F + (¢+,¢+) g-(¢~')
+ F - (¢~-,¢~-) = 2a =F-(¢0,¢o)
f
(55)
Finite elements for neutron transport
5:9
Ackroyd (1981), the expressions (51) and (52) transform the identity (37) into
Kt(* +)
=
F+ (~+o, <'t) + 1
fe_:G_l~Zaaas
[-~u(s, nsD 4=
(56)
where K +(~b +)
=
K +(¢+)
a {~.~ {¢+ (a,~)- ~+ (a, a) }} U($1rlSj)
4~
Since G and C are positive definite operators the term in parentheses in the identity (56) is nonnegative. Hence the inequality K + _< F + ( ¢ + , ¢ + ) +
1 -~ f f ¢oG-l¢'~d~'ldS
(58)
u(s,ns~) 4~r
with equality holding if ¢+ = ¢+
(59)
For then the term in parantheses in the identity (56) vanishes on account of the interface continuity condition for ¢~. On the other hand the condition
'//
K + (¢+) = F + (¢+, ¢+) + ~
¢:a-l¢'[dgldS
(60)
u(s~ns,) 4r can hold only in the limit ~ ---* oo. If it held for finite ~ the the following conditions would hold viz r + ( ¢ + - ¢+, ¢+ - ¢ + ) = 0
(61)
and
Condition (61) would imply ¢+ = ¢+, and then the continuity ¢+ at interfaces would imply that the condition (62) did not hold for ~ finite. Nevertheless maximizing K + (¢+) with A finite can bring ¢+ close to ¢+. This is shown by the geometeric interpretation of the K + (¢+) principle in section 5
580
4.2
R . T . ACKROYD
Maximum
principle
for odd-parity
transport
The K ~ - ( ¢ - ) maximum principle for odd-parity transport is obtained in the same way as the K + (¢+) maximum principle. The result is K~" ( ¢ - )
=
K- (¢-)
c-,
u(S,nSj4,)
b-
}
U(SinSj) 4 x The out of balance in this inequality is
{..n~[¢-(r~,.)-¢-(rj,~)] } F - (~b- - ¢0, ¢ - - ¢o) + A
o,s,ns,>.~, c-1{
4.3
Complementary
a--'"~[*-(~'}a-)-¢-e-b' _ c ~n--)] ;
d~dS
(64)
principles
For the classical K + (~b+) and K - ( ~ - ) maximum principles there holds the result of Davis (1968), which can be expressed as K + (¢+) - a
= F + (¢+,~b +) + F - (~bo,~b~-) = F + (¢+,~b +) - F - (¢,,~b~')
[
(66)
Corresponding to the classical complementary principles K + (¢+) - a and a - K - (~b-) there are the complementary principles K + (5+) - a' and a' - K~" (~b-) with the property
g + (¢+) - a' _
(67)
Here 2a'
=
F + (¢+,¢+) + F -
1
+ -~ f
(¢0,¢~-)
f[*+*C¢+o+¢ZG-'¢:]dadS
u ( s , o s , ) 4~ 21~'
=
F+
+
+) - r-
1
o(s,osD4. The classical complementary principles have been used for benchmarking finite element solutions in two ways, viz
Finite elements for neutron transport
58l
,, •
flux',
O(x) I
I I
',.t
fuel
', s=so l
void
/
i o=s0 1o~1
o=0
I
~'-
moderator
's=~,
,o=1 ot o~0.9 =1,
\ o~0.9
\
l i
:perfect rellector ~
O
I
I I.
0
2
/: J
: I
I
3
5
base surface
6
x
8
Figure 3: Solution of the edge cell problem 1. to place bounds on a local characteristic of a solution 2. to place bounds on a global characteristic of a solution
Bounds for a local characteristic
For the edge-cell problem of Reed (1971), Fig 3. the capture rate of neutrons in the can region has to be determined accurately for it to be used as a benchmark. The flux in the numerical solution of Galliara and Williams (1979) varies rapidly in the can. The bi-variational method of Ackroyd and Splawksi (1982) shows that the capture rate is in the range 0.7855 + 0.0036, ie. it is known to within + 0.46%. This example shows that classical complementary principles can be used to give benchmarks of worthwhile precision. Bounds for a global characteristic Many numerical methods for the disadvantage factor for a reactor cell converge slowly, as shown by table 1. This global characteristic is the ratio of the average flux in the moderator to the average flux in the fuel. The results of table 1 refer to a uranium/water slab lattice with the thicknesses of the fuel and water regions at 0.167 and 1.5 mean free paths. The juxtaposition of the uranium strong absorber and water scattering medium give rise to significant transport effects. I Method ANISN-S8 ANISN-SI6 RIPPLE DI FTRAN FELQO EK+ t FELQOE K - t
fuel]
moderator
Average Average Disadvantage 1 flux flux factor 3.00081 3.24259 1.08057 2.99733 3.29745 1.10013 3.00980 3.31400 I.I0698 2.99466 3.34050 I. 11599 2.99390 3.35246 1.11976 t29 moments, 8 quadratic elements mean=l.12087 2.99351 ] 3.35865 1.12198
Table 1: Numerical results for the average fluxes and disadvantage factor using 8 meshes
582
R.T. ACKROYD 1.6
[,
1.5
H
o
1.4 , ,
E v e n - parity
z~ O d d - parity
~~ . 3 ~ ~1.2 c21 1.1 1.0
0
5
10 15 20 25 Number of moments
J
30
Figure 4: Comparison of K + and K - flux plots using 8 quadratic elements for a uranium/water lattice cell Ackroyd and Nanneh (1988) have shown that the K + and K - principles give lower and upper bounds respectively for the disadvantage factor. The table shows that the results for the well known discrete ordinate code ANISN, the collision probability method RIPPLE , and the elegant DIFTRAN method of Nanneh and Williams (1985), (1986) have not converged fully. The K + and K - bounds for the disadvantage factor as functions of the number of Legendre functions, used to represent the directional dependence of the neutron transport, are shown in Fig 4.
5 5.1
Geometrical interpretation of the even-parity maximum principle A Hilbert
space
A geometrical interpretation of the effect of maximizing K + (~b+), with respect to the coefficients specifying the trial function ¢+, can be given using a Hilbert space that employs an even and an odd function of _~ to define a vector of the space. I f O ~ [0+,0 -] and ~ ~ [ ¢ + , ¢ - ] are a pair of arbitrary vectors of the Hilbert space associated with the K + (~b+), principle then their scalar product is defined to be
o. • = [ [(~_.vo+,G£, v¢+) + (0+,C¢+)]dV v
+//,o 6V 4.
{
]/
×
The even-parity function 0+ (_r, f/) is admissible in the K + (¢+) prinople. The odd-parity function ¢ - (_r,~) is admissible in the classical K - (¢-) principle. Whereas 0 + can be discontinuous across interfaces ¢ - is required to be continuous at interfaces for all directions crossing the interface.
Finite elements for neutron transport
5~3
The definition (69) satisfies the commutative, distributive, associative and multiplication by zero laws required of a scalar product. The commutative law holds because G and C are self-adjoint operators, and the distributive law is satisfied since G and (7 are linear operators.
5.2
Geometrical
transcription
of a variational
identity
The identity (56) can be expressed using the scalar product (69) as K + (~b+) + ((I)+ - @0+- ~o-)2 = (~_o + +~_o-) 2 with • + . . [~+,01,
~
~ [~0+,01
and
2 - ~ [o, ~:]
(70)
The results (¢__o+)2= F+ (¢+, ¢ + )
(71)
f f *z~'¢Za.dS
(2-)2=
(72)
u(s,nsO 4~
~.
2- = 0
(73)
are established readily. T h e functional K + ( ¢ + ) can be expressed in terms of ~ + , ~ +
and ¢--o- by
means of the cosine rule identity viz (*+) 2 + (~+ - [~__o++ ~__o-])2- 2~+ • (~+ - [~__o++ ~_.o-])2 = (~__o++ ~__o-)2
(74)
This identity is recast by using the result (73) as
2_~+ • ( ~ + ~ - ) - (~_+)~ + (_~+ - [ ~ + ~ - ] ) ~ = ( ~ ) ~ + ( ~ - ) 2
(7~)
Comparison of the results (70) and (75) gives K~+ (¢+) = 2 ~_+. ( ~
5.3
Geometerical
interpretation
+_~:) - (~_+)~
of the variational
(76)
principle
For a typical trim function expressed as M
¢+ (_~,_a) = ~
am + ~m + (-~, a__)
rn.~l
in terms of the basis functions ¢+ there corresponds the vector M
,__+
~ m=l
+ +
(77)
584
R.T. ACKROYD
w h e r e ~ ~ [~+, 0]. Maximizing K + (~b+) is according to the identity (70) equivalent to minimizing (~+ - ~
- ~o-) 2 with respect to the
a +. Thus (f
-*,+ -~-)
=0
(m= 1,...,M)
(78)
where --
=
z-.,m=1
m-~-m
~
f
(79)
and the~+~ are the values of the coef~cients .+ which m~ximize K~+ (~+). Summing for all m the product of ~+ and the result (78) gives
• ÷. Thus defining the target vector ~
=0
<80)
to be = ~
+ ~--o- ~ [~+, ~0]
(81)
it can be said that in seeking the optimum approximation vector _~+ *-* [~+, 0] by the minimization of (~+ - ~__o)2 the target error vector ~+ - ~
~_+
is orthogonal to ~+, i.e. -
= 0
This result is illustrated in Fig 5. The geometrical effect of seeking the optimum approximation d+ by maximizing K + (~b+) is to determine a point P of the linear M space (77) which is nearest to the unknown point Q specified by the target vector ~--o- Note that the definition (81) and the result (73) give (*0) 2 = (~_o+)2 + (~_o-)2 (83) Since (~__o-)2 is given by the expression (72) the distance QP ( ~ 2 of the target ~ vector from the exact solution vector ~+ varies inversely with A . As A --, oo the magnitude of the target vector shrinks to the magnitude of the exact solution vector, and the target vector aligns itself with the exact solution vector. Then the geometrical interpretation of the maximum principle becomes the classical interpretation, Luneberg (1969), as illustrated for the K + (~+) principle by Ackroyd (1986).
6
6.1
The
role of the
penalty
parameter
A
Effect of increasing
The effect on variational solutions for neutron transport is described here for solutions obtained with the use of the Kx (~) principle. The relevant properties of the Kx (~) functional can be
Finite elements for neutron transport
__
585
g+ - ( @^+ - % ) = 0
Error Vector
Og
Offset Vector
cg- ®o--o
Figure 5: Solution, target and optimum approximation vectors established in the way used by Ackroyd and Wilson, (1988), to obtain the corresponding properties for the K +- (~b+,~b- ) principle. The cumulative error functional Vx (~), (35) is written as
vx (¢) = i, (¢) + ~i~ (¢)
(84)
where
i,(,) = f f
V*+ *-SldadV
V 4¢
+ 2/
f
Ifi.n_l[7"-,l'dadS
(85)
tv fi.n_.< o
is a measure of the failure to satisfy the Boltzmann equation within elements and to fulfil the boundary condition, and ld (~) =
f
/
{~ "n--'i [~ ( ~ ' ~ ) -- ~ (-rJ'~)] } x
dfldS
(86)
u(s~nsi) 4,~
is a measure of interface discontinuities The basis for the trial functions ~ is assumed to contain a sub-set capable of representing classical trial functions satisfying the continuity condition for interfaces. Let ~bx be the trial function which maximizes Kx (~b) then lr (~bx) is monotone increasing and Id (~bx) is monotone decreasing with A . A s A --- oo, lr (~bx) --~ I, (~b¢) where ~c is the classical trial function which maximizes Kx (~b). Thereby ~c maximizes the functional/((~) of the maximum principle for the Boltzmann equation
586
R.T. ACKROYD
which admits only classical trial functions. Since ~, is not the maximizing trial function of KA (~) for trial functions which in general are discontinuous / i (¢c) = KA (~,) < KA ($A) The identity (37) gives KA (CA) + Ir (~bx) + Aid (~bA)= 2 a = KA (¢~) + I~ (~bc) Hence the simple inequalities I~(¢~)
<
It(Co)
<
ir (¢o)
(87)
A
indicating the trends. As )~ is increased the measure Ix (¢~) of the interface discontinuities is decreased, whereas the degree of the failure to satisfy the Boltzmann equation within elements and fulfil the boundary condition, as measured by I~ (¢~), is increased. These remarks suggest that there may be some specific values of A for which a discontinuous variational solution is superior overall to a classical continuous variational solution. In the next section a good variational solution is obtained for a specific value of A and the simplest discontinuous finite element. This specific value of A can be used as a guide in selecting an appropriate penalty parameter for a particular application.
6.2
Specific penalty
parameter
A specific value of the penalty parameter A is obtained for the K + (¢+) principle by considering the neutron transport induced in a homogeneous infinite slab of an anisotropic scattering medium by a source varying over space in a piece-wise constant fashion. For simplicity the source is taken to be isotropic and of strength S, a constant over the mth spatial element (Zm-l,Zm) of the slab. The primitive elements employed are assumed to be of constant pitch h. The shape function for a primitive element is unity within the element and zero outside. The trial function¢ + (x, p)is given, for the mth element by = am
= am
+
(88)
where cos - l p is the angle made by the direction ~_ with the x-axis. The functions am (~) maximize K~ (¢+), e is a parameter and e0m (ju) is an arbitrary variation of the trial function for the m 'h element. In taking variations of K + (¢+) the trial functions can be of the kinds ¢+ (x,/~) = ( tim (#) + tO+ (/J) tim, (/~)
for element m for all elements m' # m
(89)
because with no interface conditions to be satisfied the variations for the elements can be made independently of each other.
Finite elements for neutron transport
587
Let K+m (@+) be the contribution to K + (@+) that involves O+, then the general formulae (38) and (39) give for m an interior element of the slab 1
K+x(@+)
=
2.
- [ 6 r e + e 0 +] C [ a m + e 0 +]
dr
--1 1
-~,[am+eO~-am_,]G(-~la,,, +eO+-am_,]) } dl,
_
(90)
-1
since ~ - n / =
p for the interface Xm and ~ .n_/= - # for the interface Zm-1.
When K + (¢+) is maximized for the trial function (89) then d
+
i.e. 1
o =
h / o+,t [s~* - Cam] d. -1
-,
{
. a ( . [ a m - am+,]) -~G ( - ~
d.
(91)
[am - a~_~])
since C and O are self-adjoint operators. As 0m + can be arbitrary the condition (91) gives
_~. { c ( . [am+, a m- l ) h _ C (. [amh a~-,l) } +Cam =S+ Choosing 1
A= g
(92)
gives the central difference scheme -#
{G(p[6m+l-am]/h)-G(P[fim-am-1]/h)}+Ofim=S+
(93)
h
for comparison with the even-parity equation
for the even-parity angular flux ¢+ (x, p) in the slab with an isotropic source S +. Central difference methods have been used by Fletcher (1983) to obtain good solutions for neutron transport in systems coverable by a regular spatial mesh. In his work Fletcher has used a spherical harmonic expansion to treat the directional dependence of the even-parity angular flux. The use of the K + (@+) principle with the simplest spatial finite element, viz the primitive element, has lead to the good approximation of the central difference equation (93) for the transport equation (94). This suggests that with more refined spatial elements than the element the K + (¢+) principle has potential for producing better solutions.
588
7 7.1
R.T.
ACKROYD
N u m e r i c a l solutions K~+ (~b+) s o l u t i o n o f t h e R e e d
problem
The way in which the penalty parameter A influences the convergence of K~+ (~b+) variational solutions is shown in Fig.6 for the Reed problem. The version of the K~+ (~b+) principle used for this example has an extra term for the perfect reflecting surface, and the bare surface is accounted for by setting T to zero. The extra term is -A //(
{~" n [~b+( r , ~ ) - ~b+ (_r,fl__')]}. ×G fl) )]}
)d~dS
(95)
~V ~ 4x
where dfV~ is the additional surface presented by the perfect reflector, and 2" is the reflected direction to the incident direction fl__at r.r_. In the results of Ackroyd, Abuzid and Mirza (1992) for the Reed problem the number of Legendre functions used per region for the directional dependence of the trial function ~b+ has been chosen as follows. The flux over most of the fuel region is expected to be flat, because the source is constant over a region several means free paths thick. Two moments are probably sufficient here for ~+. With a constant flux expected for the void, a high flux induced locally by the moderator source, and a flux depression int he source-free strongly absorbing can, the great variations in the flux necessitate a fairly high Legendre expansion of 4 moments for ~+. In the far region of the moderator without a source there is a short distance from the moderator region with a source to the bare surface. Five moments are perhaps adequate for the far region of the moderator.
7.2
Averaging
results from complementary
principles
Disadvantage factors For the uranium-water lattice of section 4.3 the asymptotic value of the disadvantage factor is 1.20087. The error of the mean value of the K + and K - results for the disadvantage factor is shown in table 7.2 as a function of the number of Legendre moments employed. Thus the mean of four moment K + and K - calculations for the disadvantage factor is better then either an $16 discrete ordinates calculation or a RIPPLE result with errors at 1.78% and 1.17% respectively. The good results for the mean disadvantage factor for four or more moments is a consequence of an approximation for the flux which is provided by the mean of the K + and K fluxes. Flux profiles Problems arising in the transport of thermal radiation are concerned usually with optically thick diffusive regions. For practical reasons the spatial mesh employed has to be large compared to the mean free path for such regions. Adams (1991) has solved some typical problems using both
Finite elements for neutron transport
589
~lGo
2.4-
.o g
JXD 1'3 LL
rr
1.2,
L--~Z.'~--7,.:
<
•
09
\", 0.6.
U
\
X CELL WIDTH (CM) Figure 6: Effect of penalty parameter A on K + solutions
Number of Moments 1 2 3 4 5 10 15
Percentage Error of mean disadvantage factor 14.870 4.850 2.270 0.890 0.790 0.130 0.024
Table 2: Percentage error of mean disadvantage factor
20t •
~"
Fluxes
1.5
(rdmm~) (~) (odd) (ave=go)
8 n--o
~¢¢=¢L----.,. =~l . ° r
o.5 ---~s,~ 0.0 1 0
2oo
x
400
6oo
Figure 7: Averaging of parity fluxes (Adams 1991)
590
R . T . ACKROYD
even and odd-parity finite element transport solutions. He has considered a two region problem for an infinite slab, with the incoming flux at the left-hand face of the slab having delta function of direction about p = 0.1. The left hand region of the slab is a purely absorbing skin with a thickness of a thousanth of a mean-free path. The second region is a source-free purely scattering medium some 600 means free paths thick. A reference solution was obtained for this problem using a linear-discontinuous diseretization of the first order transport equation with a very fine spatial mesh. For the even-parity and odd-parity transport coarse mesh solutions ten elements used were for the second region. In the thin region, which is well resolved spatially by the mesh employed, the even and odd-parity solutions for the flux are well out with the reference solution. As shown in Fig.7 the average of the even and odd-parity solutions is in good accord with the reference solution. The average used does not have to be an unweighted average. Both Nanneh (1991) and Adams (1991) have used in effect
~+(z,#)=7 ~+(z,#) +(1-7) average
derived from ~b-
~+(z,#)
(96)
e v e n - p a r i t y solution
with 7 a coefficient to be determined in some way to improve ~b+ average. The even-parity derived angular flux is given in terms of if- as follows
derived
This component for the even-parity angular is in general discontinuous at interfaces The coefficients 3' can be determined, as an independent coefficient for each element, by substituting an expression of the kind (96) for each element in the Kx+ principle, with A given by the specific value given in section 6.2.
7.3
Some observations
There is the prospect of applying the complementary principles, which admit discontinuous trial functions, in the same ways as those found to be useful for the classical complementary principles, viz to obtain variational solutions for practical problems, and to provide benchmarks for test problems by constructing upper and lower bounds for both global and local characteristic. In seeking variational solutions there is one feature in which the K + principle has more flexibility than the K + principle. The optimized solution for the latter has the merit of satisfying neutron conservation for the whole system. For the former principle the trial function can be specified independently for each element. Consequently neutron conservation can be imposed on each element, ie. for each element V/
/ / {-~_" V [G~" V~b+]+ C~+ - S+ } df~dS = O V, 4~"
(98)
Finite elements for neutron transport
591
These conditions serve to impose neutron conservation for each dement, by making that coefficient in the trial function $+ which specifies the element flux dependent on the other coefficients for that element. In seeking variational solutions more efficient than the classical eontinuoua solutions attention has been focused on the method of composite solutions and the averaging of solutions; but the variational principle K + ($+) can he used as the basis of 'nodal' solutions, with neutron conservation for each element or 'node' ensured by the condition (98).
References
Ackroyd, R.T. (1978), A finite element method for neutron transport-I, some theoretical consideration. Ann. nuei. Energy, 5, 75. Ackroyd, R.T. (1981), The why and how of finite elements. Ann. nucl. Energy, 8, 539. Ackroyd, R.T. and Splawski, B.A. (1982), A finite element method for neutron transport-VI, upper and lower bounds for characteristics of solutions. Ann. nucl .Energy, 9, 315. Ackroyd, R.T. (1983a), Least squares derivation of extremum and weighted residual methods for equations of reactor physics. Ann. nucl. Energy, 10, 65. Ackroyd, I~.T. (1983h), A finite element method for neutron transport-VII, completely boundaryfree maximum principle for the first order Boltzmann equation, Ann. nucl. Energy, Lfr, 243. Ackroyd, R.T. (1986), Generalized least squares as a generator of variational principles and weighted residual methods for FEM transport methods. Prog. Nucl. Energy,18, 7. Ackroyd, R.T., Fletcher, J.K., Goddard, A.J.H.G., Issa, J., Williams, M.M.R. and Wood, J. (1987), Some recent developments in finite element methods for neutron transport theory, 381-383, Advances in Nuclear Science and Technology, 19, edited by J. Lewins and M. Becker, Plenum Press, New York. Ackroyd, R.T. and Nanneh, M.N. (1988), Upper and lower hounds for disadvantage factors as a test of an algorithm used in a synthesis method. Ann. nucl. Energy,15, 241. Ackroyd, R.T. and Wilson, W.E. (1988), Composite finite element solutions for neutron transport, Ann. nucl. Energy, 15,397. Ackroyd, P~.T. and Goddard, A.J.H.G. (1991), Evolution of a finite element method for neutron transport, 151-161. The Mathematics of Finite Element and Applications VII, edited by J. Whiteman, Academic Press, London. Ackroyd, R.T. (1992), Foundations of finite element applications to neutron transport, to appear in Prog. Nucl.Energy. ANE 19:10/t2.1~
592
R.T. ACKROYD
Ackroyd, R.T., Abuzid, O.A. and Mirza, A.M. (1992), Discontinuous even-parity solutions for neutron transport. (to be published) Adams, M.L. (1991), Even and odd-parity finite element transport solutions in the thick diffusion limit. Conference on Advances in Mathematics, Computation, and Reactor Physics, Pittsburg, PA, April 28-May 1. Axelsson, O. and Barker, V.A. (1984), A Finite Element Solution of Boundary Value Problems, Academic Press, New York. Davis, J. (1968), Transport error bounds via PN approximations Nucl.Sci.Engng, 31, 127. Fletcher, J .K. (1983), The solution of the multi-group transport equation using spherical harmonics, Nucl. Sci. Engng, 84, 33. Galliara, J. and Williams, M.M.R. (1979), A finite element~method for neutron transport-II, some practical considerations, Ann. nucl. Energy, 6,205. Goddard, A.J.H. and Williams, M.M.R. (1986) (Editors) Finite Element Methods, 18, 1-264, Prog. Nucl. Energy. Goddard, A.J.II. (1991) (Editor) Finite Element and Allied Methods in Radiation Transport, Prog. Nucl. Energy, 25, 1-305. Lewis, E.E. (1981) Finite element approximation to the even-parity transport equation, Advances in Nuclear Science and Technology, 155-225, Editors, J. Lewins and M. Becket, 13, Plenum Press, New York. Luneberg, D. (1969), Optimization by Vector Space Methods, Wiley, New York. Nanneh, M.M. and Williams, M.M.R. (1985) , A diffusion-transport theory hybrid method for calculating neutron flux distribution in slab lattices, Atom kernenergie,47, 221. Nanneh, M.M. and Williams, M.M.R. (1986), An analytical solution to the hybrid diffusiontransport equation, Annals, nucl. Energy, 13, 337. Nanneh, M.M. (1990) A Synthesis Method Based on Hybrid Principles for Finite Element Neutron Transport, Ph.D thesis, London University, 1990. Williams, M.M.R. and Goddard, A.J.H. (1981) (Editors) Finite Element Methods in Radiation Physics, Ann. nucl. Energy, 8, 539-722.
0306-4549/92 SS.OO+O.OO Copyright Q I992 Pergmon Press Ltd
FINITEELEMENTMETHODSFORNEUTRON TRANSPORTBASEDONMAXIMUMANDMINIMUM PRINCIPLESFORDISCONTINUOUSTRIAL FUNCTIONS R.T. ACKROY D Energy Systems Section, Mechanical Engineering Imperial College, Exhibition Rd, London SW7 2BX
Department
Abstract Variational principles are given for the first order Boltzmann equation for neutron transport and its equivalent second order forms for even and odd-parity transport. For these extremum principles there are no boundary conditions to be imposed, and across interfaces the trial functions can be discontinuous. Numerical solutions are presented for the even-parity principle. The scale of interface discontinuities is controlled by a penalty parameter A. As A is increased for a given basis the degree of satisfaction of the transport equation by the optimised trial function decreases to that for a classical trial function, whereas the’interface discontinuities decrease. These trends lend themselves to a simple interpretation in terms of vectors in a Hilbert space. In applications the spatial dependence of a trial function is represented by finite elements, and a spherical harmonic expansion gives the directional dependence of the trial function. The principles are applied in the method of composite solutions in which the order of the spherical harmonic expansion can be varied according to the physical needs of the regions of a system. Large values of X are used in the method of composite solutions. On the other hand specific values of X when used with the simplest finite element lead to the well-known central difference equations for neutron transport.
1
Introduction
Finite element methods for neutron transport have achieved a geometrical flexibility comparable to that of Monte Carlo methods. The progress of the finite element method has been reviewed in of Williams and Goddard (1981), Goddard and Williams (1986), and Goddard (1991). Exposition 565
566
R . T . ACKROYD
variational treatments of finite element methods for neutron transport have been given by Lewis (1981) and Ackroyd et al (1987). The former is based on the even-parity transport equation, whereas the latter treats both even and odd-parity equations and the first order of Boltzmann transport equation. In the sequel current treatments using continuous variational solutions are extended so that discontinuous variational solutions can be used to increase the efficiency of computations. Variational treatments based on both the even and odd-parity transport equations can be used to provide efficient solutions and to establish rigorous numerical benchmarks for test problems. The maximum principles for even and odd-parity transport give independent variational approximations for the fluxes. The average of these fluxes for relatively coarse approximations can be as good as or better than refined calculation for the even-parity flux. An example of Adam's work is shown in section 7.2. The complementary principles for even and odd-parity transport can be used to establish benchmarks for both local and global characteristics of a system, because the principles are of the maximum and minimum kinds. Ackroyd and Splawksi (1982) showed that close upper and lower bounds for a local characteristic, such as the capture rate in a small region, can be obtained from complementary principles. Tight bounds on a global characteristic, viz the disadvantage factor, have been obtained by Ackroyd and Nanneh (1988). Bounds for a local and a global characteristic are illustrated in section 4.3. In almost all applications of the finite element methods for neutron transport the specification of the directional dependence of the angular flux is the same for all regions. For many applications strong directional dependence of the angular flux is confined to a fraction of the system, but its effect is significant. Then a conventional treatment, whether using spherical harmonic or discrete ordinate directional representations, is wasteful computationally. The method of composite solutions is a technique of using in a variational solution directional specifications that can be chosen independently for each region according to the physical need of the region. Thus the method of composite solutions can lead to a significant reduction in the number of coefficients to be evaluated in a variational solution. Ackroyd and Wilson (1988) implemented the method with the K +- ( ¢ + , d - ) maximum principle. The even-parity trial function and the odd-parity trial function in this principle can be discontinuous across interfaces. It is this property that permits the directional specification of a trial function in the method of composite solutions to be chosen independently for each region. The efficiency of the method is improved significantly by basing it on the K +° (¢+) principle, Ackroyd (1986), Ackroyd, Mirza, Pourzand, (1992), which uses a single trial function ¢+. In the K +° (~b+) functional there is a surface integral which is a least-squares measure of the scale of the discontinuities in at interfaces. This integral involves a weighting function and a penalty parameter A. The measure of the discontinuities for the variational solution yielded by the K +° (¢+) principle varies as 1/As. In the method of composite solutions A is taken to be large relative to a specific value of A. For this specific value of A, and a natural choice for the weighting function, the variational solution given by the K +° (¢+) principle for an infinite slab of an iostropic scattering medium is the central difference equation for neutron transport. The finite elements used are of uniform pitch and the spatial dependence of the trial function over an element is a constant. This result can be extended to a medium with an arbitrary scattering law on using the K~ (¢+) maximum principle employing discontinuous trial functions for the even-parity angular flux. The concept of complementary principles can be extended to discontinuous trial functions by introducing the K~- (¢-) maximum principle using discontinuous trial functions 4 for the odd-parity angular flux.
Finite elements for neutron transport
567
A simple means of establishing the K~ ($+) and K~" ($- )principles is to establish them, as in section 4, as special cases of the KA (¢) maximum principle for the first order Boltzmann equation. The derivation is given in section 3 for discontinuous trial functions $ for the angular flux. The I~'~ ($+) and K~- ( $ - ) principles have properties of a more general kind than the classical principles for even and odd-parity transport. This difference can be illustrated geometrically with the aid of an appropriate Hilbert space, as shown for the K + (¢+) principle in section 5. The ways in which the penalty parameter affects the variational solution, with respect to the scale of discontinuities at interfaces and the degree or satisfaction of the transport equation, are reviewed in section 6. The numerical solutions of section 7.1 illustrates the effect on flux profiles of the penalty parameter for a lattice cell. The basic theory in the sequel is for neutron transport in a typical energy group of a system with impressed sources. For simplicity the boundary condition is confined to that for a surface source. The theory can be extended readily to cover bare surface, perfect reflector and albedo boundary conditions. This exposition in the sequel follows that of Ackroyd and Goddard (1991) and Ackroyd (1992).
2 2.1
A review Boltzmann
of some
transport
first order
principles
equation
The steady state angular flux ¢o (_r,~_) at r for the direction ~ is given for each energy group by an equation of the kind n_. V¢o (_~,n_) + ~¢o (_~,n__)= S(r,n_) (1) for an absorbing and scattering system of volume V with a distributed source S. The operator is given in terms of the total cross section a(_r) and the differential scattering cross section
., (~,n_ n_n_')by
r.u (_~,n_) = ~ (_~) - f~, ~s (_~,n__.n__')u (_~,_n') dn'
(2)
On the boundary 6V of V there is a source T (_r,2) for which
¢o (_r,fl__)= T (_r,~ ) ,
n.n < 0
(3)
with L~tile outward normal to 6 V at _r. Across an interface of physically distinct regions $o (_r_r,2 ) is continuous for all those directions ~ which cross the interface at _r, but it is in general discontinuous for directions lying in the tangent plane at _r to the interface.
568 2.2
R.T. ACKROYD Parity
second
order
equations
The Boltzmann equation (1) is transformed into the pair of simultaneous equations
Q. v ~ : + c¢o+ ~ . v4+ + G - ' 40
= s+ } S-
(4)
by making the Vladimirov substitutions
¢0~ (r, ~) = [¢o (r, ~) ± ¢o (r, - ~ ) ] / 2 s + (r, n_) = [ s (r, fi) + s (r, -a_)]/2 f
(5)
The odd-parity angular flux 4,o is eliminated from the parity pair (4) to yield the second order even-parity equation - _~. v [Ga_. V¢~+] + C¢+~ = s + - a_. v G S (e) The odd-parity equation
- Q_. v [c-1_~. r e 0 ] +
G-~:
= S - - ~_. V C - ~ S +
(7)
is obtained in a similar way. The operators C and G -1 are defined by
c . (r,_~) = ~ (r). (r, ~) - / 4 . ~+ (r,_a- a_') ~ (r, a__')da'
(s)
C-~,, (,-,_~) = ,, (,.),, (,-,_~) - ~ , ,,7 (,-,~ ._~') ,, (,-, n') de'
(9)
and
Expressions are available for the operators C -x, G and ~-1 in terms of a, and a. The operators ~, C and G, and their inverses, are self-adjoint. All of these operators are positive definite. Useful identities are
F,u
= Cu + + G-lu -
~-lu
= C-Xu + + Gu-
]
i
The boundary conditions for the parity angular fluxes on 6V are
(10)
Finite elements for neutron transport
569
5Z /5'
Subregions I~,l~ etc of system volume V Subregional surfaces (i) Interior to the system , , Si for ~ , (ii) A part of the surface enclosing the system Boundary surface 6V
Sj for I~. S~ for Vi
Figure h Subregions and specified surfaces of a system
~t + G [S- - ~. V4o+] = T (r,_~) } 4I+C-*[S+--~.V4~ -] =T(_r,~)
ft.n<0
(11)
and
,/,;'- - G [S- - ~. VC, t] 4,; - C -~ [S+-- a_. v¢,;']
= T (r,-_a) } -T (r, -~_)
a"" > 0
(12)
with n the outward normal to 6V. Across an interface of physically distinct regions ~ (_r,~_) are continuous for all those directions ~_ which cross the interface at _r.
570
2.3
R.T.
ACKROYD
Classical principles
Associated with the second order even and odd-parity transport equations are the classical maximum principles K + (¢+) and K - ( ¢ - ) for ~b+ a trial function for ¢+, and~b- a trial function for ¢~-. [In the sequel ~b+ (r, fl) and ¢ - (_r,fl_) are even and odd functions of _fl respectively] The g + principle is classical because (i) the Euler-Lagrange equation for the K + (¢+) functional
V 4,r
V 4~-
+
4/r
/" in"
+ (_',m
d dS
6v _q.n_<0
_/I,._.._,(::.°,,
.,,
6 V 4T
is the even-parity equation (6); (ii) ~+ has to be continuous across interfaces in the manner (c) below. A trial function ¢+ (r_r,~_) is admissible in the i f + principle if :- (a) It is a continuous function of L' and ~_ almost everywhere within the subregions ~ into which V is partitioned by the specified surfaces of Fig.(1). These surfaces comprise the boundaries of physicM regions, the edges of neutron shadows (Ackroyd 1978), and the artificial boundaries introduced for computational convenience, eg. the faces of finite elements. (b) Within each subregion 2 " V¢ + is almost everywhere a continuous function of r and l). (c) On each interface ¢+ (r, 12) is continuous at _r for all directions ~_ which cross the interface. Note that the trial function is not required to satisfy on 6V the boundary condition for a surface source. The functional K - (~b-) for the K - maximum principle is
V 4~r
-
dV V 4~r
+
4f f
6v _n.n
-
4//,~_.nl(¢-)2dfmS
(14)
6V 4~t
The conditions for ¢- to be admissible in the K - principle are of the same type as the conditions (a)-(c) above for ~b+ to be admissible in the K + principle. For finite element work a trial function ¢+ or ¢ - is expressed as a linear combination of basis functions. Each basis function is the product of a function of position _r and a function of direction l). The functions of position are the shape functions for the finite elements into which the volume
Finite elements for neutron transport
571
V has been partitioned. Either spherical harmonics or discrete ordinates are the usual means for the specification of the functions of direction. Spherical harmonics have the advantage of being the analytical companions of the Legendre functions used to specify the differential scattering cross section. Even spherical harmonics are used for ¢+ and odd spherical harmonics for ¢ - , and either requires half the number of spherical harmonics to specify the angular flux ~b= ¢+ + ¢ - . In the maximum principles g +(¢+) K - (¢-)
_
(15)
equality holds if, and only if, the trial function is identical with the corresponding parity component of ¢o. If ¢+ is the even-parity trial function that maximizes K + (¢+) for trial functions with a given basis then the corresponding odd-parity angular flux approximation is usually taken to be
(16)
:
by analogy with the exact result (4). Similarly for a maximizing trial function ¢ - the corresponding even-parity angular flux approximation is (17)
=c-'
The principles K + and K - are complementary in the following sense. It can be shown that
[[{s+c-is + +s-Gs-}d v
:
V 4~
=
2a
(say)
(18)
On defining 2/3 = K + (¢+) - K - (¢~-)
(19)
there follows for the functionals K + (¢+) - a and a - K - ( ¢ - ) the complementary principle g + (¢+) - ~ < ~ < ~ - K - ( ¢ - )
(20)
with the equalities holding if, and only if, ¢+ = ¢+ and ¢ - = ¢~-. This result is shown to be of some value in section 4.3 for the benchmarking of transport solutions.
2.4
A neo-classical
principle
Although the Euler-Lagrange method is restricted to differential equations of even order, Axelsson and Barker (1984), a maximum principle for the first order Boltzmann equation can be obtained, Ackroyd (1983 a, b), using a generalized least-squares procedure. This method is illustrated in section 3.3 by deriving a maximum principle for the first order Boltzmann equation for which discontinuous trial functions are admissible. Here the neo-classical principle f~" for the first order transport equation is stated.
572
R . T . ACKROYD
Special choices for the continuous trial function 6 yield the classical K + and K - principles. The continuity conditions for 6 are of the same kinds as those for 6 + in the K + (6 +) and K - (6-) principles. T h e / ~ (6) functional is
V 4=
-- //{(~-.~'V6)r-1~--'V6-6r6}d a d V v 4~t
6 v fi-n_<0
-
4 / / I n , a l ( 6 ) 2 dadS
(21)
6V 4~
The maximum principle is
(6) _< k (6,) = 2~
(22)
with equality holding if, and only if, 6 = 60. Substituting the special choices of trial functions
(a)
6
(b)
6 -- 6 + + 6 -
=
6 + + 6; (23)
in the/~ (6) principle leads to the K + (6 +) and K - (6-) principles respectively. This simple means of deriving complementary principles is illustrated in section 4 by deriving maximum principles for even and odd-parity transport that accept discontinuous trial functions.
2.5
P r i n c i p l e s for d i s c o n t i n u o u s
trial functions
K (6) principle A generalization of t h e / ( ( 6 ) principle is the K (6) principle
g ( 6 ) < 2~
(24)
Admissible trial functions 6 can be discontinuous across interfaces, and they do not have to satisfy the boundary condition. Equality holds if, and only if, 6 = 6o. The R"(6) functional is given by
USi 4~¢
u ( . % n S j ) 4~-
Finite elements for neutron transport
573
Jumps at r in the trial functions are ~:k ( ~ , ~ ) _ ~- (rt,fi)
~
" ~
Interface S../'IS.
Figure 2: Jumps in trial functions at r on the interface S~ FI 8j of elements I~ and where the definition (21) for/{ (~b) is extended to discontinuous trial functions. The USi is the set of all element faces which lie within the volume V, and n~ is the outward normal to the element surface Si. The set of all element interfaces is U(S~ NSj) and ~ (~.~,i~) - ~ (r~,fl) is the jump at r on the interface St N Sj in the trial function in going from the element with volume I~ to the element with volume Vj in Fig(2). The weight ~7 (£, ~) is non-negative and is usually chosen to be positive for ft. nl ~ 0, where n= is the normal to Sl. The penalty parameter A is positive. Rules for selecting ~17 and A are obtained in the same w~t.vas the rules derived in section 6 for the maximum principles for even and odd-parity transport which agimit discontinuous trial functions. K +- (~+, ~ - ) principle The K +- ( ~ + , ~ - ) principle can be used to derive all of the other principles mentioned in this section. In the principle ~+ and ~ - are trial functions for ~b+ sad ~" respectively. They can be discontinuous at interfaces, and they need not satisfy the boundary condition. For the principle
K+-
<
equality holds if, and only if, ~b+ = ~+ and ~b- - ~ .
The functional K +- (~+, ~-) is
= K÷(,+)+K-(,-)-./ f t.t$ i 4 w
u(s~ns~) 4.
with l ~ (_r,fl) non-negativeweights,usuallychosen to be positiveforQ . m # O. Putting
=W-=W
(26)
574
R . T . ACKROYD
=
+
(28)
gives the I'~ (~) principle. Using the substitutions (28) and letting A ~ oo gives t h e / ~ (~) principle. The substition = ~bo (29) transforms the principle (27) into the even-parity principle, Ackroyd (1986), r ~+° (~+) < K+ (~o+)
(30)
with equality holding if~b+ = ~b+. The way in which this principle operates can be illustrated geometrically by means of a Hilbert space, in very much the same as that give in section 5 for the related even-parity maximum principle K + (~b+). Composite solutions Ackroyd and Wilson (1988) used the K + - (~b+, ~ - ) principle for the method of composite solutions, with A large relative to a specific value for this penalty parameter. For the specific value of A , and a uniform mesh of the primitive elements of section 6.2, the optimised trial functions satisfy the central difference equations for neutron transport in an infinite slab of an isotropic scattering medium. T h e drawback of the K + - (~b+,~b- ) principle in the use of two trial functions is removed for the K +° (~b+) and K ~ - (~b') principles. The K + (~b+) principle is derived in section 4.1 It is superior to the K +° (~+) principle, because with the same kind of trial function for an infinite slab, and the same specific value of A , it gives the central difference transport equation for an arbitrary scattering medium. The K + (~+)principle is illustrated in section T with numerical solutions for a lattice cell.
3
3.1
Maximum tion Error
p r i n c i p l e for t h e first o r d e r t r a n s p o r t
equa-
measures
In constructing a maximum principle some measures are required of the errors made by a trial function ~ (_r, ~_) for the angular flux #0 (__r,~_). A measure of the error made by ¢ in not satisfying the transport everywhere within the volume V, and for all directions, is defined to be
V 4r
Since ~-1 is a positive definite operator, the volume error e~ vanishes if, and only if, ~_- V~b + E~b - S = 0
(32)
Finite elements for neutron transport
575
ie.ff is a solution of the transport equation. The failure to satisfy the boundary condition given is measured by
~,~-2/ / In.fil[~(r.,fi)-T(,',fi)12df~dS
(33)
6V fl.n_<0 where tile normal n is outward to the boundary surface 6V. To define the error incurred by the failure of ~bto satisfy the interface continuity condition, consider in Fig 2 a typical interface Si f3Sj of two spatial finite elements V~ and Vj. For a point r on S, N Sj let ~ (r~, f/) denote the value of the trial function as given by its expression for the element ~ . Similarly ~b(_rj, fl_) is the value of the trial function at r as given by the trial function's expression for I,~. The error incurred at the set U (Si N Sj) of all interfaces is specified as
¢/=A
f
{fl •n~{_~..~ [~ (r~,f~) - ~b(_r,,fl)] } ×~_1 [~(~,~)_ ~ (_.~.~)] }
f
dfldS
(34)
u(S,nS~) 4.
where ni is the outward normal at r to Si. The penalty parameter A serves to bias the interface error relative to the other errors. The cumulative error is v~ (~)
=
S])dV
/(ft. v ~ + r.~ _ s, ~-1 L~. v ~ + ~ V
6v ~-n_<0
f f
{"-"-[*("'"-)-*('-""-)1} ×r.-, {._... [~ (.., a_) _ ¢ (_._.._)] }
o(s,.s,)'For brevity the notation
f f uvdfldV = f(u, v)dV is used in the V 41r
sequel. The functionl Vx (~b) has
V
the property that Vx(if) >_ 0
(36)
because )"~--1 is a positive definite operator.
3.2
A least-squares
identity
Expansion of the least-squares error functional Vx (O~) gives the identity Kx(~b)+Vx(~b)
=
/(S,E-1S, dV + 2 / V
=
/
[f~'n[T2dfldS
6V fl.n_<0
2a
(37)
The functional Kx (d) can be expressed as
K~ (~) =
[ [2(~,s) + 2(ft.v~, r.-'s)]dv - P (~,~) V
6V fl.n<0
USi 4~r
576
R.T. ACKROYD
where the bi-linear functional F (0, ~b) is defined as ~' (0, ~b) --
/ { ( ~ . V0, ~-]-1[~. V~) "dr (0, Y]~b)} dV v
+ f/[~,nlO~bdf~dS=
F(~b,0)
(39)
6V 4x
The expression (38) is obtained by using the divergence theorem to convert the term 2f(ff, f~ • V V~)dV=f f div [Q_~b2] d~dV in the expansion of the volume component of V into surface integrals v 4~ for USi, the set of all element faces inside V, and the boundary 6V.
3.3
The
K~ (~b) m a x i m u m
principle
It follows from the definition (35) that Vx (~b,) vanishes to give
Kx (¢o) = 2~ The condition
(a)
K~ (¢0)
implies (b)
(40) = 2~ ] /
~
(41)
= ~0
as follows. Putting
x (r, n__)= ~ (_~,9_) - ~0 (r, 9_)
(42)
the condition (41a) and the identity (37) give
0= Vx(if)= f [(_Q.vx, x~-lfl__,vx)+ (x, ~cx)] dV+ f f I,,nlxZdftdS V
(43)
6V 4x
on expanding the volume integrand of Vx (~b), noting that ~-1 is self-adjoint, and applying the divergence theorem to the term f f [2. ~7(X~) dfldV. The expression (43) implies the vanishing V ,ix of X, ie.~b = ~b,, because E and its inverse are positive definite operators. The conditions (40) and (41) and the property (36) establish the maximum principle Kx (~b) _< 2~
(44)
with equality holding if, and only if, ~b= ~b. If ~b(r, _~) is taken to be continuous across each interface, for all directions D crossing that interface, then Kx (5) =/~" (~b) (45)
Finite elementsfor neutrontransport
577
where
k(¢)
=
f [2(¢,s)+2(~.v,,r='s)]d r -
P(¢,¢)
V
SV N-n_
and
k (~b)_<2a
(47)
is the maximum principle of Ackroyd (1983 a,b). For this principle admissible trial function have to satisfy the interface continuity condition.
3.4
M a x i m u m principle in t e r m s of the parity c o m p o n e n t s of the trial function
On expressing q~as the sum of the parity components ¢+ and ¢- the functionals Vx (¢) and Kx (~b) become, with the aid of the properties (10) for I[i-1 and E ,
:
+(_~. v¢+
f{
+ G-l,-
- S-.G
[ft. V, + + a - ' , - - s - l ) }
dV
V
6V fl.n_.< 0
( {ll. ~ [~+ (,~,a)- ~+ (~i,a)] } +
A
f
/
)
x G { a ' n ~ [ ¢ + ( r - / ' a ) - ¢ + (_r,,f~)]}
+ {a-~ [¢-(,~,a)-,- (,_~,a)] }
dadS
(48)
×c-1 {~. ~ [,- (~,a)_,- (,~,~_)1}
o(s,osj),,
and
+,-) =
K+(,+)+K-(,-)-2f U S i 4~"
_ ~ f
f
o(~,°~,)'-
×a{,.o_-.~ [,+(r__i, (r~,_,)-, (~_,,._)1} _i_ { fI. n_./[,;b_ a) _ ~;b_+(r_./, a) ] } ×c-' 1~,_.~_,[,- (~..,_)- ,- (~_,a)] }
dadS
(49)
578
4
4.1
R.T.
Complementary port Maximum
principles
ACKROYD
for even and odd-parity
principle for even-parity
trans-
transport
For the special choice of trial function
¢ = ¢+ + ¢:
(50)
the functions (48) and (49), taken with the parity equations (4) and the boundary conditions (11) (12), give
V~ (¢+ +¢2)
=
r + (~,+-¢+,,~,+-¢t)
+
A
f f
G{~'n~[dz+(r~'~)-ek+(£J'~)]}× { a ~ [~+ (~,~) - ~+ (_~,~_)] } dadS
(51)
u ( S , n S i ) 4~r
and Kx(¢++¢:)
-
K-(¢:)
-
[J
=K +(¢+)
( ~ {a ,~ [~+(~,a)- ~+ (_~,a)] (-~,-~)]} }
[J
×al._.n~ [~+ (.~,._)-~+ +2~. ~ . [¢+ (a,~) - ¢+ .(,~,~)] ¢: (,,~) . .
u ( S , n S ~ ) 4~r
)
d~dS (52)
Here F + (u +, v+) is the bi-linear functional
r+ (u+,v+) = [ {(~. w+,G~, w+) + (~+,Cv+)} av v
+ f f I"-.nl,,+v+d.dS
(53)
6V 4~
with the property that F + (u+,u +) is positive definite. The bilinear function F - ( u - , v - ) is
F-(u-,,~-) = [ { l ~ v u - , c - ~ v,,-)+lu-,C-~v-)}dv v
8 V 41r
with the property that F - (u-, u - ) is positive definite. In obtaining the last result from the expression (49) the integral over US/, the set of all element faces within V, has been written as an integral over the set of all interfaces U(Si I-I
St)
On making use of the results F + (¢+,¢+) g-(¢~')
+ F - (¢~-,¢~-) = 2a =F-(¢0,¢o)
f
(55)
Finite elements for neutron transport
5:9
Ackroyd (1981), the expressions (51) and (52) transform the identity (37) into
Kt(* +)
=
F+ (~+o, <'t) + 1
fe_:G_l~Zaaas
[-~u(s, nsD 4=
(56)
where K +(~b +)
=
K +(¢+)
a {~.~ {¢+ (a,~)- ~+ (a, a) }} U($1rlSj)
4~
Since G and C are positive definite operators the term in parentheses in the identity (56) is nonnegative. Hence the inequality K + _< F + ( ¢ + , ¢ + ) +
1 -~ f f ¢oG-l¢'~d~'ldS
(58)
u(s,ns~) 4~r
with equality holding if ¢+ = ¢+
(59)
For then the term in parantheses in the identity (56) vanishes on account of the interface continuity condition for ¢~. On the other hand the condition
'//
K + (¢+) = F + (¢+, ¢+) + ~
¢:a-l¢'[dgldS
(60)
u(s~ns,) 4r can hold only in the limit ~ ---* oo. If it held for finite ~ the the following conditions would hold viz r + ( ¢ + - ¢+, ¢+ - ¢ + ) = 0
(61)
and
Condition (61) would imply ¢+ = ¢+, and then the continuity ¢+ at interfaces would imply that the condition (62) did not hold for ~ finite. Nevertheless maximizing K + (¢+) with A finite can bring ¢+ close to ¢+. This is shown by the geometeric interpretation of the K + (¢+) principle in section 5
580
4.2
R . T . ACKROYD
Maximum
principle
for odd-parity
transport
The K ~ - ( ¢ - ) maximum principle for odd-parity transport is obtained in the same way as the K + (¢+) maximum principle. The result is K~" ( ¢ - )
=
K- (¢-)
c-,
u(S,nSj4,)
b-
}
U(SinSj) 4 x The out of balance in this inequality is
{..n~[¢-(r~,.)-¢-(rj,~)] } F - (~b- - ¢0, ¢ - - ¢o) + A
o,s,ns,>.~, c-1{
4.3
Complementary
a--'"~[*-(~'}a-)-¢-e-b' _ c ~n--)] ;
d~dS
(64)
principles
For the classical K + (~b+) and K - ( ~ - ) maximum principles there holds the result of Davis (1968), which can be expressed as K + (¢+) - a
= F + (¢+,~b +) + F - (~bo,~b~-) = F + (¢+,~b +) - F - (¢,,~b~')
[
(66)
Corresponding to the classical complementary principles K + (¢+) - a and a - K - (~b-) there are the complementary principles K + (5+) - a' and a' - K~" (~b-) with the property
g + (¢+) - a' _
(67)
Here 2a'
=
F + (¢+,¢+) + F -
1
+ -~ f
(¢0,¢~-)
f[*+*C¢+o+¢ZG-'¢:]dadS
u ( s , o s , ) 4~ 21~'
=
F+
+
+) - r-
1
o(s,osD4. The classical complementary principles have been used for benchmarking finite element solutions in two ways, viz
Finite elements for neutron transport
58l
,, •
flux',
O(x) I
I I
',.t
fuel
', s=so l
void
/
i o=s0 1o~1
o=0
I
~'-
moderator
's=~,
,o=1 ot o~0.9 =1,
\ o~0.9
\
l i
:perfect rellector ~
O
I
I I.
0
2
/: J
: I
I
3
5
base surface
6
x
8
Figure 3: Solution of the edge cell problem 1. to place bounds on a local characteristic of a solution 2. to place bounds on a global characteristic of a solution
Bounds for a local characteristic
For the edge-cell problem of Reed (1971), Fig 3. the capture rate of neutrons in the can region has to be determined accurately for it to be used as a benchmark. The flux in the numerical solution of Galliara and Williams (1979) varies rapidly in the can. The bi-variational method of Ackroyd and Splawksi (1982) shows that the capture rate is in the range 0.7855 + 0.0036, ie. it is known to within + 0.46%. This example shows that classical complementary principles can be used to give benchmarks of worthwhile precision. Bounds for a global characteristic Many numerical methods for the disadvantage factor for a reactor cell converge slowly, as shown by table 1. This global characteristic is the ratio of the average flux in the moderator to the average flux in the fuel. The results of table 1 refer to a uranium/water slab lattice with the thicknesses of the fuel and water regions at 0.167 and 1.5 mean free paths. The juxtaposition of the uranium strong absorber and water scattering medium give rise to significant transport effects. I Method ANISN-S8 ANISN-SI6 RIPPLE DI FTRAN FELQO EK+ t FELQOE K - t
fuel]
moderator
Average Average Disadvantage 1 flux flux factor 3.00081 3.24259 1.08057 2.99733 3.29745 1.10013 3.00980 3.31400 I.I0698 2.99466 3.34050 I. 11599 2.99390 3.35246 1.11976 t29 moments, 8 quadratic elements mean=l.12087 2.99351 ] 3.35865 1.12198
Table 1: Numerical results for the average fluxes and disadvantage factor using 8 meshes
582
R.T. ACKROYD 1.6
[,
1.5
H
o
1.4 , ,
E v e n - parity
z~ O d d - parity
~~ . 3 ~ ~1.2 c21 1.1 1.0
0
5
10 15 20 25 Number of moments
J
30
Figure 4: Comparison of K + and K - flux plots using 8 quadratic elements for a uranium/water lattice cell Ackroyd and Nanneh (1988) have shown that the K + and K - principles give lower and upper bounds respectively for the disadvantage factor. The table shows that the results for the well known discrete ordinate code ANISN, the collision probability method RIPPLE , and the elegant DIFTRAN method of Nanneh and Williams (1985), (1986) have not converged fully. The K + and K - bounds for the disadvantage factor as functions of the number of Legendre functions, used to represent the directional dependence of the neutron transport, are shown in Fig 4.
5 5.1
Geometrical interpretation of the even-parity maximum principle A Hilbert
space
A geometrical interpretation of the effect of maximizing K + (~b+), with respect to the coefficients specifying the trial function ¢+, can be given using a Hilbert space that employs an even and an odd function of _~ to define a vector of the space. I f O ~ [0+,0 -] and ~ ~ [ ¢ + , ¢ - ] are a pair of arbitrary vectors of the Hilbert space associated with the K + (~b+), principle then their scalar product is defined to be
o. • = [ [(~_.vo+,G£, v¢+) + (0+,C¢+)]dV v
+//,o 6V 4.
{
]/
×
The even-parity function 0+ (_r, f/) is admissible in the K + (¢+) prinople. The odd-parity function ¢ - (_r,~) is admissible in the classical K - (¢-) principle. Whereas 0 + can be discontinuous across interfaces ¢ - is required to be continuous at interfaces for all directions crossing the interface.
Finite elements for neutron transport
5~3
The definition (69) satisfies the commutative, distributive, associative and multiplication by zero laws required of a scalar product. The commutative law holds because G and C are self-adjoint operators, and the distributive law is satisfied since G and (7 are linear operators.
5.2
Geometrical
transcription
of a variational
identity
The identity (56) can be expressed using the scalar product (69) as K + (~b+) + ((I)+ - @0+- ~o-)2 = (~_o + +~_o-) 2 with • + . . [~+,01,
~
~ [~0+,01
and
2 - ~ [o, ~:]
(70)
The results (¢__o+)2= F+ (¢+, ¢ + )
(71)
f f *z~'¢Za.dS
(2-)2=
(72)
u(s,nsO 4~
~.
2- = 0
(73)
are established readily. T h e functional K + ( ¢ + ) can be expressed in terms of ~ + , ~ +
and ¢--o- by
means of the cosine rule identity viz (*+) 2 + (~+ - [~__o++ ~__o-])2- 2~+ • (~+ - [~__o++ ~_.o-])2 = (~__o++ ~__o-)2
(74)
This identity is recast by using the result (73) as
2_~+ • ( ~ + ~ - ) - (~_+)~ + (_~+ - [ ~ + ~ - ] ) ~ = ( ~ ) ~ + ( ~ - ) 2
(7~)
Comparison of the results (70) and (75) gives K~+ (¢+) = 2 ~_+. ( ~
5.3
Geometerical
interpretation
+_~:) - (~_+)~
of the variational
(76)
principle
For a typical trim function expressed as M
¢+ (_~,_a) = ~
am + ~m + (-~, a__)
rn.~l
in terms of the basis functions ¢+ there corresponds the vector M
,__+
~ m=l
+ +
(77)
584
R.T. ACKROYD
w h e r e ~ ~ [~+, 0]. Maximizing K + (~b+) is according to the identity (70) equivalent to minimizing (~+ - ~
- ~o-) 2 with respect to the
a +. Thus (f
-*,+ -~-)
=0
(m= 1,...,M)
(78)
where --
=
z-.,m=1
m-~-m
~
f
(79)
and the~+~ are the values of the coef~cients .+ which m~ximize K~+ (~+). Summing for all m the product of ~+ and the result (78) gives
• ÷. Thus defining the target vector ~
=0
<80)
to be = ~
+ ~--o- ~ [~+, ~0]
(81)
it can be said that in seeking the optimum approximation vector _~+ *-* [~+, 0] by the minimization of (~+ - ~__o)2 the target error vector ~+ - ~
~_+
is orthogonal to ~+, i.e. -
= 0
This result is illustrated in Fig 5. The geometrical effect of seeking the optimum approximation d+ by maximizing K + (~b+) is to determine a point P of the linear M space (77) which is nearest to the unknown point Q specified by the target vector ~--o- Note that the definition (81) and the result (73) give (*0) 2 = (~_o+)2 + (~_o-)2 (83) Since (~__o-)2 is given by the expression (72) the distance QP ( ~ 2 of the target ~ vector from the exact solution vector ~+ varies inversely with A . As A --, oo the magnitude of the target vector shrinks to the magnitude of the exact solution vector, and the target vector aligns itself with the exact solution vector. Then the geometrical interpretation of the maximum principle becomes the classical interpretation, Luneberg (1969), as illustrated for the K + (~+) principle by Ackroyd (1986).
6
6.1
The
role of the
penalty
parameter
A
Effect of increasing
The effect on variational solutions for neutron transport is described here for solutions obtained with the use of the Kx (~) principle. The relevant properties of the Kx (~) functional can be
Finite elements for neutron transport
__
585
g+ - ( @^+ - % ) = 0
Error Vector
Og
Offset Vector
cg- ®o--o
Figure 5: Solution, target and optimum approximation vectors established in the way used by Ackroyd and Wilson, (1988), to obtain the corresponding properties for the K +- (~b+,~b- ) principle. The cumulative error functional Vx (~), (35) is written as
vx (¢) = i, (¢) + ~i~ (¢)
(84)
where
i,(,) = f f
V*+ *-SldadV
V 4¢
+ 2/
f
Ifi.n_l[7"-,l'dadS
(85)
tv fi.n_.< o
is a measure of the failure to satisfy the Boltzmann equation within elements and to fulfil the boundary condition, and ld (~) =
f
/
{~ "n--'i [~ ( ~ ' ~ ) -- ~ (-rJ'~)] } x
dfldS
(86)
u(s~nsi) 4,~
is a measure of interface discontinuities The basis for the trial functions ~ is assumed to contain a sub-set capable of representing classical trial functions satisfying the continuity condition for interfaces. Let ~bx be the trial function which maximizes Kx (~b) then lr (~bx) is monotone increasing and Id (~bx) is monotone decreasing with A . A s A --- oo, lr (~bx) --~ I, (~b¢) where ~c is the classical trial function which maximizes Kx (~b). Thereby ~c maximizes the functional/((~) of the maximum principle for the Boltzmann equation
586
R.T. ACKROYD
which admits only classical trial functions. Since ~, is not the maximizing trial function of KA (~) for trial functions which in general are discontinuous / i (¢c) = KA (~,) < KA ($A) The identity (37) gives KA (CA) + Ir (~bx) + Aid (~bA)= 2 a = KA (¢~) + I~ (~bc) Hence the simple inequalities I~(¢~)
<
It(Co)
<
ir (¢o)
(87)
A
indicating the trends. As )~ is increased the measure Ix (¢~) of the interface discontinuities is decreased, whereas the degree of the failure to satisfy the Boltzmann equation within elements and fulfil the boundary condition, as measured by I~ (¢~), is increased. These remarks suggest that there may be some specific values of A for which a discontinuous variational solution is superior overall to a classical continuous variational solution. In the next section a good variational solution is obtained for a specific value of A and the simplest discontinuous finite element. This specific value of A can be used as a guide in selecting an appropriate penalty parameter for a particular application.
6.2
Specific penalty
parameter
A specific value of the penalty parameter A is obtained for the K + (¢+) principle by considering the neutron transport induced in a homogeneous infinite slab of an anisotropic scattering medium by a source varying over space in a piece-wise constant fashion. For simplicity the source is taken to be isotropic and of strength S, a constant over the mth spatial element (Zm-l,Zm) of the slab. The primitive elements employed are assumed to be of constant pitch h. The shape function for a primitive element is unity within the element and zero outside. The trial function¢ + (x, p)is given, for the mth element by = am
= am
+
(88)
where cos - l p is the angle made by the direction ~_ with the x-axis. The functions am (~) maximize K~ (¢+), e is a parameter and e0m (ju) is an arbitrary variation of the trial function for the m 'h element. In taking variations of K + (¢+) the trial functions can be of the kinds ¢+ (x,/~) = ( tim (#) + tO+ (/J) tim, (/~)
for element m for all elements m' # m
(89)
because with no interface conditions to be satisfied the variations for the elements can be made independently of each other.
Finite elements for neutron transport
587
Let K+m (@+) be the contribution to K + (@+) that involves O+, then the general formulae (38) and (39) give for m an interior element of the slab 1
K+x(@+)
=
2.
- [ 6 r e + e 0 +] C [ a m + e 0 +]
dr
--1 1
-~,[am+eO~-am_,]G(-~la,,, +eO+-am_,]) } dl,
_
(90)
-1
since ~ - n / =
p for the interface Xm and ~ .n_/= - # for the interface Zm-1.
When K + (¢+) is maximized for the trial function (89) then d
+
i.e. 1
o =
h / o+,t [s~* - Cam] d. -1
-,
{
. a ( . [ a m - am+,]) -~G ( - ~
d.
(91)
[am - a~_~])
since C and O are self-adjoint operators. As 0m + can be arbitrary the condition (91) gives
_~. { c ( . [am+, a m- l ) h _ C (. [amh a~-,l) } +Cam =S+ Choosing 1
A= g
(92)
gives the central difference scheme -#
{G(p[6m+l-am]/h)-G(P[fim-am-1]/h)}+Ofim=S+
(93)
h
for comparison with the even-parity equation
for the even-parity angular flux ¢+ (x, p) in the slab with an isotropic source S +. Central difference methods have been used by Fletcher (1983) to obtain good solutions for neutron transport in systems coverable by a regular spatial mesh. In his work Fletcher has used a spherical harmonic expansion to treat the directional dependence of the even-parity angular flux. The use of the K + (@+) principle with the simplest spatial finite element, viz the primitive element, has lead to the good approximation of the central difference equation (93) for the transport equation (94). This suggests that with more refined spatial elements than the element the K + (¢+) principle has potential for producing better solutions.
588
7 7.1
R.T.
ACKROYD
N u m e r i c a l solutions K~+ (~b+) s o l u t i o n o f t h e R e e d
problem
The way in which the penalty parameter A influences the convergence of K~+ (~b+) variational solutions is shown in Fig.6 for the Reed problem. The version of the K~+ (~b+) principle used for this example has an extra term for the perfect reflecting surface, and the bare surface is accounted for by setting T to zero. The extra term is -A //(
{~" n [~b+( r , ~ ) - ~b+ (_r,fl__')]}. ×G fl) )]}
)d~dS
(95)
~V ~ 4x
where dfV~ is the additional surface presented by the perfect reflector, and 2" is the reflected direction to the incident direction fl__at r.r_. In the results of Ackroyd, Abuzid and Mirza (1992) for the Reed problem the number of Legendre functions used per region for the directional dependence of the trial function ~b+ has been chosen as follows. The flux over most of the fuel region is expected to be flat, because the source is constant over a region several means free paths thick. Two moments are probably sufficient here for ~+. With a constant flux expected for the void, a high flux induced locally by the moderator source, and a flux depression int he source-free strongly absorbing can, the great variations in the flux necessitate a fairly high Legendre expansion of 4 moments for ~+. In the far region of the moderator without a source there is a short distance from the moderator region with a source to the bare surface. Five moments are perhaps adequate for the far region of the moderator.
7.2
Averaging
results from complementary
principles
Disadvantage factors For the uranium-water lattice of section 4.3 the asymptotic value of the disadvantage factor is 1.20087. The error of the mean value of the K + and K - results for the disadvantage factor is shown in table 7.2 as a function of the number of Legendre moments employed. Thus the mean of four moment K + and K - calculations for the disadvantage factor is better then either an $16 discrete ordinates calculation or a RIPPLE result with errors at 1.78% and 1.17% respectively. The good results for the mean disadvantage factor for four or more moments is a consequence of an approximation for the flux which is provided by the mean of the K + and K fluxes. Flux profiles Problems arising in the transport of thermal radiation are concerned usually with optically thick diffusive regions. For practical reasons the spatial mesh employed has to be large compared to the mean free path for such regions. Adams (1991) has solved some typical problems using both
Finite elements for neutron transport
589
~lGo
2.4-
.o g
JXD 1'3 LL
rr
1.2,
L--~Z.'~--7,.:
<
•
09
\", 0.6.
U
\
X CELL WIDTH (CM) Figure 6: Effect of penalty parameter A on K + solutions
Number of Moments 1 2 3 4 5 10 15
Percentage Error of mean disadvantage factor 14.870 4.850 2.270 0.890 0.790 0.130 0.024
Table 2: Percentage error of mean disadvantage factor
20t •
~"
Fluxes
1.5
(rdmm~) (~) (odd) (ave=go)
8 n--o
~¢¢=¢L----.,. =~l . ° r
o.5 ---~s,~ 0.0 1 0
2oo
x
400
6oo
Figure 7: Averaging of parity fluxes (Adams 1991)
590
R . T . ACKROYD
even and odd-parity finite element transport solutions. He has considered a two region problem for an infinite slab, with the incoming flux at the left-hand face of the slab having delta function of direction about p = 0.1. The left hand region of the slab is a purely absorbing skin with a thickness of a thousanth of a mean-free path. The second region is a source-free purely scattering medium some 600 means free paths thick. A reference solution was obtained for this problem using a linear-discontinuous diseretization of the first order transport equation with a very fine spatial mesh. For the even-parity and odd-parity transport coarse mesh solutions ten elements used were for the second region. In the thin region, which is well resolved spatially by the mesh employed, the even and odd-parity solutions for the flux are well out with the reference solution. As shown in Fig.7 the average of the even and odd-parity solutions is in good accord with the reference solution. The average used does not have to be an unweighted average. Both Nanneh (1991) and Adams (1991) have used in effect
~+(z,#)=7 ~+(z,#) +(1-7) average
derived from ~b-
~+(z,#)
(96)
e v e n - p a r i t y solution
with 7 a coefficient to be determined in some way to improve ~b+ average. The even-parity derived angular flux is given in terms of if- as follows
derived
This component for the even-parity angular is in general discontinuous at interfaces The coefficients 3' can be determined, as an independent coefficient for each element, by substituting an expression of the kind (96) for each element in the Kx+ principle, with A given by the specific value given in section 6.2.
7.3
Some observations
There is the prospect of applying the complementary principles, which admit discontinuous trial functions, in the same ways as those found to be useful for the classical complementary principles, viz to obtain variational solutions for practical problems, and to provide benchmarks for test problems by constructing upper and lower bounds for both global and local characteristic. In seeking variational solutions there is one feature in which the K + principle has more flexibility than the K + principle. The optimized solution for the latter has the merit of satisfying neutron conservation for the whole system. For the former principle the trial function can be specified independently for each element. Consequently neutron conservation can be imposed on each element, ie. for each element V/
/ / {-~_" V [G~" V~b+]+ C~+ - S+ } df~dS = O V, 4~"
(98)
Finite elements for neutron transport
591
These conditions serve to impose neutron conservation for each dement, by making that coefficient in the trial function $+ which specifies the element flux dependent on the other coefficients for that element. In seeking variational solutions more efficient than the classical eontinuoua solutions attention has been focused on the method of composite solutions and the averaging of solutions; but the variational principle K + ($+) can he used as the basis of 'nodal' solutions, with neutron conservation for each element or 'node' ensured by the condition (98).
References
Ackroyd, R.T. (1978), A finite element method for neutron transport-I, some theoretical consideration. Ann. nuei. Energy, 5, 75. Ackroyd, R.T. (1981), The why and how of finite elements. Ann. nucl. Energy, 8, 539. Ackroyd, R.T. and Splawski, B.A. (1982), A finite element method for neutron transport-VI, upper and lower bounds for characteristics of solutions. Ann. nucl .Energy, 9, 315. Ackroyd, R.T. (1983a), Least squares derivation of extremum and weighted residual methods for equations of reactor physics. Ann. nucl. Energy, 10, 65. Ackroyd, I~.T. (1983h), A finite element method for neutron transport-VII, completely boundaryfree maximum principle for the first order Boltzmann equation, Ann. nucl. Energy, Lfr, 243. Ackroyd, R.T. (1986), Generalized least squares as a generator of variational principles and weighted residual methods for FEM transport methods. Prog. Nucl. Energy,18, 7. Ackroyd, R.T., Fletcher, J.K., Goddard, A.J.H.G., Issa, J., Williams, M.M.R. and Wood, J. (1987), Some recent developments in finite element methods for neutron transport theory, 381-383, Advances in Nuclear Science and Technology, 19, edited by J. Lewins and M. Becker, Plenum Press, New York. Ackroyd, R.T. and Nanneh, M.N. (1988), Upper and lower hounds for disadvantage factors as a test of an algorithm used in a synthesis method. Ann. nucl. Energy,15, 241. Ackroyd, R.T. and Wilson, W.E. (1988), Composite finite element solutions for neutron transport, Ann. nucl. Energy, 15,397. Ackroyd, P~.T. and Goddard, A.J.H.G. (1991), Evolution of a finite element method for neutron transport, 151-161. The Mathematics of Finite Element and Applications VII, edited by J. Whiteman, Academic Press, London. Ackroyd, R.T. (1992), Foundations of finite element applications to neutron transport, to appear in Prog. Nucl.Energy. ANE 19:10/t2.1~
592
R.T. ACKROYD
Ackroyd, R.T., Abuzid, O.A. and Mirza, A.M. (1992), Discontinuous even-parity solutions for neutron transport. (to be published) Adams, M.L. (1991), Even and odd-parity finite element transport solutions in the thick diffusion limit. Conference on Advances in Mathematics, Computation, and Reactor Physics, Pittsburg, PA, April 28-May 1. Axelsson, O. and Barker, V.A. (1984), A Finite Element Solution of Boundary Value Problems, Academic Press, New York. Davis, J. (1968), Transport error bounds via PN approximations Nucl.Sci.Engng, 31, 127. Fletcher, J .K. (1983), The solution of the multi-group transport equation using spherical harmonics, Nucl. Sci. Engng, 84, 33. Galliara, J. and Williams, M.M.R. (1979), A finite element~method for neutron transport-II, some practical considerations, Ann. nucl. Energy, 6,205. Goddard, A.J.H. and Williams, M.M.R. (1986) (Editors) Finite Element Methods, 18, 1-264, Prog. Nucl. Energy. Goddard, A.J.II. (1991) (Editor) Finite Element and Allied Methods in Radiation Transport, Prog. Nucl. Energy, 25, 1-305. Lewis, E.E. (1981) Finite element approximation to the even-parity transport equation, Advances in Nuclear Science and Technology, 155-225, Editors, J. Lewins and M. Becket, 13, Plenum Press, New York. Luneberg, D. (1969), Optimization by Vector Space Methods, Wiley, New York. Nanneh, M.M. and Williams, M.M.R. (1985) , A diffusion-transport theory hybrid method for calculating neutron flux distribution in slab lattices, Atom kernenergie,47, 221. Nanneh, M.M. and Williams, M.M.R. (1986), An analytical solution to the hybrid diffusiontransport equation, Annals, nucl. Energy, 13, 337. Nanneh, M.M. (1990) A Synthesis Method Based on Hybrid Principles for Finite Element Neutron Transport, Ph.D thesis, London University, 1990. Williams, M.M.R. and Goddard, A.J.H. (1981) (Editors) Finite Element Methods in Radiation Physics, Ann. nucl. Energy, 8, 539-722.