- Email: [email protected]
On the calculation of flux in slightly subcritical reactors with external neutron sources
Prottre~s in '~'uclear Fnerqy, Vol. 23. No. 3. pp. 1g1-189. 1990. Printed in Great Britain. All rights reserved,
i
0149~1970 90 S 0 , 0 1 ) - 5 0 1991 Pergamon Press pie
ON THE CALCULATION OF FLUX IN SLIGHTLY SUBCRITICAL REACTORS WITH EXTERNAL NEUTRON SOURCES J. PLANCHARD Electricit6 de France, Etudes et Recherches. 1, Avenue du G~n6ral de Gaulle. 92141 Clamart, France (Receired 21 March 1990)
Abstract--The behaviour of the flux inside a subcritical reactor in the presence of external neutron sources is examined. It is shown, in particular, that the flux can be approximated by the flux resulting from eigenvalue calculation as the reactor approaches its critical state. A method based on the perturbation technique is described. allowing an estimation of spatial effects on the flux by the local sources.
The critical equations are then conventionally:
1. INTRODUCTION It is a well known fact that the determination of the power distribution inside a nuclear reactor core is equivalent to solving an eigenvalue problem; this search process is usually denominated as a critical calculation. This approach presupposes that no source of neutrons other than fissions is present inside the core. In reality, there are always sources in reactors:
AtOl(x) = -div(dt(xJgrad Ot(x))+aa~(x)Ol(x) va :~ {x)Ot (x)+ wr::(x)O2(x) 2 AJ92(x) =- - div(d2(x)grad O:(x)) + a~:(x)O2(x)
--- a ~ t ( x ) ,
(1)
in which: - - t h e natural radioactivity and the residual products of fissions create neutrons, - - p e r m a n e n t radioactive sources (curium, americium, etc.) are placed intentionally in cores to trigger the chain reaction during start-up operations, - - t h e cosmic rays and the interaction of ~ and y-rays with oxygen and hydrogen also produce neutrons, etc.
--d~, ao., are respectively the diffusion coefficient and the macroscopic total absorption cross-section of group k, - - v a : , is the fission cross-section of group k, multiplied by the average number of neutrons emitted by a fission, m a R is the slowing-down cross-section "of the fast group, - - 2 is the eigenvalue (or multiplication factor).
The intensity of the different sources is of course very low, so they are neglected in calculations. The purpose of this paper is to show, in a mathematical way, that the neutron flux resulting from sources is virtually proportional to the fundamental mode of critical calculation; the proportionality factor is then determined by the source intensity and the antireactivity of the core. A method is also proposed to investigate the spatial effects of point sources (for start-up operating) on the fundamental mode.
Boundary conditions on the frontier of the core must be added to equations (1); to fix the ideas, zero flux conditions are chosen (zero current or albedo conditions can be also considered). Equations (1) can be rewritten in a more condensed form: AO = FO/>.,
(2)
in which:
2. NOTATIONS AND SOME REMINDERS
For the sake of simplicity and to clarify the ideas, only two energy groups are considered; the cases with larger numbers of groups may be extended with no difficulty.
For the chosen boundary conditions, the differential operator A is invertible and equation (2) is equivalent to 181
182
J. PLANCHARD
A - IFO = 20. It is shown (see Habetler and Martino, 1961) that A - tF has eigenvalues in numerable infinity 2,, (n = 0 , 1, 2, etc.) which converge towards zero as n tends to infinity; this stems from the fact that A - , is a compact operator for the L ~ - - n o r m . The n u m b e r zero is an eigenvalue only when F has a kernel not reduced to
3. BEHAVlOUR OF A SUBCRITICAL REACTOR IN THE PRESENCE OF NEUTRON SOURCES
3.1. Let us consider a subcritical reactor that is characterized by the operators A and F as defined in the previous section. We thus have
{0}. In addition, according to the K r e i n - R u t m a n theorem (Krein and Rutman, 1962), the first eigenvalue 2 o is positive, simple and equal to the radius spectral of A - iF:
2o=r(A-'F)-max n
12.{,
(3)
r(A - I F ) < 1.
(7)
Let
/s,(x)'~ s'x'=~,s2(x) ) .
,
be a neutron source inside f2. The flux ~ resulting from the source must satisfy the equation
n = 0 , 1, 2, etc.,
(,4-F) O=s,
and the associated eigenvector (referred as the fundamental mode) is positive throughout the domain f~ occupied by the reactor, i.e.:
or equivalently ( I - - A - 1 F ) 0 = A - is; it follows from equation (7) that this last equation is invertible since the spectrum ofA - t F d o e s not contain the number 1. Moreover, we have:
¢O,k(X)>0, X~f2 and x ~ r , F = b o u n d a r y off~, (k is group index).
(I--A-1F)-t=I+A-tF+(A-tF)E+
Moreover, Oo is the only positive eigenvector of A - 1F; it therefore corresponds to the physical solution of the problem and the other eigenmodes are referred as the 'harmonics' in neutronic literature. The adjoint problem is associated with equation (2): A*O* = F*@*/2*,
(4)
where A*, F* are respectively the adjoint operators of A and F. In a similar way, (4) has a numerable sequence of eigenvalues 2* n = 0, 1, 2, etc., and
2~ =r(A*- XF*)=).o,
(5)
associated with a positive eigenvector ~ (which is the only positive eigenvector of equation (4)); Oo is called the adjoint flux (while Oo is the direct flux). There is an orthogonality relation between the O. and O~ ( F ¢ . , ¢*)--
1 if r e = n ; 0 if m:~n;
+(A - ~ F ) " + . . . .
in equation (6), the parentheses define the standard inner product on the functional space LZ(f~):
(u'v)= ~nu'(x)6'(x)dx'i=, ~ where the bar denotes the complex conjugate. The associated norm is Iloll = Io, All these properties are true whatever the n u m b e r of groups. The reactor is said to be critical when 2 o = I (it is subcritical if 2 o < I and supercritical if 2 0 > 1).
"'" (9)
in which the series converge* because the spectral radius of A - t F is smaller than unity. In addition, each Ak- t is a strictly positive operator: in fact it is an integral operator with a strictly positive kernel (see Protter and Weinberger, 1967), and the same is consequently true for A -~. But A-~F is only a non-negative operator, because of the presence of a reflector. In general, the source s is physically nonnegative, but A - Xs is strictly positive; it thus results from equation (9) that the flux 0 is also strictly positive inside the domain f~ occupied by the reactor, Another consequence of the inequality (7) is that the flux ~ can be calculated in an iterative manner; starting from an arbitrary ~OI°~, one defines the sequence 0 ~"~ by:
A~k~"+l~=F~"~+ s, (6)
(8)
(10)
and ~,1"~--*~kas m--* + ~ , since it is well known fact that the iterative process is convergent if, and only if, r(A-IF)< 1 (see Varga, 1982). However, the convergence is very slow when the spectral radius is close to one (the convergence can be markedly accelerated by using the Chebychev polynomial technique). Note that, when s = 0 (no source), the algorithm is equivalent to a critical computation: in this case, ~l,.~ tends to ¢o in direction.
* The convergence is not uniform, except if the norm of A - 1F is smaller than one.
Calculation of flux in subcritical reactors 3.2. About the expression of the reactivity The reactivity is defined by the ratio
)-o- I P=
20
(ll)
Let s(x) be a source and ~k the corresponding flux. Taking the inner product of equation (8) with the adjoint flux, we get:
t83
calculation of the reactivity, since it necessitates knowledge of both the adjoint flux and the solution ¢, o f ( A - F ) f f = s . Besides, i f ~ is known, )-0 and p are automatically determined. The main point of interest is to lead to a physical interpretation of the adjoint flux ~ . Let ~0_ t be the solution of ( A - F ) ~k=s with sk(x) = - - 1 , k = 1 and 2. Then p is written: 2
E I.
( ( A - F)ff, (l)J)= (s, (I)~)= (O, ( A - F)*(/)~),
k=t
P=
(qJ_~, Y * ~ )
Let us set dpk(x ) =
whence another expression of the reactivity:
(s, o*0) P=
(0, F*(b~)"
(12)
This formula was derived by Greenspan (1975).
a)j.dx) dx (ff _ ,, F*(I)~ )'
and renormalize (I)~ such that (~h_ t, F*(I)J) = 1. Then we have: P= Z
dpk(x)"
(14)
k=l
3.2.1. Consequences of relation (12). (a) Relation (12) remains valid even if the function s(x) takes both positive and negative values, in which case ¢, is generally no longer positive. When the reactor is critical, p is zero and (8) has a solution if, and only if, the source s is orthogonal to the adjoint flux (I)~, due to the Fredholm Alternative Theorem (see for instance Taylor, 1961). The relation (12) is still true when the reactor is supercritical, providing, however, that the number 1 is not an eigenvalue of A - ~F (according to the alternative theorem); to be more accurate, if 2, = 1 for a certain index n, (s, (I)*) must be equal to zero and qJ is then determined to within an additive term ~(1), (in which the constant ~ is arbitrary), but p is uniquely determined by (12) since, whatever the value of ~, we have:
-(s, ,~'g) (O+xO., F*O'~)
-(s, ,v~) (~k, F*O~) - p '
owing to the orthogonality relation (6). In the case of supercritical reactors, p is necessarily positive and if, s is non-negative, we have:
(s, ~ )
gk, F * ( 1 ) ~ ) = - - - < 0 , P
(13)
from the positivity of (I)~, Since F*(I)~ is non-negative, the inequality (13) shows that ff must be negative in some regions of the reactor and therefore cannot correspond to a physical solution of (8); that means there is no steady solution: in fact. the chain reaction starts racing, a well-known phenomenon in supercritical reactors. (bl Relation (12) is of no practical interest in the
(l)g.k(x) thus 'measures' the contribution to reactivity of the neutrons of group k located at point x in the reactor. Note that qJ_ t is negative subcritical reactor and leads to change the sign of adjoint flux in this case.
the the for the
4. THE CASE OF A SLIGHTLY SUBCRITICAL REACTOR IN THE PRESENCE OF SOURCES WITH LOW INTENSITY 4.1. It was seen in the previous section how a steady flux ~0 exists in the presence of neutron sources if, and only if, the reactor is subcritical. A question immediately arises: what happens to ~, when the system approaches criticality? In this situation it has already been observed that the calculation of the flux becomes difficult because the convergence of the algorithm is extremely low; this fact is confirmed by numerical experiments (see for instance West, 19871. Moreover. these computations show as it is seen on Fig. 1 that tends to be proportional to the flux ~o obtained by a standard critical calculation (Palmiotti and Salvatores, 1982). The proportionality factor z~ obviously depends on the source intensity and antireactivity. The purpose of this section is to study the behaviour of :~ with respect to these parameters. Let a critical reactor be characterized by the operators A and F; we have then the criticality relations: ( A - F)Oo =0, ( A * - F * ) ~ * = O , ).o = r(A - 1F) = 1, (:I)o , (I)~ are positive.
(161
Now the same reactor is considered in a slightly
184
J. PLANCHARD
subcritical state. It is defined by the operators A and F , = r / F , in which r/is a real parameter such that: 0 < ~ / < 1, r/close to 1.
(16)
We suppose, in order to simplify the presentation, that we are able to act uniformily on the fuel enrichment of the reactor. Clearly we have: r(A - I F~)=r/ < l
in which ,ui=,~.[ ~ satisfies tSLjl> 1 for j > O , since the reactor is subcritical. It is then clear that
~,=(s.
~)~o
+ O(l -r/).
(18)
in which the term O ( 1 - r / ) is the order of the antireactivity 1 - q.'t Thus the flux ft, corresponding to the source s T, is close to the flux obtained by critical calculation, if s, is of the order of the antireactivity.
and
(*)
[/(A - G ) - ' [I---,+ :c as r/~ 1.
Let ~b~(x) be the solution of (A - r/F)0,(x) = s~(x) -= (1 -
r/)s(x).
R e m a r k 4. I. Assuming that, instead of s T= ( I - r/)s, s T= ( 1 - r/Vs, then one has the following properties:
- - i f p > 1, then ¢,--*0 as r/--* 1, in which s(x) is non-negative: s, is therefore a 'small" source whose amplitude is of the order of the antireactivity 1 - r / (s describes the spatial shape of the source and 1 - r / i s a rescaling factor). Since s Ttends to zero as r/---,l, t p , = ( A - F ~ ) - l s , may be expected to converge toward something. It will be effectively shown that ~ , ~ ~ o when i - r / i s small enough. Let us write 0r as
i=0
where the ~ ' s are the eigenvectors of A - ~ F (one implicitly assumes that they form a complete basis). Then: (A - F~)~Oq= Z :~dA - r/F)O, = ~ ~,(/.t, - r/)~, = sT, i
i
where/t~ = 2 ( ~, i = 0, 1, 2, etc. By taking the inner product of this relation with the adjoint harmonics:
- - i f p < 1, then
110.11-'+2
(more exactly Or = 0(( l - r/)v - t )).
4.2.
In reality, the sources are given and fixed; they are of small amplitude and one cannot act on them. The problem which arises is as follows: how can the reactor power be brought to a given level P? The reactor power is an integral of the type:
where ~t is a proportionality constant. The desired power P is linked to the corresponding antireactivity I - r / ( i . e . necessary to obtain the power P). Indeed, according to equation (18), we have: P~-(s, ¢b'~)(h, ~o).
:q(lh - r/) (F**, (b* ) = ( 1 - r/) (s, dO*),
Let us normalize *o such that (h, ~ o ) = 1. Therefore
i
whence, from the orthogonality condition (6):
P = (s, *J).
However, s(x) is unknown; what is known (at least theoretically) is the real source: s,,~2= s,. Then
=i = (s, q)j*) 1 - r/, J = 0, 1, 2, etc.
~tj
- -
(19)
r~
In particular, since/~o = 1,
P=(s.,~)=
% = (s, ~ ) .
Thus:
l-r/
'
or
~b,=(s, ~J)q~o+ i l - r / (s, q~*)q~j, ~=~ # ~ - q
(17)
(sf.~t. ~ ) 1 -,1 = - -
P
* The norm of a bounded operator B is defined by i
t:, = H_w~SUp 0 ,,B~ Hwhere'l~ll=(k~ f n t ~ k ( x ) ' 2 d x ) .
(20)
This is the necessary antireactivity for the supplied power to be equal to P. Obviously, this antireactivity is practically equal to zero at nominal power. t O(x) means a quantity of the order of :~.
Calculation of flux in subcritical reactors
Remark 4.2. In a reactor, the adequate control rod motions enable the desired power P to be reached. We then have:
From equation (22) A'(z)=Aa~ and using equation (6). the antireactivity may be written: ~(~)= ~(Aa,¢ o, ¢o)+O(~-). *
=
where p denotes the control parameter; p is chosen so that the antireactivity satisfies: (21 )
in which It(p) is the inverse of the spectral radius of A(p)-~F. This point is examined in the following subsections in more detail.
4.3. In practice, the power of the reactor is controlled by control rod displacements or boron concentration (for PWR), leading to changes in the absorption crosssections. For instance, in the event of variation of boron concentration, we set
A(~)= A + ~ Aa,(x),
1 =
_
~(e,) -
(s, 0~)(1)o
e(ha~¢o, q)~)
+ 0(1),
(24)
where O(1) contains the contribution of all the harmonics q)~. It is clear from equation (23) and equation (24), that the power of the reactor, for small e, is inversely proportional to the antireactivity.
Remark 4.4. Suppose that the variation in antireactivity is due to control rod motion and denote by S Othe surface of the bottom of this control rod at the position for which the reactor is critical. ~ being the displacement, with respect to S o, of the control rod, it is not difficult to see that
(22)
with Aaa(x)= diag(Aaak(x)),Aaakbeing the absorption cross-section variation of group k corresponding to a reference boron concentration; the parameter measures the change of this boron concentration. The reactor defined by A and F is assumed to be critical when ~=0, and we denote by ;.(~), ~(~) the fundamental eigensolution corresponding to e. Obviously ~o = • (0) and 20 =).(0). The antireactivity for e # 0 is ~(~) -= _ p ( ~ )
(23)
If q~(~) is the flux in the presence of a source six), clearly:
(A(p)- F)O(p) s:~,t
1 -- ~t(p) = (s~,~t, O ~ ) / P
185
A'(~)l,=o=a~r° ,Sso(X), where Aao is the control rod absorption cross-section and 3so(X) is the Dirac distribution whose support is the surface S o. Then formulae (23) and (24) are still valid, but the factor (Aaa¢ 0, ~ ) must be replaced by
Z Aa°kfso¢o
dSo.
- ;.(~) ).(~) _
and using a limited expansion of).(e): -~).'(0) ~(~) - - 4- O ( ~ 2 ) , 1 + t).'(0)
where 2'(0) is the derivative of 2(~) with respect to ~ at ~=0. In order to estimate Z(0) the critical calculation equation is differentiated with respect to parameter so that the derivative ¢'(e) satisfies:
( A(e) - ).~))dO'(~)= - ( A'(~) + ~)"(~)F)~(~); because the derivative ~'(~) exists, the right-side member of the above equality must be orthogonal to the adjoint flux ¢*(e), implying
,:.'(~)
(,4 '(~)~(~), ~*(~))
/.(~)2
(F~(~), ~*(~))
5. COMPUTATION OF THE FLUX BY A PERTURBATION METHOD 5.1. It is always assumed that the reactor defined by A, F is critical and we are interested in the calculation o f ~ , the solution of ( A - (1 - ~ ) F ) 0 , =
es0( -= s),
(25)
where eso is a small non-negative source and e is the antireactivity of the reactor characterized by A and (1 - t ) F (here 1 - e plays the role of q in Section 4). It was previously observed that the direct computation of W, is difficult when e is too small, because of the wrong convergence of the iterative process (I0), but it was seen that 0, = (So, ~ ) ¢ o + O(~). However one may be interested by the term O(e); for example, one wants to study the disturbances of the
186
J.
PLANCHARD
2(~)
50
,
25(
,
:
__
o
" ~
3(~)
Criticalcalculauon
),
--,
Flux with source
\\ \ '\ ~(t~). --
I 4fXXX) --
_,
:i
//i,I
i
),
Radius (cm) 0
5(}
I~X)
150
200,
")50.
3(~)
Fig. I. Comparison of the flux obtained by critical calculation and the flux resulting from a source in a slightly subcritical reactor (from Palmiotti and Salvatores, 1982).
flux caused by a local source and these spatial effects are just described by O(~). This term is expressed with the harmonics q)i, but these functions may be complex (except for one-group problems) and they are not easy to obtain through numerical methods. To overcome this difficulty, the source and eF are considered as small perturbations of the critical reactor and °d~ is written as an asymptotic expansion of e, namely:
0,(x) = :%% (x) -t- ~1~/1 i X ) +
8202
(X) Jr-'''
-1" EPI~p(X)-t-
....
(26) in which ~p~, ~2, etc., have to be determined. For this purpose, the expansion equation (26) is set into equation (25) and the different terms in e ~' are identified. Doing that, we implicitly assume that equation (26) is valid whatever the values of e. The following equalities are successively obtained:
{
~o(A - F)~ o = O, (A -
F)O,
=
- :%F~ o + so ,
(A -- K)¢, = - F O r , (A - F ) ¢ 3 = - F ¢ , ,
etc.
(27)
It is noted that the first equation (27)implies that the leading term of the asymptotic expansion necessarily contains D o. Because A - F is singular, the right member of equation (27)., must be orthogonal to ~ and then (So, 0 ~ ) ~0 -
(F~o
' q:)~) •
which we already knew! Denoting by s' I ( x ) the right-side member of equation (27)2 corresponding to this value of :%, we search the solution 01,0 of (A - F I e 1.o = s', ,
(28)
with the condition: (~,.o, ~o) =0'
(29)
The Fredholm Alternative Theorem guarantees the existence of a unique solution for equation (28) and equation (29). the computation of ~,.o can be carried out by the Gandini-Usachev iterative process (Gandini, 1967; Usachev, 1964) with acceleration by a Chebychev polynomial technique (see Gomit and Planchard, 1985). So the general solution for equation (27) l is obviously ofthe form tk t ( x ) = ~ t . o ( X ) + : t t r b o ( X )
Calculation of flux in subcritical reactors in which the coefficient ~t has to be defined. Because it is intended that ~bz exists, the right-side of equation [27)., must be necessarily orthogonal to O~, leading to (Fff t.o, O~) :¢~ = (FO o. O~) ' In a similar manner, ~z, ~a, etc. can be calculated and the process may be indefinitely pursued.
Remark 5.1.
187
0 ( ~ - ~); more precisely there exists a certain positive constant CO such that:
c°
when ~ is small enough. (See the appendix for the value of Co.) Now, let us consider the exact solution "P~ of equation (25) and suppose the different terms ~'t, ~.,, etc. of equation (26) have been calculated. Taking the first terms up to zz~bz (to fix the ideas), we then define:
It is possible to write the source in the form
s(x) = eSo(X) + e2s I (x) + • - • + eksk_ ~(x) + . . . ; then the different terms sk would be present in the right-sides of equations (27). In particular, if So=0, one has ~ o = 0 and then e t ( x ) = ~ l O o ( X ) with (s,, ~t -
O~)
(FOo ' @~)
Consequently, the flux is of the order of e : ~ , ( x ) = ~ ~¢o(X) + 0(~2). More generally, ifs o = s i = • • - = s~_ ~= 0 and s~ ~: O, then ~ ( x ) = ~£~@o(X) + 0(~:~+ ~) with (s~,
~u -
O~)
(FO o , O*)
Remark 5.2. If s(x)=So(X)+tst(x)+ . . . . it is necessary to set $ , ( x ) = ~ o , - t O o ( X ) + $l(x)+e~b2(x)+ . . . . and then (s o,
~o =
O~)
(F¢o '
¢~)-
Remark 5.3. If the perturbation is caused by e Aa, (instead ofeF), all the calculations remain unchanged except that F must be replaced by Aa~ in the expressions of coefficients ~o, ~t, ~2, etc.
5.2. A posteriori justification of the asymptotic expansion of ~ , The above calculations are purely formal since it is not sure that the series ~o¢o+e~b t +e2~b2 + . . . converges. However, this technique is currently used by physicists in q u a n t u m theory even in cases where the series does not converge, but they are only interested in the sum of the first terms (see F e y n m a n n , 1965; Rellich, 1969). Our objective is now to give a justification of this point. Firstly, it is observed that I[( A - ( 1 - e ) F ) - ' l l =
(301
£
0~ = ¢ , -
(~o¢o + ~ t + ~0,)-
Clearly we have, after some handling: (A - (1 - e.)F)0~ = -
PFO,,
and it results directly from equation (30) that
I[o, ll <-Co fll More generally, + ~kOk), then
if
211.
0~=t,b~-(:%Oo+:q~k I + - "
I1¢11-
(3 1
Thus, the flux ~ can be approached by a sum of the first terms of the asymptotic expansion, even when the series is not convergent; the inequality equation (31) gives a bound of the error caused by this approximation: it depends on the last term ~,~ calculated in this way. This a posteriori justification is widespread in perturbation theory methods and many other applications can be seen in Lions (1973). The method presented in this section necessitates functions Ok(x) to be determined and only a few functions are required if ~ is small enough; because they are corrective terms, it is not necessary to compute them with a high accuracy and only coarse discretisation meshes can be used.
6. CONCLUSION Thus, the problem of neutron diffusion in the presence of sources always admits a steady positive solution in the case of subcritical reactors. In the event of low intensity of sources and if the reactor is close to criticality, this steady solution can be approximately achieved by standard critical eigenvalue calculation; the flux obtained in this way is determined to within a multiplicative factor, and it was shown that this factor is fixed by the source intensity and the antireactivity. A method, based on asymptotic expansion of the flux, is proposed for the study of the spatial effects of sources on the fundamental eigenmode obtained by critical calculation. We hope in a future paper (by J. P. West and the author) to present numerical results.
188
J. PLANCHARD
Acknowledgments--The author thanks M. Salvatores and G. Palmiotti who have kindly communicated to him the curve of Fig. 1.
where Q = A - 'F. Obviously, the Oi are the eigenfunctions of Q, and y~ can be expressed as Y,= ~
i=O
REFERENCES
Blanchon F., Ha-Duong T. and Planchard J. (1988) Numerical methods for solving the reactor kinetic equations. Proy. Nucl. Energy 22, 2, 173-180. Feynmann R. (1965) La nature de la Physique. Editions du Seuil, Paris (French translation). Gandini A. (1967) A generalized perturbation method for bilinear functionals of the real and adjoint fluxes. J. Nucl. Energy 21, 735. Gomit J. M. and Planchard J. 0985) Perturbations conservant la criticit6 et fonction Importance. Bulletin des Etudes et Recherches, Electricit6 de France, S6rie A, 1, 17-22. Greenspan E. (1975) A source multiplication reactivity. Nucl. Sci. Engng 55, 103-105. Habetler G. J. and Martino M. A. (1961) Existence theorems and spectral theory for the multigroup diffusion model. Proc. Syrup. Appl. Math., I 1, Nuclear Reactor Theory, Am. Math. Soc., Providence (R.I), pp. 127-139. Krein M. G. and Rutman M. A. (1962) Linear operators leaving invariant a cone in a Banach space. Am. Math. Soc., Translation No. I0, 199-325. Lions J. L. (1973) Perturbations singuli~res dans les probldmes aux limites et en contr61e optimal. Lecture Notes in Mathematics 323, Springer Vertag, Berlin. Palmiotti G. and Salvatores, M. (1982) (private communication). Planchard J. and West J. P. (in preparation). Protter M. H. and Weinberger H. F. (1967) Maximum Principles in Differential Equations, Prentice Hall, New York. Rellich F. (1969) Perturbation Theory of Eigenvalue Problems, Gordon and Breach, New York. Riesz F. and Nagy B. (1955) Le?ons d'analysefonctionnelle, Acad~mie des Sciences de Hongrie, 3rd Ed., Budapest. English translation: Functional Analysis (1965) Ungar Publishing Company, New York. Taylor A. E. (1961) Functional Analysis, Wiley, New York. Usachev L. N. (1964) Perturbation theory for the breeding ratio and for other number ratios pertaining to various reactor processes. J. Nucl. Energy, 18, 571. Varga R. S. (1962) Matrix lterative Analysis, Prentice Hall, New York. West J. P. Performance des Modules de Diffusion des Syst~mes CARDIF (code SNAP-3D et CCRR--Calculs Homoo~nes et InhomogOnes d 2 Dimensions, Electricit~ de France, Internal Report No. HT-11/8718. APPENDIX
Estimate of the constant CO Our purpose is to obtain an upper bound of IIEA-(1WJ-'11 This norm is defined by: ~,~ Ily, lllllfll, where y, is the solution of [A -- (1 -- 8)F-Jy~ = f or equivalently:
[1 - ( 1 - e ) Q ] y , = A - 'f=-9,
(A1)
~i(~i
with
(g, uJ) ~o =
(~, u*) , cq=
1 -- (1 - ~))'i
, i>0
where u~' are the eigenfunctions of Q*, with the orthonormalization condition (~i, u*)=(}ii- Because • ~, u~' form a biorthonormal sequence, the numbers (9, u~*) are the components of g with respect to this sequence. Consequently, a positive constant C exists such that (see Riesz and Nagy, 1955):
.?)l
"< cIIgll,
where C does not depend on 9Clearly
Ily
'
max i>o
II-l(1-e)2i
gll
The numerical calculations show that practically the second eigenvalue ;-t is positive and close to 20( = 1), so that one gets
Ily, II-< IIA-' IIIISlI. Then
Co =
C
JIA-'II- CllA-'ll ,
The difficulty is that we are unable to exhibit the value of C, unless the ~ and u~' are explicitly known. However, that can be achieved if A and F are symmetric operators: this is the case in one-group problems. In this situation, A being symmetric and positive definite, the square root A r of A(A = (A½)2) can be defined (see Riesz and Nagy, 1955) and equation (A1) is equivalent to: [1 - ( I - e ) 0 ] L = f
(A2)
where 1
!
1
1
Q = A -~FA -~, A - ~ = (A ~) - 1
-'-7 is symmetric positive definite and has the same eigenvalues as Q. It ensues that the norm of the inverse RE of I - ( 1 - ~ ) Q is:
Calculation of flux in subcritical reactors 1
where d, is the distance from 1 to the spectrum of (1-~)0, i.e.
d,=ll -(1-
189
and the minimum is reached for a certain function u 0 (u o is the eigenfunction associated with ~,). From the maximum principle (Protter and Weinberger, 1967), uo is positive inside fL We have ~, = ( ( A - F + eF)u o, Uo ) + e( Fuo , Uo )
~)X,I = e.
This is due to the fact that the norm of a symmetric operator is equal to the modulus of its greatest eigenvalue. Thus
II;:11
I1 1.
on the other hand, IIfGI_
- 11 II;,11.
IlY
IIA Since the norm of A - r is strictly equal to the square root of the greatest eigenvalue of A - t, or equivalently, of the inverse of the lowest eigenvalue fl of A, we then have l
Ily.ll-< llSII, and consequently Co = 1/ft.
(A3)
For PWR cores, fl is of the order of the absorption cross-section, i.e. 10 -2 (see Blanchon et al., 1988); more precisely
> min ( ( A - F ) u , u ) + e ( F u o, Uo), It,ll = t
and m i n ( ( A - F ) u , u) is zero because the reactor II~t = t
characterized by A, F is critical and A is positive definite. Thus '~, > e( Fu o, %). It is noted that u o implicitly depends on e, but Uo~-~o where ~o is the fundamental mode of the reactor, and then
I[*o[l:
CO= (FOo ' Oo)"
(A4)
The estimates (A3), (A4) are rigorously valid for one-group problems. Now, it is well-known that multigroup systems may be approximated by onegroup equations. Then, these estimates remain valid with F and fl computed for the equivalent one-group problem. It is important to note that the exact value of may be easily obtained by inverse iteration method (a few iterations are sufficient).
/~> min a.(x). Remark
Another estimate of C Ocan be derived. Indeed, let us come back to (At). Clearly
When the perturbation is due to variation of try, eF must be replaced by e A% in equation (Al), so that equation (A4) becomes:
1
N ( A - - f +eF)-tll =-~t E
where ~ is the smallest eigenvalue of A - F + eF: ~t,= rain ( ( A - F + e F ) u u), i1,!1= 1
Co =
ll*oll (Aa.~o, ~ o ) '
(g4 bis)
In (A4 bis), Aa, is the reference absorption crosssection variation for the equivalent one-group problem.
i
0149~1970 90 S 0 , 0 1 ) - 5 0 1991 Pergamon Press pie
ON THE CALCULATION OF FLUX IN SLIGHTLY SUBCRITICAL REACTORS WITH EXTERNAL NEUTRON SOURCES J. PLANCHARD Electricit6 de France, Etudes et Recherches. 1, Avenue du G~n6ral de Gaulle. 92141 Clamart, France (Receired 21 March 1990)
Abstract--The behaviour of the flux inside a subcritical reactor in the presence of external neutron sources is examined. It is shown, in particular, that the flux can be approximated by the flux resulting from eigenvalue calculation as the reactor approaches its critical state. A method based on the perturbation technique is described. allowing an estimation of spatial effects on the flux by the local sources.
The critical equations are then conventionally:
1. INTRODUCTION It is a well known fact that the determination of the power distribution inside a nuclear reactor core is equivalent to solving an eigenvalue problem; this search process is usually denominated as a critical calculation. This approach presupposes that no source of neutrons other than fissions is present inside the core. In reality, there are always sources in reactors:
AtOl(x) = -div(dt(xJgrad Ot(x))+aa~(x)Ol(x) va :~ {x)Ot (x)+ wr::(x)O2(x) 2 AJ92(x) =- - div(d2(x)grad O:(x)) + a~:(x)O2(x)
--- a ~ t ( x ) ,
(1)
in which: - - t h e natural radioactivity and the residual products of fissions create neutrons, - - p e r m a n e n t radioactive sources (curium, americium, etc.) are placed intentionally in cores to trigger the chain reaction during start-up operations, - - t h e cosmic rays and the interaction of ~ and y-rays with oxygen and hydrogen also produce neutrons, etc.
--d~, ao., are respectively the diffusion coefficient and the macroscopic total absorption cross-section of group k, - - v a : , is the fission cross-section of group k, multiplied by the average number of neutrons emitted by a fission, m a R is the slowing-down cross-section "of the fast group, - - 2 is the eigenvalue (or multiplication factor).
The intensity of the different sources is of course very low, so they are neglected in calculations. The purpose of this paper is to show, in a mathematical way, that the neutron flux resulting from sources is virtually proportional to the fundamental mode of critical calculation; the proportionality factor is then determined by the source intensity and the antireactivity of the core. A method is also proposed to investigate the spatial effects of point sources (for start-up operating) on the fundamental mode.
Boundary conditions on the frontier of the core must be added to equations (1); to fix the ideas, zero flux conditions are chosen (zero current or albedo conditions can be also considered). Equations (1) can be rewritten in a more condensed form: AO = FO/>.,
(2)
in which:
2. NOTATIONS AND SOME REMINDERS
For the sake of simplicity and to clarify the ideas, only two energy groups are considered; the cases with larger numbers of groups may be extended with no difficulty.
For the chosen boundary conditions, the differential operator A is invertible and equation (2) is equivalent to 181
182
J. PLANCHARD
A - IFO = 20. It is shown (see Habetler and Martino, 1961) that A - tF has eigenvalues in numerable infinity 2,, (n = 0 , 1, 2, etc.) which converge towards zero as n tends to infinity; this stems from the fact that A - , is a compact operator for the L ~ - - n o r m . The n u m b e r zero is an eigenvalue only when F has a kernel not reduced to
3. BEHAVlOUR OF A SUBCRITICAL REACTOR IN THE PRESENCE OF NEUTRON SOURCES
3.1. Let us consider a subcritical reactor that is characterized by the operators A and F as defined in the previous section. We thus have
{0}. In addition, according to the K r e i n - R u t m a n theorem (Krein and Rutman, 1962), the first eigenvalue 2 o is positive, simple and equal to the radius spectral of A - iF:
2o=r(A-'F)-max n
12.{,
(3)
r(A - I F ) < 1.
(7)
Let
/s,(x)'~ s'x'=~,s2(x) ) .
,
be a neutron source inside f2. The flux ~ resulting from the source must satisfy the equation
n = 0 , 1, 2, etc.,
(,4-F) O=s,
and the associated eigenvector (referred as the fundamental mode) is positive throughout the domain f~ occupied by the reactor, i.e.:
or equivalently ( I - - A - 1 F ) 0 = A - is; it follows from equation (7) that this last equation is invertible since the spectrum ofA - t F d o e s not contain the number 1. Moreover, we have:
¢O,k(X)>0, X~f2 and x ~ r , F = b o u n d a r y off~, (k is group index).
(I--A-1F)-t=I+A-tF+(A-tF)E+
Moreover, Oo is the only positive eigenvector of A - 1F; it therefore corresponds to the physical solution of the problem and the other eigenmodes are referred as the 'harmonics' in neutronic literature. The adjoint problem is associated with equation (2): A*O* = F*@*/2*,
(4)
where A*, F* are respectively the adjoint operators of A and F. In a similar way, (4) has a numerable sequence of eigenvalues 2* n = 0, 1, 2, etc., and
2~ =r(A*- XF*)=).o,
(5)
associated with a positive eigenvector ~ (which is the only positive eigenvector of equation (4)); Oo is called the adjoint flux (while Oo is the direct flux). There is an orthogonality relation between the O. and O~ ( F ¢ . , ¢*)--
1 if r e = n ; 0 if m:~n;
+(A - ~ F ) " + . . . .
in equation (6), the parentheses define the standard inner product on the functional space LZ(f~):
(u'v)= ~nu'(x)6'(x)dx'i=, ~ where the bar denotes the complex conjugate. The associated norm is Iloll = Io, All these properties are true whatever the n u m b e r of groups. The reactor is said to be critical when 2 o = I (it is subcritical if 2 o < I and supercritical if 2 0 > 1).
"'" (9)
in which the series converge* because the spectral radius of A - t F is smaller than unity. In addition, each Ak- t is a strictly positive operator: in fact it is an integral operator with a strictly positive kernel (see Protter and Weinberger, 1967), and the same is consequently true for A -~. But A-~F is only a non-negative operator, because of the presence of a reflector. In general, the source s is physically nonnegative, but A - Xs is strictly positive; it thus results from equation (9) that the flux 0 is also strictly positive inside the domain f~ occupied by the reactor, Another consequence of the inequality (7) is that the flux ~ can be calculated in an iterative manner; starting from an arbitrary ~OI°~, one defines the sequence 0 ~"~ by:
A~k~"+l~=F~"~+ s, (6)
(8)
(10)
and ~,1"~--*~kas m--* + ~ , since it is well known fact that the iterative process is convergent if, and only if, r(A-IF)< 1 (see Varga, 1982). However, the convergence is very slow when the spectral radius is close to one (the convergence can be markedly accelerated by using the Chebychev polynomial technique). Note that, when s = 0 (no source), the algorithm is equivalent to a critical computation: in this case, ~l,.~ tends to ¢o in direction.
* The convergence is not uniform, except if the norm of A - 1F is smaller than one.
Calculation of flux in subcritical reactors 3.2. About the expression of the reactivity The reactivity is defined by the ratio
)-o- I P=
20
(ll)
Let s(x) be a source and ~k the corresponding flux. Taking the inner product of equation (8) with the adjoint flux, we get:
t83
calculation of the reactivity, since it necessitates knowledge of both the adjoint flux and the solution ¢, o f ( A - F ) f f = s . Besides, i f ~ is known, )-0 and p are automatically determined. The main point of interest is to lead to a physical interpretation of the adjoint flux ~ . Let ~0_ t be the solution of ( A - F ) ~k=s with sk(x) = - - 1 , k = 1 and 2. Then p is written: 2
E I.
( ( A - F)ff, (l)J)= (s, (I)~)= (O, ( A - F)*(/)~),
k=t
P=
(qJ_~, Y * ~ )
Let us set dpk(x ) =
whence another expression of the reactivity:
(s, o*0) P=
(0, F*(b~)"
(12)
This formula was derived by Greenspan (1975).
a)j.dx) dx (ff _ ,, F*(I)~ )'
and renormalize (I)~ such that (~h_ t, F*(I)J) = 1. Then we have: P= Z
dpk(x)"
(14)
k=l
3.2.1. Consequences of relation (12). (a) Relation (12) remains valid even if the function s(x) takes both positive and negative values, in which case ¢, is generally no longer positive. When the reactor is critical, p is zero and (8) has a solution if, and only if, the source s is orthogonal to the adjoint flux (I)~, due to the Fredholm Alternative Theorem (see for instance Taylor, 1961). The relation (12) is still true when the reactor is supercritical, providing, however, that the number 1 is not an eigenvalue of A - ~F (according to the alternative theorem); to be more accurate, if 2, = 1 for a certain index n, (s, (I)*) must be equal to zero and qJ is then determined to within an additive term ~(1), (in which the constant ~ is arbitrary), but p is uniquely determined by (12) since, whatever the value of ~, we have:
-(s, ,~'g) (O+xO., F*O'~)
-(s, ,v~) (~k, F*O~) - p '
owing to the orthogonality relation (6). In the case of supercritical reactors, p is necessarily positive and if, s is non-negative, we have:
(s, ~ )
gk, F * ( 1 ) ~ ) = - - - < 0 , P
(13)
from the positivity of (I)~, Since F*(I)~ is non-negative, the inequality (13) shows that ff must be negative in some regions of the reactor and therefore cannot correspond to a physical solution of (8); that means there is no steady solution: in fact. the chain reaction starts racing, a well-known phenomenon in supercritical reactors. (bl Relation (12) is of no practical interest in the
(l)g.k(x) thus 'measures' the contribution to reactivity of the neutrons of group k located at point x in the reactor. Note that qJ_ t is negative subcritical reactor and leads to change the sign of adjoint flux in this case.
the the for the
4. THE CASE OF A SLIGHTLY SUBCRITICAL REACTOR IN THE PRESENCE OF SOURCES WITH LOW INTENSITY 4.1. It was seen in the previous section how a steady flux ~0 exists in the presence of neutron sources if, and only if, the reactor is subcritical. A question immediately arises: what happens to ~, when the system approaches criticality? In this situation it has already been observed that the calculation of the flux becomes difficult because the convergence of the algorithm is extremely low; this fact is confirmed by numerical experiments (see for instance West, 19871. Moreover. these computations show as it is seen on Fig. 1 that tends to be proportional to the flux ~o obtained by a standard critical calculation (Palmiotti and Salvatores, 1982). The proportionality factor z~ obviously depends on the source intensity and antireactivity. The purpose of this section is to study the behaviour of :~ with respect to these parameters. Let a critical reactor be characterized by the operators A and F; we have then the criticality relations: ( A - F)Oo =0, ( A * - F * ) ~ * = O , ).o = r(A - 1F) = 1, (:I)o , (I)~ are positive.
(161
Now the same reactor is considered in a slightly
184
J. PLANCHARD
subcritical state. It is defined by the operators A and F , = r / F , in which r/is a real parameter such that: 0 < ~ / < 1, r/close to 1.
(16)
We suppose, in order to simplify the presentation, that we are able to act uniformily on the fuel enrichment of the reactor. Clearly we have: r(A - I F~)=r/ < l
in which ,ui=,~.[ ~ satisfies tSLjl> 1 for j > O , since the reactor is subcritical. It is then clear that
~,=(s.
~)~o
+ O(l -r/).
(18)
in which the term O ( 1 - r / ) is the order of the antireactivity 1 - q.'t Thus the flux ft, corresponding to the source s T, is close to the flux obtained by critical calculation, if s, is of the order of the antireactivity.
and
(*)
[/(A - G ) - ' [I---,+ :c as r/~ 1.
Let ~b~(x) be the solution of (A - r/F)0,(x) = s~(x) -= (1 -
r/)s(x).
R e m a r k 4. I. Assuming that, instead of s T= ( I - r/)s, s T= ( 1 - r/Vs, then one has the following properties:
- - i f p > 1, then ¢,--*0 as r/--* 1, in which s(x) is non-negative: s, is therefore a 'small" source whose amplitude is of the order of the antireactivity 1 - r / (s describes the spatial shape of the source and 1 - r / i s a rescaling factor). Since s Ttends to zero as r/---,l, t p , = ( A - F ~ ) - l s , may be expected to converge toward something. It will be effectively shown that ~ , ~ ~ o when i - r / i s small enough. Let us write 0r as
i=0
where the ~ ' s are the eigenvectors of A - ~ F (one implicitly assumes that they form a complete basis). Then: (A - F~)~Oq= Z :~dA - r/F)O, = ~ ~,(/.t, - r/)~, = sT, i
i
where/t~ = 2 ( ~, i = 0, 1, 2, etc. By taking the inner product of this relation with the adjoint harmonics:
- - i f p < 1, then
110.11-'+2
(more exactly Or = 0(( l - r/)v - t )).
4.2.
In reality, the sources are given and fixed; they are of small amplitude and one cannot act on them. The problem which arises is as follows: how can the reactor power be brought to a given level P? The reactor power is an integral of the type:
where ~t is a proportionality constant. The desired power P is linked to the corresponding antireactivity I - r / ( i . e . necessary to obtain the power P). Indeed, according to equation (18), we have: P~-(s, ¢b'~)(h, ~o).
:q(lh - r/) (F**, (b* ) = ( 1 - r/) (s, dO*),
Let us normalize *o such that (h, ~ o ) = 1. Therefore
i
whence, from the orthogonality condition (6):
P = (s, *J).
However, s(x) is unknown; what is known (at least theoretically) is the real source: s,,~2= s,. Then
=i = (s, q)j*) 1 - r/, J = 0, 1, 2, etc.
~tj
- -
(19)
r~
In particular, since/~o = 1,
P=(s.,~)=
% = (s, ~ ) .
Thus:
l-r/
'
or
~b,=(s, ~J)q~o+ i l - r / (s, q~*)q~j, ~=~ # ~ - q
(17)
(sf.~t. ~ ) 1 -,1 = - -
P
* The norm of a bounded operator B is defined by i
t:, = H_w~SUp 0 ,,B~ Hwhere'l~ll=(k~ f n t ~ k ( x ) ' 2 d x ) .
(20)
This is the necessary antireactivity for the supplied power to be equal to P. Obviously, this antireactivity is practically equal to zero at nominal power. t O(x) means a quantity of the order of :~.
Calculation of flux in subcritical reactors
Remark 4.2. In a reactor, the adequate control rod motions enable the desired power P to be reached. We then have:
From equation (22) A'(z)=Aa~ and using equation (6). the antireactivity may be written: ~(~)= ~(Aa,¢ o, ¢o)+O(~-). *
=
where p denotes the control parameter; p is chosen so that the antireactivity satisfies: (21 )
in which It(p) is the inverse of the spectral radius of A(p)-~F. This point is examined in the following subsections in more detail.
4.3. In practice, the power of the reactor is controlled by control rod displacements or boron concentration (for PWR), leading to changes in the absorption crosssections. For instance, in the event of variation of boron concentration, we set
A(~)= A + ~ Aa,(x),
1 =
_
~(e,) -
(s, 0~)(1)o
e(ha~¢o, q)~)
+ 0(1),
(24)
where O(1) contains the contribution of all the harmonics q)~. It is clear from equation (23) and equation (24), that the power of the reactor, for small e, is inversely proportional to the antireactivity.
Remark 4.4. Suppose that the variation in antireactivity is due to control rod motion and denote by S Othe surface of the bottom of this control rod at the position for which the reactor is critical. ~ being the displacement, with respect to S o, of the control rod, it is not difficult to see that
(22)
with Aaa(x)= diag(Aaak(x)),Aaakbeing the absorption cross-section variation of group k corresponding to a reference boron concentration; the parameter measures the change of this boron concentration. The reactor defined by A and F is assumed to be critical when ~=0, and we denote by ;.(~), ~(~) the fundamental eigensolution corresponding to e. Obviously ~o = • (0) and 20 =).(0). The antireactivity for e # 0 is ~(~) -= _ p ( ~ )
(23)
If q~(~) is the flux in the presence of a source six), clearly:
(A(p)- F)O(p) s:~,t
1 -- ~t(p) = (s~,~t, O ~ ) / P
185
A'(~)l,=o=a~r° ,Sso(X), where Aao is the control rod absorption cross-section and 3so(X) is the Dirac distribution whose support is the surface S o. Then formulae (23) and (24) are still valid, but the factor (Aaa¢ 0, ~ ) must be replaced by
Z Aa°kfso¢o
dSo.
- ;.(~) ).(~) _
and using a limited expansion of).(e): -~).'(0) ~(~) - - 4- O ( ~ 2 ) , 1 + t).'(0)
where 2'(0) is the derivative of 2(~) with respect to ~ at ~=0. In order to estimate Z(0) the critical calculation equation is differentiated with respect to parameter so that the derivative ¢'(e) satisfies:
( A(e) - ).~))dO'(~)= - ( A'(~) + ~)"(~)F)~(~); because the derivative ~'(~) exists, the right-side member of the above equality must be orthogonal to the adjoint flux ¢*(e), implying
,:.'(~)
(,4 '(~)~(~), ~*(~))
/.(~)2
(F~(~), ~*(~))
5. COMPUTATION OF THE FLUX BY A PERTURBATION METHOD 5.1. It is always assumed that the reactor defined by A, F is critical and we are interested in the calculation o f ~ , the solution of ( A - (1 - ~ ) F ) 0 , =
es0( -= s),
(25)
where eso is a small non-negative source and e is the antireactivity of the reactor characterized by A and (1 - t ) F (here 1 - e plays the role of q in Section 4). It was previously observed that the direct computation of W, is difficult when e is too small, because of the wrong convergence of the iterative process (I0), but it was seen that 0, = (So, ~ ) ¢ o + O(~). However one may be interested by the term O(e); for example, one wants to study the disturbances of the
186
J.
PLANCHARD
2(~)
50
,
25(
,
:
__
o
" ~
3(~)
Criticalcalculauon
),
--,
Flux with source
\\ \ '\ ~(t~). --
I 4fXXX) --
_,
:i
//i,I
i
),
Radius (cm) 0
5(}
I~X)
150
200,
")50.
3(~)
Fig. I. Comparison of the flux obtained by critical calculation and the flux resulting from a source in a slightly subcritical reactor (from Palmiotti and Salvatores, 1982).
flux caused by a local source and these spatial effects are just described by O(~). This term is expressed with the harmonics q)i, but these functions may be complex (except for one-group problems) and they are not easy to obtain through numerical methods. To overcome this difficulty, the source and eF are considered as small perturbations of the critical reactor and °d~ is written as an asymptotic expansion of e, namely:
0,(x) = :%% (x) -t- ~1~/1 i X ) +
8202
(X) Jr-'''
-1" EPI~p(X)-t-
....
(26) in which ~p~, ~2, etc., have to be determined. For this purpose, the expansion equation (26) is set into equation (25) and the different terms in e ~' are identified. Doing that, we implicitly assume that equation (26) is valid whatever the values of e. The following equalities are successively obtained:
{
~o(A - F)~ o = O, (A -
F)O,
=
- :%F~ o + so ,
(A -- K)¢, = - F O r , (A - F ) ¢ 3 = - F ¢ , ,
etc.
(27)
It is noted that the first equation (27)implies that the leading term of the asymptotic expansion necessarily contains D o. Because A - F is singular, the right member of equation (27)., must be orthogonal to ~ and then (So, 0 ~ ) ~0 -
(F~o
' q:)~) •
which we already knew! Denoting by s' I ( x ) the right-side member of equation (27)2 corresponding to this value of :%, we search the solution 01,0 of (A - F I e 1.o = s', ,
(28)
with the condition: (~,.o, ~o) =0'
(29)
The Fredholm Alternative Theorem guarantees the existence of a unique solution for equation (28) and equation (29). the computation of ~,.o can be carried out by the Gandini-Usachev iterative process (Gandini, 1967; Usachev, 1964) with acceleration by a Chebychev polynomial technique (see Gomit and Planchard, 1985). So the general solution for equation (27) l is obviously ofthe form tk t ( x ) = ~ t . o ( X ) + : t t r b o ( X )
Calculation of flux in subcritical reactors in which the coefficient ~t has to be defined. Because it is intended that ~bz exists, the right-side of equation [27)., must be necessarily orthogonal to O~, leading to (Fff t.o, O~) :¢~ = (FO o. O~) ' In a similar manner, ~z, ~a, etc. can be calculated and the process may be indefinitely pursued.
Remark 5.1.
187
0 ( ~ - ~); more precisely there exists a certain positive constant CO such that:
c°
when ~ is small enough. (See the appendix for the value of Co.) Now, let us consider the exact solution "P~ of equation (25) and suppose the different terms ~'t, ~.,, etc. of equation (26) have been calculated. Taking the first terms up to zz~bz (to fix the ideas), we then define:
It is possible to write the source in the form
s(x) = eSo(X) + e2s I (x) + • - • + eksk_ ~(x) + . . . ; then the different terms sk would be present in the right-sides of equations (27). In particular, if So=0, one has ~ o = 0 and then e t ( x ) = ~ l O o ( X ) with (s,, ~t -
O~)
(FOo ' @~)
Consequently, the flux is of the order of e : ~ , ( x ) = ~ ~¢o(X) + 0(~2). More generally, ifs o = s i = • • - = s~_ ~= 0 and s~ ~: O, then ~ ( x ) = ~£~@o(X) + 0(~:~+ ~) with (s~,
~u -
O~)
(FO o , O*)
Remark 5.2. If s(x)=So(X)+tst(x)+ . . . . it is necessary to set $ , ( x ) = ~ o , - t O o ( X ) + $l(x)+e~b2(x)+ . . . . and then (s o,
~o =
O~)
(F¢o '
¢~)-
Remark 5.3. If the perturbation is caused by e Aa, (instead ofeF), all the calculations remain unchanged except that F must be replaced by Aa~ in the expressions of coefficients ~o, ~t, ~2, etc.
5.2. A posteriori justification of the asymptotic expansion of ~ , The above calculations are purely formal since it is not sure that the series ~o¢o+e~b t +e2~b2 + . . . converges. However, this technique is currently used by physicists in q u a n t u m theory even in cases where the series does not converge, but they are only interested in the sum of the first terms (see F e y n m a n n , 1965; Rellich, 1969). Our objective is now to give a justification of this point. Firstly, it is observed that I[( A - ( 1 - e ) F ) - ' l l =
(301
£
0~ = ¢ , -
(~o¢o + ~ t + ~0,)-
Clearly we have, after some handling: (A - (1 - e.)F)0~ = -
PFO,,
and it results directly from equation (30) that
I[o, ll <-Co fll More generally, + ~kOk), then
if
211.
0~=t,b~-(:%Oo+:q~k I + - "
I1¢11-
(3 1
Thus, the flux ~ can be approached by a sum of the first terms of the asymptotic expansion, even when the series is not convergent; the inequality equation (31) gives a bound of the error caused by this approximation: it depends on the last term ~,~ calculated in this way. This a posteriori justification is widespread in perturbation theory methods and many other applications can be seen in Lions (1973). The method presented in this section necessitates functions Ok(x) to be determined and only a few functions are required if ~ is small enough; because they are corrective terms, it is not necessary to compute them with a high accuracy and only coarse discretisation meshes can be used.
6. CONCLUSION Thus, the problem of neutron diffusion in the presence of sources always admits a steady positive solution in the case of subcritical reactors. In the event of low intensity of sources and if the reactor is close to criticality, this steady solution can be approximately achieved by standard critical eigenvalue calculation; the flux obtained in this way is determined to within a multiplicative factor, and it was shown that this factor is fixed by the source intensity and the antireactivity. A method, based on asymptotic expansion of the flux, is proposed for the study of the spatial effects of sources on the fundamental eigenmode obtained by critical calculation. We hope in a future paper (by J. P. West and the author) to present numerical results.
188
J. PLANCHARD
Acknowledgments--The author thanks M. Salvatores and G. Palmiotti who have kindly communicated to him the curve of Fig. 1.
where Q = A - 'F. Obviously, the Oi are the eigenfunctions of Q, and y~ can be expressed as Y,= ~
i=O
REFERENCES
Blanchon F., Ha-Duong T. and Planchard J. (1988) Numerical methods for solving the reactor kinetic equations. Proy. Nucl. Energy 22, 2, 173-180. Feynmann R. (1965) La nature de la Physique. Editions du Seuil, Paris (French translation). Gandini A. (1967) A generalized perturbation method for bilinear functionals of the real and adjoint fluxes. J. Nucl. Energy 21, 735. Gomit J. M. and Planchard J. 0985) Perturbations conservant la criticit6 et fonction Importance. Bulletin des Etudes et Recherches, Electricit6 de France, S6rie A, 1, 17-22. Greenspan E. (1975) A source multiplication reactivity. Nucl. Sci. Engng 55, 103-105. Habetler G. J. and Martino M. A. (1961) Existence theorems and spectral theory for the multigroup diffusion model. Proc. Syrup. Appl. Math., I 1, Nuclear Reactor Theory, Am. Math. Soc., Providence (R.I), pp. 127-139. Krein M. G. and Rutman M. A. (1962) Linear operators leaving invariant a cone in a Banach space. Am. Math. Soc., Translation No. I0, 199-325. Lions J. L. (1973) Perturbations singuli~res dans les probldmes aux limites et en contr61e optimal. Lecture Notes in Mathematics 323, Springer Vertag, Berlin. Palmiotti G. and Salvatores, M. (1982) (private communication). Planchard J. and West J. P. (in preparation). Protter M. H. and Weinberger H. F. (1967) Maximum Principles in Differential Equations, Prentice Hall, New York. Rellich F. (1969) Perturbation Theory of Eigenvalue Problems, Gordon and Breach, New York. Riesz F. and Nagy B. (1955) Le?ons d'analysefonctionnelle, Acad~mie des Sciences de Hongrie, 3rd Ed., Budapest. English translation: Functional Analysis (1965) Ungar Publishing Company, New York. Taylor A. E. (1961) Functional Analysis, Wiley, New York. Usachev L. N. (1964) Perturbation theory for the breeding ratio and for other number ratios pertaining to various reactor processes. J. Nucl. Energy, 18, 571. Varga R. S. (1962) Matrix lterative Analysis, Prentice Hall, New York. West J. P. Performance des Modules de Diffusion des Syst~mes CARDIF (code SNAP-3D et CCRR--Calculs Homoo~nes et InhomogOnes d 2 Dimensions, Electricit~ de France, Internal Report No. HT-11/8718. APPENDIX
Estimate of the constant CO Our purpose is to obtain an upper bound of IIEA-(1WJ-'11 This norm is defined by: ~,~ Ily, lllllfll, where y, is the solution of [A -- (1 -- 8)F-Jy~ = f or equivalently:
[1 - ( 1 - e ) Q ] y , = A - 'f=-9,
(A1)
~i(~i
with
(g, uJ) ~o =
(~, u*) , cq=
1 -- (1 - ~))'i
, i>0
where u~' are the eigenfunctions of Q*, with the orthonormalization condition (~i, u*)=(}ii- Because • ~, u~' form a biorthonormal sequence, the numbers (9, u~*) are the components of g with respect to this sequence. Consequently, a positive constant C exists such that (see Riesz and Nagy, 1955):
.?)l
"< cIIgll,
where C does not depend on 9Clearly
Ily
'
max i>o
II-l(1-e)2i
gll
The numerical calculations show that practically the second eigenvalue ;-t is positive and close to 20( = 1), so that one gets
Ily, II-< IIA-' IIIISlI. Then
Co =
C
JIA-'II- CllA-'ll ,
The difficulty is that we are unable to exhibit the value of C, unless the ~ and u~' are explicitly known. However, that can be achieved if A and F are symmetric operators: this is the case in one-group problems. In this situation, A being symmetric and positive definite, the square root A r of A(A = (A½)2) can be defined (see Riesz and Nagy, 1955) and equation (A1) is equivalent to: [1 - ( I - e ) 0 ] L = f
(A2)
where 1
!
1
1
Q = A -~FA -~, A - ~ = (A ~) - 1
-'-7 is symmetric positive definite and has the same eigenvalues as Q. It ensues that the norm of the inverse RE of I - ( 1 - ~ ) Q is:
Calculation of flux in subcritical reactors 1
where d, is the distance from 1 to the spectrum of (1-~)0, i.e.
d,=ll -(1-
189
and the minimum is reached for a certain function u 0 (u o is the eigenfunction associated with ~,). From the maximum principle (Protter and Weinberger, 1967), uo is positive inside fL We have ~, = ( ( A - F + eF)u o, Uo ) + e( Fuo , Uo )
~)X,I = e.
This is due to the fact that the norm of a symmetric operator is equal to the modulus of its greatest eigenvalue. Thus
II;:11
I1 1.
on the other hand, IIfGI_
- 11 II;,11.
IlY
IIA Since the norm of A - r is strictly equal to the square root of the greatest eigenvalue of A - t, or equivalently, of the inverse of the lowest eigenvalue fl of A, we then have l
Ily.ll-< llSII, and consequently Co = 1/ft.
(A3)
For PWR cores, fl is of the order of the absorption cross-section, i.e. 10 -2 (see Blanchon et al., 1988); more precisely
> min ( ( A - F ) u , u ) + e ( F u o, Uo), It,ll = t
and m i n ( ( A - F ) u , u) is zero because the reactor II~t = t
characterized by A, F is critical and A is positive definite. Thus '~, > e( Fu o, %). It is noted that u o implicitly depends on e, but Uo~-~o where ~o is the fundamental mode of the reactor, and then
I[*o[l:
CO= (FOo ' Oo)"
(A4)
The estimates (A3), (A4) are rigorously valid for one-group problems. Now, it is well-known that multigroup systems may be approximated by onegroup equations. Then, these estimates remain valid with F and fl computed for the equivalent one-group problem. It is important to note that the exact value of may be easily obtained by inverse iteration method (a few iterations are sufficient).
/~> min a.(x). Remark
Another estimate of C Ocan be derived. Indeed, let us come back to (At). Clearly
When the perturbation is due to variation of try, eF must be replaced by e A% in equation (Al), so that equation (A4) becomes:
1
N ( A - - f +eF)-tll =-~t E
where ~ is the smallest eigenvalue of A - F + eF: ~t,= rain ( ( A - F + e F ) u u), i1,!1= 1
Co =
ll*oll (Aa.~o, ~ o ) '
(g4 bis)
In (A4 bis), Aa, is the reference absorption crosssection variation for the equivalent one-group problem.