- Email: [email protected]
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems
Engineering Analysis with Boundary Elements 11 (1993) 39-45
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems Masafumi Itagaki Japan Atomic Energy Research Institute, 2-4 Shirakota-Shirane, Tokai, Ibaraki 319-11, Japan
&
Carlos A. Brebbia Wessex Institute of Technology~University of Portsmouth, Ashurst Lodge, Ashurst, Southampton SO4 2AA, UK (Received 12 November 1991; accepted 29 January 1992) The multiple reciprocity method has been applied to one-group fission neutron source iteration problems. The domain integral related to the fission neutron source in the mth source iteration stage is transformed into a series of m - 1 boundary integrals. These integrals require the zero order to the m - 1 order fundamental solutions and the boundary neutron fluxes and currents calculated in the previous m - 1 source iterations. Results obtained in a simple test calculation indicate that the source iteration based on the present technique converges rapidly to the desired eigenvalue as long as a certain condition is satisfied.
Key words: boundary element method, multiple reciprocity method, higher order fundamental solution, neutron diffusion equation, fission source iteration, modified Helmholtz equation. 1 INTRODUCTION
reciprocity method (MRM)' which is known as another elegant technique to transform domain integrals into boundary ones. This method is based on a repeated application of the reciprocity theorem using a sequence of higher order fundamental solutions. This method does not require any approximated source functions, and enables one to receive accurate solutions. The present paper discusses the application o f this method to fission neutron source iterative problems in nuclear reactor analyses. Applications in this paper are restricted to one-energy-group problems for simplicity. In the application of the M R M to neutron diffusion problems, one uses the higher order fundamental solution of the modified Helmholtz equation, and the series of transformed boundary integrals consists of a finite number of terms.
The boundary element method (BEM) 1 has been applied in many engineering fields since 1978. The first application of the BEM to a neutron diffusion problem was made by the first author and further research is now under development. 2-7 The main difficulty when applying boundary element methods to neutron diffusion problems has been the computations of the domain integrals resulting from the presence of volume sources. In a previous paper, the authors have shown that the domain integral related to a fission neutron source can be transformed into an equivalent boundary one by the use of the dual reciprocity method. 7 The accuracy in the results obtained by the method, however, depends on approximated source functions, and the boundary integrations caused by the functions are rather complicated. Recently, Nowak & Brebbia 8 proposed the 'multiple
2 N E U T R O N D I F F U S I O N E Q U A T I O N AND FISSION S O U R C E I T E R A T I O N S C H E M E
Engineering Analysis with Boundary Elements 0955-7997/93/$05.00 © 1993 Elsevier Science Publishers Ltd.
The one-energy-group neutron diffusion equation is 39
40
M. Itagaki, C.A. Brebbia
written in the form
is denoted as
-DV2t~(r) + Ea~(r) =
(vEf/kgf)¢(r)
(1)
where D, Ea, vEf, keff and 4~(r) are called diffusion coefficient, absorption cross section, neutron emission cross section, effective multiplication factor (eigenvalue) and neutron flux. Equation (1) can be reduced to the Helmholtz equation V2,(r) + B2m0(r) = 0
(2)
if one defines the quantity B2m which is called 'material buckling' =
v ilkefr- ~'~a D
(3)
In a previous paper one of the authors demonstrated that eqn (2) can be solved using the BEM. 2 In this case one has to use the determinant search technique to find the eigenvalue. For reactor physics analysis one wishes to find only the eigenvalue of largest magnitude. The Helmholtz equation has many eigenvalues and it is known that the determinant search is troublesome for finding the maximum eigenvalue. In the field of nuclear reactor analyses, eqn (1) is usually considered as a modified Helmholtz equation with a source term, rather than the simple Helmholtz equation. That is,
V2O(r)-
k2$(r) +
S(r)/D
= 0
(4)
where the quantity k 2 is defined by k2 = r
o/o
vEfO(r)/keff
f So = 1"0
-ll
(m = 1)
(7a)
(m => 2)
(7b)
A unit fission source distribution S ( r ) = So = 1.0 is assumed as the initial estimate S 0). Obviously in this case a domain integral defined by Vc= I n S(')(r)df~ is equal to the total volume or area of the domain f~ under consideration. Here, another quantity is defined
vEf(b(m)(r)
G(m)(r) =
(m >=2)
(8)
from the calculated neutron flux distribution. The new estimate of the eigenvalue is computed as follows:
k(m) eft =
JG(m)df~/Js(m)df~ fl
= Jn
fl
o(m)df~/{fn dP(m-1)df~/k~f-1)}
(9)
Using this new eigenvalue, the source distribution is re-normalized as S (re+l) =
G(m)/k~ )
(10)
in such a way that the quantity (5)
and the fission neutron source distribution is given as S(r) =
_ DV2~b(m) + ~adp(m) = S (m)
(6)
Equation (6) shows that the fission source distribution is related to the unknown neutron flux distribution $(r). One can here introduce the 'power method' or the 'fission source iteration technique', as it is called in the field of nuclear reactor analyses, to find the maximum eigenvalue keff. The power method is guaranteed to converge automatically to the eigenvalue of largest magnitude. 9 This is the advantage of the power method over the determinant search for the standard Helmholtz equation. In addition, nuclear analysts should often solve the multi-energy-group neutron diffusion equations, i.e. simultaneous differential equations where each source term contains unknown fluxes of different energy groups. In this case each of the group's diffusion equation cannot be considered as a Helmholtz equation any more. It should be noted that the source iteration technique is also applicable to such multi-group problems. An outline of the source iteration technique is described as follows. At a general ruth stage, the neutron diffusion equation
Pc = Jfl S('+I)(I~) d~ =
IflV~f+ (m)d " l k ~ r )
(]l)
is constant through the iteration. This normalized source is used again in the neutron flux computations, and the source iteration is repeated until a given convergence criterion is satisfied. The proof of the convergence with this method is described in textbooks of numerical analyses. 9 The validity of the method was also demonstrated for the BEM with a domain source integral 2 and for the dual reciprocity BEM 7 in the authors' previous papers.
3 BOUNDARY INTEGRAL REPRESENTATION OF THE CRITICAL EIGENVALUE The multiplication factor or eigenvalue, kdr is originally defined by a ratio of two different domain integrals, as shown in eqn (9). However, these domain integrals can be transformed into boundary-only integrals. From eqn (7), one can write,
Io
1 •(m) dfl = ~a
Io
°I
s(m) d n + ~ a n V2t~(m) d fl
(12)
for a homogeneous domain ~. Using Gauss' theorem,
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems
boundary integral expression for the domain integral defined by eqn (19). Here the authors define the following domain integral for convenience:
eqn (12) can be rewritten as
[flfb(m'dQ:~-alfls(m) dl-l-~--aaJFJ(m)dF 1
"
J
J(m)dF
"'eft
(13) where the quantity J(") = - D [O¢(")/On] is the normal neutron current. The source distribution S0)(= So) is uniform when m = 1, hence for this case one obtains In gb(l)d F t = ~S~°aIn df~ - ~l-alr J(l) dF _ So . a -
1 _Iv j(I) dF
(14)
where A = J'n d r is the total volume or area. In this manner the calculation of the eigenvalue in eqn (9) which is, in principle, involved in domain integrals, can now be expressed in terms of boundary integrals only.
TO NEUTRON
DIFFUSION
6 i -m--0
(L = 0)
The authors also define the normal neutron current for the mth iteration stage
-D[OC ~/On]
(17)
which corresponds to the Lth order fundamental solutions ¢~(L}. The boundary integral equation corresponding to eqns (7a) and (7b) for ruth source iteration can be written in the form 2'3
hci~)lra) : JI (#(re'J;(°)dr --If"J(m)~9;(O)OF + QIm, (18)
where the choice for parameter ~'~is 0/(2~r) for a boundary point at a corner with an angle of 0. The domain source integral Q~m) is defined as
Qlm) = J¢lsCm)¢~(0)dr
k~_,) Jn ~(m-')(V2q~7(L+O - k2~bT(L+'))df~
-- /.(m-l)V~f (Ift q~(L+l)(v2o(m-l)-k2o(m-l))d" n-eft [ ~b(m_l)0 ¢~*(L+I) +Jr
~nn ;
dr
--IV gb~(L+l)0) O~b(m-l) n all' _ /Sb-(m-1)V~f{Is(m_l)dp~(L+l)d~[2
:',
= ~t.(m_l)
{
+'>JIm-"dr} c Z i ( m - 1, L + 1)
(21a) When m = 1 and L == 1, the uniform source SO)(= So) can be taken outside the domain integral, so one obtains Zi(1, L) = Ja S (1)~b*(L)d r
(16)
and the following quantity S ~(~)=
.r:f
Zi(m, L) =
(15a)
V2~b~(L) - k2~b~(L) + ~b~(L-I) = 0 (L = 1,2,...) (15b)
j (m) = _D[O~)(m)lOn]
(20)
When m > 2 and L > 0, the quantity Zi(m, L) can be rewritten as
PROBLEMS
Here the authors introduce a sequence of higher order fundamental solutions o~(L) which are defined by the following subsidiary equations, each of which corresponds to eqns (7a) and (7b). V2q~*(°) - k2~b~(°) +
Zi( m , L) = In S (m)~b~(L)dr
+J
4 MULTIPLE RECIPROCITY T E C H N I Q U E APPLIED
41
-- Jo + : ' _Sol 0
=
+
s,l,o:-,,
'
Ir j;(L) dr + ~-fZi(1,L-
1)
(21b)
When m = 1 and L = 0, the quantity Zi(1,0) is equal to the integral defined by eqn (19): Zi(1, O) = .In
S(1)~b~(°)dr-
QI l)
(21c)
The above domain integral QI l) for a uniform source can be transformed into a boundary-only integral
(19)
The aim in the present paper is to find an equivalent
using Gauss' theorem. The derivation to obtain eqn (22) is found in the authors' previou~ papers. 2'3
M. Itagaki, C.A. Brebbia
42
Using the above three types of recurrence relations (21a)-(21c), the domain integral Q}m) can be transformed into a boundary-only integral as follows. The source integral QI 1) for the first iteration stage can be expressed by the following boundary integral:
--- /]~f (m-i){Ir (~(m-1)J:O)dF _ Jr ~;(l)j(m-I) Dk eff 4
/./Ef
(m-l)
Dkeff
dr}
/]~f (J ~(m-2)j~(2) dF (m-2)
Dk eff
r
O}') = Jfl S(')~b; (°) df~ = (23a)
.rs
4
Using eqn (21a), the source integral for the second stage is expressed by
• Zi(m - 2, 2)
Of 2) = If~ S (2)t~(0) d~ __ /~b (l)#']zf (Zi(1,
Finally one obtains the following boundary-only integral expression which is equivalent to the domain integral (19) for the fission neutron source in the mth source iteration.
l)+ Ir ~b(l)J;(')dr
m-l~ L
uE
\
:"
The quantity Zi(1, 1) in eqn (23) is given as
_ So jT/l) dr)
Z;(I' 1)= ~2 { Q}') -DJr
where Q}I) is known from eqn (23a). Similarly + 013) = Ia S(3)~b;(°)dO
,_-2-~--£) • Zi(1, m - 1) kZ=l Dkeff J
(m >_-2) (24a)
- nz./2)vEf{Z/(2, 1)+ Jr
q~(2)j;(l)dF
where
Zi(1,m-1) =-~{Zi(1,m-2)-~JrJ*(m-l) dr } (m > 2)
(24b)
(m = l)
(24c)
and
: FIb(2) + n~(2) " nl.(O
VEf
vY],f
+ Fib(2) -"-''eft " Dk~
. z'(l'2)
(23c)
where Zi(1,2) = ff1 { Z i ( 1 , 1 )-S °D[ Jr J*(2)dF} ' In general, the source integral Q}") for the ruth iteration (m > 2) can be, therefore, expressed in the following form: Q}m) : f s(m)(b~(O) df~ = Zi(m, O) d _
uEf
4 D ~ I ) Zi(m - 1, 1)
There are no unknown variables in eqns (24a)-(24c). These expressions are applicable to both two- and three-dimensional problems•
5 HIGHERORDER FUNDAMENTALSOLUTIONS FOR TWO-DIMENSIONALPROBLEMS Higher order fundamental solutions are required for the multiple reciprocity computations. In the case of a twodimensional, one-group neutron diffusion problem, the fundamental solution of zero order is given by
(~;(o) = Ko(kr)/(2rc)
(25)
where the quantity k is given as k 2 = Za/D. The Lth order fundamental solution to eqn (15b) has the form
dp~(L) = Az(kr)LK£(kr)
(26)
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems where the coefficients AL are obtained recurrently from the relationship
AL = AL_I/(2Lk 2)
for L = 1,2, 3,...
1.2
~
10 °
/
(27)
The derivation of eqns (26) and (27) are found in another paper by the authors. 1° It is noted that the function o~(L-)has no singularity for L > 1, so that the integration of these functions does not require any special technique.
I
43
1.1
lO-Z ~0
v
¢0 :>
'
1.0
eigenvalue
C 0
k e f f = 1.0
I 0 - = .,~
g
O a
bo / x
%J
10 `3
0.9 ¸ to
6 RESULTS OF A S I M P L E TEST CALCULATION
g
8 LU 0
A simple test calculation has been made. A homogeneous square domain with dimensions of 50crn by 50cm has been chosen for the example. A zero-flux boundary condition is set for two sides and a zerocurrent condition for the other two sides of the outer boundary, as shown in Fig. 1. In this case the analytic flux solution is given by
(~(x,y) cx sin (B~x) . sin (Bc y)
(28)
with the geometrical buckling Baa = 2(1r/100) 2. When the neutron source iteration converges, the neutron diffusion eqn (1) should be equivalent to the Helmholtz equation V2q~ +/~Gq~ = 0
(29)
with the relationship =
D
- B2m
(30)
The quantity B~ is called the material buckling. For convenience one assumes here that the nuclear con-
Y
(cm) 5O
10-*
0.7 0
I I 20 Number
f I 40 of source
I I 60 iterations
~ 80
I
10 " s 100
Fig. 2. Convergence behaviour of source iteration. stants are given as: D=l.0(cm),
~ a - - 0 " 0 1 ( c m -1)
and
v~f = 0.011 9739 (cm -1) so that the geometrical buckling is equal to the material buckling when kerr = 1"0. For the BEM calculation, each side of the square domain is divided into five equal 'linear' boundary elements. The source iteration started with a flat source distribution S0)(r) = 1.0 as an initial estimate. Figure 2 shows the convergence behaviour of the eigenvalue ken' from the first to the 100th iteration, where the higher order fundamental solutions were generated up to the 100th order. The dotted line represents the variation in the eigenvalue deviation at the ruth iteration defined by
~=0 e --
uZf=O.O119739 T.a =0.01 J=0
0.8
cm-'
cm-t
D = 1.0 cm
J= 0
¢~=0
50 (cm)
Fig. 1. Problem specifications for the test calculation.
X
/.(m-l)
(31)
During the early stage of iterations, the calculated eigenvalue rapidly approached the true value (ken' = 1-0): the 7th solution was kerr= 0.99843 and the 8th was kerr = I'00070. This rapid convergence suggests that there is no serious mistake in the present formulation. Note that acceptable results are obtained if one stops the iteration here. However, after the 8th iteration the calculated eigenvalue gradually goes away from the true value. The deviation e has the minimum value e = 1.18 x 10-3 at the 12th iteration, but the deviation increases again as the iteration proceeds. It is interesting to point out that the eigenvalue seems to converge finally to the value of koo = UEf/Ea(= 1"19739). The quantity ko~ is called the 'infinite multiplication factor' for an infinite domain or a domain with no neutron leakage. In fact, the calcu-
44
M. Itagaki, C.A. Brebbia
lated neutron currents along the boundary represent very small values at the 100th iteration.
7 DISCUSSION According to the authors' further investigation, the convergence of the Lth order term in eqn (24a) at the mth source iteration depends on whether the following condition vEf • t-(m-Z ) r / = ~a %ff
L~
1
< 1"0
(32)
is satisfied or not. The calculated eigenvalue approaches the true value efficiently as long as the quantity r/is less than unity. The reason for this can be described as follows. As shown in eqn (24a), the Lth order fundamental solution $i*(L) given by eqn (26), always associates with the factor m-I
V~f
H i)b(m-L) L=1 ~'~eff then one can make another definition of the Lth order solution as
~(L) = JiL(kr)LKL(kr)
(33)
with the relationship of the new coefficients
-
AL-I
Um_t )
(34)
2L. ~efr
In order to investigate the convergence of the newly defined fundamental solutions, one estimates the following quantity
ff,(L)
u~flEa (kr)LKL(kr) = lim TM = lim lr(m_L ) • (kr)L_lKL_l(kr) r-~O~
uEf/E a
2 L-I" (L - 1)!
_
r/= 2L. k ~ -L) " 2L-2" (L - 2)! -
uEf/Ea
L - 1
b (m-L)
L
(36) The value of the function -4L xLKL(x) is sequentially calculated using the recurrence formula
~,+1xL+IKL+I(X) = X2{~L_lXL-IKL_I(X)} • ('4~+l/A~-~)
+ 2I~(~,xLK~(x)} • (,L+I//L-~)
(37)
This recurrence scheme starts with the values of the
functions K0(x) and xKl (x), which are calculated using polynomial approximations for the modified Bessel functions found in the literature) l The polynomial approximations originally contain some round-off errors. It is understandable that the errors are magnified if the condition r} < 1.0 is not satisfied during the recurrence calculations. This is supposedly the reason why the source iteration does not converge to the true value. Why did the method seemed to converge finally to the false eigenvalue
k(m-L) VEf eft ~ k ~ = Ea in the test calculation? Because of the round-off error, the source integral QI m) was not calculated properly. The calculated eigenvalue has no physical meaning any more and gradually increased. Finally it had a value very close to k~ and the quantity r/ happened to become less than unity (r1 ~ (L - 1)/L < 1.0). Consequently the iteration apparently converged because of no addition of round-off error accumulation. Nowak 12 pointed out that a similar phenomenon may occur in his M R M computations for solving Poisson-type problems. He also recommended that one should scale and reformulate the problem, i.e. all dimensions should be divided by the maximum dimension. However, such a scaling cannot be applied to the authors' modified Helmholtz problems. One should seek another way of removing the numerical instability.
8 CONCLUSIONS The multiple reciprocity method has been applied to one-group fission neutron source iteration problems. The domain integral related to the fission neutron source in the ruth source iteration is transformed into a series of m - 1 boundary integrals. These integrals require the zero order to the r n - 1 order fundamental solutions and the boundary neutron fluxes and currents calculated in the previous m - 1 source iterations. The boundary-only integral expression proposed here for the iterative fission source is applicable to both twoand three-dimensional criticality problems. Any information in a domain is not required in this formulation, so the advantages of the BEM can be perfectly preserved. In addition, there is no need to assume any approximated source functions like ones introduced in the dual reciprocity method. Results obtained in a simple test calculation indicate that the source iteration based on the present techniqu e converges rapidly to the desired eigenvalue as long as the condition given by eqn (32), r / < 1'0, is satisfied. Unless the condition is satisfied, a numerical error is accumulated as the iteration proceeds. In order to avoid such situations, further
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems consideration should be incorporated into the present technique.
ACKNOWLEDGMENT The authors wish to express their gratitude to Dr. A.J. Nowak o f Silesian Technical University, Poland, for his valuable comments on the work and manuscript.
REFERENCES 1. Brebbia, C.A. The Boundary Element Method for Engineers, Pentech Press, London, 1978. 2. Itagaki, M. Boundary element methods applied to twodimensional neutron diffusion problems, Journal of Nuclear Science and Technology, 1985, 22(6), 565. 3. Itagaki, M. Boundary element techniques for two-dimensional nuclear reactor calculations, Engineering Analysis, 1987, 4(4), 190. 4. Itagaki, M. & Brebbia, C.A. Boundary element method applied to neutron diffusion problems, Proc. 10th Int.
5.
6.
7.
8.
9. 10.
11. 12.
45
Conf. BEM, Southampton University, UK, SpringerVerlag, 1988. Itagaki, M. & Brebbia, C.A. Boundary element method applied to neutron diffusion problems of a thin layer sandwiched between two zones, Proc. 12th Int. Conf. BEM, Hokkaido University, Japan, Springer-Verlag, 1990. Itagaki, M. & Brcbbia, C.A. Space-dependent core/ reflector boundary conditions generated by the boundary clement method for pressurized water reactors, Nucl. Sci. Eng., 1991, 107, 246. Itagaki, M. & Brcbbia, C.A. Boundary clement formulation of fission neutron source problems using only boundary integrals, Engineering Analysis with Boundary Elements, 1991, 8(5), 239. Nowak, A.J. & Brebbia, C.A. The multiple-reciprocity method. A new approach for transforming BEM domain integral to the boundary, Engineering Analysis with Boundary Elements, 1989, 6(3), 164-168. Wachspress, E.L. Iterative Solution of Elliptic Systems, Prentice-Hall, Inc., 1966. Itagaki, M. & Brebbia, C.A. Generation of higher order fundamental solutions to the two-dimensional modified Helmholtz equation, Engineering Analysis with Boundary Elements, 1993, 11(1), 87-90. Abramowitz, M. & Stcgun, I.A. Handbook of Mathematical Functions, Dover, 1965. Nowak, A.J. Private Communication, 1991.
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems Masafumi Itagaki Japan Atomic Energy Research Institute, 2-4 Shirakota-Shirane, Tokai, Ibaraki 319-11, Japan
&
Carlos A. Brebbia Wessex Institute of Technology~University of Portsmouth, Ashurst Lodge, Ashurst, Southampton SO4 2AA, UK (Received 12 November 1991; accepted 29 January 1992) The multiple reciprocity method has been applied to one-group fission neutron source iteration problems. The domain integral related to the fission neutron source in the mth source iteration stage is transformed into a series of m - 1 boundary integrals. These integrals require the zero order to the m - 1 order fundamental solutions and the boundary neutron fluxes and currents calculated in the previous m - 1 source iterations. Results obtained in a simple test calculation indicate that the source iteration based on the present technique converges rapidly to the desired eigenvalue as long as a certain condition is satisfied.
Key words: boundary element method, multiple reciprocity method, higher order fundamental solution, neutron diffusion equation, fission source iteration, modified Helmholtz equation. 1 INTRODUCTION
reciprocity method (MRM)' which is known as another elegant technique to transform domain integrals into boundary ones. This method is based on a repeated application of the reciprocity theorem using a sequence of higher order fundamental solutions. This method does not require any approximated source functions, and enables one to receive accurate solutions. The present paper discusses the application o f this method to fission neutron source iterative problems in nuclear reactor analyses. Applications in this paper are restricted to one-energy-group problems for simplicity. In the application of the M R M to neutron diffusion problems, one uses the higher order fundamental solution of the modified Helmholtz equation, and the series of transformed boundary integrals consists of a finite number of terms.
The boundary element method (BEM) 1 has been applied in many engineering fields since 1978. The first application of the BEM to a neutron diffusion problem was made by the first author and further research is now under development. 2-7 The main difficulty when applying boundary element methods to neutron diffusion problems has been the computations of the domain integrals resulting from the presence of volume sources. In a previous paper, the authors have shown that the domain integral related to a fission neutron source can be transformed into an equivalent boundary one by the use of the dual reciprocity method. 7 The accuracy in the results obtained by the method, however, depends on approximated source functions, and the boundary integrations caused by the functions are rather complicated. Recently, Nowak & Brebbia 8 proposed the 'multiple
2 N E U T R O N D I F F U S I O N E Q U A T I O N AND FISSION S O U R C E I T E R A T I O N S C H E M E
Engineering Analysis with Boundary Elements 0955-7997/93/$05.00 © 1993 Elsevier Science Publishers Ltd.
The one-energy-group neutron diffusion equation is 39
40
M. Itagaki, C.A. Brebbia
written in the form
is denoted as
-DV2t~(r) + Ea~(r) =
(vEf/kgf)¢(r)
(1)
where D, Ea, vEf, keff and 4~(r) are called diffusion coefficient, absorption cross section, neutron emission cross section, effective multiplication factor (eigenvalue) and neutron flux. Equation (1) can be reduced to the Helmholtz equation V2,(r) + B2m0(r) = 0
(2)
if one defines the quantity B2m which is called 'material buckling' =
v ilkefr- ~'~a D
(3)
In a previous paper one of the authors demonstrated that eqn (2) can be solved using the BEM. 2 In this case one has to use the determinant search technique to find the eigenvalue. For reactor physics analysis one wishes to find only the eigenvalue of largest magnitude. The Helmholtz equation has many eigenvalues and it is known that the determinant search is troublesome for finding the maximum eigenvalue. In the field of nuclear reactor analyses, eqn (1) is usually considered as a modified Helmholtz equation with a source term, rather than the simple Helmholtz equation. That is,
V2O(r)-
k2$(r) +
S(r)/D
= 0
(4)
where the quantity k 2 is defined by k2 = r
o/o
vEfO(r)/keff
f So = 1"0
-ll
(m = 1)
(7a)
(m => 2)
(7b)
A unit fission source distribution S ( r ) = So = 1.0 is assumed as the initial estimate S 0). Obviously in this case a domain integral defined by Vc= I n S(')(r)df~ is equal to the total volume or area of the domain f~ under consideration. Here, another quantity is defined
vEf(b(m)(r)
G(m)(r) =
(m >=2)
(8)
from the calculated neutron flux distribution. The new estimate of the eigenvalue is computed as follows:
k(m) eft =
JG(m)df~/Js(m)df~ fl
= Jn
fl
o(m)df~/{fn dP(m-1)df~/k~f-1)}
(9)
Using this new eigenvalue, the source distribution is re-normalized as S (re+l) =
G(m)/k~ )
(10)
in such a way that the quantity (5)
and the fission neutron source distribution is given as S(r) =
_ DV2~b(m) + ~adp(m) = S (m)
(6)
Equation (6) shows that the fission source distribution is related to the unknown neutron flux distribution $(r). One can here introduce the 'power method' or the 'fission source iteration technique', as it is called in the field of nuclear reactor analyses, to find the maximum eigenvalue keff. The power method is guaranteed to converge automatically to the eigenvalue of largest magnitude. 9 This is the advantage of the power method over the determinant search for the standard Helmholtz equation. In addition, nuclear analysts should often solve the multi-energy-group neutron diffusion equations, i.e. simultaneous differential equations where each source term contains unknown fluxes of different energy groups. In this case each of the group's diffusion equation cannot be considered as a Helmholtz equation any more. It should be noted that the source iteration technique is also applicable to such multi-group problems. An outline of the source iteration technique is described as follows. At a general ruth stage, the neutron diffusion equation
Pc = Jfl S('+I)(I~) d~ =
IflV~f+ (m)d " l k ~ r )
(]l)
is constant through the iteration. This normalized source is used again in the neutron flux computations, and the source iteration is repeated until a given convergence criterion is satisfied. The proof of the convergence with this method is described in textbooks of numerical analyses. 9 The validity of the method was also demonstrated for the BEM with a domain source integral 2 and for the dual reciprocity BEM 7 in the authors' previous papers.
3 BOUNDARY INTEGRAL REPRESENTATION OF THE CRITICAL EIGENVALUE The multiplication factor or eigenvalue, kdr is originally defined by a ratio of two different domain integrals, as shown in eqn (9). However, these domain integrals can be transformed into boundary-only integrals. From eqn (7), one can write,
Io
1 •(m) dfl = ~a
Io
°I
s(m) d n + ~ a n V2t~(m) d fl
(12)
for a homogeneous domain ~. Using Gauss' theorem,
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems
boundary integral expression for the domain integral defined by eqn (19). Here the authors define the following domain integral for convenience:
eqn (12) can be rewritten as
[flfb(m'dQ:~-alfls(m) dl-l-~--aaJFJ(m)dF 1
"
J
J(m)dF
"'eft
(13) where the quantity J(") = - D [O¢(")/On] is the normal neutron current. The source distribution S0)(= So) is uniform when m = 1, hence for this case one obtains In gb(l)d F t = ~S~°aIn df~ - ~l-alr J(l) dF _ So . a -
1 _Iv j(I) dF
(14)
where A = J'n d r is the total volume or area. In this manner the calculation of the eigenvalue in eqn (9) which is, in principle, involved in domain integrals, can now be expressed in terms of boundary integrals only.
TO NEUTRON
DIFFUSION
6 i -m--0
(L = 0)
The authors also define the normal neutron current for the mth iteration stage
-D[OC ~/On]
(17)
which corresponds to the Lth order fundamental solutions ¢~(L}. The boundary integral equation corresponding to eqns (7a) and (7b) for ruth source iteration can be written in the form 2'3
hci~)lra) : JI (#(re'J;(°)dr --If"J(m)~9;(O)OF + QIm, (18)
where the choice for parameter ~'~is 0/(2~r) for a boundary point at a corner with an angle of 0. The domain source integral Q~m) is defined as
Qlm) = J¢lsCm)¢~(0)dr
k~_,) Jn ~(m-')(V2q~7(L+O - k2~bT(L+'))df~
-- /.(m-l)V~f (Ift q~(L+l)(v2o(m-l)-k2o(m-l))d" n-eft [ ~b(m_l)0 ¢~*(L+I) +Jr
~nn ;
dr
--IV gb~(L+l)0) O~b(m-l) n all' _ /Sb-(m-1)V~f{Is(m_l)dp~(L+l)d~[2
:',
= ~t.(m_l)
{
+'>JIm-"dr} c Z i ( m - 1, L + 1)
(21a) When m = 1 and L == 1, the uniform source SO)(= So) can be taken outside the domain integral, so one obtains Zi(1, L) = Ja S (1)~b*(L)d r
(16)
and the following quantity S ~(~)=
.r:f
Zi(m, L) =
(15a)
V2~b~(L) - k2~b~(L) + ~b~(L-I) = 0 (L = 1,2,...) (15b)
j (m) = _D[O~)(m)lOn]
(20)
When m > 2 and L > 0, the quantity Zi(m, L) can be rewritten as
PROBLEMS
Here the authors introduce a sequence of higher order fundamental solutions o~(L) which are defined by the following subsidiary equations, each of which corresponds to eqns (7a) and (7b). V2q~*(°) - k2~b~(°) +
Zi( m , L) = In S (m)~b~(L)dr
+J
4 MULTIPLE RECIPROCITY T E C H N I Q U E APPLIED
41
-- Jo + : ' _Sol 0
=
+
s,l,o:-,,
'
Ir j;(L) dr + ~-fZi(1,L-
1)
(21b)
When m = 1 and L = 0, the quantity Zi(1,0) is equal to the integral defined by eqn (19): Zi(1, O) = .In
S(1)~b~(°)dr-
QI l)
(21c)
The above domain integral QI l) for a uniform source can be transformed into a boundary-only integral
(19)
The aim in the present paper is to find an equivalent
using Gauss' theorem. The derivation to obtain eqn (22) is found in the authors' previou~ papers. 2'3
M. Itagaki, C.A. Brebbia
42
Using the above three types of recurrence relations (21a)-(21c), the domain integral Q}m) can be transformed into a boundary-only integral as follows. The source integral QI 1) for the first iteration stage can be expressed by the following boundary integral:
--- /]~f (m-i){Ir (~(m-1)J:O)dF _ Jr ~;(l)j(m-I) Dk eff 4
/./Ef
(m-l)
Dkeff
dr}
/]~f (J ~(m-2)j~(2) dF (m-2)
Dk eff
r
O}') = Jfl S(')~b; (°) df~ = (23a)
.rs
4
Using eqn (21a), the source integral for the second stage is expressed by
• Zi(m - 2, 2)
Of 2) = If~ S (2)t~(0) d~ __ /~b (l)#']zf (Zi(1,
Finally one obtains the following boundary-only integral expression which is equivalent to the domain integral (19) for the fission neutron source in the mth source iteration.
l)+ Ir ~b(l)J;(')dr
m-l~ L
uE
\
:"
The quantity Zi(1, 1) in eqn (23) is given as
_ So jT/l) dr)
Z;(I' 1)= ~2 { Q}') -DJr
where Q}I) is known from eqn (23a). Similarly + 013) = Ia S(3)~b;(°)dO
,_-2-~--£) • Zi(1, m - 1) kZ=l Dkeff J
(m >_-2) (24a)
- nz./2)vEf{Z/(2, 1)+ Jr
q~(2)j;(l)dF
where
Zi(1,m-1) =-~{Zi(1,m-2)-~JrJ*(m-l) dr } (m > 2)
(24b)
(m = l)
(24c)
and
: FIb(2) + n~(2) " nl.(O
VEf
vY],f
+ Fib(2) -"-''eft " Dk~
. z'(l'2)
(23c)
where Zi(1,2) = ff1 { Z i ( 1 , 1 )-S °D[ Jr J*(2)dF} ' In general, the source integral Q}") for the ruth iteration (m > 2) can be, therefore, expressed in the following form: Q}m) : f s(m)(b~(O) df~ = Zi(m, O) d _
uEf
4 D ~ I ) Zi(m - 1, 1)
There are no unknown variables in eqns (24a)-(24c). These expressions are applicable to both two- and three-dimensional problems•
5 HIGHERORDER FUNDAMENTALSOLUTIONS FOR TWO-DIMENSIONALPROBLEMS Higher order fundamental solutions are required for the multiple reciprocity computations. In the case of a twodimensional, one-group neutron diffusion problem, the fundamental solution of zero order is given by
(~;(o) = Ko(kr)/(2rc)
(25)
where the quantity k is given as k 2 = Za/D. The Lth order fundamental solution to eqn (15b) has the form
dp~(L) = Az(kr)LK£(kr)
(26)
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems where the coefficients AL are obtained recurrently from the relationship
AL = AL_I/(2Lk 2)
for L = 1,2, 3,...
1.2
~
10 °
/
(27)
The derivation of eqns (26) and (27) are found in another paper by the authors. 1° It is noted that the function o~(L-)has no singularity for L > 1, so that the integration of these functions does not require any special technique.
I
43
1.1
lO-Z ~0
v
¢0 :>
'
1.0
eigenvalue
C 0
k e f f = 1.0
I 0 - = .,~
g
O a
bo / x
%J
10 `3
0.9 ¸ to
6 RESULTS OF A S I M P L E TEST CALCULATION
g
8 LU 0
A simple test calculation has been made. A homogeneous square domain with dimensions of 50crn by 50cm has been chosen for the example. A zero-flux boundary condition is set for two sides and a zerocurrent condition for the other two sides of the outer boundary, as shown in Fig. 1. In this case the analytic flux solution is given by
(~(x,y) cx sin (B~x) . sin (Bc y)
(28)
with the geometrical buckling Baa = 2(1r/100) 2. When the neutron source iteration converges, the neutron diffusion eqn (1) should be equivalent to the Helmholtz equation V2q~ +/~Gq~ = 0
(29)
with the relationship =
D
- B2m
(30)
The quantity B~ is called the material buckling. For convenience one assumes here that the nuclear con-
Y
(cm) 5O
10-*
0.7 0
I I 20 Number
f I 40 of source
I I 60 iterations
~ 80
I
10 " s 100
Fig. 2. Convergence behaviour of source iteration. stants are given as: D=l.0(cm),
~ a - - 0 " 0 1 ( c m -1)
and
v~f = 0.011 9739 (cm -1) so that the geometrical buckling is equal to the material buckling when kerr = 1"0. For the BEM calculation, each side of the square domain is divided into five equal 'linear' boundary elements. The source iteration started with a flat source distribution S0)(r) = 1.0 as an initial estimate. Figure 2 shows the convergence behaviour of the eigenvalue ken' from the first to the 100th iteration, where the higher order fundamental solutions were generated up to the 100th order. The dotted line represents the variation in the eigenvalue deviation at the ruth iteration defined by
~=0 e --
uZf=O.O119739 T.a =0.01 J=0
0.8
cm-'
cm-t
D = 1.0 cm
J= 0
¢~=0
50 (cm)
Fig. 1. Problem specifications for the test calculation.
X
/.(m-l)
(31)
During the early stage of iterations, the calculated eigenvalue rapidly approached the true value (ken' = 1-0): the 7th solution was kerr= 0.99843 and the 8th was kerr = I'00070. This rapid convergence suggests that there is no serious mistake in the present formulation. Note that acceptable results are obtained if one stops the iteration here. However, after the 8th iteration the calculated eigenvalue gradually goes away from the true value. The deviation e has the minimum value e = 1.18 x 10-3 at the 12th iteration, but the deviation increases again as the iteration proceeds. It is interesting to point out that the eigenvalue seems to converge finally to the value of koo = UEf/Ea(= 1"19739). The quantity ko~ is called the 'infinite multiplication factor' for an infinite domain or a domain with no neutron leakage. In fact, the calcu-
44
M. Itagaki, C.A. Brebbia
lated neutron currents along the boundary represent very small values at the 100th iteration.
7 DISCUSSION According to the authors' further investigation, the convergence of the Lth order term in eqn (24a) at the mth source iteration depends on whether the following condition vEf • t-(m-Z ) r / = ~a %ff
L~
1
< 1"0
(32)
is satisfied or not. The calculated eigenvalue approaches the true value efficiently as long as the quantity r/is less than unity. The reason for this can be described as follows. As shown in eqn (24a), the Lth order fundamental solution $i*(L) given by eqn (26), always associates with the factor m-I
V~f
H i)b(m-L) L=1 ~'~eff then one can make another definition of the Lth order solution as
~(L) = JiL(kr)LKL(kr)
(33)
with the relationship of the new coefficients
-
AL-I
Um_t )
(34)
2L. ~efr
In order to investigate the convergence of the newly defined fundamental solutions, one estimates the following quantity
ff,(L)
u~flEa (kr)LKL(kr) = lim TM = lim lr(m_L ) • (kr)L_lKL_l(kr) r-~O~
uEf/E a
2 L-I" (L - 1)!
_
r/= 2L. k ~ -L) " 2L-2" (L - 2)! -
uEf/Ea
L - 1
b (m-L)
L
(36) The value of the function -4L xLKL(x) is sequentially calculated using the recurrence formula
~,+1xL+IKL+I(X) = X2{~L_lXL-IKL_I(X)} • ('4~+l/A~-~)
+ 2I~(~,xLK~(x)} • (,L+I//L-~)
(37)
This recurrence scheme starts with the values of the
functions K0(x) and xKl (x), which are calculated using polynomial approximations for the modified Bessel functions found in the literature) l The polynomial approximations originally contain some round-off errors. It is understandable that the errors are magnified if the condition r} < 1.0 is not satisfied during the recurrence calculations. This is supposedly the reason why the source iteration does not converge to the true value. Why did the method seemed to converge finally to the false eigenvalue
k(m-L) VEf eft ~ k ~ = Ea in the test calculation? Because of the round-off error, the source integral QI m) was not calculated properly. The calculated eigenvalue has no physical meaning any more and gradually increased. Finally it had a value very close to k~ and the quantity r/ happened to become less than unity (r1 ~ (L - 1)/L < 1.0). Consequently the iteration apparently converged because of no addition of round-off error accumulation. Nowak 12 pointed out that a similar phenomenon may occur in his M R M computations for solving Poisson-type problems. He also recommended that one should scale and reformulate the problem, i.e. all dimensions should be divided by the maximum dimension. However, such a scaling cannot be applied to the authors' modified Helmholtz problems. One should seek another way of removing the numerical instability.
8 CONCLUSIONS The multiple reciprocity method has been applied to one-group fission neutron source iteration problems. The domain integral related to the fission neutron source in the ruth source iteration is transformed into a series of m - 1 boundary integrals. These integrals require the zero order to the r n - 1 order fundamental solutions and the boundary neutron fluxes and currents calculated in the previous m - 1 source iterations. The boundary-only integral expression proposed here for the iterative fission source is applicable to both twoand three-dimensional criticality problems. Any information in a domain is not required in this formulation, so the advantages of the BEM can be perfectly preserved. In addition, there is no need to assume any approximated source functions like ones introduced in the dual reciprocity method. Results obtained in a simple test calculation indicate that the source iteration based on the present techniqu e converges rapidly to the desired eigenvalue as long as the condition given by eqn (32), r / < 1'0, is satisfied. Unless the condition is satisfied, a numerical error is accumulated as the iteration proceeds. In order to avoid such situations, further
Multiple reciprocity boundary element formulation for one-group fission neutron source iteration problems consideration should be incorporated into the present technique.
ACKNOWLEDGMENT The authors wish to express their gratitude to Dr. A.J. Nowak o f Silesian Technical University, Poland, for his valuable comments on the work and manuscript.
REFERENCES 1. Brebbia, C.A. The Boundary Element Method for Engineers, Pentech Press, London, 1978. 2. Itagaki, M. Boundary element methods applied to twodimensional neutron diffusion problems, Journal of Nuclear Science and Technology, 1985, 22(6), 565. 3. Itagaki, M. Boundary element techniques for two-dimensional nuclear reactor calculations, Engineering Analysis, 1987, 4(4), 190. 4. Itagaki, M. & Brebbia, C.A. Boundary element method applied to neutron diffusion problems, Proc. 10th Int.
5.
6.
7.
8.
9. 10.
11. 12.
45
Conf. BEM, Southampton University, UK, SpringerVerlag, 1988. Itagaki, M. & Brebbia, C.A. Boundary element method applied to neutron diffusion problems of a thin layer sandwiched between two zones, Proc. 12th Int. Conf. BEM, Hokkaido University, Japan, Springer-Verlag, 1990. Itagaki, M. & Brcbbia, C.A. Space-dependent core/ reflector boundary conditions generated by the boundary clement method for pressurized water reactors, Nucl. Sci. Eng., 1991, 107, 246. Itagaki, M. & Brcbbia, C.A. Boundary clement formulation of fission neutron source problems using only boundary integrals, Engineering Analysis with Boundary Elements, 1991, 8(5), 239. Nowak, A.J. & Brebbia, C.A. The multiple-reciprocity method. A new approach for transforming BEM domain integral to the boundary, Engineering Analysis with Boundary Elements, 1989, 6(3), 164-168. Wachspress, E.L. Iterative Solution of Elliptic Systems, Prentice-Hall, Inc., 1966. Itagaki, M. & Brebbia, C.A. Generation of higher order fundamental solutions to the two-dimensional modified Helmholtz equation, Engineering Analysis with Boundary Elements, 1993, 11(1), 87-90. Abramowitz, M. & Stcgun, I.A. Handbook of Mathematical Functions, Dover, 1965. Nowak, A.J. Private Communication, 1991.