- Email: [email protected]
Synthesis of conservative discontinuous solutions for diffusion theory
Ann. Nucl. Energy Vol. 23, No. 17, pp. 1381-1399, 1996 Copyright 0 1996 Elsevier science Ltd Printed in Great Britain. All rights reserved PII: s0306-4549(96)ooo17-5 03064549196 $15.00+0.00
Pergamon
SYNTHESIS
OF CONSERVATIVE DISCONTINUOUS DIFFUSION THEORY
R.T. Ackroyd, A. M. Gashut,
0.
SOLUTIONS
FOR
A. Abuzid
Mechanical Engineering Department, Imperial College of Science, Technology and Medicine, Exhibition Road, London, SW7 2BX, U.K. (e-mail: [email protected])
(Received 14 April 1995) ABSTRACT-A finite difference scheme for diffusion theory is presented which is based on the finite element method, and which is as efficient and accurate as the nodal method using a quadratic approximation in the local coordinates for the flux. A series of test problems in one and two dimensional geometries show that the method produces accurate solutions without the need for the several iterations of Copyright 0 1996 Elsevier Science Ltd the nodal method.
1. INTRODUCTION This paper deals with the development of a single and double quadratic non-iterative nodal method based on a discontinuous finite element solution obtained by the K’ variational principle of Ackroyd (1986a). To start the solution, a piecewise constant discontinuous finite element solution is found. For brevity these elements are called primitive elements. The K’ principle with discontinuous finite elements aims at cutting out as far as possible iterations. Usually two applications of the K’ principle are sufficient to improve upon the primitive element solution which is in effect the solution of the well known central difference equation of the diffusion theory. The following steps are considered : 1. Determine a piecewise constant discontinuous primitive elements with the K’ principle.
finite
2. Use the variational solution profiles say, 41 and Br.
of 4 to form two different
3. Combine
solutions
the extrapolated
4. Optimize 5. Repeat
(~2, varying
4s by maximizing if desired
to obtain
extrapolated
solution
(4) by the use of
solutions
with quadratic
in the form
4s = dl with the coefficient
element
+ (1 -
from element
a)
01
to element.
K’ (&) with respect to oz. +,,+I, (say). 1381
(1)
1382
R. T. Ackroyd et al.
1. MATHEMATICAL
FORMULATION
The I(* variational principle The diffusion equation is
(2) The K* (4) principle states that for a volume V partitioned by a surface Sr,
I(* (4)
=
/ {2&S’- C42 - D (V+)2} dV ”
-
J w[+(cto)-$(r-o)12dS Sl
-
I
F($o&o) t /,, (‘*rdJ2dS
(3)
where
F(+oo,40)
=
/ {D(V+)2 t C+2}dV ”
and Sb denoting a bare surface. The penalty parameter w is positive. It can be chosen so that for rectangular mesh and primitive elements the variational solution is the same as the central difference equation for diffusion. Formulae for this w are given in section (4.2). The (p •l-1)th trial function &,+rl(z, y) is given in terms of &(z, y) and an associated function 8,(z, y). For example over the element V”” of Figure 1, with local coordinate (t, q), A+1 (L 7) = o::iA
(L 11)+ (1 - or-r) 0, (& 77)
(5)
where &, at p=l is the flux obtained from the discontinuous finite element, 0, is its associated function. The coefficients crp”+;,(p = 1, . . . . .P - 1) are determined by maximising K* (&,+I) with respect to the coefficients ap”+; for all m and n of the mesh. The flux +,, (t, YJ)is defined over the element V”” for -1 5 t 5 1, -1 < q 5 1, in terms of the Legendre polynomials of e and q
4p(t, 8) = 4”” (P)PO(0 fi (rl)
t
$P;“”(P>- 4Y WI
Pl(O PO(17)
t ; WY (P)- 44mn (PII PO(0 Pl (VI
where :
+
f [4Y
(?+-I-~;"
t
; [4;" (PI t4T"
(P> - 24""
(P,
b)lPz
- W”” bwm
(0
w?)
p2(rl)
(f-9
1383
Diffusion theory
Figure 1: Typical element (m,n) for non-iterative anticlockwise (1,2,3,4). t, q are as defined before
7*..*P, 4mn (PI 4T” (P)
nodal; method.
The surfaces are numbered
p=l
is the average of &, over V”“. is the average of +P taken over the face i,(i=1,2,3,4) of the element.
The associated function flP(E, q) for V”” is expressed in the following from; 0, (6 V)
=
0”” (P) PO(0 PO(rl)
+
f PY
(P) - OF
(P)l Pl
[@Y
(P>
-
or
(P)J
(0 PO
PO (0
(77)
Pl
(7))
+
-5
+
;
[@T-”
(P)
+
@Y
(PI
-
2@“”
(PN
p2
(0
PO
(7)
+
f
[W
(P)
+
Or;”
(P)
-
z3””
(P)]
PO (0
p2
(71)
(7)
where : 0”‘” (p) is the average of 6, over V”” Or” (p) is the average of 8, taken over the face i, (i=1,2,3,4)
of the V”“.
The neutron conservation condition for 8, (t, q) over the element V”“, /_: /_:
[-Dm”v2e, + Yv,
-
sq =0
is used to express 0”” (p) in t erms of the 07” (p), which are determined as follows;
(8)
1384
R. T. Ackroydet
Construction
al.
of q51(t, 11) and 81 (E, 7)
(a) The double quadratic scheme To start the calculations the functions &([, q) and 81 (I, 7) are constructed from the results of the simplest possible finite element calculations based on the K* (4) principle, i.e the calculations based on the primitive elements. For primitive elements the trial function 4 = amn, where the am” are constants. Construction of $~r(t, q) from the am” dl (6, II) is defined by equation (6) with the #” (1) defined by the average interface values
1 [am+l*n +
amn]
C#J~(1) = i [am-l*n +
amn] amn]
4;” (1) = @” (1) =
1[am*n+l + amn]
4”” (1) = f [am~n-l +
(9)
The average element flux 4”‘” (1) is specified by the conservation condition (10) Construction of 8i([, q) from amn 0 ([,q) is defined by equation (7) with the ey (1) are the average values of 8i (E, 11) over the surfaces (i = 1,2,3,4) and 0”” (1) is the average of 81 ([,v) over Pn. These averages are determined in terms of the following procedures; Fit a quadratic to the amn at < = 0, the mid-plane of V”“; amtlln at t = 2 mid-plane of Vm+l+; and am-l+ at t = -2 the mid-plane of Vm-l*n. i.e. if(t)
= [am-l*nt (t - 2) - 2amn (< - 2) (< t 2) t amtltnt (e + 2)]
(11)
Set; @y (1) = f(1) = L [3amt1+ + 6a*n @y” (1) = f(l) = i [3a*-l+ + 6a*n
_ a*-ltn] _
a*+l+]
(12)
Similarly @r @r
(I) = f(I) (1) = f(I)
= 1[3a”‘*“+r + 6am” _ am*n-i] = i [3am.n-r + 6e”‘n _ e?n+r]
(13)
Also 0”” (1) is got from the conservation condition;
J-l/_: [-D”“v2elt zmnel- smn]d&iv =o
(14)
(6) The single quadratic scheme In this scheme the flux $r([, II) is taken to be a constant over an element , i.e for V”” &(<,q)
=4”“(l)
= f#~r(l) =amn(i
= 1,..,4)
The construction of the component function 8, (c, 9) is based on dP (t, 11)in the following way : Let ST”(p) be th e average value of &p(<,q) (t, q) along DB of Figure 1.
(15)
Diffusion theory
1385
A quadratic P(E) is fitted in Figure 1 to the average values ,7-l*” the mid-element lines GH,DB and EF respectively, p (0 = ; [QY-l,n (P) E (E - 2) - 2T”
(p), \k‘y” (p) @y+‘” (p) for
(t - 2) (< + 3) + Q:+‘*” (P) E (t + 2,]
(16)
The functions Oy” (p) and Oyn (p) are specified as @Y (P) = f(l) W(P) = U-1)
(17)
1
i.e o;ln
(p) = i 3q+‘+
@Tn (p) = ;
t
_q+l.n
(P) + f3’QY (P) - Q?”
(P)]
(P) + 6’K’” (P) t 39:-l’”
(18)
(P)] I
Similarly by defining STn (p) as the average value of &, ([, 0) along CA in Figure 1, i.e
(19) The above fitting procedures gives @yn (p) = i [307’“+’ (p) + SST 0:”
(P) = 5 [-%““+1
(p) - @;Vn-l (p)]
(p) t SKI’? (p) + 3!P;*“-’ (p)]
(20)
Also 0”” (p) satisfies the conservation condition for V””
Thus all the terms of eP ([, q) are known. Hence it can be determined using expression (7). Element Conservation for the Double Quadratic Scheme In the sequel 6,, p 2 2 is constructed for the double quadratic so that it satisfies neutron conservation condition for every element. In particular for V”“. l$s=cz~~t[1-cX~]&
(22)
satisfies the neutron conservation condition because 42 and 02 both satisfy the condition. Continuing the procedures show that &+I, p > _ 1 satisfies the neutron conservation condition, Ackroyd (1991). Determination
of ~$2(I, q)
The trial function 42 (5,~) is a linear combination for each element of +l(t,q)
and e1 (E, 71)
i.e 42 (L 77) = 4” 41 (<,‘I)
-I- (1
-
e)
b(k
7)
(23)
Wliere the coefficients a?” for the element V”” are found by maximising K’ (&) with respect to cy?n. Since $1 is an approximate solution of the diffusion equation for the impressed source,S, and 01 is another approximate solution for the same source S, then & as defined above can be
1386
R. T. Ackroyd et al.
regarded as an approximate solution for the given source. Note that 4s satisfies the conservation condition for V”” + Cmn4& - S”“‘] d(dq = 0
I_: J_: [-D”“V242 because its component for V”“.
functions
3. DOUBLE Having as
found
41 and $2 satisfy
QUADRATIC
by construction
PROCEDURES
& using steps (1) and (2) on calculates
(24)
the conserva.tion
FOR
condition
p 2 2
for p 2 2 &+I (t, q) which is defined
: htdb?)
with the function function Let @y”
=
$3A4777)
+
[I - $$$ad
(25)
coefficients c$!!i to be found by maximizing IC with respect to orti. The component r$p (t, ~7) is k nown from previous calculations. The construction of the component & (I, q) is based on &, (t, 77) in the following way : (p) be th e average value of 4(p) (E, v) along DB of Figure 1.
‘K’YP)
=
;J_:9,(Wdrl
=
; [6Y
(P) - &F’ (~11
(P) - d?
A quadratic P(t) is fitted in Figure 1 to the average values ST-‘” the mid-element lines GH,DB and EF respectively, P (E) = ; [*:-r7n The functions
(p), ar”
(26) (p) @y+‘” (p) for
(p) t ([ - 2) - 297” (E - 2) (E + 3) + rlY+“” (p) t (5 + 2)]
Oy” (p) and Oyn (i) are specified
(27)
as
@Y(P) = f(l) wn (P) = f (-1)
(28)
i.e @yn (p) = ; $$y”” @y Similarly
by defining
(p) = ;
t
(P) •t @J’~” (P) - %“-lTn (P)] (P) + @T
~IJ?” (p) as the average
ep (P) . The above fitting
_q+‘*”
procedures @r
(29)
(P) + 3912-“” (P)]
value of &, (t, 0) along CA in Figure
1, i.e
= =
; [Wmn (P) - $7” (P) - d4mnb)]
(30)
gives
(p) = f [3’l?;‘“+’ (p) + SST” (p) - ,y”-’
@Tn (P) = i [-Q?“+
(p) + 6@yn (p) + 397”-’
(p)] (p)]
(31)
Diffusion theory
Also 0”” (p) satisfies
the conservation
condition
1387
for V””
using expression (7). In Thus all the terms of 0 (t, 7) are known. Hence it can be determined so that it satisfies neutron conservation condition for every the sequel B,, p 1 2 is constructed element. In particular for V”“.
satisfies the neutron conservation condition because 4s and 0s both satisfy the condition. Continuing the procedures show that 4 r+r, p > 1 satisfies the neutron conservation condition.
4. NUMERICAL
EXAMPLES
In order to investigate the usefulness of the method discussed above, a computer code, DNMFD, (Direct Nodal Method For Diffusion) was developed in 1-D (slab) and 2-D x-y geometry, three sample problems in one and two dimensions were used to test the method for accuracy of solution. One-dimensional
problems
Two problems were presented. The first is a test problem consists of two regions. It is used to check the accuracy of the non- iterative nodal method solution against the analytical solution. The geometry of the problem is given in Figure 2 and Table 1 presents the material properties. The boundary conditions used are perfect reflection at x=0.0 cm and at x=2.0 cm. These boundary conditions are no more than symmetry conditions. Consequently the tests are only of the capabilities of the schemes for solving the differential equation for diffusion. The flux distributions obtained with the discontinuous finite element using primitive element, double quadratic non-iterative nodal method and the classical (iterative) nodal method for 20 nodes, together with the analytical solution are presented in Figure 3 for comparison. Twelve iterations were used for the nodal solution. The percentage error in these solutions, as calculated in reference to the analytical solution is presented in Figures 4, 5, 6 and 7. A maximum error of 0.16% in the primitive finite element solution and 0.12% in the single quadratic non-iterative nodal solution is seen Figure 4. The single quadratic solution shows marginally better agreement with the exact solution than does the primitive finite element solution. The error in the noniterative nodal solution reduce to 0.045% when the double quadratic scheme is employed for the same number of nodes as seen in Figure 5. This show that the non-iterative nodal results are significantly more accurate than the primitive finite element results relative to the analytical solution. To assess the double quadratic non-iterative nodal solution against the iterative nodal solution, Figures 6 and 7 present the error distributions in both solutions for 20 nodes. A maximum error of 0.4212% in the iterative nodal solution with 6 iterations and 0.045% in the double quadratic non-iterative nodal solution is seen in Figure 6. The iterative nodal error reduced to 0.046% when 12 iterations are involved as seen in Figure 7. This show that non-iterative nodal method. gives results as good as the iterative nodal method but much faster. The second problem is a slab lattice cell consisting of one mean free path thickness of natural uranium as fuel and three mean free path of graphite as a moderator with a unit source. The problem has been solved by Nanneh(1990) for transport calculations. The geometry of the
R. T. Ackroydet
1388
al.
1.0
0.0
b
2.0
Wcm)
Figure 2: Geometry of the test problem
Table 1: Group cross sections and source density for Problem 1 Cross Setions
Region 1
0 Source Density 1
1.0
Region 2
;I 1
0.0
0
: . Iterative 1.25- 8 Primitive element _ 0 Non-iterative
l.O-
Figure 3: Analytical, primitive element, double quadratic solutions for Problem 1
non-iterative
and iterative
nodal
1389
Diffusion theory
9 %
on
0.5
0.7S
1.0
,.I
1.5
1.n
Di5mcehllhox-(m)
Figure 4: Error in central difference for 20 nodes for Problem 1
(K’)
and single quadratic
non- iterative
nodal
solutions
nodal
solutions
Primitive element
0:s
on
IR
125
I.5
1.78
Dk?tMoehlheXdlrecabn(an)
Figure 5: Error in central difference for 20 nodes for Problem 1
(I(*) and double
quadratic
non- iterative
1390
R. T.Ackroyd et al.
Dls?ancehIhsXdiredlW,(an)
Figure 6: Error in double quadratic non- iterative and iterative nodal solutions for 20 nodes nodes after 6 iterations for Problem 1
Figure 7: Error in double quadratic non- iterative and iterative nodal solutions for 20 nodes after 12 iterations for Problem 1
1391
Diffusion theory
Figure
8: Geometry
of the repeating
lattice
cell problem
problem is given by Figure 8 while the material properties and the source density in each region are presented in Table 2. This problem is also a test of the capability of DNMFD calculations to give accurate results as compared to the finite element continuous solution presented by a diffusion theory solution given by the K+ principle encoded on EVENT, de Oliveira (1994). The solution obtained from DNMFD, primitive finte element method and iterative nodal method are presented in Figure 9, together with the continuous finite element solution for comparison. Figure 10 presents the error distribution of the primitive finite element and the double quadratic non-iterative nodal solutions relative to the continuous finite element solution. A localised maximum error of -1.3% in the primitive finite element solution and -0.91% in the double quadratic non-iterative nodal solution with 3 nodes in the fuel region and 10 nodes in the moderator region can be seen in the figure. To check the double quadratic non-iterative nodal solution against the iterative nodal solution, Figures 11 and 12 show the percentage error in both solutions relative to the continuous finite element solution based on K+ solution. A maximum error of -0.72% in the iterative nodal solution after 6 iterations and -0191% in the double quadratic non- iterative nodal solution is seen in Figure 11. The error in the iterative nodal solution reduced to 0.5% when 10 iterations are involved. The non-iterative nodal method shows once again that it is a promising method in comparison to the continuous finite element method and that it overall provides a solution nearly as accurate as the iterative nodal method solution but much faster. To assess the single quadratic non-iterative nodal method in obtaining accurate solutions for lattice problems in comparison to the contiimous finite element solution, Figure 13 presents the error distribution in the primitive finite element solution and the single quadratic non-iterative nodal solutions relative to the continuous finite element solution (K+). As it is seen from the figure, the single quadratic non-iterative nodal method gives solution of comparable accuracy This primitive solution is the same to the use of the K* principle with primitive elements. as that given by the classical difference equation for diffusion. In the sequel only the double quadratic scheme is considered as a non-iterative nodal scheme, as it is superior in accuracy to the single quadratic scheme, and it does not involve extra computation.
R. T. Ackroyd et al.
1392
Table 2: Group
cross sections
Cross Setions
and source density
for Problem
Region
2
1
Region
2
111 1 Source Density
1
1
1.0
0.0
fl
Iterative Primitive element Non-iterative Finite element (K+)
P
IO
1,
12
la
14
olstance in the x ckaion (cm)
Figure
9: Neutron
flux profile for the one group, one dimensional
lattice
cell problem
Diffusiontheory
1393
4.!26-
Primitive element Non-iterative
Q.S-
B
Ii
is
i 47S-
.o-
-1
-l.zs-
1 P
10
11
12
II
14
13
19
1,
OiaanmhtheXd~n(an) Figure 10: Error in primitive finite element and double quadratic non-iterative solutions for 13 nodes for Problem 2
nodal method
0 Non-iterative
Figure 11: Error in double quadratic non-iterative for 13 nodes after 6 iterations for Problem 2
nodal and iterative nodal methods solutions
1394
R. T. Ackroyd et
al.
Iterative 0 Non-iterative
0
Dismnca IntheX dlreabn(an) Figure 12: Error in double quadratic non-iterative nodal and iterative nodal methods solutions for 13 nodes after 10 iterations for Problem 2
LI Primitive element 0 Non-iterative
Figure 13: Error in primitive finite element and single quadratic non-iterative solutions for 13 nodes for Problem 2
nodal method
1395
Diffusion theory Two-dimensional
The penalty
Problems
parameter
w discussed
in section
(2) is given by
(34)
One case study is considered in 2-D (x-y) geometry. The problem is used by Fletcher for transport test calculations. The geometry and the neutronic parameters are given if Figure 14 and Table 3 respectively. It has been selected to check the validity of the non- iterative nodal method in obtaining solution of comparable accuracy to the continuous finite element solution As it can be seen based on K+ principle, de Oliveira (1986,1994), and using quadratic elements. from the figure, the problem is a square source problem with a constant source in a (1.2c97~)~ conditions for the sides adjacent to the source region are perfect reflector, region. Boundary while that for the other sides are bare surface conditions. This problem has been solved with the x and y pitches being 0.2 cm. The double quadratic non-iterative nodal, the primitive finite element , the K+ and the iterative nodal solutions are presented in Figure 15 at y=O.l cm and Table 4 gives these results for comparison. From both the table and the figure, it is clear that all the methods provide solutions near to the continuous finite element solution based on the K+ principle. Figures 16, 17 and 18 give the error distributions in these solutions relative to K+. A maximumerror of 3.1% in the primitive finite element solution and 2.4% in the double quadratic non-iterative nodal solution can be seen in Figure 16. Comparison between the errors in the double quadratic non- iterative nodal method and the iterative nodal method (see Figures 17 and 18) show that the non-iterative nodal method gives solutions overall as accurate as the iterative nodal but much faster.
6.0
“IOU I.2
Figure
rczz
14: Geometry
6.0
of Fletcher
X(m)
problem
1396
R. T. Ackroydet al. Table 3: Cross sections and source data for Fletcher’s problem Cross Section cm-’
Region 1
Region 2
Et E,
1.0 0.0
1.0 0.0
Source Density
0.69444
0.0
Table 4: Continuous finite element (I(+), iterative nodal solutions at y=O.l cm
non-iterative nodal, primitive finite element and
Distance along x axis (cm)
Continuous Finite Element Method (Pi)
Iterative Nodal Method (18 Iterations)
Primitive Finite Element Method
Non-Iterative Nodal Method
0.1000E+00 0.3000E+OO 0.5000E+OO 0.7000E+OO 0.9000E+OO O.llOOE+Ol O.l300E+Ol O.l500E=Ol O.l700E+Ol O.l900E+Ol 0.2100E+Ol 0.2300E+Ol 0.2500E+Ol 0.2700E+Ol 0.2900E+Ol 0.3100E+Ol 0.3300E+Ol 0.3500E+Ol 0.3700EtOl 0.3900E+Ol
0.541233+00 0.531813+00 0.511593+00 0.477623+00 0.425043+00 0.346493+00 0.252533+00 0.173023+00 0.118383+00 0.809543-01 0.553703-01 0.37896EOl 0.25961EOl 0.178033-01 0.12223EOl 0.840123-02 0.578073-02 0.398193-02 0.274553-02 0.189493-02
0.5427443+00 0.5334193+00 0.5133983+00 0.4797673+00 0.4277093+00 0.3499363+00 0.2500263+00 0.1713033+00 0.1172073+00 0.801531EOl 0.5482273-01 0.3752153-01 0.2570403-01 0.1762733-01 0.121021EO1 0.831813E02 0:572359E02 0.39425OEO2 0.271839E02 0.1876153-02
0.542363+00 0.533243+00 0.513683+00 0.480923+00 0.430353+00 0.354993+00 0.244713+00 0.168393+00 0.115773+00 0.795673-01 0.546983-01 0.376223-01 0.258973-01 0.178423-01 0.123043-01 0.849343-02 0.586863-02 0.405893-02 0.280983-02 0.194683-02
0.5423933+00 0.5331763+00 0.5133733+00 0.4801603+00 0.4289093+00 0.3525843+00 0.2464913+00 0.1693983+00 O.l16300E+OO 0.798152EOl 0.547815EOl 0.3761513-01 0.258436EOl 0.1776883-01 0.122263EOl 0.8419043-02 0.5801593-02 0.400061EO2 0.2760393-02 0.1905653-02
Note that the primitive, non-iterative and iterative nodal solutions have a jump in the error at the interface of the source and source-free regions. Both the non-iterative and the iterative nodal solutions recover from the error jump. 5. CONCLUSIONS Based on the results of the numerical examples presented in the previous section, we can make a number of general inferences about the potential of the non-iterative nodal method
1397
Diffusion theory 0.6
0.5
-
NOtI- h-auw !soluuoll Duf. sohdon
cmtra
wuuan conunwur (Pl) !soluuon
1 nwahw 0 0.4
0.3
0.2
0.7
1.5
Distance
2.0
2.5
3.0
3.5
in the X direction (cm)
15: FEMFD, DNMFD and INMFD fluxes in comparison to the K+ finite element solution at Y=O.l cm
Figure
0.6
1 .o
1.5
2.0
2.5
3.0
3.5
Distance in the X direction (cm)
Figure 16: Primitive element and double quadratic non-iterative nodal method error in reference to K+ solution at Y=O.l cm
1398
R. T. Ackroyd
et al.
Figure 17: Iterative nodal and double quadratic non-iterative nodal methods error in reference to K+ solution after 2 iterations at Y=O.l cm
-2o!
,,,,,,,,,,,,I
,,,,I
0.5
1.0
1.5
,,,,,,,,1,,1,,,,,,,
2.0
2.5
3.0
5.5
Distance in the X direction (cm)
Figure 18: Iterative nodal and double quadratic non-iterative nodal methods error in reference to K+ solution after 18 iterations at Y=O.l cm
1399
Diffusion theory
as new numerical technique. The results show fairly good agreement with the analytical and results. Due to the iterative nature of the classical nodal the continuous finite element (K+) expansion method, it involves more calculations in converging towards the true solution than do the direct approach of non-iterative, double quadratic nodal method. The non-iterative nodal method shows promise as a simple method since its use does not require the specification of convergence criteria and the initial guess of the classical nodal method. It starts with the central difference approximation given by the use of primitive elements in the K* principle and replaces that solution with a quadratic profile for each element and with neutron conservation holding for each element. Iterations can be implemented as indicated in section 3 for the double quadratic scheme if required. In developing a finite element’method (Non-iterative Nodal Method) which is akin to conventional nodal methods two new techniques have been implemented with discontinuous finite element method : i. Ensuring neutron balance the flux in each element.
element
by element
by the addition
ii. Constructing a polynomial profile for the flux over an element, in a flux calculation only one unknown per element.
of constant
components
for
by using at any one step
method for neutron The technique (i) is also used by Umpleby (1994) in his discontinuous component of the angular flux is transport in a generalized way. Each spherical harmonic made to satisfy neutron conservation element by element. The direct nodal method and his differ in the treatment of discontinuity. In the direct nodal the trial function is completely discontinuous between elements. In his scheme the trial function have both continuous and discontniuous components. The discontinuities are introduced in the non-iterative method by the choice of a penalty parameter which ensures with the simplest of elements, primitive elements, the central difference equation for diffusion theory. The interface term in the K* principle, as weighted by the penalty In Umpleby’s scheme the element parameter, serves to couple the elements into a system. coupling is achieved via the continuous component of the trial function. The techniques (i) and (“) 11 are also applicable to transport theory, and their investigation for transport calculations is one of the topics under study by Mirza (1994). The use of primitive elements (constant trial function) reduces the number of unknowns, since there is only one node per element. However, higher order trial functions can be used with the I(’ principle. REFERENCES 0. A. (1994) D iscontinuous Finite Element Solutions for Neutron Diffusion PhD Thesis, University of London. Ackroyd, R.T. (1986a) Prog. Nucl. Energy, 18,45. Ackroyd, R.T. (1986b) Prog. Nucl. Energy, 18, 7. Ackroyd, R.T. (1992) Ann. Nucl. Energy, 19, 565. Ackroyd, R.T., Abuzid, 0. A. and Gashut, A. M. (1995) Ann. Nucl. Energy, in press. de Oliveira, C.R.E. (1986) Ann. Nucl. Energy, 13, 663. de Oliveira, C.R.E. (1994) Private communication. Nanneh, M. M. (1990) A Synthesis Method Bzsed on Hybrid Principles for Finite Element Neutron Transport, Phd Thesis, University of London. Mirza, A. M. (1994) Private communication. Umpleby, A. P. (1994) Private communication. Abuzid,
and Transport,
Pergamon
SYNTHESIS
OF CONSERVATIVE DISCONTINUOUS DIFFUSION THEORY
R.T. Ackroyd, A. M. Gashut,
0.
SOLUTIONS
FOR
A. Abuzid
Mechanical Engineering Department, Imperial College of Science, Technology and Medicine, Exhibition Road, London, SW7 2BX, U.K. (e-mail: [email protected])
(Received 14 April 1995) ABSTRACT-A finite difference scheme for diffusion theory is presented which is based on the finite element method, and which is as efficient and accurate as the nodal method using a quadratic approximation in the local coordinates for the flux. A series of test problems in one and two dimensional geometries show that the method produces accurate solutions without the need for the several iterations of Copyright 0 1996 Elsevier Science Ltd the nodal method.
1. INTRODUCTION This paper deals with the development of a single and double quadratic non-iterative nodal method based on a discontinuous finite element solution obtained by the K’ variational principle of Ackroyd (1986a). To start the solution, a piecewise constant discontinuous finite element solution is found. For brevity these elements are called primitive elements. The K’ principle with discontinuous finite elements aims at cutting out as far as possible iterations. Usually two applications of the K’ principle are sufficient to improve upon the primitive element solution which is in effect the solution of the well known central difference equation of the diffusion theory. The following steps are considered : 1. Determine a piecewise constant discontinuous primitive elements with the K’ principle.
finite
2. Use the variational solution profiles say, 41 and Br.
of 4 to form two different
3. Combine
solutions
the extrapolated
4. Optimize 5. Repeat
(~2, varying
4s by maximizing if desired
to obtain
extrapolated
solution
(4) by the use of
solutions
with quadratic
in the form
4s = dl with the coefficient
element
+ (1 -
from element
a)
01
to element.
K’ (&) with respect to oz. +,,+I, (say). 1381
(1)
1382
R. T. Ackroyd et al.
1. MATHEMATICAL
FORMULATION
The I(* variational principle The diffusion equation is
(2) The K* (4) principle states that for a volume V partitioned by a surface Sr,
I(* (4)
=
/ {2&S’- C42 - D (V+)2} dV ”
-
J w[+(cto)-$(r-o)12dS Sl
-
I
F($o&o) t /,, (‘*rdJ2dS
(3)
where
F(+oo,40)
=
/ {D(V+)2 t C+2}dV ”
and Sb denoting a bare surface. The penalty parameter w is positive. It can be chosen so that for rectangular mesh and primitive elements the variational solution is the same as the central difference equation for diffusion. Formulae for this w are given in section (4.2). The (p •l-1)th trial function &,+rl(z, y) is given in terms of &(z, y) and an associated function 8,(z, y). For example over the element V”” of Figure 1, with local coordinate (t, q), A+1 (L 7) = o::iA
(L 11)+ (1 - or-r) 0, (& 77)
(5)
where &, at p=l is the flux obtained from the discontinuous finite element, 0, is its associated function. The coefficients crp”+;,(p = 1, . . . . .P - 1) are determined by maximising K* (&,+I) with respect to the coefficients ap”+; for all m and n of the mesh. The flux +,, (t, YJ)is defined over the element V”” for -1 5 t 5 1, -1 < q 5 1, in terms of the Legendre polynomials of e and q
4p(t, 8) = 4”” (P)PO(0 fi (rl)
t
$P;“”(P>- 4Y WI
Pl(O PO(17)
t ; WY (P)- 44mn (PII PO(0 Pl (VI
where :
+
f [4Y
(?+-I-~;"
t
; [4;" (PI t4T"
(P> - 24""
(P,
b)lPz
- W”” bwm
(0
w?)
p2(rl)
(f-9
1383
Diffusion theory
Figure 1: Typical element (m,n) for non-iterative anticlockwise (1,2,3,4). t, q are as defined before
7*..*P, 4mn (PI 4T” (P)
nodal; method.
The surfaces are numbered
p=l
is the average of &, over V”“. is the average of +P taken over the face i,(i=1,2,3,4) of the element.
The associated function flP(E, q) for V”” is expressed in the following from; 0, (6 V)
=
0”” (P) PO(0 PO(rl)
+
f PY
(P) - OF
(P)l Pl
[@Y
(P>
-
or
(P)J
(0 PO
PO (0
(77)
Pl
(7))
+
-5
+
;
[@T-”
(P)
+
@Y
(PI
-
2@“”
(PN
p2
(0
PO
(7)
+
f
[W
(P)
+
Or;”
(P)
-
z3””
(P)]
PO (0
p2
(71)
(7)
where : 0”‘” (p) is the average of 6, over V”” Or” (p) is the average of 8, taken over the face i, (i=1,2,3,4)
of the V”“.
The neutron conservation condition for 8, (t, q) over the element V”“, /_: /_:
[-Dm”v2e, + Yv,
-
sq =0
is used to express 0”” (p) in t erms of the 07” (p), which are determined as follows;
(8)
1384
R. T. Ackroydet
Construction
al.
of q51(t, 11) and 81 (E, 7)
(a) The double quadratic scheme To start the calculations the functions &([, q) and 81 (I, 7) are constructed from the results of the simplest possible finite element calculations based on the K* (4) principle, i.e the calculations based on the primitive elements. For primitive elements the trial function 4 = amn, where the am” are constants. Construction of $~r(t, q) from the am” dl (6, II) is defined by equation (6) with the #” (1) defined by the average interface values
1 [am+l*n +
amn]
C#J~(1) = i [am-l*n +
amn] amn]
4;” (1) = @” (1) =
1[am*n+l + amn]
4”” (1) = f [am~n-l +
(9)
The average element flux 4”‘” (1) is specified by the conservation condition (10) Construction of 8i([, q) from amn 0 ([,q) is defined by equation (7) with the ey (1) are the average values of 8i (E, 11) over the surfaces (i = 1,2,3,4) and 0”” (1) is the average of 81 ([,v) over Pn. These averages are determined in terms of the following procedures; Fit a quadratic to the amn at < = 0, the mid-plane of V”“; amtlln at t = 2 mid-plane of Vm+l+; and am-l+ at t = -2 the mid-plane of Vm-l*n. i.e. if(t)
= [am-l*nt (t - 2) - 2amn (< - 2) (< t 2) t amtltnt (e + 2)]
(11)
Set; @y (1) = f(1) = L [3amt1+ + 6a*n @y” (1) = f(l) = i [3a*-l+ + 6a*n
_ a*-ltn] _
a*+l+]
(12)
Similarly @r @r
(I) = f(I) (1) = f(I)
= 1[3a”‘*“+r + 6am” _ am*n-i] = i [3am.n-r + 6e”‘n _ e?n+r]
(13)
Also 0”” (1) is got from the conservation condition;
J-l/_: [-D”“v2elt zmnel- smn]d&iv =o
(14)
(6) The single quadratic scheme In this scheme the flux $r([, II) is taken to be a constant over an element , i.e for V”” &(<,q)
=4”“(l)
= f#~r(l) =amn(i
= 1,..,4)
The construction of the component function 8, (c, 9) is based on dP (t, 11)in the following way : Let ST”(p) be th e average value of &p(<,q) (t, q) along DB of Figure 1.
(15)
Diffusion theory
1385
A quadratic P(E) is fitted in Figure 1 to the average values ,7-l*” the mid-element lines GH,DB and EF respectively, p (0 = ; [QY-l,n (P) E (E - 2) - 2T”
(p), \k‘y” (p) @y+‘” (p) for
(t - 2) (< + 3) + Q:+‘*” (P) E (t + 2,]
(16)
The functions Oy” (p) and Oyn (p) are specified as @Y (P) = f(l) W(P) = U-1)
(17)
1
i.e o;ln
(p) = i 3q+‘+
@Tn (p) = ;
t
_q+l.n
(P) + f3’QY (P) - Q?”
(P)]
(P) + 6’K’” (P) t 39:-l’”
(18)
(P)] I
Similarly by defining STn (p) as the average value of &, ([, 0) along CA in Figure 1, i.e
(19) The above fitting procedures gives @yn (p) = i [307’“+’ (p) + SST 0:”
(P) = 5 [-%““+1
(p) - @;Vn-l (p)]
(p) t SKI’? (p) + 3!P;*“-’ (p)]
(20)
Also 0”” (p) satisfies the conservation condition for V””
Thus all the terms of eP ([, q) are known. Hence it can be determined using expression (7). Element Conservation for the Double Quadratic Scheme In the sequel 6,, p 2 2 is constructed for the double quadratic so that it satisfies neutron conservation condition for every element. In particular for V”“. l$s=cz~~t[1-cX~]&
(22)
satisfies the neutron conservation condition because 42 and 02 both satisfy the condition. Continuing the procedures show that &+I, p > _ 1 satisfies the neutron conservation condition, Ackroyd (1991). Determination
of ~$2(I, q)
The trial function 42 (5,~) is a linear combination for each element of +l(t,q)
and e1 (E, 71)
i.e 42 (L 77) = 4” 41 (<,‘I)
-I- (1
-
e)
b(k
7)
(23)
Wliere the coefficients a?” for the element V”” are found by maximising K’ (&) with respect to cy?n. Since $1 is an approximate solution of the diffusion equation for the impressed source,S, and 01 is another approximate solution for the same source S, then & as defined above can be
1386
R. T. Ackroyd et al.
regarded as an approximate solution for the given source. Note that 4s satisfies the conservation condition for V”” + Cmn4& - S”“‘] d(dq = 0
I_: J_: [-D”“V242 because its component for V”“.
functions
3. DOUBLE Having as
found
41 and $2 satisfy
QUADRATIC
by construction
PROCEDURES
& using steps (1) and (2) on calculates
(24)
the conserva.tion
FOR
condition
p 2 2
for p 2 2 &+I (t, q) which is defined
: htdb?)
with the function function Let @y”
=
$3A4777)
+
[I - $$$ad
(25)
coefficients c$!!i to be found by maximizing IC with respect to orti. The component r$p (t, ~7) is k nown from previous calculations. The construction of the component & (I, q) is based on &, (t, 77) in the following way : (p) be th e average value of 4(p) (E, v) along DB of Figure 1.
‘K’YP)
=
;J_:9,(Wdrl
=
; [6Y
(P) - &F’ (~11
(P) - d?
A quadratic P(t) is fitted in Figure 1 to the average values ST-‘” the mid-element lines GH,DB and EF respectively, P (E) = ; [*:-r7n The functions
(p), ar”
(26) (p) @y+‘” (p) for
(p) t ([ - 2) - 297” (E - 2) (E + 3) + rlY+“” (p) t (5 + 2)]
Oy” (p) and Oyn (i) are specified
(27)
as
@Y(P) = f(l) wn (P) = f (-1)
(28)
i.e @yn (p) = ; $$y”” @y Similarly
by defining
(p) = ;
t
(P) •t @J’~” (P) - %“-lTn (P)] (P) + @T
~IJ?” (p) as the average
ep (P) . The above fitting
_q+‘*”
procedures @r
(29)
(P) + 3912-“” (P)]
value of &, (t, 0) along CA in Figure
1, i.e
= =
; [Wmn (P) - $7” (P) - d4mnb)]
(30)
gives
(p) = f [3’l?;‘“+’ (p) + SST” (p) - ,y”-’
@Tn (P) = i [-Q?“+
(p) + 6@yn (p) + 397”-’
(p)] (p)]
(31)
Diffusion theory
Also 0”” (p) satisfies
the conservation
condition
1387
for V””
using expression (7). In Thus all the terms of 0 (t, 7) are known. Hence it can be determined so that it satisfies neutron conservation condition for every the sequel B,, p 1 2 is constructed element. In particular for V”“.
satisfies the neutron conservation condition because 4s and 0s both satisfy the condition. Continuing the procedures show that 4 r+r, p > 1 satisfies the neutron conservation condition.
4. NUMERICAL
EXAMPLES
In order to investigate the usefulness of the method discussed above, a computer code, DNMFD, (Direct Nodal Method For Diffusion) was developed in 1-D (slab) and 2-D x-y geometry, three sample problems in one and two dimensions were used to test the method for accuracy of solution. One-dimensional
problems
Two problems were presented. The first is a test problem consists of two regions. It is used to check the accuracy of the non- iterative nodal method solution against the analytical solution. The geometry of the problem is given in Figure 2 and Table 1 presents the material properties. The boundary conditions used are perfect reflection at x=0.0 cm and at x=2.0 cm. These boundary conditions are no more than symmetry conditions. Consequently the tests are only of the capabilities of the schemes for solving the differential equation for diffusion. The flux distributions obtained with the discontinuous finite element using primitive element, double quadratic non-iterative nodal method and the classical (iterative) nodal method for 20 nodes, together with the analytical solution are presented in Figure 3 for comparison. Twelve iterations were used for the nodal solution. The percentage error in these solutions, as calculated in reference to the analytical solution is presented in Figures 4, 5, 6 and 7. A maximum error of 0.16% in the primitive finite element solution and 0.12% in the single quadratic non-iterative nodal solution is seen Figure 4. The single quadratic solution shows marginally better agreement with the exact solution than does the primitive finite element solution. The error in the noniterative nodal solution reduce to 0.045% when the double quadratic scheme is employed for the same number of nodes as seen in Figure 5. This show that the non-iterative nodal results are significantly more accurate than the primitive finite element results relative to the analytical solution. To assess the double quadratic non-iterative nodal solution against the iterative nodal solution, Figures 6 and 7 present the error distributions in both solutions for 20 nodes. A maximum error of 0.4212% in the iterative nodal solution with 6 iterations and 0.045% in the double quadratic non-iterative nodal solution is seen in Figure 6. The iterative nodal error reduced to 0.046% when 12 iterations are involved as seen in Figure 7. This show that non-iterative nodal method. gives results as good as the iterative nodal method but much faster. The second problem is a slab lattice cell consisting of one mean free path thickness of natural uranium as fuel and three mean free path of graphite as a moderator with a unit source. The problem has been solved by Nanneh(1990) for transport calculations. The geometry of the
R. T. Ackroydet
1388
al.
1.0
0.0
b
2.0
Wcm)
Figure 2: Geometry of the test problem
Table 1: Group cross sections and source density for Problem 1 Cross Setions
Region 1
0 Source Density 1
1.0
Region 2
;I 1
0.0
0
: . Iterative 1.25- 8 Primitive element _ 0 Non-iterative
l.O-
Figure 3: Analytical, primitive element, double quadratic solutions for Problem 1
non-iterative
and iterative
nodal
1389
Diffusion theory
9 %
on
0.5
0.7S
1.0
,.I
1.5
1.n
Di5mcehllhox-(m)
Figure 4: Error in central difference for 20 nodes for Problem 1
(K’)
and single quadratic
non- iterative
nodal
solutions
nodal
solutions
Primitive element
0:s
on
IR
125
I.5
1.78
Dk?tMoehlheXdlrecabn(an)
Figure 5: Error in central difference for 20 nodes for Problem 1
(I(*) and double
quadratic
non- iterative
1390
R. T.Ackroyd et al.
Dls?ancehIhsXdiredlW,(an)
Figure 6: Error in double quadratic non- iterative and iterative nodal solutions for 20 nodes nodes after 6 iterations for Problem 1
Figure 7: Error in double quadratic non- iterative and iterative nodal solutions for 20 nodes after 12 iterations for Problem 1
1391
Diffusion theory
Figure
8: Geometry
of the repeating
lattice
cell problem
problem is given by Figure 8 while the material properties and the source density in each region are presented in Table 2. This problem is also a test of the capability of DNMFD calculations to give accurate results as compared to the finite element continuous solution presented by a diffusion theory solution given by the K+ principle encoded on EVENT, de Oliveira (1994). The solution obtained from DNMFD, primitive finte element method and iterative nodal method are presented in Figure 9, together with the continuous finite element solution for comparison. Figure 10 presents the error distribution of the primitive finite element and the double quadratic non-iterative nodal solutions relative to the continuous finite element solution. A localised maximum error of -1.3% in the primitive finite element solution and -0.91% in the double quadratic non-iterative nodal solution with 3 nodes in the fuel region and 10 nodes in the moderator region can be seen in the figure. To check the double quadratic non-iterative nodal solution against the iterative nodal solution, Figures 11 and 12 show the percentage error in both solutions relative to the continuous finite element solution based on K+ solution. A maximum error of -0.72% in the iterative nodal solution after 6 iterations and -0191% in the double quadratic non- iterative nodal solution is seen in Figure 11. The error in the iterative nodal solution reduced to 0.5% when 10 iterations are involved. The non-iterative nodal method shows once again that it is a promising method in comparison to the continuous finite element method and that it overall provides a solution nearly as accurate as the iterative nodal method solution but much faster. To assess the single quadratic non-iterative nodal method in obtaining accurate solutions for lattice problems in comparison to the contiimous finite element solution, Figure 13 presents the error distribution in the primitive finite element solution and the single quadratic non-iterative nodal solutions relative to the continuous finite element solution (K+). As it is seen from the figure, the single quadratic non-iterative nodal method gives solution of comparable accuracy This primitive solution is the same to the use of the K* principle with primitive elements. as that given by the classical difference equation for diffusion. In the sequel only the double quadratic scheme is considered as a non-iterative nodal scheme, as it is superior in accuracy to the single quadratic scheme, and it does not involve extra computation.
R. T. Ackroyd et al.
1392
Table 2: Group
cross sections
Cross Setions
and source density
for Problem
Region
2
1
Region
2
111 1 Source Density
1
1
1.0
0.0
fl
Iterative Primitive element Non-iterative Finite element (K+)
P
IO
1,
12
la
14
olstance in the x ckaion (cm)
Figure
9: Neutron
flux profile for the one group, one dimensional
lattice
cell problem
Diffusiontheory
1393
4.!26-
Primitive element Non-iterative
Q.S-
B
Ii
is
i 47S-
.o-
-1
-l.zs-
1 P
10
11
12
II
14
13
19
1,
OiaanmhtheXd~n(an) Figure 10: Error in primitive finite element and double quadratic non-iterative solutions for 13 nodes for Problem 2
nodal method
0 Non-iterative
Figure 11: Error in double quadratic non-iterative for 13 nodes after 6 iterations for Problem 2
nodal and iterative nodal methods solutions
1394
R. T. Ackroyd et
al.
Iterative 0 Non-iterative
0
Dismnca IntheX dlreabn(an) Figure 12: Error in double quadratic non-iterative nodal and iterative nodal methods solutions for 13 nodes after 10 iterations for Problem 2
LI Primitive element 0 Non-iterative
Figure 13: Error in primitive finite element and single quadratic non-iterative solutions for 13 nodes for Problem 2
nodal method
1395
Diffusion theory Two-dimensional
The penalty
Problems
parameter
w discussed
in section
(2) is given by
(34)
One case study is considered in 2-D (x-y) geometry. The problem is used by Fletcher for transport test calculations. The geometry and the neutronic parameters are given if Figure 14 and Table 3 respectively. It has been selected to check the validity of the non- iterative nodal method in obtaining solution of comparable accuracy to the continuous finite element solution As it can be seen based on K+ principle, de Oliveira (1986,1994), and using quadratic elements. from the figure, the problem is a square source problem with a constant source in a (1.2c97~)~ conditions for the sides adjacent to the source region are perfect reflector, region. Boundary while that for the other sides are bare surface conditions. This problem has been solved with the x and y pitches being 0.2 cm. The double quadratic non-iterative nodal, the primitive finite element , the K+ and the iterative nodal solutions are presented in Figure 15 at y=O.l cm and Table 4 gives these results for comparison. From both the table and the figure, it is clear that all the methods provide solutions near to the continuous finite element solution based on the K+ principle. Figures 16, 17 and 18 give the error distributions in these solutions relative to K+. A maximumerror of 3.1% in the primitive finite element solution and 2.4% in the double quadratic non-iterative nodal solution can be seen in Figure 16. Comparison between the errors in the double quadratic non- iterative nodal method and the iterative nodal method (see Figures 17 and 18) show that the non-iterative nodal method gives solutions overall as accurate as the iterative nodal but much faster.
6.0
“IOU I.2
Figure
rczz
14: Geometry
6.0
of Fletcher
X(m)
problem
1396
R. T. Ackroydet al. Table 3: Cross sections and source data for Fletcher’s problem Cross Section cm-’
Region 1
Region 2
Et E,
1.0 0.0
1.0 0.0
Source Density
0.69444
0.0
Table 4: Continuous finite element (I(+), iterative nodal solutions at y=O.l cm
non-iterative nodal, primitive finite element and
Distance along x axis (cm)
Continuous Finite Element Method (Pi)
Iterative Nodal Method (18 Iterations)
Primitive Finite Element Method
Non-Iterative Nodal Method
0.1000E+00 0.3000E+OO 0.5000E+OO 0.7000E+OO 0.9000E+OO O.llOOE+Ol O.l300E+Ol O.l500E=Ol O.l700E+Ol O.l900E+Ol 0.2100E+Ol 0.2300E+Ol 0.2500E+Ol 0.2700E+Ol 0.2900E+Ol 0.3100E+Ol 0.3300E+Ol 0.3500E+Ol 0.3700EtOl 0.3900E+Ol
0.541233+00 0.531813+00 0.511593+00 0.477623+00 0.425043+00 0.346493+00 0.252533+00 0.173023+00 0.118383+00 0.809543-01 0.553703-01 0.37896EOl 0.25961EOl 0.178033-01 0.12223EOl 0.840123-02 0.578073-02 0.398193-02 0.274553-02 0.189493-02
0.5427443+00 0.5334193+00 0.5133983+00 0.4797673+00 0.4277093+00 0.3499363+00 0.2500263+00 0.1713033+00 0.1172073+00 0.801531EOl 0.5482273-01 0.3752153-01 0.2570403-01 0.1762733-01 0.121021EO1 0.831813E02 0:572359E02 0.39425OEO2 0.271839E02 0.1876153-02
0.542363+00 0.533243+00 0.513683+00 0.480923+00 0.430353+00 0.354993+00 0.244713+00 0.168393+00 0.115773+00 0.795673-01 0.546983-01 0.376223-01 0.258973-01 0.178423-01 0.123043-01 0.849343-02 0.586863-02 0.405893-02 0.280983-02 0.194683-02
0.5423933+00 0.5331763+00 0.5133733+00 0.4801603+00 0.4289093+00 0.3525843+00 0.2464913+00 0.1693983+00 O.l16300E+OO 0.798152EOl 0.547815EOl 0.3761513-01 0.258436EOl 0.1776883-01 0.122263EOl 0.8419043-02 0.5801593-02 0.400061EO2 0.2760393-02 0.1905653-02
Note that the primitive, non-iterative and iterative nodal solutions have a jump in the error at the interface of the source and source-free regions. Both the non-iterative and the iterative nodal solutions recover from the error jump. 5. CONCLUSIONS Based on the results of the numerical examples presented in the previous section, we can make a number of general inferences about the potential of the non-iterative nodal method
1397
Diffusion theory 0.6
0.5
-
NOtI- h-auw !soluuoll Duf. sohdon
cmtra
wuuan conunwur (Pl) !soluuon
1 nwahw 0 0.4
0.3
0.2
0.7
1.5
Distance
2.0
2.5
3.0
3.5
in the X direction (cm)
15: FEMFD, DNMFD and INMFD fluxes in comparison to the K+ finite element solution at Y=O.l cm
Figure
0.6
1 .o
1.5
2.0
2.5
3.0
3.5
Distance in the X direction (cm)
Figure 16: Primitive element and double quadratic non-iterative nodal method error in reference to K+ solution at Y=O.l cm
1398
R. T. Ackroyd
et al.
Figure 17: Iterative nodal and double quadratic non-iterative nodal methods error in reference to K+ solution after 2 iterations at Y=O.l cm
-2o!
,,,,,,,,,,,,I
,,,,I
0.5
1.0
1.5
,,,,,,,,1,,1,,,,,,,
2.0
2.5
3.0
5.5
Distance in the X direction (cm)
Figure 18: Iterative nodal and double quadratic non-iterative nodal methods error in reference to K+ solution after 18 iterations at Y=O.l cm
1399
Diffusion theory
as new numerical technique. The results show fairly good agreement with the analytical and results. Due to the iterative nature of the classical nodal the continuous finite element (K+) expansion method, it involves more calculations in converging towards the true solution than do the direct approach of non-iterative, double quadratic nodal method. The non-iterative nodal method shows promise as a simple method since its use does not require the specification of convergence criteria and the initial guess of the classical nodal method. It starts with the central difference approximation given by the use of primitive elements in the K* principle and replaces that solution with a quadratic profile for each element and with neutron conservation holding for each element. Iterations can be implemented as indicated in section 3 for the double quadratic scheme if required. In developing a finite element’method (Non-iterative Nodal Method) which is akin to conventional nodal methods two new techniques have been implemented with discontinuous finite element method : i. Ensuring neutron balance the flux in each element.
element
by element
by the addition
ii. Constructing a polynomial profile for the flux over an element, in a flux calculation only one unknown per element.
of constant
components
for
by using at any one step
method for neutron The technique (i) is also used by Umpleby (1994) in his discontinuous component of the angular flux is transport in a generalized way. Each spherical harmonic made to satisfy neutron conservation element by element. The direct nodal method and his differ in the treatment of discontinuity. In the direct nodal the trial function is completely discontinuous between elements. In his scheme the trial function have both continuous and discontniuous components. The discontinuities are introduced in the non-iterative method by the choice of a penalty parameter which ensures with the simplest of elements, primitive elements, the central difference equation for diffusion theory. The interface term in the K* principle, as weighted by the penalty In Umpleby’s scheme the element parameter, serves to couple the elements into a system. coupling is achieved via the continuous component of the trial function. The techniques (i) and (“) 11 are also applicable to transport theory, and their investigation for transport calculations is one of the topics under study by Mirza (1994). The use of primitive elements (constant trial function) reduces the number of unknowns, since there is only one node per element. However, higher order trial functions can be used with the I(’ principle. REFERENCES 0. A. (1994) D iscontinuous Finite Element Solutions for Neutron Diffusion PhD Thesis, University of London. Ackroyd, R.T. (1986a) Prog. Nucl. Energy, 18,45. Ackroyd, R.T. (1986b) Prog. Nucl. Energy, 18, 7. Ackroyd, R.T. (1992) Ann. Nucl. Energy, 19, 565. Ackroyd, R.T., Abuzid, 0. A. and Gashut, A. M. (1995) Ann. Nucl. Energy, in press. de Oliveira, C.R.E. (1986) Ann. Nucl. Energy, 13, 663. de Oliveira, C.R.E. (1994) Private communication. Nanneh, M. M. (1990) A Synthesis Method Bzsed on Hybrid Principles for Finite Element Neutron Transport, Phd Thesis, University of London. Mirza, A. M. (1994) Private communication. Umpleby, A. P. (1994) Private communication. Abuzid,
and Transport,