- Email: [email protected]
Solutions of the diffusion equation in cones and wedges
Ann.nucl.Energy,Vol. 15, No. 4, pp. 191-199, 1988
0306-4549/88 $3.00+0.00 Copyright © 1988 Pergamon Press plc
Printed in Great Britain. All rights reserved
SOLUTIONS OF THE DIFFUSION EQUATION IN CONES A N D WEDGES E. M. RODRIGUEZ VENTURA, l M. M. R. WILLIAMS 2 and J. W O O D l Nuclear Engineering Department, Queen Mary College, University of London, Mile End Road, London E1 4NS U.K. and 2Nuclear Engineering Department, The University of Michigan, Cooley Building, N. Campus, Ann Arbor, MI 48109-2104, U.S.A.
(Received l August 1987; in revisedform 29 September 1987) Abstract--The neutron diffusion equation has been solved for cones and finite wedges. Criticality conditions are derived linking the dimensions of the cones and wedges to the nuclear properties of the system. The fundamental mode flux is evaluated numerically for a number of cases. We have also discovered an interesting minimum in the leakage rate as the cone or wedge angle passes through a certain value for a fixed volume. Some practical implications of this result are pointed out. The exact solutions derived have been used as benchmarks to assess the accuracy of the finite element code TRIPAC in R-Z and X-Y-Z geometries for a four-group criticality problem in cones and wedges. Excellent agreement is found between the analytical results and the code for both eigenvalues and fluxes.
l. INTRODUCTION
(called the diffusion equation in reactor physics) in
Solutions of the diffusion equation in nonstandard geometries are ofintrinsic interest but also ofpractical value since the results may be used to check the accuracy of computer codes. With this in mind we present some solutions of the equation (Glastone and Edlund, 1953): V2~b + B 2~b --- 0,
(1)
subject to zero flux on the boundary for nonreentrant cones and wedges. Our solutions are given analytically and we present the flux distribution and the critical condition in terms o f the physical dimensions of the cone or wedge. These results are employed to obtain a comparison with results from the finite element code T R I P A C and lend further support to its accuracy, In addition to obtaining analytical solutions, we are able to show that, for a given buckling, the critical cone volume and corresponding surface area exhibit a minimum for a certain value of the cone angle, Similarly, for a fixed critical volume, the buckling has a minimum at a particular cone angle. Analogous behaviour arises for wedges, 2. SOLUTION OF THE DIFFUSION EQUATION FOR A BARE CONE
2.1. General solution and eiyenvalues Following the usual procedure (Morse and Feshbach, 1953), we can write the Helmholtz equation
spherical coordinates as : 1 O [ 2g~b'~
1
~
sm0~0
+B ~b=0,
r ~ r ~ r ~r)+r2sinO~O (2) where we note that ~b = e~(r,O), 0 being the angle measured with respect to the z-axis. There is no dependence on the azimuthal angle due to symmetry. The boundary conditions are : (a) the flux is everywhere finite,
(3)
(b) q~(r, 00) = 0, (c) qS(Ro,O)= 0,
(4) (5)
where R0 is the cone radius and 00 is its half-angle. We restrict 00 < n/2 because we wish to avoid reentrant conditions. However, values of 00 > n/2 could have physical significance if it is assumed that the volume contained within the reentrant line-of-sight region is black to neutrons. Using the standard procedure of separation of variables, we can write the solution to equation (2) which satisfies condition (3) as : A
~b(r, 0) = ~ J ~ + 1.2(Br)Pv(cos0),
,/r
(6)
where A is an arbitrary constant and v is a parameter to be determined. J,+ t/2(Z) is a Bessel function and P~(cos0) is a Legendre polynomial. 191
E . M . RODR1GUEZ VENTURA et
192 B o u n d a r y c o n d i t i o n (4) leads to :
Table la. Zeros of Bessel function of first kind of order 3~/40o
P~(COS0o) = 0; all v/> 0.
(7)
The roots of this e q u a t i o n lead to the a p p r o p r i a t e value o f v. In fact unless v is a n integer there is a n infinite n u m b e r o f solutions, b u t we are interested only in the smallest root since this c o r r e s p o n d s to the f u n d a m e n t a l mode. T h e precise roots of e q u a t i o n (7) require some effort to o b t a i n b u t a good first approxim a r i o n arises from the following asymptotic representation, viz. ( A b r a m o w i t z a n d Stegun, 1965) : 2
ev(cos0)
~
L~3
762 F(Vq'- 1)
r~
x sin
(v+~)O+ ~ + O ( v -3'2)
(8)
where v ~> 1 a n d 0 is confined to the interval (e, i t - e) where e is a small n u m b e r . T h e greater the value of 0, the more accurate is e q u a t i o n (8). Neglecting the term o f O(v 3/2) a n d applying condition (7) we find n (v + ½)00 + ~ = nn, n -- 0, 1,2 . . . . (9) T h e f u n d a m e n t a l m o d e arises w h e n n = 1 (n = 0 corresponds to negative v) a n d thus : 3rt v - 40o
1 2"
(10)
By c o m p a r i s o n with an exact calculation, it m a y be s h o w n t h a t this formula is accurate to better t h a n 2 % even for 00 = 1.8 °. It becomes progressively more accurate as 0o increases being exact for 0 o = n/2. T h e flux c a n n o w be written as :
A dp(r,O)-=-~rJ3~(Br)e3~ 400
1(cos 0). 400 5
3tt/40o
Oo
3/2 2 5/2 3 7/2 9/2 11/2 6 7 8 25/2 39/2
~z/2 3~/8 3n/lO 7r/4 3n/14 r~/6 3rr/22
n/8
3g/28 3z/32 3n/50 3r~/78
(BRo)exact (BRo)approx. 4.493409 5.135620 5.763459 6.380160 6.987932 8.182561 9.355812 9.936110 11.08637 12.22509 17.250455 24.878005
4.516 5.141 5.756 6.365 6.967 8.158 9.333 9.915 11.072 12.220 17.291 24.979
The accuracy of this relationship is better t h a n 0.5% t h r o u g h o u t the range of values o f 00 s h o w n in Table 1. As an indicator of the overall accuracy of equations (10) a n d (13) we note t h a t a n accurate d e t e r m i n a t i o n of BR o for 00 = 1.8 ° leads to BRo = 80.53, whilst the a p p r o x i m a t i o n yields BRo = 82.39, a n error of 2.3%. As 0o increases, this error progressively becomes smaller, e.g. at 0o = 13.22 °, the error is 2.1% a n d at 00 = 54.73 ° it is - 0 . 7 % . Thus o u r empirical result would a p p e a r to be a useful a p p r o x i m a t i o n for bare cones. In some calculations, it is m o r e practical to select a n integer value of v a n d find the c o r r e s p o n d i n g value of 0o from the e q u a t i o n : P~(cos 00) = 0. Such roots are simply the G a u s s - L e g e n d r e quadrature p o i n t s a n d are k n o w n very accurately. W i t h a n integer value of v it is easy to find the roots of
Jv+t/2(BRo), (11)
Finally, using e q u a t i o n (5), we have :
J3n (BR0) = 0,
al.
(12)
from the tables in A b r a m o w i t z a n d Stegun (1965). We shall a d o p t this procedure in our c o m p a r i s o n with the finite element m e t h o d to be described below a n d give some illustrative results in Table 1b.
40~ which yields the critical c o n d i t i o n between B a n d the physical dimensions of the cone. T h e roots of e q u a t i o n (12) have been f o u n d from tables of Bessel functions ( A b r a m o w i t z a n d Stegun, 1965) a n d are listed in Table la. There is a n infinite n u m b e r of such roots for a given 00 b u t we are only c o n c e r n e d with the one t h a t c o r r e s p o n d s to a n o n negative eigenfunction. A useful empirical relation between BRo a n d 00 is f o u n d to be :
BRo = 3.07250ff°'9415 +2.5076.
(13)
Table lb. Exact roots and half-anglesfor cone v 1 2 3 4 5 6 7 8
BRo 4.493409 5.763459 6,987932 7,58834 8.77148 9.93611 11.08637 12.22509
00 90 54.73561 39.23152 30.55559 25.01733 21.17690 18.35785 16.20078
Diffusion equation in cones and wedges
2.2. Optimum size o f cone
3o
Let us assume that we have a cone of fixed buckling, i.e. fixed net leakage, and require to know how the volume and surface area vary with 0o. The volume of a cone is given by the expression:
25
2nR 3
v = - ~ - - ( 1 - c o s 00).
(14)
Using the empirical relationship (13) for Ro as a function of 0o (remember that B is fixed), we find that the volume is : 2~
'
193
3
V(Oo) = 3 ~ ( ~ 0 f f " + ~ ) (I --cos0o)
% 2o
(15) 15
where ~, fl and 7 are given in (13). It is readily shown that V(Oo) has a minimum at 00 = 0.5259 (30.132°). Similarly, the surface area of the cone which is given by :
S(Oo) = 2~zR2(l --cos 0o) +TzR ~ sin 0o
(16)
also has a minimum at 00 = 0.559 (32°). Both V(Oo) and S(Oo) are shown in Fig. 1 as functions of 0o. Let us now consider a cone of constant volume and examine the buckling as a function of 00. Using the earlier formulae we readily find that : x~2/3 B2(0o) = \ ~ j
(~0J+7)2(l-cos0o)
2'3. (171
I 1o
2oI
I I I 30 4o 5o 8otdegr~sl
I 6o
I 70
1 8o
Fig. 2. The variation of the buckling with cone half-angle, 00, for a fixed volume; the ordinate is in arbitrary units.
2.3. Flux distribution in a cone To illustrate the flux shape in the cone, we consider the following parameters: 0o = ~/6 (30 °) and hence BRo = 8.182561, from Table 1. We also note that v = 4 and thus the flux distribution is from equation (6):
This is illustrated graphically in Fig. 2 and exhibits a minimum at 00 ~- n/6. Clearly there is a critical value of cone angle which minimizes the leakage. The possibility of using such a phenomenon as a reactor control device should be noted.
220 21o 200
I I
-
VtOo)
-----
S(eol
I
"( n "y,,2 ~b(r, 0) = \2B-rrJ Jg/2(Br)P4(cos 0).
The flux is given in Table 2 for 0 = 0 for various values of the ratio r/Ro. We note that a maximum arises at r = 0.685R0. It is also interesting to examine the angular variation of the flux which follows the function P4(cos 0). The values o f P 4 ( c o s 0) for a range of
Table 2. Flux along central axis of cone,
v/Ro
4~(r,O)
0.9t7 0.856 0.733 0.685 0.672 0.611 0.489 0.367 0.244 0.122 0.0611
0.0793 0.1327 0.1968 0.2015 0.2008 0.1870 0.1249 0.05615 0.01408 1.011 x 10 3 6.539 x 10 ~
180 1
170
~ \
/ U \
160
/
~k
" ~ /
I 20
I 30
J
150 I 10
J 40
I 50
I 60
J 70
I 80
I 90
Degrees
Fig. I. The variation of the critical volume and surface area with cone half-angle 00 ; the ordinate is in arbitrary units.
(18)
0.0
0
0.0
E . M . RODRIGUEZ VENTURA et
194
Table 3. Angular flux across cone 0
of radius R0 a n d height H. The zero flux b o u n d a r y conditions are :
P4(cos 0)
o.o
l.O
5
0.9623
20
0.4750
10
0.8532
25
0.2465
30
0.0234
30.55559
al.
0.0
0 are s h o w n in T a b l e 3. T h e flux is n o t precisely zero at 0 = 30 ° due to the a p p r o x i m a t e f o r m u l a for v of e q u a t i o n (10). In fact, it is zero at the smallest root of P 4 ( c o s 0 ) = 0, i.e. 0 = 30.55559 °. The fluxes are illustrated graphically in Fig. 3.
(a)
qS(r, O, z) = O,
(20)
(b)
q~(r, 00, z) = 0,
(21)
(c)
q~(R0, 0, z) = 0,
(22)
(d)
qS(0, 0, z) = 0,
(23)
(e)
~b(r, 0, 0) = 0,
(24)
(f)
0(r, 0, H ) = 0.
(25)
U s i n g separation of variables, we find t h a t the solution can be written :
dp(r,O,z) = CJu,(Brr)sin(plO)sin(B,z),
(26)
P l = X/0o
(27)
where
3. S O L U T I O N OF T H E DIFFUSION EQUATION F O R A
and
BARE, FINITE WEDGE
3.1. Generalsolution andeigenvalues
(H)2
T h e diffusion e q u a t i o n for a wedge of finite height
B2 =
= B2
B 2.
(28)
is:
1 8 [" 8dA
~r~)+~
1 82¢
02,~
a02 +7;~ +B2~=0
B 2 is the radial buckling a n d is f o u n d from c o n d i t i o n
(19)
(23) viz.
J~oo(BrRo) = where ~b = ~b(r, 0, z). W e assume t h a t one face o f the wedge is in the plane 0 = 0 and the other in the plane 0 -- 00. The wedge is
r
0.
(29)
Thus we have a relation between the physical dimensions of the wedge a n d its nuclear properties, B 2. F r o m tables o f zeros of Bessel functions we can o b t a i n Table 4. Again, for convenience in calculation, an empirical relationship between BrRo a n d 00 has been o b t a i n e d thus :
10
BrRo =
4.048800 0.9378 _.~2.50.
(30)
Table 4. Zeros of Besssel function of first kind of order, x/0o .
Fig. 3. Flux distribution inside a cone with 0o = x/6, R0 = 10. The curve on the right hand side of the figure represents the radial flux distribution and that on the left the angular distribution ; actual magnitudes are arbitrary.
rt/Oo
0o
(BrRo) exact
(B, Ro) approx.
1 3.2 2 5/2 3 7/2 4 9/2 5 11/2 6 13/2 7
r~ 2x/3 x/2 2rt/5 x/3 2x/7 x/14 2~/9 re/5 2~z/l 1 ~z/6 2n/13 ~/7
3.83171 4.493409 5.135620 5.763459 6.380160 6.987932 7.58834 8.182561 8.77148 9.355812 9.936110 10.512835 11.08637
3.884 4.524 5.151 5.768 6.377 6.981 7.578 8.171 8.760 9.346 9.928 10.507 11.083
15/2 8
2n/15 r~/8
11.657032 12.22509 17.250455 24.878005
11.657 12.228 17.284 24.933
25/2 39/2
2r¢/25
2rt/39
Diffusion equation in cones and wedges The accuracy of the expression may be judged from Table 4. For 0o < 120 ° it is better than 0.7%. 3.2. Optimum size o f wedge As in the case of the cone, we look for an angle which minimizes the buckling for a given wedge volume. Using equation (30) and the volume of the wedge, we get: H 0 ob +c)-, B~(Oo) = 2~Oo(a
(31)
where a, b and e are defined by equation (30). As Fig. 4 shows, the radial buckling has a m i n i m u m and by differentiation of equation (31) we find: -
-
= 83.15 °.
(32)
3.3. Flux distribution in a wedge For 00 = ~z/4, we find that the flux distribution can be written: rcz) (o(r, O, z) = J4(Brr) sin (40) sin ~ (33) where B_ = 7t/Hand BrRo = 7.58834. The maximum value of q~ along the radial direction arises when: d J4 (Brr) _ O, dr
(34)
which leads to Brrmax = 5.31755, i.e. rmax ~---0.70075R0. Of course, in the axial direction, the maximum is at z = 11/2 and in the 0-direction at 0o/2.
15e
~
-,
100
195
4. NUMERICAL VERIFICATIONOF CODE AND THEORY TO verify the theoretical relationships derived in the previous sections, particular examples of the cone and wedge configurations are solved by means of a multigroup finite element transport code called T R I P A C . This code is a development of two earlier codes, namely F E L I C I T (Wood and Williams, 1984; Wood, 1986) and E V E N T (De Oliveira, 1986, 1987). A feature of T R I P A C , as of its precursors, is its ability to handle complex geometries, and the two configurations considered in this paper provide an interesting test of this aspect of the code's capabilities. The numerical solutions that we give also provide useful benchmarks for proving other techniques and c o m puter codes for which the claim of geometrical flexibility is also advanced. In this paper we are concerned with diffusion theory, i.e. P,/P~ solutions, but of course T R I P A C is also capable of providing higherorder transport solutions, namely, P,/P/. Two cases of the cone are considered, and two of the wedge. Four-group linear anisotropic scattering data is used for the homogeneous systems. The data is based on the simple H T G R core data used by White and Frank (1987). Our data is listed in Appendix A ; the corresponding value of materials buckling is Bm2ax= 1.9701604 × 10 2. In the computations, iterations are continued until the eigenvalue (kc~) and eigenvector have converged, respectively, to 10- 6 and 10- 5. Also the eigenvector is normalized to unit fission neutrons produced in the system, i.e. 1 ['~' k ~ J J d E d ~ ' ~ Z j ~b(E, r) = 1. In the case of the wedge, the eigenvalue is actually normalized to 1/'4 of the total volume--this is explained below. It is perhaps worth mentioning that this particular data has the fundamental and next lower eigenvalue close together (White and Frank, 1987). Thus, unless a power iteration convergence acceleration scheme is employed in the code, a large number of iterations may be necessary to produce the accuracy we quote in T R I P A C , outer iteration acceleration is achieved by Chebyshev extrapolation. 4.1. The cone
50 I t I I I t L I I I I t I I I J ~o 40 6o 80 leo lzo 14o 16o 8 0 (degrees)
Fig. 4. Variation of radial buckling B,2 with wedge angle 0o for constant volume; the ordinate is in arbitrary units,
The cone is solved in R - Z geometry using 81 quadratic elements. The element mesh employed is illustrated in Fig. 5. The exact and T R I P A C results are compared in Tables 5 and 6. Clearly, the agreement is excellent. The flux distribution along the z-axis of the smaller cone is shown in Fig. 6. Figure 7 shows the variation of critical volume of a cone for fixed buckling as a function of the cone
196
E . M . RODRIGUEZ VENTURAet al.
Z!.
13.7.87
Cone 4-Grp
.~'-;~-.-.~'.i
R-Z Diff. AppCone 21 =
" ---
x
~ -'2
¢qoo
. /
f::
:: '.?
i" i 'Z/'
:."f)
'"
"'}:"
~ c~ ,::; :
, ii"
•
!
x ,~.
rq ~ o¢~
::
8
..?'
::= ~
~_
~
~- o~
emoo
i./.
i] ,:::
i - - . . . - : ; -.- - •.. !!:
"! e! o o
i.
/11 ii?-:; !,::::.:.:::ii:::::"::..... '
X i..:::" ii ,:-
~?
~
•
x
o o ~
x
~
~,
,~o,~ e~ r~
?,
"].._........i:
I '
R
0
Progro,,'n P o s t - - p r o c e s s ~ F i g . 5. F l u x
contours Contour
Mo~'ch ~9 , 9 7
for group increment
4 neutrons i n 1 x 10
in
21 ° c o n e .
4.
~rq
half-angle. This is a four-group calculation using the data in Appendix A ; it confirms the analytical result in Section 2.2 and indicates an optimum half-angle of about 30 °.
~ ~
~ '~ ~. "~'t~
4.2. The wedge This problem is solved in X - Y - Z g e o m e t r y . T h e basic element in the X - Y plane is an arbitrarily ori-
(~
~ .~
Diffusion equation in cones and wedges
197
Table 6. Comparison of critical flux distributions in 54 ° cone, group 4 neutrons
r
z
TRIPAC
0.0 0.0 0.0 0.0 0.0 0.0 0.0 95.777 126.706 152.234
390.082 307.960 287.429 246.368 184.776 123.184 61.592 248.827 160.352 336.739
x
105
EXACT
0.1324 0.6317 0.7155 0.7991 0.7031 0.4173 0.1234 0.6213 0.3215 0.1960
105
x
~.
0.1294 0.6271 0.7116 0.7970 0.7020 0.4173 0.1231 0.6208 0.3215 0.1969
~_
Note: cone is solved by TRIPAC in R Z geometry.
ented triangle. Again, a quadratic triangle is used. The subdivision of the X-Yplane is shown in Fig. 8. These triangles form the cross sections of prisms whose axes are parallel to the z-axis. Each prism so formed is further subdivided into 3 tetrahedra. Thus the basic element in space is, in this case, a quadratic tetrahedron (with 10 nodes). Because of symmetry, the problem need only be solved in a 1/4 wedge by taking the origin to be at the mid-height of the wedge and exploiting symmetry boundary conditions along appropriate axes. As with the cone, the boundary condition on the surface of the wedge is zero flux. The exact and TRIPAC results are compared in Tables 7 and 8. In TRIPAC, 735 quadratic tetrahedra are used for each wedge. An isometric view of the flux in a critical wedge is shown in Fig. 9. For the wedge, as for the cone, the agreement between theory and TRIPAC results is excellent. Appendix B indicates how the volume-averaged fluxes are calculated for cones and wedges.
0.9
,,Dt~ ~
tt'5
,.~
~Le~ m m:~ .--: ~ ~
~ × ~=
~ ~ eqtt5
~ ~"
~5
-
~2
~;
~s
=
~
~
"~ "~ .~ _~ .~
~, ~ ~~
"~
b-
13" 7 " 8 7 Cone 4 - G r p R - Z Diff. App Cone 21 ° Tota I. flux or dose rate
0.8 0.7 0.6
o ×
~=
05 0.4
E
0.3
0.2 0.1
''
i 0
I00
200
300
400
500
600
700
I 800
rj
Z
Program Post-processor March 1 9 8 7
Fig. 6. Variation of flux along z-axis for 21': cone (group 4 neutrons).
~
~ x~
198
E . M . RODRIGUEZ VENTURA et al. Critical volume versus cone ongte 7.0
20'7.87 Wedge 4-Grp X - Y - Z
6.5
~~
~
Diff. Approx.90"
6.0
5,5 I 40
5.(
I 20
I 30
I 40
I 50
i
60
I 70
I 80
YI/~X/~//~//~
I 90
'
/
\
/
'
/
'
/
'
0
Fig. 7. Variation of critical volume of cone with half-angle 00 for fixed buckling using four-group data.
IRefteetion
X
GEM Vr 1.9.86 Fig. 8. Subdivision of X Y plane into quadratic triangles for 90 ° wedge. Because of symmetry, only 1/2 of X- Y plane need be solved and 1/4 of volume of wedge.
5. CONCLUSIONS AND SUMMARY
Some exact solutions of the diffusion equation in cones and wedges have been obtained leading to critical conditions and fundamental mode flux distributions. We have observed that the volume and surface area of cones and wedges go through a m i n i m u m at a certain value of cone or wedge angle for a fixed buckling..Similarly, the buckling has a m i n i m u m for
14.7-87 Wedge 4 - G r p X - Y - Z Scalar flux group 4
fixed volume at a critical angle. The latter condition implies a m i n i m u m leakage condition and has implications for reactor control and storage of fissile material. We have also used these exact solutions to assess the accuracy of the T R I P A C finite element code. A
Diff.Approx.45 °
2.75
[- 2.75
225~
FV 225
1.75~
7 1.75
1.25 ~
!
0.75 --
['-
0.25
~- 0.25 0.0
D
1.0
2.25
2.0
1.75 3.0
1.25 4.0
102 Z axis x 102 Vertical axis x 10 -5 R
aXiS
X
0,75
5.0
0.25
6.0
R
Program Post- processor March 1987
Fig, 9. Isometric view of flux in 45' wedge, for the plane Y = 0 (group 4 neutrons).
1.25
Diffusion equation in cones and wedges
199
Table 8. C o m p a r i s o n o f critical flux distributions in 90 ° wedge, g r o u p 4 neutrons x
y
85.714 171.429 257.143 342.857 428.571 514.286 190.536 414.445 355.528 300.917
z
0.0 0.0 0.0 0.0 0.0 0.0 97.0872 99.740 308.062 160.844
T R I P A C x 105
E X A C T x 105
0.3324 1.1782 2.0833 2.5521 2.2798 1.3032 0.6853 1.1970 0.2642 0.9993
0.3370 1.1737 2.0712 2.5390 2.2733 1.3072 0.6826 1.1958 0.2666 0.9964
0.0 0.0 0.0 0.0 0.0 0.0 70.6011 84.7213 0.0 70.6011
N o t e : wedge is solved by T R I P A C in X - Y - Z geometry.
f o u r - g r o u p b a r e r e a c t o r criticality p r o b l e m is solved for a variety o f w e d g e s a n d cones. B o t h eigenvalues
0.28248588E+00 0.121900E-02 0.132000E-02 O.OOO00000E+OO0.00000000E+00 0.25816688E+00 0.23100000E-01
a n d critical fluxes are s h o w n to be in excellent agreem e n t w i t h the analytical results. C o n e s a n d w e d g e s are a sensitive test o f a n y n u m e r i c a l p r o c e d u r e a n d lend f u r t h e r c o n f i d e n c e in t h e a c c u r a c y o f t h e T R I -
0.28985507E+OO 0.a54000E 02 O.OOOOOOO0E+OO0.00000000E+00 0.0 0.0 Normalizedfissionspectrum 0.9675
0.0
0.0
0.025816688E +00 0.0
0.67S000E-02 0.00000000E+00 0.28531507E+00 0.0 0.028531507E+00 0.0325 0.0 0.0
P A C code.
APPENDIX B Average Fluxes
REFERENCES Abramowitz M. and Stegun I. A. (1965)Handbook of Mathematical Functions. Dover, New York. Glasstone S. and Edlund M. (1953) Nuclear Reactor Theory, Van Nostrand, New York. Morse P. M. and Feshbach H. (1953) Methods of Theoretical Physics. McGraw-Hill, New York. White J. R. and Frank B. R. (1987) International Topical
In order to obtain average fluxes it is necessary to calculate :
1 (" ~ = V Jv dr q~(r). For the cone with :
Meeting on Advances in Reactor Physics, Mathematics and Computation, Vol. 3, p. 1385. American Nuclear Society. Wood J. and Williams M. M. R. (1984) Prog. Nucl. Energy 14, 2 I. Wood J. (1986) Ann. nucl. Energy 13, 91. De Oliveira C. R. E. (1986) Ann. nucL Energy 13, 227. De Oliveira C. R. E. (1987) Ph.D. Thesis, University of London.
1 (p(r, O) = ~
Jv+ I/2(B,)Pr( cOs 0),
x/r this expression becomes: 3 /'R q~ - R3(I - c o s 00) J0 drr3/2J~+j/2(B,) x ('1°° dO sin 0 P~ cos 0). ,)0
APPENDIX A E
]E~
E~o(E, ~ El)
E,o(E, --* E2)
For integer values of v, it is easy to carry out the 0 integration and although the integration over the Bessel function can be represented by an infinite sum, it is more convenient to evaluate this directly by quadratures. For wedges we have:
E~, (E~ ~ El)
Y's~(Ei --* E2)
~p(r, O,z) = J~,, (Brr) sin (p. 0) sin (Bzz),
For each energy group i: v-'Zl
Note: diffusion coefficient D = l/3Z,r where E,, = Z - E s l . whence Four-group data
8
0.14184397E+00 0.844000E-04 0.12945957E+00 0.12300000E 01 0.012945957E +00 0.0
0.784000E 04 0.00000000E+00 0.0
0.00000000E+00 0.0
0,27777778E + 00 0.145100E-02 0.00000000E+00 0.27255678E+00 0.0 0.027255678E+00
0.488000E- 03 0.37700000E-02 0.0
0.00000000E+00 0.0
repro
(~--7~2(/~-D0)2[~r'" Jo |
dxxJ~(x),
where t/i = 71"/0o.
0306-4549/88 $3.00+0.00 Copyright © 1988 Pergamon Press plc
Printed in Great Britain. All rights reserved
SOLUTIONS OF THE DIFFUSION EQUATION IN CONES A N D WEDGES E. M. RODRIGUEZ VENTURA, l M. M. R. WILLIAMS 2 and J. W O O D l Nuclear Engineering Department, Queen Mary College, University of London, Mile End Road, London E1 4NS U.K. and 2Nuclear Engineering Department, The University of Michigan, Cooley Building, N. Campus, Ann Arbor, MI 48109-2104, U.S.A.
(Received l August 1987; in revisedform 29 September 1987) Abstract--The neutron diffusion equation has been solved for cones and finite wedges. Criticality conditions are derived linking the dimensions of the cones and wedges to the nuclear properties of the system. The fundamental mode flux is evaluated numerically for a number of cases. We have also discovered an interesting minimum in the leakage rate as the cone or wedge angle passes through a certain value for a fixed volume. Some practical implications of this result are pointed out. The exact solutions derived have been used as benchmarks to assess the accuracy of the finite element code TRIPAC in R-Z and X-Y-Z geometries for a four-group criticality problem in cones and wedges. Excellent agreement is found between the analytical results and the code for both eigenvalues and fluxes.
l. INTRODUCTION
(called the diffusion equation in reactor physics) in
Solutions of the diffusion equation in nonstandard geometries are ofintrinsic interest but also ofpractical value since the results may be used to check the accuracy of computer codes. With this in mind we present some solutions of the equation (Glastone and Edlund, 1953): V2~b + B 2~b --- 0,
(1)
subject to zero flux on the boundary for nonreentrant cones and wedges. Our solutions are given analytically and we present the flux distribution and the critical condition in terms o f the physical dimensions of the cone or wedge. These results are employed to obtain a comparison with results from the finite element code T R I P A C and lend further support to its accuracy, In addition to obtaining analytical solutions, we are able to show that, for a given buckling, the critical cone volume and corresponding surface area exhibit a minimum for a certain value of the cone angle, Similarly, for a fixed critical volume, the buckling has a minimum at a particular cone angle. Analogous behaviour arises for wedges, 2. SOLUTION OF THE DIFFUSION EQUATION FOR A BARE CONE
2.1. General solution and eiyenvalues Following the usual procedure (Morse and Feshbach, 1953), we can write the Helmholtz equation
spherical coordinates as : 1 O [ 2g~b'~
1
~
sm0~0
+B ~b=0,
r ~ r ~ r ~r)+r2sinO~O (2) where we note that ~b = e~(r,O), 0 being the angle measured with respect to the z-axis. There is no dependence on the azimuthal angle due to symmetry. The boundary conditions are : (a) the flux is everywhere finite,
(3)
(b) q~(r, 00) = 0, (c) qS(Ro,O)= 0,
(4) (5)
where R0 is the cone radius and 00 is its half-angle. We restrict 00 < n/2 because we wish to avoid reentrant conditions. However, values of 00 > n/2 could have physical significance if it is assumed that the volume contained within the reentrant line-of-sight region is black to neutrons. Using the standard procedure of separation of variables, we can write the solution to equation (2) which satisfies condition (3) as : A
~b(r, 0) = ~ J ~ + 1.2(Br)Pv(cos0),
,/r
(6)
where A is an arbitrary constant and v is a parameter to be determined. J,+ t/2(Z) is a Bessel function and P~(cos0) is a Legendre polynomial. 191
E . M . RODR1GUEZ VENTURA et
192 B o u n d a r y c o n d i t i o n (4) leads to :
Table la. Zeros of Bessel function of first kind of order 3~/40o
P~(COS0o) = 0; all v/> 0.
(7)
The roots of this e q u a t i o n lead to the a p p r o p r i a t e value o f v. In fact unless v is a n integer there is a n infinite n u m b e r o f solutions, b u t we are interested only in the smallest root since this c o r r e s p o n d s to the f u n d a m e n t a l mode. T h e precise roots of e q u a t i o n (7) require some effort to o b t a i n b u t a good first approxim a r i o n arises from the following asymptotic representation, viz. ( A b r a m o w i t z a n d Stegun, 1965) : 2
ev(cos0)
~
L~3
762 F(Vq'- 1)
r~
x sin
(v+~)O+ ~ + O ( v -3'2)
(8)
where v ~> 1 a n d 0 is confined to the interval (e, i t - e) where e is a small n u m b e r . T h e greater the value of 0, the more accurate is e q u a t i o n (8). Neglecting the term o f O(v 3/2) a n d applying condition (7) we find n (v + ½)00 + ~ = nn, n -- 0, 1,2 . . . . (9) T h e f u n d a m e n t a l m o d e arises w h e n n = 1 (n = 0 corresponds to negative v) a n d thus : 3rt v - 40o
1 2"
(10)
By c o m p a r i s o n with an exact calculation, it m a y be s h o w n t h a t this formula is accurate to better t h a n 2 % even for 00 = 1.8 °. It becomes progressively more accurate as 0o increases being exact for 0 o = n/2. T h e flux c a n n o w be written as :
A dp(r,O)-=-~rJ3~(Br)e3~ 400
1(cos 0). 400 5
3tt/40o
Oo
3/2 2 5/2 3 7/2 9/2 11/2 6 7 8 25/2 39/2
~z/2 3~/8 3n/lO 7r/4 3n/14 r~/6 3rr/22
n/8
3g/28 3z/32 3n/50 3r~/78
(BRo)exact (BRo)approx. 4.493409 5.135620 5.763459 6.380160 6.987932 8.182561 9.355812 9.936110 11.08637 12.22509 17.250455 24.878005
4.516 5.141 5.756 6.365 6.967 8.158 9.333 9.915 11.072 12.220 17.291 24.979
The accuracy of this relationship is better t h a n 0.5% t h r o u g h o u t the range of values o f 00 s h o w n in Table 1. As an indicator of the overall accuracy of equations (10) a n d (13) we note t h a t a n accurate d e t e r m i n a t i o n of BR o for 00 = 1.8 ° leads to BRo = 80.53, whilst the a p p r o x i m a t i o n yields BRo = 82.39, a n error of 2.3%. As 0o increases, this error progressively becomes smaller, e.g. at 0o = 13.22 °, the error is 2.1% a n d at 00 = 54.73 ° it is - 0 . 7 % . Thus o u r empirical result would a p p e a r to be a useful a p p r o x i m a t i o n for bare cones. In some calculations, it is m o r e practical to select a n integer value of v a n d find the c o r r e s p o n d i n g value of 0o from the e q u a t i o n : P~(cos 00) = 0. Such roots are simply the G a u s s - L e g e n d r e quadrature p o i n t s a n d are k n o w n very accurately. W i t h a n integer value of v it is easy to find the roots of
Jv+t/2(BRo), (11)
Finally, using e q u a t i o n (5), we have :
J3n (BR0) = 0,
al.
(12)
from the tables in A b r a m o w i t z a n d Stegun (1965). We shall a d o p t this procedure in our c o m p a r i s o n with the finite element m e t h o d to be described below a n d give some illustrative results in Table 1b.
40~ which yields the critical c o n d i t i o n between B a n d the physical dimensions of the cone. T h e roots of e q u a t i o n (12) have been f o u n d from tables of Bessel functions ( A b r a m o w i t z a n d Stegun, 1965) a n d are listed in Table la. There is a n infinite n u m b e r of such roots for a given 00 b u t we are only c o n c e r n e d with the one t h a t c o r r e s p o n d s to a n o n negative eigenfunction. A useful empirical relation between BRo a n d 00 is f o u n d to be :
BRo = 3.07250ff°'9415 +2.5076.
(13)
Table lb. Exact roots and half-anglesfor cone v 1 2 3 4 5 6 7 8
BRo 4.493409 5.763459 6,987932 7,58834 8.77148 9.93611 11.08637 12.22509
00 90 54.73561 39.23152 30.55559 25.01733 21.17690 18.35785 16.20078
Diffusion equation in cones and wedges
2.2. Optimum size o f cone
3o
Let us assume that we have a cone of fixed buckling, i.e. fixed net leakage, and require to know how the volume and surface area vary with 0o. The volume of a cone is given by the expression:
25
2nR 3
v = - ~ - - ( 1 - c o s 00).
(14)
Using the empirical relationship (13) for Ro as a function of 0o (remember that B is fixed), we find that the volume is : 2~
'
193
3
V(Oo) = 3 ~ ( ~ 0 f f " + ~ ) (I --cos0o)
% 2o
(15) 15
where ~, fl and 7 are given in (13). It is readily shown that V(Oo) has a minimum at 00 = 0.5259 (30.132°). Similarly, the surface area of the cone which is given by :
S(Oo) = 2~zR2(l --cos 0o) +TzR ~ sin 0o
(16)
also has a minimum at 00 = 0.559 (32°). Both V(Oo) and S(Oo) are shown in Fig. 1 as functions of 0o. Let us now consider a cone of constant volume and examine the buckling as a function of 00. Using the earlier formulae we readily find that : x~2/3 B2(0o) = \ ~ j
(~0J+7)2(l-cos0o)
2'3. (171
I 1o
2oI
I I I 30 4o 5o 8otdegr~sl
I 6o
I 70
1 8o
Fig. 2. The variation of the buckling with cone half-angle, 00, for a fixed volume; the ordinate is in arbitrary units.
2.3. Flux distribution in a cone To illustrate the flux shape in the cone, we consider the following parameters: 0o = ~/6 (30 °) and hence BRo = 8.182561, from Table 1. We also note that v = 4 and thus the flux distribution is from equation (6):
This is illustrated graphically in Fig. 2 and exhibits a minimum at 00 ~- n/6. Clearly there is a critical value of cone angle which minimizes the leakage. The possibility of using such a phenomenon as a reactor control device should be noted.
220 21o 200
I I
-
VtOo)
-----
S(eol
I
"( n "y,,2 ~b(r, 0) = \2B-rrJ Jg/2(Br)P4(cos 0).
The flux is given in Table 2 for 0 = 0 for various values of the ratio r/Ro. We note that a maximum arises at r = 0.685R0. It is also interesting to examine the angular variation of the flux which follows the function P4(cos 0). The values o f P 4 ( c o s 0) for a range of
Table 2. Flux along central axis of cone,
v/Ro
4~(r,O)
0.9t7 0.856 0.733 0.685 0.672 0.611 0.489 0.367 0.244 0.122 0.0611
0.0793 0.1327 0.1968 0.2015 0.2008 0.1870 0.1249 0.05615 0.01408 1.011 x 10 3 6.539 x 10 ~
180 1
170
~ \
/ U \
160
/
~k
" ~ /
I 20
I 30
J
150 I 10
J 40
I 50
I 60
J 70
I 80
I 90
Degrees
Fig. I. The variation of the critical volume and surface area with cone half-angle 00 ; the ordinate is in arbitrary units.
(18)
0.0
0
0.0
E . M . RODRIGUEZ VENTURA et
194
Table 3. Angular flux across cone 0
of radius R0 a n d height H. The zero flux b o u n d a r y conditions are :
P4(cos 0)
o.o
l.O
5
0.9623
20
0.4750
10
0.8532
25
0.2465
30
0.0234
30.55559
al.
0.0
0 are s h o w n in T a b l e 3. T h e flux is n o t precisely zero at 0 = 30 ° due to the a p p r o x i m a t e f o r m u l a for v of e q u a t i o n (10). In fact, it is zero at the smallest root of P 4 ( c o s 0 ) = 0, i.e. 0 = 30.55559 °. The fluxes are illustrated graphically in Fig. 3.
(a)
qS(r, O, z) = O,
(20)
(b)
q~(r, 00, z) = 0,
(21)
(c)
q~(R0, 0, z) = 0,
(22)
(d)
qS(0, 0, z) = 0,
(23)
(e)
~b(r, 0, 0) = 0,
(24)
(f)
0(r, 0, H ) = 0.
(25)
U s i n g separation of variables, we find t h a t the solution can be written :
dp(r,O,z) = CJu,(Brr)sin(plO)sin(B,z),
(26)
P l = X/0o
(27)
where
3. S O L U T I O N OF T H E DIFFUSION EQUATION F O R A
and
BARE, FINITE WEDGE
3.1. Generalsolution andeigenvalues
(H)2
T h e diffusion e q u a t i o n for a wedge of finite height
B2 =
= B2
B 2.
(28)
is:
1 8 [" 8dA
~r~)+~
1 82¢
02,~
a02 +7;~ +B2~=0
B 2 is the radial buckling a n d is f o u n d from c o n d i t i o n
(19)
(23) viz.
J~oo(BrRo) = where ~b = ~b(r, 0, z). W e assume t h a t one face o f the wedge is in the plane 0 = 0 and the other in the plane 0 -- 00. The wedge is
r
0.
(29)
Thus we have a relation between the physical dimensions of the wedge a n d its nuclear properties, B 2. F r o m tables o f zeros of Bessel functions we can o b t a i n Table 4. Again, for convenience in calculation, an empirical relationship between BrRo a n d 00 has been o b t a i n e d thus :
10
BrRo =
4.048800 0.9378 _.~2.50.
(30)
Table 4. Zeros of Besssel function of first kind of order, x/0o .
Fig. 3. Flux distribution inside a cone with 0o = x/6, R0 = 10. The curve on the right hand side of the figure represents the radial flux distribution and that on the left the angular distribution ; actual magnitudes are arbitrary.
rt/Oo
0o
(BrRo) exact
(B, Ro) approx.
1 3.2 2 5/2 3 7/2 4 9/2 5 11/2 6 13/2 7
r~ 2x/3 x/2 2rt/5 x/3 2x/7 x/14 2~/9 re/5 2~z/l 1 ~z/6 2n/13 ~/7
3.83171 4.493409 5.135620 5.763459 6.380160 6.987932 7.58834 8.182561 8.77148 9.355812 9.936110 10.512835 11.08637
3.884 4.524 5.151 5.768 6.377 6.981 7.578 8.171 8.760 9.346 9.928 10.507 11.083
15/2 8
2n/15 r~/8
11.657032 12.22509 17.250455 24.878005
11.657 12.228 17.284 24.933
25/2 39/2
2r¢/25
2rt/39
Diffusion equation in cones and wedges The accuracy of the expression may be judged from Table 4. For 0o < 120 ° it is better than 0.7%. 3.2. Optimum size o f wedge As in the case of the cone, we look for an angle which minimizes the buckling for a given wedge volume. Using equation (30) and the volume of the wedge, we get: H 0 ob +c)-, B~(Oo) = 2~Oo(a
(31)
where a, b and e are defined by equation (30). As Fig. 4 shows, the radial buckling has a m i n i m u m and by differentiation of equation (31) we find: -
-
= 83.15 °.
(32)
3.3. Flux distribution in a wedge For 00 = ~z/4, we find that the flux distribution can be written: rcz) (o(r, O, z) = J4(Brr) sin (40) sin ~ (33) where B_ = 7t/Hand BrRo = 7.58834. The maximum value of q~ along the radial direction arises when: d J4 (Brr) _ O, dr
(34)
which leads to Brrmax = 5.31755, i.e. rmax ~---0.70075R0. Of course, in the axial direction, the maximum is at z = 11/2 and in the 0-direction at 0o/2.
15e
~
-,
100
195
4. NUMERICAL VERIFICATIONOF CODE AND THEORY TO verify the theoretical relationships derived in the previous sections, particular examples of the cone and wedge configurations are solved by means of a multigroup finite element transport code called T R I P A C . This code is a development of two earlier codes, namely F E L I C I T (Wood and Williams, 1984; Wood, 1986) and E V E N T (De Oliveira, 1986, 1987). A feature of T R I P A C , as of its precursors, is its ability to handle complex geometries, and the two configurations considered in this paper provide an interesting test of this aspect of the code's capabilities. The numerical solutions that we give also provide useful benchmarks for proving other techniques and c o m puter codes for which the claim of geometrical flexibility is also advanced. In this paper we are concerned with diffusion theory, i.e. P,/P~ solutions, but of course T R I P A C is also capable of providing higherorder transport solutions, namely, P,/P/. Two cases of the cone are considered, and two of the wedge. Four-group linear anisotropic scattering data is used for the homogeneous systems. The data is based on the simple H T G R core data used by White and Frank (1987). Our data is listed in Appendix A ; the corresponding value of materials buckling is Bm2ax= 1.9701604 × 10 2. In the computations, iterations are continued until the eigenvalue (kc~) and eigenvector have converged, respectively, to 10- 6 and 10- 5. Also the eigenvector is normalized to unit fission neutrons produced in the system, i.e. 1 ['~' k ~ J J d E d ~ ' ~ Z j ~b(E, r) = 1. In the case of the wedge, the eigenvalue is actually normalized to 1/'4 of the total volume--this is explained below. It is perhaps worth mentioning that this particular data has the fundamental and next lower eigenvalue close together (White and Frank, 1987). Thus, unless a power iteration convergence acceleration scheme is employed in the code, a large number of iterations may be necessary to produce the accuracy we quote in T R I P A C , outer iteration acceleration is achieved by Chebyshev extrapolation. 4.1. The cone
50 I t I I I t L I I I I t I I I J ~o 40 6o 80 leo lzo 14o 16o 8 0 (degrees)
Fig. 4. Variation of radial buckling B,2 with wedge angle 0o for constant volume; the ordinate is in arbitrary units,
The cone is solved in R - Z geometry using 81 quadratic elements. The element mesh employed is illustrated in Fig. 5. The exact and T R I P A C results are compared in Tables 5 and 6. Clearly, the agreement is excellent. The flux distribution along the z-axis of the smaller cone is shown in Fig. 6. Figure 7 shows the variation of critical volume of a cone for fixed buckling as a function of the cone
196
E . M . RODRIGUEZ VENTURAet al.
Z!.
13.7.87
Cone 4-Grp
.~'-;~-.-.~'.i
R-Z Diff. AppCone 21 =
" ---
x
~ -'2
¢qoo
. /
f::
:: '.?
i" i 'Z/'
:."f)
'"
"'}:"
~ c~ ,::; :
, ii"
•
!
x ,~.
rq ~ o¢~
::
8
..?'
::= ~
~_
~
~- o~
emoo
i./.
i] ,:::
i - - . . . - : ; -.- - •.. !!:
"! e! o o
i.
/11 ii?-:; !,::::.:.:::ii:::::"::..... '
X i..:::" ii ,:-
~?
~
•
x
o o ~
x
~
~,
,~o,~ e~ r~
?,
"].._........i:
I '
R
0
Progro,,'n P o s t - - p r o c e s s ~ F i g . 5. F l u x
contours Contour
Mo~'ch ~9 , 9 7
for group increment
4 neutrons i n 1 x 10
in
21 ° c o n e .
4.
~rq
half-angle. This is a four-group calculation using the data in Appendix A ; it confirms the analytical result in Section 2.2 and indicates an optimum half-angle of about 30 °.
~ ~
~ '~ ~. "~'t~
4.2. The wedge This problem is solved in X - Y - Z g e o m e t r y . T h e basic element in the X - Y plane is an arbitrarily ori-
(~
~ .~
Diffusion equation in cones and wedges
197
Table 6. Comparison of critical flux distributions in 54 ° cone, group 4 neutrons
r
z
TRIPAC
0.0 0.0 0.0 0.0 0.0 0.0 0.0 95.777 126.706 152.234
390.082 307.960 287.429 246.368 184.776 123.184 61.592 248.827 160.352 336.739
x
105
EXACT
0.1324 0.6317 0.7155 0.7991 0.7031 0.4173 0.1234 0.6213 0.3215 0.1960
105
x
~.
0.1294 0.6271 0.7116 0.7970 0.7020 0.4173 0.1231 0.6208 0.3215 0.1969
~_
Note: cone is solved by TRIPAC in R Z geometry.
ented triangle. Again, a quadratic triangle is used. The subdivision of the X-Yplane is shown in Fig. 8. These triangles form the cross sections of prisms whose axes are parallel to the z-axis. Each prism so formed is further subdivided into 3 tetrahedra. Thus the basic element in space is, in this case, a quadratic tetrahedron (with 10 nodes). Because of symmetry, the problem need only be solved in a 1/4 wedge by taking the origin to be at the mid-height of the wedge and exploiting symmetry boundary conditions along appropriate axes. As with the cone, the boundary condition on the surface of the wedge is zero flux. The exact and TRIPAC results are compared in Tables 7 and 8. In TRIPAC, 735 quadratic tetrahedra are used for each wedge. An isometric view of the flux in a critical wedge is shown in Fig. 9. For the wedge, as for the cone, the agreement between theory and TRIPAC results is excellent. Appendix B indicates how the volume-averaged fluxes are calculated for cones and wedges.
0.9
,,Dt~ ~
tt'5
,.~
~Le~ m m:~ .--: ~ ~
~ × ~=
~ ~ eqtt5
~ ~"
~5
-
~2
~;
~s
=
~
~
"~ "~ .~ _~ .~
~, ~ ~~
"~
b-
13" 7 " 8 7 Cone 4 - G r p R - Z Diff. App Cone 21 ° Tota I. flux or dose rate
0.8 0.7 0.6
o ×
~=
05 0.4
E
0.3
0.2 0.1
''
i 0
I00
200
300
400
500
600
700
I 800
rj
Z
Program Post-processor March 1 9 8 7
Fig. 6. Variation of flux along z-axis for 21': cone (group 4 neutrons).
~
~ x~
198
E . M . RODRIGUEZ VENTURA et al. Critical volume versus cone ongte 7.0
20'7.87 Wedge 4-Grp X - Y - Z
6.5
~~
~
Diff. Approx.90"
6.0
5,5 I 40
5.(
I 20
I 30
I 40
I 50
i
60
I 70
I 80
YI/~X/~//~//~
I 90
'
/
\
/
'
/
'
/
'
0
Fig. 7. Variation of critical volume of cone with half-angle 00 for fixed buckling using four-group data.
IRefteetion
X
GEM Vr 1.9.86 Fig. 8. Subdivision of X Y plane into quadratic triangles for 90 ° wedge. Because of symmetry, only 1/2 of X- Y plane need be solved and 1/4 of volume of wedge.
5. CONCLUSIONS AND SUMMARY
Some exact solutions of the diffusion equation in cones and wedges have been obtained leading to critical conditions and fundamental mode flux distributions. We have observed that the volume and surface area of cones and wedges go through a m i n i m u m at a certain value of cone or wedge angle for a fixed buckling..Similarly, the buckling has a m i n i m u m for
14.7-87 Wedge 4 - G r p X - Y - Z Scalar flux group 4
fixed volume at a critical angle. The latter condition implies a m i n i m u m leakage condition and has implications for reactor control and storage of fissile material. We have also used these exact solutions to assess the accuracy of the T R I P A C finite element code. A
Diff.Approx.45 °
2.75
[- 2.75
225~
FV 225
1.75~
7 1.75
1.25 ~
!
0.75 --
['-
0.25
~- 0.25 0.0
D
1.0
2.25
2.0
1.75 3.0
1.25 4.0
102 Z axis x 102 Vertical axis x 10 -5 R
aXiS
X
0,75
5.0
0.25
6.0
R
Program Post- processor March 1987
Fig, 9. Isometric view of flux in 45' wedge, for the plane Y = 0 (group 4 neutrons).
1.25
Diffusion equation in cones and wedges
199
Table 8. C o m p a r i s o n o f critical flux distributions in 90 ° wedge, g r o u p 4 neutrons x
y
85.714 171.429 257.143 342.857 428.571 514.286 190.536 414.445 355.528 300.917
z
0.0 0.0 0.0 0.0 0.0 0.0 97.0872 99.740 308.062 160.844
T R I P A C x 105
E X A C T x 105
0.3324 1.1782 2.0833 2.5521 2.2798 1.3032 0.6853 1.1970 0.2642 0.9993
0.3370 1.1737 2.0712 2.5390 2.2733 1.3072 0.6826 1.1958 0.2666 0.9964
0.0 0.0 0.0 0.0 0.0 0.0 70.6011 84.7213 0.0 70.6011
N o t e : wedge is solved by T R I P A C in X - Y - Z geometry.
f o u r - g r o u p b a r e r e a c t o r criticality p r o b l e m is solved for a variety o f w e d g e s a n d cones. B o t h eigenvalues
0.28248588E+00 0.121900E-02 0.132000E-02 O.OOO00000E+OO0.00000000E+00 0.25816688E+00 0.23100000E-01
a n d critical fluxes are s h o w n to be in excellent agreem e n t w i t h the analytical results. C o n e s a n d w e d g e s are a sensitive test o f a n y n u m e r i c a l p r o c e d u r e a n d lend f u r t h e r c o n f i d e n c e in t h e a c c u r a c y o f t h e T R I -
0.28985507E+OO 0.a54000E 02 O.OOOOOOO0E+OO0.00000000E+00 0.0 0.0 Normalizedfissionspectrum 0.9675
0.0
0.0
0.025816688E +00 0.0
0.67S000E-02 0.00000000E+00 0.28531507E+00 0.0 0.028531507E+00 0.0325 0.0 0.0
P A C code.
APPENDIX B Average Fluxes
REFERENCES Abramowitz M. and Stegun I. A. (1965)Handbook of Mathematical Functions. Dover, New York. Glasstone S. and Edlund M. (1953) Nuclear Reactor Theory, Van Nostrand, New York. Morse P. M. and Feshbach H. (1953) Methods of Theoretical Physics. McGraw-Hill, New York. White J. R. and Frank B. R. (1987) International Topical
In order to obtain average fluxes it is necessary to calculate :
1 (" ~ = V Jv dr q~(r). For the cone with :
Meeting on Advances in Reactor Physics, Mathematics and Computation, Vol. 3, p. 1385. American Nuclear Society. Wood J. and Williams M. M. R. (1984) Prog. Nucl. Energy 14, 2 I. Wood J. (1986) Ann. nucl. Energy 13, 91. De Oliveira C. R. E. (1986) Ann. nucL Energy 13, 227. De Oliveira C. R. E. (1987) Ph.D. Thesis, University of London.
1 (p(r, O) = ~
Jv+ I/2(B,)Pr( cOs 0),
x/r this expression becomes: 3 /'R q~ - R3(I - c o s 00) J0 drr3/2J~+j/2(B,) x ('1°° dO sin 0 P~ cos 0). ,)0
APPENDIX A E
]E~
E~o(E, ~ El)
E,o(E, --* E2)
For integer values of v, it is easy to carry out the 0 integration and although the integration over the Bessel function can be represented by an infinite sum, it is more convenient to evaluate this directly by quadratures. For wedges we have:
E~, (E~ ~ El)
Y's~(Ei --* E2)
~p(r, O,z) = J~,, (Brr) sin (p. 0) sin (Bzz),
For each energy group i: v-'Zl
Note: diffusion coefficient D = l/3Z,r where E,, = Z - E s l . whence Four-group data
8
0.14184397E+00 0.844000E-04 0.12945957E+00 0.12300000E 01 0.012945957E +00 0.0
0.784000E 04 0.00000000E+00 0.0
0.00000000E+00 0.0
0,27777778E + 00 0.145100E-02 0.00000000E+00 0.27255678E+00 0.0 0.027255678E+00
0.488000E- 03 0.37700000E-02 0.0
0.00000000E+00 0.0
repro
(~--7~2(/~-D0)2[~r'" Jo |
dxxJ~(x),
where t/i = 71"/0o.