Composite finite element solutions for neutron transport

Ann. nuel. Energy, Vol. 15, No. 8, pp. 397~,19, 1988 Printed in Great Britain. All rights reserved

0306-4549/88 $3.00+0.00 Copyright © 1988 Pergamon Press plc

COMPOSITE FINITE E L E M E N T SOLUTIONS FOR NEUTRON TRANSPORT R. T. ACKROYD and W. E. WILSON Nuclear Power Group, Department of Engineering, Queen Mary College, University of London

(Received 22 December 1987) Abstract--Composite solutions for neutron transport employ a finite element representation for a spatial dependence of the angular flux and spherical-harmonic expansions for its directional dependence. The orders of the spherical harmonic expansions are varied from region to region, according to the physical characteristics of the regions, so that the massive calculations inherent with high-order expansions can be confined to the regions for which they are essential. The transitions from low- to high-order expansions are made at the interfaces of elements. The K ÷ - boundary-free maximum principle provides a convenient means of optimizing a composite approximate solution for the angular flux. The way in which penalty weighting of interfaces is used to optimize composite solutions is investigated theoretically and numerically. The method is illustrated by solutions of two benchmark problems. An estirdate for penalty weightings is obtained by deriving a central difference analogue of the even-parity transport equation. This provides also a demonstration of the use of admissible, wildly nonconforming elements in a maximum principle.

1. INTRODUCTION

1.1. Discontinuous approximations for neutron transport In this section, the development of variational principles and weighted residual methods which can admit discontinuous functions as approximate solutions is reviewed. The exact solution for the angular flux in the transport equation is a continuous function of position r and direction l'~ for points within a physical region. At a point on the interface of different media, the angular flux is continuous except for those directions lying in the tangent plane to the interface. Thus, in 1-D slab geometry, the flux is continuous in angle except at /~ = 0. Martin and Duderstadt (1977) have shown that when solving the first-order transport equation for problems with strong heterogeneities, there can be substantial discontinuities at # = 0 on some of the interfaces. An accurate treatment of such discontinuities is undertaken conveniently by employing, for example, a double P N approximation, or its equivalent in Walsh functions as shown by Splawski (1981). In seeking good and economical approximate solutions of neutron transport problems, discontinuities in the approximate solutions are introduced deliberately for the following practical reasons : (i) the problem domain may contain regions of greatly differing material properties. Thus it may be advantageous to allow discontinuities in the angular flux at interfaces so that appropriate basis functions can be applied in each region to provide a close approximation to the solutions of the transport equation. The price paid is a discontinuity on the boundary of the region. F o r example, diffusion theory is usually accurate enough for those parts of a reactor core which are comprised of fuel assemblies of uniform composition, since fission tends to keep the flux almost isotropic. In the neighbourhood of control rods, the core reflector boundary, and between parts of the core with dissimilar enrichments, transport effects are important. In shields, the neutron angular flux becomes increasingly biased forward with increasing penetration, i.e. more and more anisotropic, and so at great depths, transport theory should be used. If large parts of a system can be treated by diffusion theory and the full transport method is reserved for those regions which need it, some savings in computation can be effected ; (ii) the gradient of the angular flux may be very steep in places, e.g. at the interface of a source region with a strong absorber. Rather than trying to obtain this result with continuous elements, it may be simpler to allow a discontinuity at the interface ; (iii) discontinuous trial functions may result in the solution converging faster. 397 ANE

1 ~,,'8.- [

398

R.T. ACKROVDand W. E. WILSON

The motivation for the construction of variational principles and weighted residual methods which admit spatially-discontinuous trial functions originally came from the synthesis method in which a given problem is reduced to one of lower dimension, e.g. plane geometry calculations at various axial zones in a reactor. Thus, discontinuities in the axial trial functions would allow appropriate trial functions to be used in each axial zone. Pomraning (1966) first applied the Roussopolous stationary variational principle to transport theory and allowed trial functions which were discontinuous, and so did not obey the boundary conditions. Nelson (1975) then formulated a generalized Roussopolous principle from which other stationary variational principles can be derived. The principle admits discontinuous trial functions so that different basis functions can be applied on either side of an interface. Pitk/iranta (1977) has considered variational principles admitting discontinuous trial functions for the diffusion and transport equations. He analyzed the stability and convergence of such solutions. Duderstadt and Martin (1979) have applied discontinuous phase-space elements in the integral law (or weak) form of the transport equation. Discontinuous elements are employed in the phase-space element code ZEPHYR of Mordant (1975) and also in several of the well-known finite element/discrete ordinates codes, e.g. ONETRAN and TRIDENT. In the sequel, a maximum principle is employed which admits discontinuous trial functions. As it is an extremum principle, increasing the basis for the trial functions automatically improves the approximate solution. The principle also involves penalty weights for interfaces and perfect reflectors. By increasing the weights and using appropriate bases for the trial functions, discontinuities can be eliminated practically, if desired. The principle shows that when close satisfaction of the perfect reflector and interface conditions are imposed, the satisfaction of the transport equation by the approximate solution is impaired. In obtaining a good approximation there are therefore these two conflicting trends. In the applications, large weightings are used to enforce near satisfaction of the interface and reflector conditions. This provides a convenient means of employing different orders of spherical harmonics for the directional specification of the angular flux in the various regions of a system. The same kind of spatial approximation is used, in the examples for all regions, but this can be varied if required. This method of composite solutions is based on an extremum principle using free trial functions, viz trial functions satisfying no interface or boundary conditions.

1.2. A variational principle employing free trial functions The K + principle (Ackroyd, 1982) admits two free trial functions ~b+(r, ~ ) for the even- and odd-parity angular fluxes q~+. Neither of the trial functions needs to satisfy any boundary or interface conditions. (Hence the appellation free.) The principle is a fundamental one, because special cases give rise to a family of variational principles (Ackroyd, 1983). Some of these are complementary principles for the second-order equations in q~ff, another is a maximum principle for the first-order equation. The method also leads to a family of weightedresidual equations of the Galerkin and Petrov-Galerkin kinds. Weighted-residual methods are obtained in this way for both the first- and second-order transport equations. These developments are summarized by Ackroyd (1986). Here the K + (~b+, ~b-) principle is applied to the numerical determination of composite approximate solutions. The K+(q5+, 66 ) principle can be stated in the following way:

K + (dp+,dp)<~Iv{(C-'S+,S+)+(GS-,S

)}dV+2fsfn s

.~'nlT2df~dS=2c~(say),

(1)

.n,~ 0

for a system with volume V (Fig. 1). Here (., .) denotes Sn" x" dfL The even- and odd-parity components of the volume source are S ±. The surface source is T. Equality holds if, and only if q5+ = 4)+ and ~b- = q~o. The operators C - ~and G occur in the following way in the even-parity equation : - l ' ~ . V[G~-Vq~ +] + Cdp~ = S + - ~ . VGS .

(2)

The functional K + is given by :

K + (c~+,dP )=K+(dP+)+K-(49-)+L(d?+,c~ )+M(~b+,~b-)--Ir(~b+,O )--Ipr(~b+,~b-).

(3)

Here K+(~b+) and K-(~b-) are the functionals specifying the well known maximum principles, K + for ~b~-and K for q~o, of the second-order forms of the transport equation. Whereas the trial functions in these classical principles have to satisfy a continuity condition at interfaces, and a reflection condition for perfect reflectors, there are no such restrictions on the trial functions in the K + - principle. The functionals L and M are defined

Finite element solutions for neutron transport Surface with

399

IBare surface

a so

Perfect reflecting surface Spr

Surface with a source

Bore surface Legends Vi • Vj typical regions/elements

rr,,--

Sj surfoce of Vj enclosed by regions/elements

- -

Si

"

Vi

Vi in the surfoce of the system = := = Sb bore surfacer Ss surfaces with imposed sources × ~ × × Spr perfect reflecting surfoce

......

S~

"

Fig. 1. Specification of system. by :

M(O+,O-)=-2fsf~'nO+O-df~dS,

(4)

L(~+,O-)=-2f~ ~'niO+O-dDdS

(5)

pr

for the perfect reflector Spr and

$1

Here Si is that part of the surface of the element Vi which lies entirely inside V, and the normal ni is outward to Si. Instead of using interface continuity conditions on O + to link the elements of V together into a system, the positive weighted surface integral : [ w+(r,f~)[0+(r+0,f~)-0+(r-0,f~)] 2 )

If(o+,qS-)=f~
(6)

is used in the way outlined below. This integral is taken over all the interfaces of the elements, i.e. over The two sides of an interface Si c~ Sj are specified by the position vectors r___0, and n is normal to the interface. In applications :

u (S~c~Ss).

w + = Aft + (r, ~ ) , w- = Aft- (r, ~ ) ,

(7)

with 2 a positive parameter and w±(r, ~ ) are specified positive weights. In the sequel it is shown that If(b +, 0 - ) steadily decreases as 2 is increased, i.e. increasing 2 tends to enforce continuity on 0 + and 0 at interfaces. A similar enforcing role is undertaken by the functional : Ipr = 2 { . {_ I f l ' n l f f ( r , n ) { [ O + ( r , n ) - - O + ( r , I I * ) ] 2 + t O - ( r , fl)--O-(r,l'~*)]z} dDdS,

(8)

O% r

where fl* denotes the reflected direction to the incident direction f~ on the perfect reflector apr. The functional I~r is used to enforce the reflection condition by increasing 2 with the directional dependencies of 0 + and 0-

400

R.T.

ACKROYD a n d W. E. WILSON

represented by an even and an odd spherical-harmonic expansion, respectively. The functional If is used in a more general way than Ipr, because its role in the method of composite solutions is to join up elements so that different orders of spherical-harmonic expansions can be employed for the directional dependence of a trial function on the two sides of an interface. Another way of employing trial functions with spherical-harmonic expansions which vary from region to region, is to make the orders of the expansions change across the border elements of regions and use classical variational principles. The proposed method of composite trial functions avoids the complications of the alternative procedure. The illustrations given are for slab systems for which the functionals L, M, If and Ipr simplify to

if'

L(~b+,gb ) = 2 . = ,

-~ [(a+(x~+O,l~)C~-(xi+O,#)-(o+(xi-O,#)r~

M(~b+, q~- ) = 2(1--/.)

~)+(a,#)4)-(a,#)#d#--2(1--lb) 1

(x~ - 0, #)l # d#,

(9)

r~+(b,#)4)-(b,#)#d#,

(lO)

1

[#l{w+(xi,#)[~+(xi-t-O,#)--([)+(xi-O,~)]2 Ww (Xi, l~)[~)-(xi"l-O,I-t )

If((b+, ~b-) = i=1

1

- ~ - ( x , - 0 , ~ ) ] 2 d#, Ipr(~b+,~-)=

[~[{w+(d.~)[&(d.~)-¢+(d.-~)12+w i=1

(d,,~)[4-(d.~)-~

(d,, -~)12} d#.

(11) (12)

1

The notation is : x = Xl, x 2. . . . , x N = the interfaces, x = d~ and d 2 = perfect reflectors, w -+ = the positive weighting functions, bare surface or surface with a source, perfect reflector,

atxb,,-- {;

bare surface or surface with a source, perfect reflector,

2. SIMPLE THEORETICAL BOUNDS FOR THE EFFECTS OF WEIGHT VARIATIONS

2.1. Least-squares error of an approximate solution

In this section, bounds as a function of weighting are derived for the component functionals of the leastsquares error W + - (~b+, 4)-) of the approximations ~b+ and ~b- for ~b~-and 4)0 respectively. This least-squares error is defined by Ackroyd (1982) and it satisfies : t" f ( c K+-(qS+'~b-)+W+-(qS+'O

Is+,s+))

>dV+2Jst" Jn s

]~'nlTZd, dS=2~(say).

(13)

-n<0

Numerical confirmation of these bounds is given in Section 3. However, for large weights these bounds are crude, but the effects of large weights can be established by asymptotic analysis, as in Section 4. The functional W + (~b+, ~ - ) can be written in the following way : W + - (q~+, ~b-) = /r(~b + , ~b-) -t-,~Id (~b + , q~-),

where I, and Ia are given in terms of the usual functionals defining W + given by :

(14)

as follows. Thus /r(~b+, ~b-) is

/r(~b+, q~-) = I(q~+, 4)-) +J(~b +, ~b-) + l,(~b + , ~b ) + Ib(q~+ , ~b-),

(lSa)

Finite element solutions for neutron transport

401

where :

I(dp+,q5 ) = fv(G[ll'V4J+ +G-'dp--S J(~+, ¢ ) = f < c - ' [ n - v ¢ -

],l-l-V4b++G

+c~ + -S+l,n.v¢

'c~--S-)dV,

(15b)

- +c¢ + -s +) dr,

(15c)

Ib(¢~+,4~-)=f~b{f,..>oln.nl[4~+-4~-]2dn+f,..
(15d)

I'(c~+'¢-)=fs,{fn..>oll 'nl[dP+-dP -T(r'-l'l)]2+fn .... .ll'nl[+++q$---T(r,l-D]2df~}dS.(15e) Also 2Ia(¢ +, 4b ) is defined by : 2Id(q$+,q$ - ) = lr(q~+,~b

)+Ior(dp+,¢ -) =

2

~'n

~+~(r,n)[q~_(r+0,n)_q~_(r_0,n)12)d~taa

.

(16)

+f,foI~'nwtr'")~+[¢-(r,~)--c~-(r,f~*)lZ]d~dS ....

(

[¢+(r, n ) - ¢+(r, n*)Y~ ,~

The component lr(q~+, q~ ) of the least-squares error arises from the errors made in not satisfying the transport equation within the elements, on the bare surface and the surface with a source. The component Ia(q~+, q~-) arises from the mismatch of the trial functions across interfaces and the failure of the trial functions to satisfy the conditions for a perfect reflector. Thus the component la(¢ +, ~b ) has contributions from those surfaces for which essential boundary conditions are imposed for trial functions admissible in the classical K ÷ and K principles. With the definition : K+-(~b+,~b ) = K + ( q $ + ) + K ( 0 5 ) + L ( q $ + , q $ ) + M ( ~ b + , ¢ - )

(17)

and the expression (3) the variational principle (1) becomes

K+-(4, +, 4>-) = K : - ( 4 ,

+ , ,~-)-~&(4,

+ , ¢-)

< 2~.

(18)

In the special case of ~+ and q$- satisfying the interface conditions and the reflection conditions/d(~b +, 4b-) vanishes, and the variational principle (1) reduces to : K~+ (~b+,q5- ) ~< 2~.

(19)

2.2. Variational inequalities In this analysis, the trial functions are linear combinations of appropriate basis functions of r and gL even functions of g~ for ~+ and odd functions of f~ for 4~ • The bases are assumed to have the property that they are capable of representing at least one even-parity trial function and one odd-parity trial function of the kinds which satisfy the interface continuity condition and the perfect reflector conditions. The bases are insufficient for the exact representation of q~+ and ~b . In the expressions : M ±

q5-+(r, ~ ) = ~ am ~ q~m e (r, l'~),

(20)

m--I

the components 4)~+ and q~ are the even and odd basis functions, respectively. The basis functions are linearly independent and they are defined locally for each element. Thus a basis function as a function of r vanishes outside the element for which it is specified by a polynomial in the local coordinates of the element. In general the coefficients a~, are independent. For 2 > 2' >/ 1 let q$~ be the trial functions which maximize K + - (q$+, q~-) -- 2Id (¢ +, q$-) and let 4b~ be the trial functions which maximize K + -[q$+, q ~ - ] - 2'Id(4b + , ~b-). Specifically:

402

R . T . ACKROYD and W. E. WILSON M ±

CZ

Z

M ±

am~q~,, ± ± and¢~

m=l

Since C f are the unique pair maximizing

K + - (~+,

~

=

+ ~ m+ " am2"

(21)

m=l

q~-) - 2Id(¢ + , q~-)

the inequality

:

K~+ - ( e l , qS;) -- 2Id (q~f, 0 ; ) < K~+ - (q~-, q~j-) -- 2Ia (qS~-, q~-),

(22)

K + (~b~-,~bZ) - 2'Io(~b + , ~bj-) < K + - (q~+, $ ; ) -2'Ia(~b + , ~bj:).

(23)

holds. Similarly :

Since2>2'>tl: K + - (tk+ , ~b~:) - 2/0 (¢+, ~b~:) < K + - (q~+, q~-) - 2Id (q~+, ~b~-) < K + - (~b+ , qSj-) - 2'Ia (q~-, ~b~-)

< K+-(4+,¢;)-2"I~(¢+,¢;).

(24)

The sum of the first and third entries is less than the sum of the second and fourth entries. Hence : 0 < [2-2'][I~@~,, $ ; ) - I ~ ( ~ +, ¢ ; ) ] , i.e.

or

O<~Id(¢+'ck;) < Id(¢;, ¢;:) t" lid(¢?, ¢;) -I~(¢ +, ¢;)]/[,~-,V] < 0

(25)

Thus the measure Ia(~b~-, q~j-) of interface and reflector errors decreases with increased weight 2. The third and fourth entries of the chain (24) give with (25) : K + - (~b~-, q~;) < K + - ( q ~ , ~b~;).

(26)

Hence both K + - (q~-, q~j-) and I a (q~-, ~bj-) are monotonic decreasing with increasing 2. The identity (13) with the relations (14) and (18) give for arbitrary q~± : K + - (¢+, ~b-) -q-Ir ((~ + , ~ - )

=

2tZ.

(27)

Since K + - ( ¢ +, ~bj-) decreases with increasing 2, I~(¢j ~, ¢~-) is an increasing function of 2. This result may seem unwelcome, b e c a u s e / ~ ( ¢ + , ~bj-) is the component of the least-squares error arising from the failure of ¢ ~ to satisfy the transport equation within the elements and the natural boundary conditions for the bare surface and the surface with a source. Nevertheless, the monotonic increasing/r(~b +, ~bj-) is bounded for the following reason. Consider the set of classical trial functions, i.e. functions which satisfy the interface-continuous condition and the perfect reflection condition. F o r these functions Id vanishes and the functional K + - reduces to the functional K + - . F o r the set of classical trial functions let ~b~ be the pair that maximizes K + - . In the set of all possible trial functions q~ is the pair which maximizes K + - (q~+, q~- ) - 21d (~b+, ~b- ). Therefore : K~+- (q~+, ~b~-) ~< K~+-(~b~- , ~bj-)-).Ia(~b~, ~bj-).

(28)

Subsequently in Section 4.1, it is shown that equality is approached in this result as 2 ~ ~ . Hence K¢+ - (q~-, ~bj-) is bounded b e l o w / ~ + - (~b+ , q5¢-) and since the identity (27) holds, it follows that Ir(~b~-, ¢~-) is bounded above b y / r (O+ , ¢¢- ). This upper bound is approached as 2 --* oo. The above results can be presented in a simple way when the trial functions ~b+ are specified on the bare surface and the surface with a source, i.e. the external surfaces of the physical system.? As the relative weighting 2 is increased the least-squares error Id (O + , ~b£) arising from the interfaces and the perfect reflector is reduced, whereas the contribution to the least-squares error /r(O~-, ¢ £ ) resulting from failure to satisfy the parity transport equations within the elements is increased. Briefly enforcing the essential boundary conditions impairs the degree of satisfaction of the transport equation within the elements. t Perfect reflectors occur in a physical system as an internal surface dividing a part of a system from its mirror image. For computational purposes the perfect reflector is taken as an external surface on part of the system.

Finite element solutions for neutron transport

403

The overall error W+-(~b+,q~ -) is given by expression (14) as I~(qS~-,q~-)+)`Id(q~-,qS;) i.e. 2e +).lo(qS[, q~a-)--K + -(4~-, q~-) on using the results (13) and (18). Furthermore, using the inequality (28) the overall error becomes : W + - (4~+ , ~b~) ~< 2c~--K + - (qS~+ , ~bc = / r ((]~c+ , [#c-),

(29)

i.e. the greatest overall least-squares error is made by the classical approximations q5+ . 2.3. Some simple bounds f o r errors

(i) Id(qS~-, q~-). Adding the inequalities (23) and (28) gives : /'d (~b+ , ~b~-) < [K + - (~b~,, qS~) -)`'Id(q~', qS;,) - K + (q9+ , q5~-)]/[2 -)`'].

(30)

An alternative upper bound is given by the inequalities (26) and (28), viz. to (ok;, 4 ; ) < [K+ - (4~+ , 4~; ) - Kc+ - (ok+ , 4V )]/2.

(3 l)

2* = [K~+ - (~b~, $~,) - K ~ - (~b~+ , qSd-)]/Id(q~,, q~j:).

(32)

These bounds intersect at :

The inequality : K~+ - (ckJ, ck~-) < K [

(dp],, (o~:) -2'I~(q9}, ~b~),

(33)

holds, as a special case of inequality (28). Hence 2* > 2', and for 2 > 2", the upper bound (30) is lower than the bound (31) : (ii) K~+ (q~-,~b~-). Eliminating Id(~b+, ~b£) from the inequalities (23) and (28) gives: K~

)/ (ek~;,~Z) < K~+-(dp+,dpz.)--2"Id(dp~,ck£)+ ~ [ K c

+

(qS~+,,~b~)-)`'Id(~b~,,qS~:)-K+-(~b~+,~bj)].

(34)

Another upper bound for K~+ - (~b~, q5£) is given by the inequalities (23) and (31), viz. 2' K~+- (~b~-,qS;) < K~-(~b~, q~a,)--2'Id(q~,,~b~,)+ ~[K[-((o~,,dp~,)-K~+-(O~+,d);)].

(35)

The bound (34) is lower than the bound (35) if). > 2*: (iii) L(q~J-, 4~-)Lower bounds for Ir(~b~-, ~b;7) can be obtained from inequalities (34) and (35) with the use of the identity (27). Inequality (34) gives rise to : )`' + L(~b~,~bg) > /r(q52:', 4~') + 2 ' I d ( 4 ~ , 4~:)-- 2 Z ~ [ L ( ~ b ~ , 4~-)-L(~b~,, ~b~)-2'Id(~bf, ~b~,)].

(36)

The inequality (33) and the identity (27) show that : Ir(~ + , ~b~-) > /r(qS~', ~b~7)+)`'Id(q~ + , q5~7).

(37)

Hence the inequality (36) gives the lower bound as an increasing function of 2. The asymptotic value /r(4~, 4~') +).'Id(q~, qS~,) of L(~b~-, ~bZ) is also an increasing function of 2'. The inequality (35) leads to : )`' + /~ (~b~-, ~b; ) > /r (~b~, q~z:) + ).'Is (qSf, ~bz,) - ~ [It (q~o, qS~-) - I~(~b+ , ~bz,)].

(38)

Ir(4~) , ~bj) ~> I~(4~-, ~ - ) ,

(39)

The upper bound :

404

R.T. ACKROYD and W. E. WILSON

1..1 Asymptotic result(39)

~

h

\. N.

/

~%,

Jpper bound (39) for ~(e~,(~.)

/

bound(381 -for I r ( ~ ~) mseswithincreasingh' Retofive weighting X

#

--x

Upper bound (30)for Id (~,<~)

I\,// ". I \ 1

Upper bound (311 x t \

for Id(~,~)

I".,,...~..

X

X~ Retafive weighting X

Fig. 2. Trends of bounds with increased weighting.

has been established as a consequence of the result (28) and the identity (27). For 2 > 2", the lower bound (36) is higher than the lower bound (38). The trends of the above bounds for ld(~b2-,~b~-), K + (~b~-,~b~-) and/~(~b~-, ~b~-) as a function of the relative weighting are sketched in Fig. 2. 3. EXAMPLES OF COMPOSITE SOLUTIONS

3.1. Bounds on the least-squares errors f o r a lattice cell problern

The results described below were obtained from the DUPLEX code of Wilson (1985). The directional representation of the trial functions ~b± for slab systems is by expansions in Legendre polynomials. The trial functions ~b+ of the K + principle are represented spatially as piecewise-continuous functions with their discontinuities at the interfapes of the elements. Thus for a linear interpolation between the ends of an element, there are two coefficients for that element to be specified independently of the coefficients for the other elements. Bounds on Id(~b~-,~b~-)and Ir(~b~-,~b~-)are shown in Figs 3 and 4, respectively for a reactor cell. This problem for a slab cell with uranium fuel and a graphite moderator has been treated analytically by Weinberg and Wigner (1958). Here one example of the treatment of this problem by the composite method is given. The results are for a P0 approximation for ~b+ and a P~ approximation for ~b-. The mesh had twelve elements of the same length. The weighting functions ~-+ were taken to be unity, and the relative weighting 2 was varied from 1 to 10 9. The upper bounds of Fig. 3 for ld(~b+, ~bz) do not differ significantly for the present problem. Lines A and B refer to results (30) and (31), respectively. Although the bounds on Id (~b~-,~b£) are crude, they are of practical value because the reduction in Id induced by a change in the relative weighting 2 is underestimated by the bounds. A similar remark applies to the lower bound of Fig. 4 for/r(~b~-, ~b[). Whereas the elementary bounds for Id(~b~-,~b~-) vary as 1/2, the calculated values of Ia(q~-,~bz) vary asymptotically as 1/2 z. The elementary lower bounds for Ir(~b~-,~b~) also fall short of the upper bound

Finite element solutions for neutron transport

405

1010-~i

15j-1

Upper bound (39)

10-31 40-51

Idl®\,~llo_~i I0-91

~(~+%.'*X)50 / / . 1

. . . .

Lowerbound(38)

10-111

10-17 Retofive weighting

-5 1 10

~.

Fig. 3. Bounds on Id (~b~-, ~b~-) for a lattice cell problem.

10 10 10 Relative weighting %.

10

Fig. 4. Bounds on/r(~ +, ~b).) for a lattice cell problem.

/~(~b~+,q~). More precise asymptotic behaviour is established theoretically in Section 4.5 for Id, and the corresponding asymptotic behaviour for Ir is given in Section 4.6. In applying the K ÷ (~b+, q~-) principle, the practical problems of choosing ~-+(r, fZ) and the relative weighting 2 have to be resolved. A number of choices for #-+ have been considered, but the simplest ~-+(r, f4) = 1 is satisfactory for numerical work. The analysis of Section 6 provides a first estimate of ,L

3.2. Reed's edge-cell problem Reed (1971) used the problem of Fig. 5 to provide a severe test for various difference schemes associated with the discrete ordinates method. The data of Table 1 show the regions of the edge-cell have very different physical properties. The bare surface in a small system also increases the difficulty of obtaining a good solution. Solutions of this problem which employ a weighted residual method with a discrete ordinates representation of the angular flux can suffer from spurious oscillations due to the ray effect. Some examples are reported by Ackroyd and Splawski (1982). The presence of strong absorbers in a cell with a bare surface can make the angular flux strongly anisotropic in places without isotropic sources. Consequently, Galliara and Williams (1979) had to use high-order spherical-harmonic expansion everywhere to obtain an accurate solution with their precise finite-element method based on the K ÷ (~b+) principle. By contrast, high-order spherical-harmonic expansions can be confined in the method of composite solutions to the regions where they are necessary. In the fuel region, the low-order P] expansion can be used for the angular flux, because neutron streaming is negligible except in the neighbourhood of the interface of the fuel and clad. The neutron streaming is negligible in the fuel region with its fiat isotropic source, and the absorption of a neutron very close to its point of emission. A P, expansion is specified therefore for 0 < x < 1.9, the greater part of the fuel region. A P7 expansion is used in the remaining 5 m.f.ps to the interface with the clad to allow for the interface transient in the angular

Table 1. Material properties of the edge-cell

ANE

151~- D

Region

Width (cm)

Absorption cross-section (cm 1)

Total cross-section (cm 1)

Source (ncm Zs 1)

1 2 3 4 5

2.0 1.0 2.0 1.0 2.0

50.0 5.0 10-s~0.0 0.1 0.1

50.0 5.0 10-s~0.0 1.0 1.0

50.0 0.0 0.0 1.0 0.0

406

R.T. ACKROYDand W. E. WILSON

Region

1

2

Fuet

Ctod

I I

•

i

Uniform isotropic source in n strong ,absorber

0

I

3 Air-gap

I I

t,.

I

Moderator

I I

5

-~1Isofropi~ scatterer I and weak absorber

mltd I absorber I

=

IUniform I I sofropic I

I I I I I

source_ _ I

No source

I

I

5

3

II Perfect reflector

Bare surface

x(cm)

Typical Flux plot for o Variational solution (after 5alliar'a and Williams)

Ftux

I

I I

--~X

Weighted residual sotufion reported by Ackroyd and Sptowski Flux

Fig. 5. Reed's edge-cell problem.

flux. For the purely-absorbing clad, the air-gap and the source region 4, a P7 is used, because the surface x = 5 is effectively bare for a neutron travelling leftwards across it from the source into the air-gap. Ahead of such a neutron, is eventually a considerable thickness of a pure absorber. Region 5 absorption is small compared with isotropic scattering. In the absence of significant streaming, a P1 expansion for the angular flux would be adequate. There is significant streaming in the direction of x increasing because of the presence of the bare surface at x = 8. A P3 approximation for the angular flux is assumed therefore for region 5. The composite P1,3,7 solution closely follows both the standard P7 finite-element solution of Fig. 6, which is obtained using the classical K + (~b+) principle, and the F N exact solution given by Garcia and Sievert (1980). The errors in both flux approximations are shown in Fig. 7 to be effectively confined to the neighbourhoods of the interfaces of regions and the bare surface. The mesh used for both the P~,3,7 and P7 calculations is given in Table 2. The mesh is refined at the interfaces of the dissimilar regions and in narrow source region 4. The latter refinement ensures accurate shedding of neutrons to the left and right of the flux peak.

3.3. A multi-layered shieM problem In the edge-cell problem the fluxes varied across the small system by about a decade. The range of flux variation for the shield problem is about seven decades. The cross-sections varied strongly from region to

Finite element solutions for neutron transport

407

2"210-

2"ff 81.8-

6-

1-6-

k-

1./+-

Comparisons with FN method

2-

1.2-

/

P

,7" 1"0.G

composite "6-

b_ 2N_~-

das~cat

-6-

.2o

-10

xcc.sl

I_ Pl

-'

I_

P7

-'-

P1

P7 Clossica[ sotution Plj3,7 £omposite solution

P7

LI

P3

"

X

P3

Fig. 6. Composite and classical finite element solutions for the edge-cell flux.

zx 0 ....

Fig. 7. Errors of composite and classical solutions of the edge-cell problem.

region of the edge-cell, but the variation for the shield problem is small. The isotropic sources in the shield are remote from the bare-surface, where the neutron flux has to be known fairly accurately. The shield is specified by Ackroyd and Williams (1986). The media are taken to be pure absorbers to make the angular flux strongly anisotropic with increasing depth of penetration into the shield. A natural selection for the composite angular approximation as a function of penetration depth x is shown in Fig. 8. The composite Pi,3,7 finite-element solution is compared with P3 and P7 classical finite element solutions in Fig. 8. A t the bare surface the Pi,3,7 and P7 fluxes are in good agreement, but the P3 flux is lower by a factor of about 5. A P ~5 classical finite-element solution provides the benchmark used in assessing the accuracies of the composite Pi,3,7 solution with the P7 classical finite-element solution, and the solution by well-established A N I S N finite-difference discrete ordinate code. The percentage errors made by the Pi,3,7 spherical-harmonic methods and the discrete ordinate method A N I S N $6 are shown in Fig. 9. All the methods incur small errors in the vicinity of the places where the source is changed in magnitude. As the bare surface is approached, there is small and increasing error made by the P7 calculation. The increasing error made by the P1,3,7 calculation is larger than that lbr the P7 calculation, because of the additional errors induced by the sudden changes in the angular expansions at x = 100 and 170. Nevertheless, at deep penetration, the PI,3,7 composite calculation is more accurate than the A N I S N $6 calculation. 4. ASYMPTOTIC ANALYSIS 4.1. Some limits The elementary bounds (30) and (31) show that l i m ~ I d ( ~ b ~ - , ~ ) = 0 , but do not show that l i m ~ Ir(~b~-, ~b~) = Ir(~b~+, ~b~-) and l i m ~ K~+ - (~b+, ~b~-) = K + - (~b~+ , ~p~-) as indicated by numerical Table 2. Mesh for edge-cellproblem Region

Range of x

1

0.0-1.9

2 3 4 5

1.9-2.0 2.0-2.1 2.1-2.9 2.9 3.0 3.0-5.0 5.0-6.0 6.0-8.0

Number of linear elements 2

10 10 8 I0 2 20 10

408

R . T . ACKROYDand W. E. WILSON 10~10./_

10610s-

\.

O-

~0~-

o

~ 1 0 ~-

I.>

-2-

\

_~.. -6-

10-

,~,~ and Pv

1-

lO-Z "

\\

-8-

\\

-10-

10-'" Rnngesfor Pt~,7method ,PI---'i'-P~, = I 50

100

P3~

,

,

~ P7

150 200 250 Penetration x (cm)

,

7t

,

~00

350

.3.v

-12_% = o

\ P,I~["-P'

so

, ~]~

P'~

~oo ~so z~o ~o

\

ANISN $6

\,'J -I-'

~6o ~o

,60

Penefrofion depfh x(cm)

Fig. 8. Flux in a thick shield as a function of penetration depth x.

Fig. 9. Errors of shield solutions relative to P~ s benchmark.

examples. The limits for L(q~-, 0 ; ) and K~+ (~b~, 4)2-) are established below, and asymptotic formulae are given for the ways in which the limits for I d ( 0 ; , q ~ ; ) , /r(q~-, ~b~-) and K~+-(0~-,(b;) are approached as l i m a ~ q~f = 0~+. With 2' fixed, the expressions (30) and (31) give: lim Id(q~- q~;) = 0.

(40)

Since Id (0~-, 0 £ ) is a positive-definite functional, this limit implies: lim [q~f (r + 0, ~ ) - q~ (r - 0,1~)] = 0,

(41)

lim [q~±(r, ~ ) - 4 -+(r, fl*)] = 0,

(42)

on the interfaces, and

on the perfect reflector. Let q~+ = l i m ~ ~bf, then ~b~ + satisfies the interface and perfect reflector conditions of the classical K ÷ (0 +) principle; and hence it is admissible in that principle. Similarly 4% is admissible in the classical K ( 0 - ) principle. Thus : lim K+ (~b~-) = K+ (q~+),

lim K - (~b£) = K - (q~).

(43)

The property (26) shows that K +-(4~-, ~ b ; ) - K + (q~+,~b~-) is monotone decreasing with increasing 2. The inequality (28) shows that K + - ( 0 +, q ~ ; ) - K + - ( 0 +, q~-) is bounded below by zero, and hence it has a limit as 2 ~ ~ , i.e. lim K +

J.~oo

(~b~-, 4~-) ~> K + - ( 0 + , 0 U ) = K + ( 0 + ) + K -

(qg~).

(44)

K + (~b+ , ~b~-) resolves into components K + (~b+) and K~ (~b~-), because terms L(q9+ , q~d) and M(q9 + , ~b~-) for definition (17) vanish. Definition (17) and limits (43) give lim K + - (~b~-,~b~-) = lim

[K+(c],~)+K-(c#;)+L(c~ +, ~b;) + M ( q ~ - , ~b;)] = lim [K+ (~b~-) +K-(q~;)] = K+(q~2)+K-(0~).

(45)

This limit and inequality (44) give :

K+((a+)+K (0~) >~K+(c#+)+K ( ~ - ) .

(46)

Finite element solutions for neutron transport

409

The greater-than sign is inadmissible, because if the greater-than sign held then at least one of the inequalities : K+(~b~+) > K+ ($~+),

K-(Ooo) > K - (~b~-),

would have to hold. This is impossible, because q~+ and ~b~ are by definition the maximizers of K+(~b +) and K-(~b-), respectively. Hence the inequality (44) reduces to: lira K ~ +

- (~b~-, ~bj-) = K + - (tk~+ , qS~-).

(47)

The identity (27) then gives : lira I~ (the-,

2~oo

~ b £ ) = 2 ~ - K J - (~b~+ , ~b; = I~ (~b~+ ,

~bc ).

(48)

4.2. Representationof trialfunctionsfor large2 Let : M ±

=

ama~bm(r, ~ ) ,

(49)

m=l

be the maximizing pair of K + -(~b + , q S - ) - 2Id(q~+, q~-) for trial functions q~± of kind (20). The basis functions q~ are allocated to the elements in independent sets. The basis functions allocated to a given element vanish identically beyond the bounding surface of the element. The optimized classical approximations ~b+ can be written in the way : M ±

m~l

but some of the coefficients are related to coefficients of the same parity. For example, those coefficients c~ associated with nodes either on an interface, or a perfect reflector, are related because ~b~ + has to satisfy both the interface-continuity and the perfect reflection conditions. Some of the coefficients c,~ are zero because they are associated with nonconforming basis functions which cannot be used in representation of classical trial functions. Since the basis functions can represent classical approximations : M +

~b-+ = qSf + ~ d~4~,~,

(51)

m=l

represents an arbitrary pair of trial functions. The maximizing pair (49) can be written as : M ±

=

d,,~b~.

(52)

m~l

4.3. Optimizationof trialfunctions Having in principle calculated the coefficients c~ for the optimized classical approximation, the optimized approximation pair ~b~ can be found by maximizing K~*- (th + , q g - ) - 2Id(~b+, 0 - ) with respect to the d,~ of expression (51). In expression (17) for K~+-(q~+,$-) the functionals K+(~b +) and K-(qS-) are given by Ackroyd (1981) as: K+ (qg+) =

fvfn

{2¢h+S+ + 2 ( ~ ' V ~ b + ) G S - -n'VqS+(Gf~'V~b+)-~b+C~b+} d f l d V

+4 fo

.<0

;a'nl(tk+)2dl)d5 'g

'g~'nlO+Tdf~dS-[so~b~s,

(53)

{2dp-S- +2(~~.V49-)C-'S+-n'Vdp-(C-lf'~.gq~ )-tip

G-Iq~-}dndV

410

R . T . ACKROYDand W. E. WmsoN

The other terms are given by the expressions (4) and (5). The functional 21d(~b+, ~b ) is defined explicitly by expression (16). The result of maximizing K + (~b+, q~ )-2Id(q5 +, ~b-) with respect to d + can be written, for the maximizing values of coefficients d+~ as : M +

M +

M-

2 ~ A+kd~ = ~ B+mkdk+ ~ C , , ~ d k + D + ( S + , G S - , T ) k--I

k--1

(m = 1,2 . . . . . M + ) .

(54)

k=l

The A,,+k terms in these equations come from either the interface or the perfect reflector component of 2Id(~b+ , ~b ). The B,,+ terms arise from the second degree terms in q5+ and 1"~.Vq~+ in the expression for the K+(~b +) component of K + (~+,4~-). The interface components L(q~+,~b ) and M(~b +, q~-) of K + - (~b+, ~b-) give rise to the C~k terms. The constants D + (S +, G S - , T) arise from the first degree terms in ~+ and ~'Vq~ + of the K+(q~ +) component of K+-(q~+,~b-). Similarly, the maximization of K + (~b+, q5 )-21~(~b +, ~b-) with respect to dm gives equations of the kind: M-

M

M +

2 ~. Am~dea = ~ B,.,d~+ ~ C , . J ~ + D m ( S - , C k--1

k=l

'S+,T)

( m = 1,2 . . . . . M ).

(55)

k-I

4.4. Asymptotic behaviour Of ld (dp+ , dp£ ) For 2 >> 2* > 1 bounds of the kind :

km+

and ]dma[ <

k.;

(56)

must hold for some constants k + and k=. Otherwise b o u n d (31) for the positive-definite quadratic /d(~b+, ~bZ) in d+a would be violated. Hence for 2 >> 2* > 1 :

Z

+ + ~ D+(S +, GS-, T)/2, Am~dk~

(m

1,2 . . . . . M +)

k=l

(57)

M-

Amkdka,,~DT~(S , C - ' S + , T ) / 2 ,

( m = 1,2 . . . . . M - )

k=l

Thus for 2 >> 2* > 2' > 1 the maximizing coefficients satisfy : d ~ ~ h,~/2

(58)

for some constants h,~ determined by sources S +, S - and T. Consequently:

ld(dp~, (o£) ~ A2/22,

(59)

for a constant A z depending on the sources. The observed behaviour of Id(q~, ~b£) for the example of the uranium-graphite lattice cell is in accordance with the above asymptotic analysis.

4.5. Asymptotic behaviour Of Ir(dp+ , c~£ ) andK~+ - (dp+ , ~ ) Put, in accordance with the specification (52) : ,p~ = ~g + 0~ with M ~+

o~ = Z d~+~, m--I

to evaluate/

r(~b).+ , q~ ). The full definition of/~(q5 +, q~ ) is given by Ackroyd (1982), viz.

(60)

Finite element solutions for neutron transport

Ir(0 +, O ) =

;v{(G[D.'vO++G-~O -

+ +G-'O-

-S-],~'V0

-S-)+(.C

411

I[~'~'V0--t-C0+ --S+]. ~~*V0-

+CO+-S+)}dV+fs~{f...>olD.'nl[O+-O--T(r,-D.)]2d~+fn..
- T(r,")]2 d ~ } d S

+fsb{f..n>o,.'n,[O+-O-]Zd.-q-f...
(61,

Thus :

I

2(G[~'V~b++G

'0£-S

],

~'VOf+G

+ I ~ +(G["'v°~ +G 1°;]' "'v°~ +G-'°~) Ir(0],0~ -) = Ir(0~+,0~ -)

'0;-) 1

[( dV

.,Vl+z1 I+
~.vo;+co~>

3

+ ;~{fn-n> 0 ~L'n L[-2[0+--0~--T(r'--D')][0~--0Z]] d [0~ o- OZ]2+ + J.['..
[-2[0+ + 0~- -- T(r, ~ ) ] [ 0 + + 0;]-] . _ ] , . +[0: + o ; ] 2

J°~'~

~"n L

(62) Consequently for the asympotic approximations (58) one has : B

D2

(63)

/~(4:'~, 0; ) ~ L(O~+,O; ) - 7 + T ~where B depends on 4)+ and d,~z but D 2 depends only on d,~. The coefficient B has to be positive because :

L(0 + , Oc) i> L(0 +, 0~-) + 2 ~ ( 0 2 , 0 ; ) .

(64)

This inequality is established in the same way as the inequality (37). The asympotic result : K~+ ( 0 + , 0 i - ) ~ 2 ~ -

[

B

D~3

Ir(O$,O2)--~+-ff,

(65)

follows from the identity (27) and the result (63).

5. WEIGHTING FOR MESH REFINEMENT

5.1. Mesh refinement with constant weighting The mesh is refined by partitioning elements of the original mesh into smaller elements. The generalized least-squares error tends to zero as the mesh is progressively refined, but its interface component may diverge before converging. An example of this effect is shown in Fig. 11 for the uranium-graphite cell of Section 3.2. Let ~k± be a typical pair of trial functions for the refined mesh. The bases for the refined mesh are chosen sufficiently large so that they can act as bases for trial function pairs 0 -+ of the kind appropriate to the original mesh. Trial functions of the kinds 0 ± are continuous across the additional interfaces created in refining the mesh.

412

R . T . ACKROYD and W. E. WILSON OP=P OP'= P'

Y

OP" P"

I I

~

R)C=(Ir(*~Ab~),Id(*~,*~) )

-X Ir(dJ~,~b~)

:T

pvTzRY7

T _[ -I

Fig. 10. Effects of mesh refinement and increased weighting. ,3072

16'

,768

12'

lO~Id(®~,®-X.)

,12288 $2/+576

t,-

f,192

Numberof eLemenfsvaried ReLafiveweightingk fixed

,48

,12

'r(%%' Fig. I 1. Effect of mesh refinement at a fixed weighting.

Let q$~ be the maximizing pair of K+-(q$+, ¢ ) for the original mesh and a relative weighting 2'. The q$~ are admissible trial functions in the K + -(qs +, 4s-) principle for the refined mesh. Let ~s~, be the maximizing pair of K + -(qs +, qx ) for the refined mesh with relative weighting 2', then: K + - ( ¢ + , ¢ ~ : ) = K + ( ¢ ~ , , ¢ a : ) - - 2 ' I d ( ¢ + , ¢ ~ : ) ~< K[ (~O~,~Oz)--2'Ia(~O~,~) = K+-(~D~,,~D~:).

(66)

This result becomes : < r ( ¢ z , ¢2:) + 2 ' I d ( ¢ +, ¢ ; ) ,/1 + (z) 2

Geometrically this means that the point : R~, = [Ir(~/]~7,,~7), Id(~b+ , ffZ)],

2

(67)

Finite element solutions for neutron transport

413

generally lies inside the triangle OA'B" of Fig. 10, while the point : 0k, = [Ir(~b~, ff~:), Id(~b~,,~b;)], lies on the line A'B' with the equation in perpendicular form :

X+2Y

- p'.

(68)

~/1 + (2') 2 As the mesh is refined the point Ra, is shown to be forced nearer the origin by its confinement to the triangle OAB'. Consider among the trial functions ~-+ for the refined mesh the subset which satisfies the interfacecontinuity condition and the perfect reflector condition. Let ~ be the maximizing pair of K + -(~+, ~ - ) for this subset, then : K~+ - ($~+,~k¢) ~< K~+- ($~,, $;)--2'I,~(~0~,,$.e), or equivalently :

Ir(q'~',q';)+2"Id(q'+,¢;) ~

2

~<

Ir(¢+,qJ;) ~ 1 + (2') 2

- p" (say).

(69)

Thus the point Ra, lies generally inside the triangle OAB", where the line AB" has the equation :

X+2Y _

_

-

p".

(70)

The intercepts are :

OA = Ir(t~+,~kg-)

andOB" = Ir(~k+,~kc)/2 '.

(71)

Now K+ - (O~+, 0 J ) = K+ (0 +) + K - (OJ),

(72)

because the terms L(~,~+ , qJj) and M(O~+, 0 J ) vanish for the definition (17) as the consequence of the interface and boundary conditions imposed on ~O~. A property of the classical principles K+(~O÷) and K ( ~ ) is that K÷(O~+) and K - ( O j ) can be made arbitrarily close to K+(q~-) and K-(~bo) respectively by progressively increasing the bases for ~+ and ~,-. Consequently K~+-(~k~+,~kj) can be made arbitrarily close to K ÷ (~b~-)+K-(q~o) = 2~. The identity (27) then implies Ir(O~+, O j ) can be made arbitrarily small. Thus for the fixed weighting 2', the line AB" can be made to approach the origin, keeping parallel with the line A'B'. In this way, with a constant relative weighting 2', the least-squares errors Ir(O~', O;) and Ia(O +, ~Oj:) can be made arbitrarily small by progressively refining the mesh. The use of mesh refinement with a constant weighting can be an inefficient procedure, as the example given in Fig. 11 for the uranium-graphite lattice shows.

5.2. Mesh refinement with constant weighting and fixed interface error When the mesh is refined the increase in the number of interfaces where discontinuities in the trial functions can occur provides an enhanced facility for better satisfaction of the parity equations of neutron transport. The increase in the number of interfaces provides potential freedom for the least-squares interface error to increase. Here the cases are examined for which the least-squares interface error does not decrease. Consider trial functions Z -+ for the refined mesh which are constrained so that for the constant relative weighting 2', there is no increase in the interface error, i.e. /d(Z + , •- ) = /d(~ + , ~ )

= /d(~+ , Z~7),

(73)

where Xff are the optimized trial functions for relative weighting )~'. The optimized trial functions q~ for the original mesh can be regarded as a trial function pair for the refined mesh. They are a particularly fixed pair of the pairs ~±. The pairs ~± are constrained members of the ~O± for the refined mesh. Hence the inequalities :

K~+ (dp~,,dp;)-2"Id(q~,,dp;) <~Kc+-()~,,Z;)--)JId(z+,Z;) < Kc+-(~k+,~k;)-2'Id(~k~,,~;).

(74)

414

R.T. ACKROYDand W. E. WILSON

These inequalities, the constraint (73) and the identity (27) give : lr(%~,)~:) ~< Ir(~b],, q~z)

)

/~(O+ , •2;) + ~-'Id(~ + , O;;) < Ir ((%[, Xa') + 2'Ia (4)+ , qga,)~

(75)

The first of these inequalities shows, other things being held constant, that mesh refinement reduces the leastsquares error component arising from failure of the trial function pair to satisfy the parity transport equations. If the mesh refinement leads to an increase in the interface component of the least-squares error, i.e. if: Id(O+, ~/,~,) > Id(qg~,,~b~,),

(76)

I r ( ~ , , ~Oz) < Ir(x +, Zz) ~< lr(qg~, qg~:).

(77)

as in the example of Fig. 11, then :

Mesh refinement with constant weighting in this case is achieving an opposite effect to that desirable effect observed for a fixed mesh and increased weighting, viz. /r and Id as increasing and decreasing functions respectively of the relative weighting. 5.3. Mesh refinement with increased weighting Let the relative weighting for the refined mesh of Section 5.1 be increased from 2' to 2. In place of inequality (69) there is the inequality : ~< 41+22

P (say),

(78)

41+22

i.e. the point R~ = [Ir($ + , S f ) , Id ( $ f , S f )] lies generally within the triangle OAB of Fig. 10. The intercepts of the line AB are Ir(~'~+, $ c )

and

Ir($~+, $£)/2.

The point R~ can be located approximately within the triangle OAB as follows. The asymptotic analysis of Sections 4.4 and 4.5 shows that for sufficiently large 2 : B

A2 Id(qJ~,~Of) ~ 2y. The point R~ is therefore close to the parabola CA of Fig. 10 with the equation : B2 [X--/r(Oc+, ~k£)] 2 = ,42 Y.

(79)

(80)

As 2 is increased, the intercept OB decreases and Ra moves towards A in a path close to the parabola, CA. In practice for slab systems, it has been found that mesh refinement leads to improved approximations if the relative weighting is increased proportionately to the increase to the number of interfaces. 6. APPLICATIONS OF NONCONFORMING ELEMENTS

6.1. Uses o f nonconforming elements The K+-(q~ +, ~b-) maximum principle has been applied above for trial functions ¢+ represented by bases which admit trial functions with the potential to satisfy the interface-continuity condition and the perfect reflection condition. Amongst the elements used, there had to be a subset of conforming elements. As the relative weighting 2 --, oo the optimized trial functions ¢ f --* ¢~ where ¢+ maximizes K+(~b +) and ~b£ maximizes K-(~b-) for ~b+ satisfying the interface-continuity condition and the perfect reflection condition. The elements used for the representation of ~b-+ are consequently conforming. Here, by means of a specific application, the "legality" of the use of nonconforming elements with K + (~b+, q~-) principle is illustrated. The use of nonconforming elements with classical variational principles constitutes a "variational crime" in the sense of Strang and Fix (1973). If the nonconforming elements satisfy

Finite element solutions for neutron transport

415

the Iron's patch test (Wait and Mitchell, 1985), useful results can be obtained. The use below of nonconforming elements with K + (4)+, qS-) show that with a particular choice of weights one obtains the well-known difference equations of Fletcher (1974) for the even-order spherical-harmonic moments of the angular flux. These equations were obtained originally using the calculus of finite differences. The relative weighting obtained for the nonconforming elements provides a trial value of 2 when the K + -(~b +, q~-) principle is used with elements which include a conforming subset.

6.2. Model problem for a homogeneous slab Consider a system covered with the linear mesh of Fig. 12. Instead of maximizing K+-(q5 +, ~b ) the equivalent minimization of/r(q5 +, q~--)+21d(q5 +, q5 ) is undertaken. The equivalence stems from identity (13) and expression (14). The homogeneous slab has bare surfaces at its ends, so 2Id(~b+, ~b ) reduces to the functional If(~b+, ~b ) given by expression (11). The parity equations for the homogeneous slab are :

d4,o

# - d x + C¢+ = S + '

(a)

l (81)

dq~+

~x+G-'q~o=S-,

(b)J

for the exact angular fluxes ~bff (x, #). Since the slab is homogeneous, C and G - ~ are independent of position. The definition (15) of the functional I,(~b+, ~b ) gives:

xG[Y~x

+G '~--S-]}dydx+fo!P{[dP+(O,P)+q 5 (O,p)]2+[~b+(a,#)-~b (a,p)]2}dp

- t "1 o ~{[¢+(0'~)-q~ (0'")12+[q~+(a'~)+¢-(a'")]2}d"' (82) and Ir(q~+, qS-) is defined for slabs by expression (11).

6.3. Trial functions and a weightin9 scheme The trial function qS+(x, y) is taken to be of the extreme kind: q~+(x,p) = a,(#)

¢*(X,I-L)

for t h e / t h element

(xl, xi+ j),

(83)

m ,i]{p.) =~" nimP2m(l.i) m=O !

I

F

I ! I

:-o

i

x~

I

I÷!,,.., :

I I I

I I l

: ! I ! I

! Oitl~l !

I

! l I

,L

1

.L

~

-

-

zi

xi\,

.L -

¢+(x,p.) for IJ. fixed Bare surfaces

! I

,

I

.k --

I

!

.t v

,

,k.

X'N

1 v

-

X

XN÷ I =n

A -I

Fig. 12. Specification of legitimate nonconforming elements for the K+-(0 +, ¢ ) principle.

416

R . T . ACKROYD and W. E. WILSON

where ai are given by Legendre expansions : M

ai(p) = ~ aimP2,,,(P), (i = 1..... N)

(84)

m=0

with arbitrary coefficients

aim.At

the slab ends the conditions imposed are : ~b+ (0, p) = a, (p),

c~+(a, p) = aN(P).

(85)

The trial function q~-(x, p) is specified for the element (x~, xi+ l) by:

bi(p)+G-lq5 - = S-, (i= 1..... N).

(86)

The function bi(p) is regarded as an approximation for p (d~b0/dx)(x, #) in the ith element. The specification (86) for q~- is suggested by the form of the second of the parity equations (81). A t the ends of the slab, the specification is defined for intervals [(x~, xz) and (XN, XU+ 0]'

6.4. Variationallyderiveddifference equations A two-stage process is used to determine the functions ai(p) and bi(p). In the first stage, the coefficients a~m in expansion (84) are determined with all bi(p), save bl(P) and bu(p), taken to be independent of a~(p). The pairs at(p) with b l(p) and aN(p) with bN(p) are chosen so that the angular fluxes ~b+ (0, p ) + q~-(0, p) and q~+ (a, p) + gb-(a, p) satisfy bare surface boundary conditions. The functions a~(p) are then determined by minimization of the functional Ir(q~+, ~b-) + 2Id(gb +, gb-) with respect to coefficients at,,, having specified the relative weighting 2 and ff+(xi, p). The weighting scheme proposed for an isotropic scattering medium is : 1

ah ~ + = Ipr.

Having determined aim, and thus ai(p) as local p (dc~o/dx), are calculated in the second stage.

Ppr

i.e. w + (x, p) = a h '

(87)

approximations for q~+(x, p), the bi(p) as approximations for

F o r the first stage the bare surface boundary conditions are satisfied exactly by defining bt(p) and bu(p) as follows. Put :

bl(P)= {S-_~S-G lal(p), p0 ] and

bN(p) =

~S-_-G-lau(p), +G-lax(p),

,

(88)

•

(89)

p>0 p < 0

to give ~b+(0,p)+~b-(0,p) = 0,

p > 0

(a)

gb+ (0, p) -- ~b- (0, p) = O,

p < 0

(b)

~ + ( a , p ) - ~ b - ( a , p ) = 0,

p > 0

(c)

fb+(a,,u)+~cp-(a,p) =

p < 0

(d)

0,

Results (89a) and (89d) show that the trial functions ~b-+ satisfy at the ends of the slab the bare boundary conditions for inward directions. Minimization of/~(~b ÷, q~ )+2Id(q~ ÷, ~b ) with respect to the a ~ involves, for the element i, the interval (xl, xi+ 1) and the interfaces xi and xi+ l, provided the element does not abut a bare surface.

Finite element solutions for neutron transport

417

The quantity to be minimized is therefore :

I_l #2

M

2

~{I~=o(ai,,,--ai

M

2

l.m)P2mOt)1 JrI~__o(ai+l,m--aim)P2,,,(#)]}d# Jr

~ ~fl~,ax,

0=

1

~ - e2k(P)

m=0

1 I

aimCP2m(#) - S + L

X CMinimization with respect to

-}-

d-I

=0

J

J

'[-IdGS-+(~a,.,Cp2.,(#)l_S+]+[b,]G[b,])d#dx" L# ~ - x

~=0

/

(90)

j

aik gives :

(a,+ l,m-- 2aimjr a,_ l,,.)P2m (#) +

I # ~ + _ i

P~(#)dx d#

( k = 0 ..... M),

(91)

( k = 0 . . . . . M).

(92)

m=0

because C is a self-adjoint operator. This result can be written as : G [• + (Xi+ 1, #) -- 2(j~) +

0 =

(Xi, #) Jr 1~+ (Xi_l, 1 ~",+,r

h

J.L

Jr C~) q- (Xi, #)

#)]

dGS-] "] S+(x,#)-#~(x,#)Jdx~d#

Since the even-parity equation for ~b~-,Ackroyd (1981), can be written for a homogeneous slab as : d 2a'+ dS-#2G~d---xZ° + C ~ b ~ - - S + + # G d~x = 0,

(93)

a result of the kind (92) can be obtained alternatively by integrating the finite difference approximation :

_#2

+

+

~ - G[q~0 (xi+,,#)-Zq%

+ (xi,#)+q~o (xl

+ + dGS(xi,#)-S (x~,#)+#~d~-x (X,,#) ~ O,

1,#)]+Cq~0

(94)

weighted by P2k(#), This alternative method was first used by Fletcher (1974) in the finite difference version of the sphericalharmonic treatment of neutron transport in the M A R C code.

6.5. Boundary conditions for the difference equations Boundary conditions for the difference equations (92) are obtained as follows. For the bare surface boundary x, = 0 the coefficients a,,, of aft#) occur in general in the surface terms of Id(q~+, q~-) for x, and Xz, and in the volume term of Ir(q~+, 4)-). The choice of b l(#), (88), leads to the conditions (89a) and (89b), which make the contribution of the surface x~ = 0 to Ir(q~+, ~b-) vanish. The condition (88) implies with equation (86) for i=1:

~-a,(~), #>0

@-=(a,(#), #<0

and0 ~< x ~< x2.

(95)

Consequently, the functional to be minimized with respect to a~ i s : f,#2

_ a,,,)p2,,(#)t2 d#

+

I:21__

I f [m~=oa lrnfe2m(#) l - S + C- Im~O a lmf e2m(#) - S+ l l

'\+ {G-'[-a,(#) sgn#l-S-}G{G-'[-a,(#) sgn#]-S }/

d# dx,

(96)

418

R.T. ACKROYDand W. E. WILSON

where it will be recalled that M

a,(t~) = ~ almP2m(~). m=0

A result of the same kind as the boundary condition (96) can be obtained in the same way for the boundary x=a.

These boundary conditions at x = 0 and x = a differ from those used by Fletcher (1974), where the treatment of bare surfaces is approximate. In the present treatment, a condition is imposed on ~b- so that the trial function ~b+ + q~ , for the angular flux ~b0~ + ~bo satisfies the bare boundary conditions exactly. For practical applications the scheme (87) can be replaced by the estimate : 1

ah' if+ = 1.

(97)

6.6. The odd-parity approximation The odd-parity ~b could be determined by assigning a weighting to the interfaces, and minimizing the functional Ir(~b+, q~-) + 2/0 (4)+, q~-) with respect to coefficients b#, (say) specifying the bi(#) by means of oddorder Legendre expansions. A more physically appropriate way for determining b~(#) is perhaps as follows. In introducing equation (86) b~(/~) were regarded as an approximation for/~ (d~b~-/dx) (x, t~). Approximations of the kind : bi(/~) =/~[ai+, (p) -a~(#)]/[xi+, - x i ]

for (i = 2 . . . . . N - 1),

(98)

are suggested for further investigation. 7. DISCUSSION

7.1. Pros and cons o f the composite method The examples given to illustrate the method of composite solutions show that the method has practical potential for reactor physics and shielding problems. The method provides a convenient means of varying, region by region, the order of the spherical-harmonic expansion of an accurate finite-element approximation for the angular flux. Application of the composite method is suggested for the reactor and shield problems described below. Transport treatments of control rods in which a high order approximation for the angular flux is used everywhere in the reactor, whether of the spherical-harmonic or discrete ordinate kind, involve a lot of unnecessary calculations. Transport effects within the reactor core are confined to zones surrounding control rods, bands straddling the core-reflector interface, and possibly the interface of different enrichment zones. Outside the above zones and bands the diffusion theory approximation holds sufficiently well. The calculation of neutron streaming along ducts in shields using finite-elementmethods has been undertaken iteratively by Ziver et al. (1986) and Ackroyd et al. (1986) have used iterative and extrapolation methods. For both the iterative and the extrapolation methods for neutron streaming problems, two or three transport calculations have to be undertaken. There are zones surrounding the ducts where the effects of neutron streaming in the ducts has its effect. Outside of these zones the attenuation of the neutron flow through the shield is practically independent of the presence of the ducts. A P3 calculation is often adequate for the bulk shield, but a P5 or P7 calculation may be required for the duct zone. The bulk shield P3 calculation in (X, Y) geometry in the above example requires the determination of 4 moments per node per energy group, whereas the P7 calculation for the zone needs 16 moments per node per energy group. The composite method would permit a P3,7 approximation to be used for the duct problem instead of a P7 approximation throughout the shield. A practical drawback to the present version of the composite method is the need for two trial functions. One way of overcoming this drawback is to make the particular choice q5 = ~bo, the exact odd-parity angular flux. The result is the K~ 0(q~+) principle, as briefly described by Ackroyd (1986). The trial function q~+ does not have to satisfy any interface or boundary conditions. Another possible treatment using one trial function in the composite method is described below. It involves the use of the smoothing operator used by Ackroyd

Finite element solutions for neutron transport

419

et al. (1987), in a survey of projection and conservation methods in a variational treatment of neutron transport. 7.2. Alternative composite method using a single trial function An alternative to the use of K2-°(q5 +) is sketched below. It has the merit that the bare surface boundary condition can be satisfied exactly. Consider those nodes of a system which do not lie on a bare surface. Using this set of nodes the trial function q~ is made dependent on ¢+ by choosing: d? = $ G [ S - - ~ .

Vq~+],

a choice prompted by the form of the parity equation :

~o =G[S - ~ . V ~ ] , in the notation of Ackroyd (1981). The symbol $ denotes a smoothing operator, which ensures that ~b- satisfies the continuity condition for interfaces and the perfect reflector condition. Consider a node r, shared by m elements which employ spherical-harmonic expansions for q5+ of the same order. F o r each element sharing the node r, a value of G [ S - ( r , , ~ ) - f~-V~b + (r,, ~)] can be calculated. Take th-(r,, 1"~) to be the average of these m values. When the elements sharing r, employ different sphericalharmonic expansions for q~+, the values of the G[S (r,, ~ ) - ~ • V4~+ (r,, ~ ) ] are calculated taking into account those orders of spherical harmonics which are c o m m o n to all the nodes. Nodes on the base surface qS- can be defined in terms of q5+ by generalizing the scheme outlined for slab systems by equations (88) and (89). 7.3. Mesh refinement with increased weighting The way in which weighting influences approximate solutions has been explained for a fixed mesh, and the findings have been confirmed for representative problems. The derivation of difference equations in Sections 6.2-6.4 provides an estimate of the relative weighting (97) for slab systems. Further work is required to justify rules of the kind (97) for 2- and 3-D systems. REFERENCES

Ackroyd R. T. (1981) Ann. nucl. Energy 8, 539. Ackroyd R. T. (1982) Ann. nucl. Energy 9, 95. Ackroyd R. T. (1983) Ann. nucl. Energy 10, 65. Ackroyd R. T. (1986) Prog. nucl. Energy 18, 45. Ackroyd R. T. and Splawski A.B. (1982) Ann. nuN. Energy 8, 717. Ackroyd R. T. and Wilson W. E. (1986) Pro9. nucl. Energy 18, 39. Ackroyd R. T., Issa, J. G. and Riyait N. S. (1986) Prog. nucl. Energy 18, 85. Ackroyd R. T., Issa J. G., Riyait N. S. and Nanneh M. (1987) Projection and Conservation Methods for Neutron Transport. MAFELAP, Brunel University. Dunderstadt J. J. and Martin W. R. (1979) Transport Theory. Wiley-Interscience, New York. Fletcher J. K. (1974) TRG Report, 2849(R), UKAEA. Galliara J. and Williams M. M. R. (1979) Ann. nucl. Energy 6, 205. Garcia R. D. M. and Sievert C. E. (1980) Nucl. Sci. Engng 53. Martin W. R. and Duderstadt J. J. (1977) Nucl. Sci. Engng 62, 371. Mordant M. (1975) IAEA Mtg, Bologna, Italy, Nov. 3-5. Nelson P. (1975) Nucl. Sci. Engng 56, 340. Pitkg_ranta J. (1977) Trans. Theory Statist. Phys. 6, 169. Pomraning G. C. (1966) J. nucl. Energy A/B 20, 617. Reed W. H. (1971) Nucl. Sci. Engng 46, 309. Splawski A. B. (1981) Finite element methods for neutron transport calculations. Ph.D. Thesis, University of London. Strang G. and Fix G. J. (1973) An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, New Jersey. Wait R. and Mitchell A. R. (1985) Finite Element Analysis and Applications. Wiley, New York. Weinberg A. M. and Wigner E. P. (1958) The Physical Theory of Neutron Chain Reactors. University of Chicago Press. Wilson W. E. (1985) Discontinuous finite elements in neutron transport analysis. Ph.D. Thesis, University of London. Ziver A. K., Issa J., Riyait N. and Watmough M. H. (1986) Prog. nucl. Energy 18, 215.