The pseudo-harmonics perturbation method: Application to a PWR

Ann. nucl. Eneroy, Vol. 14, No. 2, pp. 99-106, 1987

0306-4549/87 $3.00+0.00 Copyright © 1987 Pergamon Journals Ltd

Printed in Great Britain. All rights reserved

THE PSEUDO-HARMONICS PERTURBATION METHOD: A P P L I C A T I O N TO A P W R F. C. DA SILVA,t S. WAINTRAUBand Z. D. THOM~ COPPE, Nuclear Engineering Programme, Federal University of Rio de Janeiro, P.O. Box 68501, 21 944 Rio de Janeiro, Brazil A. C. M. ALVIM~ CNEN, 22294 Rio de Janeiro, Brazil (Received 30 June 1986) Abstract--A modified version of the pseudo-harmonics method for perturbed neutron flux reconstruction is described together with a number of applications relevant to a PWR-type nuclear reactor, with three enrichment zones and with material constants representative of those of the ANGRA-I power plant. Global and localized perturbations were simulated. The eigenvalue convergence criterion showed that only a few orders of approximation were needed to calculate the neutron flux alterations, indicating the potentiality of the method investigated. The number of pseudo-harmonics used depended in each case on the type and intensity of the perturbation considered.

fission operator which are associated with k = 0 have to be included in the expansion. Determination of the coefficients associated with these eigenvectors may be a difficult task. An example of this method, developed by Gandini (1981), was tested by Palmiotti (1983) for two different fast-reactor configurations, using a onegroup, 1-D model. The s-modes (Filippone and Conversano, 1983), pertaining to the asymptotic problem operator, do not exhibit the same type of problems, since, on physical grounds (Henry, 1975), for a given configuration, an c~= 0 value, should it occur, is expected to be unique. It should be mentioned, however, that in multigroup problems, since both k- and s-mode-associated operators are not self-adjoint, it may well be possible for complex and/or degenerate eigenvalues to occur. In this case, both methods may require special routines to deal with complex numbers and/or compute degenerate eigenvectors to form the basis. The pseudo-harmonics method (Gomit et al., 1982) expands the perturbed flux in eigenfunctions of an operator which is only a part of the reference (unperturbed) problem operator. This 'partial operator' is chosen so as to decouple the multigroup equations, when determining the harmonics. In this manner, each harmonic of a given group is 'isolated' from the contributions of other groups. The group operators so chosen are self-adjoint, which from the viewpoint of perturbation theory has two practical benefits : (1) it is not necessary to corn-

1. INTRODUCTION

As is well-known, perturbation methods in nuclear reactor physics are widely used in the analysis of nuclear reactors, due to the computational economy which they allow. In particular, perturbation methods have been considered for the reconstruction of the neutron flux after alteration by system parameter changes (Gandini, 1978, 1981; Gomit et al., 1982; Palmiotti, 1983). By these methods, the altered flux is generally expanded in terms of eigenfunctions calculated at unperturbed conditions. One of the methods for neutron flux reconstruction is the so-called 'implicit method' (Gandini, 1978) using a recurrent scheme for the auxiliary functions calculation. Since a new calculation of these auxiliary functions has to be run for each new perturbation, this method has limited practical applications. To avoid this difficulty, methods for solving the system explicitly with respect to the perturbation have been developed (Gandini, 1978, 1981 ; Gomit et al., 1982), which use the expansion of the perturbed flux in a basis of eigenfunctions (modes) pertaining to a chosen operator, related in some sense to the problem. The k-modes (Filippone and Conversano, 1983), corresponding to the unperturbed problem, present difficulties related to the singularity of the fission operator. Due to this singularity, the eigenvectors of the tSPresent addresses: tENEA-CRE, Casaccia, Rome, Italy; :~COPPE, Federal University of Rio de Janeiro, Brazil. 99

100

F.C. DA SILVAet al.

pute the unperturbed adjoint group fluxes; and (2) there are neither complex eigenvalues and eigenvectors nor eigenvector degeneracy.

Expanding the perturbed flux, ~', and the perturbed eigenvalue 2', one has ~t=

2. THE PSEUDO-HARMONICS METHOD

~ ~ ( k ) = {~(0)_{_ ~.~ ~(k) k=0 k=l

(5)

and

The pseudo-harmonics method, developed by Gomit et al. (1982), is discussed using a multigroup formalism. The diffusion equation for a critical system can be expressed by ( A - & F ) ¢ o = O,

(1)

where A = B + S, 20 is the fundamental eigenvalue, F is the fission operator, B represents the leakage plus removal operator and S is the scattering operator. The unperturbed neutron flux ~0, in the fundamental mode, is represented by

2'= ~2(k)=2(0)+~2 k=0

(6)

where the index k in the above equations means order of approximation in the expansion, and obviously ~(0) = ~0 and ;t(°) = 20, with ~(k) satisfying the following system of equations :

{

(A--20F)~ (°) = 0 (A - 2oF)~ (k) = s (k),

where ~1,0

(k),

k=l

k = 1,2, 3. . . . ;

(7a) (7b)

k

s (k) = - 6 A d p ( k - ')+ ~. 2(°F~ (k-0 I=1

¢o =

~b~.o

,

k--l + ~ 2(t)aF~ (k-t- ').

(2)

(8)

I=0

~G,0

where G indicates the number of energy groups. If a perturbation occurs in a nuclear reactor, the k

,~(k) =

I

k-I

l=1

1=0

perturbed neutron flux is the solution of the following equation : (A' - 2 ' F ' ) ~ ' = 0

Considering that the operator (A - 2 o F ) is singular, the kth equation (7b) will have a solution if, and only if, (~o, s (k)) = 0. This condition implies that

(3)

where

(9)

where ( ) indicates integration in a volume as well as matrix operations, which are also used. The adjoint flux in the fundamental mode, ~*, satisfies the following equation : (A*--XoF*)~* = 0,

(10)

A' = A+fA,

where A* and F* are the adjoint operators of A and F, respectively, and

B" = B + f B , S" = S + 6 S , F" = F + 6F,

.4,*,1o

2' = 20+62

4go [ .

and

(11)

4~,o J ,

=

(4)

One observes that the general solution of the differential equation (7b) is not unique, since all functions of the form ~(k) -±.~(k) - - ,(k)a~ ~ 1/10T Wpart"

(12)

The pseudo-harmonics perturbation method satisfy the system (7b), where c (k) are arbitrary constants of the homogeneous solutiont and ~b(k) is the particular solution of the kth equation of system (7b). The solutions 4~(p~rt.,for all k, belong to a space void of ~b0. This space, denominated H0, is spanned by the eigenvectors of the higher modes (called harmonics) 4~ (i = 1, 2, 3, . . .) which satisfy the following equations : (A -- ).,POq~, = 0, These harmonics relations :

i = 1,2,3 . . . .

satisfy

the

f

101

a)e.i(r)Oggd(r) dr = 3ij

0 i

(13)

"

0

orthogonality

,-- gth row j = 1,2,3 . . . . .

O.)g,j

We,j =

fori¢n;W,n=O,

I:

1,2 . . . . ; (14)

i

where $* (n = l, 2 . . . . ) are the adjoint eigenvectors of the higher modes, which satisfy ( A * - - 2 , F * ) ¢ * = 0,

n = 1,2 . . . . .

(15)

The vector 4 ~ . can be represented as a linear combination of the vectors ~'e.J, the set of which form a basis for H0 : ~(k) part.

----

~

~ ~(k)~ . J,e g,J"

The pseudo-harmonics are the eigenfunctions of the operator Bg, representing leakage ( - V" DgV) plus removal (Z~), for each energy group, where g = 1,2 . . . . . G,

(17)

and

The vectors Wea belong to a space spanned by the harmonics ~ ( i = 1, 2, 3 . . . . ) and the flux 4~0 (fundamental mode). So, the vectors wej can be expressed as wga = ~ b~4~i,

The operator B e has real and positive eigenvalues ?gj, for each group g, and the associated eigenfunctions (pseudo-harmonics) cog.j are the solutions of the following equations : j=

1,2,3 . . . . .

(21)

i=0

(18)

The pseudo-harmonics obey the same b o u n d a r y conditions which are imposed on the unperturbed neutron flux. The orthogonality relation among them is expressed by

(22)

( gp*FWed)

Defining a space Ho, void of 4~0, one can construct a basis for this space, i.e. the set of vectors qte,j which result from we,j, after the contributions of the ~0 components are removed ::~
X~ = Z ,e + D e B z2 +Z~9~e' ,. B~ = axialbuckling.

BefDg,j=~g,jfDg,j,

•

where the coefficients b~ are determined by the orthogonality relation (14) as

3. DETERMINATION OF T H E BASIS VECTORS ~ , j

Be= -V'DgV+X~.,

(20)

0

(16)

j~lg=l

(19)

F r o m the pseudo-harmonics c%j, one can construct G vectors wad, having the following form :

0

(~*F~b~)=0,

dr.

~ke,j = We,j

(~o,F~o) ~0.

(23)

Having constructed this basis, the next step is to calculate the coefficients ~t),ka ) of the expansion (16).

4. DETERMINATION OF T H E COEFFICIENTS g~g)

Substituting expression (I 6) into equation (12) and using the resulting expression in the kth equation of system (7b), one has = s k,

(24)

j~lg=l

t Explicitly assuming solution ~(k)as the sum of an homogeneous part (q~0) and a particular one (4~t.), void of the fundamental flux mode, allows us to univocally determine the coefficients c~k),following a given flux normalization criteflon [see equations (33) and (34)], once the coefficients a)k) have been univocally determined. ~:In contrast to the derivation by Gomit et al. (1982), in the present method the fundamental flux (~0) mode contamination is removed from the eigenvectors Wgd.

F r o m equations (23) and (24),

j=lg=l

where (4~*Fwej)

(26)

F.C. DASILVAet al.

102

Using equation (1) and considering that A = B + S, where

HI,I

gl I , N B1

0

U2,1

Bo

B=

U2,N

U=

Bo

UG, 1

it follows that

HG,N G

~),~(B+S--zoF)wy.j= s .

(27)

and

j=lg--I

=

However, since Bwg,yg,j, one has

~, o~}~)u[7o•,wo•y+(S-).oF)wv,,]=s{k'. j--I

(28)

rl]

.

•

tin

12 rtl

...

•

I 1 rNN

r~]

..

22 rl 1

"

12 riN

• "

IG rll

•

12 rNN

..

rf~ ... r'~(~

" • •

22 rlN

• -

2G rll

,..

2G FIN

• • •

22 rNN

• • .

2G rN 1

• • •

2G rNN

.

r l~N

...

r~

.. •

r la aN

G2 rNN

...

GG rNi

...

GG rNN

...

IG rlN

:

.q--I

Truncating the series (28) in the Nth vector w.. for each group g, multiplying it on the left by %.~ and integrating over the whole volume, one obtains: N

G

{ g,k )~. E E r•~"vn!~) ~ .,(k) Hg, i(Zi, t i, g ~ /--Ig-.

(29)

21 rl I

..,

21 rl N

~21 1N I

"• "

r21 NN

22 rNI

•

rl N

r~

G1 rNN

G2 rNi

I

where

f

u¢., = y~,,J%,,(r)%,.~(r)

GI rN[

dr,

r~"u = (w~,,i(S- 2oF)w..j). v ( k ) = ( w ¢ i , s (k)) i.q ""

~ i = 1,2 . . . . . N (g 1,2 . . . . . G.'

Expressions (29) constitute a system of N x G linear equations, which can be expressed in matrix notation as

(U + ~)~¢k~ = v~k),

(30)

where (k)

~(k) i

~(k)

,

V (k) :

:Z~)2

,{k) l l,(J

ll( k N,G

...

,..

•

The components ~}~ of the vector a(k) are the coefficients of expansion (16)• The U and ~ matrix elements are calculated by using only the parameters of the reference reactor (unperturbed system) and so, they are calculated once, for any perturbation caused in the system• It should be noted that the whole effect of the perturbation is concentrated in the Vk) vector components• Thus, for a certain number of energy groups and having fixed the number of pseudo-harmonics, the matrix (U + E) ~ can be stored, for later use in the computation of the ~~k) ~ j , ~ . From a practical point of view, this is extremely useful, since for a certain order of approximation k and for a given perturbation, the ~)~ are easily calculated by the woduct of the matrix (U + ~) 1, which is computed only once, and the column vector v(k), which depends on the order of the approximation and on the perturbation imposed on the system•

5. D E T E R M I N A T I O N

OF THE

CONSTANTS

c tkl

To determine the constants c(k) of expression (12), one needs to impose the condition that the total fission

The pseudo-harmonics perturbation method rate be equal to a constant value P, before and after the perturbation. That is,

103

Table 1. Reference reactor characteristics Region

(31)

(Y,r~bo) = ( ] ~ b ; ) = P, where

W i d t h (cm) Enrichment (%)

Z;f=

Zfo

and

I~;=

kfG

~;.g

.

k;oj

Condition (31) should be valid for any order of approximation considered. Substituting ~', in condition (31), by the expansion

r

~(o

0+

+= p.

Recalling that expression (31) is valid for any order of approximation, (]B~b') "-~ P+(Z~b(k)> = P.

(33)

Substituting equation (12) into equation (13) and considering the expansion (16) with only N pseudoharmonics, one has _j.gn~ k) (Y,'f~g.j)

](

(l~,'rtko)

)

72 2.1

29 2.6

22 3.1

Z "~(')

(32)

=lj=l

3

l 0

I

c (k) = --

2

geometry and axial buckling of 7.3678 x 10 5 cm 2) and two energy groups. The order of approximation that is required in the flux expansion, according to expression (5), was determined by the following eigenvalue convergence criterion : ~(k) x ~< e,

(5) truncated in a certain order k, one has (I2;4~'> "~

1

.

(34)

where k is the required order when the above criterion is satisfied. In all cases presented in this paper, a value o f ~ = 1 0 4 w a s adopted. The results which follow show a comparison between the pseudo-harmonics method and direct calculation of perturbations in the reference reactor. In all cases presented, the effects of the perturbations on the fast-group flux were less pronounced, due to the type of perturbations used. For this reason, the fast flux is not presented here. However, it should be noted that the perturbed fast flux was very accurately predicted using the same order of approximation and the same number of pseudo-harmonics which were used, in each case, to reproduce the perturbed thermal flux. Case 1

The particular solution 4¢pka~.of the kth equation of system (7b) is determined by using expansion (16) after having calculated the coefficients ~k) --j,g The homogeneous solution of the kth equation of system (7b) can be determined by using the constants c(k) already calculated. In this way, the general solution expressed by equation (12) becomes known and therefore, one can construct the perturbed flux, using expansion (5), calculated up to the kth order of approximation.

6. R E S U L T S

Several perturbations were simulated in a thermal reactor to verify the potentiality and accuracy of the methodology. The reactor was a 1876 MW PWR, with three enrichment zones, similar to the ANGRA-I reactor. The parameters used are shown in Tables 1 and 2. Calculations were done in I-D (cylindrical

The case considered is a perturbation simulating an absorber (with 4.8 cm radius) in the centre of the reactor, characterized by a macroscopic absorption cross-section 30% higher than the unperturbed one of the same region, for both energy groups. Calculations were done with 1, 5, 11 and 22 pseudoharmonics. The flux was correctly reproduced with 22 pseudo-harmonics and eigenvalue convergence was achieved with the third-order approximation. The perturbed eigenvalue 2', calculated by both the direct method and the pseudo-harmonics method, agreed up to 5 significant figures, with a value of 1.0010. In Fig. 1, the results obtained for the perturbed flux are shown for each number of pseudo-harmonics considered. For comparison, the perturbed flux obtained by direct calculation is also shown. Case 2

This perturbation represents a 20% increase in the macroscopic cross-sections of all enrichment zones.

104

F. C. DA SILVAet al. Table 2. Material constants for the referencereactor m

g

D~, a

~,m

V~m

~m

~,~l~g

~,m,B

Xm~

1

1 2

1.3096E+ 0 4.7339E-- 1

1.1687E--2 8.9091E-2

6.9560E-- 3 1.3378E- 1

2.7370E-- 3 5.5055E-2

1.0004E--2

1.5198E 4 7.5144E-3

1.0 0.0

2

1 2

1.3094E+0 4.7152E-1

1.2053E 1.0126E

7.8450E-3 1.6018E-1

3.0950E 6.5919E

9.5280E-3

1.4609E-4 7.3985E-3

1.0 0.0

3

1 2

1.3096E+0 4.6920E 1

1.2410E- 2 1.1287E- 1

8.6910E 1.8527E

3.4370E- 3 7.6244E- 2

9.1180E

1.4101E-4 7.3007E- 3

1.0 0.0

2 1

3 1

10¸

3

10

/

-~

- -

~'~/ 7~

b

3 2

Direct Calculation I Harmonic

.... -----

6

~-~ /- -

~'~

5 Harmonics 11 Harmonics ......... 22 Harmonics

"~ ~

''

- ....

Direct CalculatiOn 1 Harmonic

.....

5Har.

ics

.....

~ ~

PerturbGtion Zone

ImJ~J Perturbation Zorze

4

N

N 2

R~ION I O,

'

.

2 .

.

.

i

', 3 \ 1 ~

I REGIONI ~ \

i

i

~

I

,

I

k~O

,

,

4O

80

,

120

I

I

t60

I

200

RADIUS ( cm )

RADIUS ( cm )

Fig. 1. Perturbed flux calculated with third- order approximation. Type of perturbation: 30% increase in Z=.

Fig. 2. Perturbed flux calculated with fourth-order approximation. Type of perturbation: 20% increase in E,.

Eigenvalue convergence was achieved with the fourthorder approximation :

bers of harmonics. Figure 3 shows the perturbed fluxes, calculated with 5 and 25 harmonics. This shows the effect of the number of pseudo-harmonics which are used to construct the perturbed flux shape, mainly in the perturbation region. In this case, it was necessary to consider the fourthorder approximation to get convergence in eigenvalue %. The values obtained are

2' (direct calculation) = 1.2332 and 2' (pseudo-harmonics) = 1.2330. Figure 2 shows the results obtained for the perturbed thermal flux with 1 and 5 harmonics. However, it should be noted that the perturbed fluxes obtained by direct calculation and the pseudo-harmonics method with 10 harmonics, are practically coincident.

1' (direct calculation) = 1.0116 and 2' (pseudo-harmonics) = 1.0120.

Case 3

In this case, the perturbation represents the inclusion of an extra absorber (of 3.6 cm width) in the centre of the second enrichment zone, having a macroscopic absorption cross-section 30% higher than the unperturbed one in the same region, for both energy groups. Several calculations were done for different num-

Case 4

In this case, the perturbation consisted of the inclusion, in the central part of the second region, of a new enrichment zone (of 3.6 cm width), thus locally modifying all the material constants. The new material constants were chosen identical to those of the third region.

The pseudo-harmonics perturbation method

105

The eigenvalues obtained are

20- -

Direct

.... .....

5 25

2' (direct calculation) = 0.9984

Colculation

Hormonics Harmonics

Perturbation

and

Zone

c~

2' (pseudo-harmonics) = 0.9983.

x

7. CONCLUSIONS

!~6~ [ [ , ~ntll ', 31

~-6t0N1 c~

4o

8o

,2o

RADIUS

,~o

( cm )

Fig. 3. Perturbed flux calculated with fourth-order approximation. Type of perturbation: 30% increase in Ea.

10;

8

~

- -

Direct Colculation 25 Harmonics

M

Perturbotion Zone

i,o

X

~

~

! ~,oN 2 ! ~:~3 : r

;,

;

; 40

;

i

; 80

RADIUS

I 120

I 160

I 200

( cm )

Fig. 4. Perturbed flux calculated with second-order approximation. Type of perturbation: inclusion of a new 3.1% enrichment region. (All perturbed parameters are identical to those of region 3.)

Eigenvalue convergence was achieved with the second-order approximation, but 25 harmonics were needed in order to reproduce the perturbed flux shape. Figure 4 shows the perturbed thermal flux for comparison between the direct calculation and pseudoharmonics results.

One of the main advantages of the pseudo-harmonics method lies in the choice of the basis used to represent the perturbed fluxes. This basis is formed by the eigenfunctions of an eigenvalue problem, whose operator is the self-adjoint part of the original problem operator. These eigenfunctions are real, orthogonal and form a complete set. From the numerical point of view, the fact that the operator is self-adjoint makes it easy to obtain the eigenfunctions. The analysis of the cases studied shows that, for a given level of intensity, the more the perturbation is localized, the larger the number of the pseudoharmonics required is, in order to correctly reproduce the altered flux shape. This effect is well-evidenced in Cases 1, 3 and 4. Given a convergence criterion on the eigenvalue, depending on the desired accuracy of the calculation, a few orders of approximation have been sufficient to reproduce flux and fundamental value changes, allowing a rapid convergence. In conclusion, this explicit method has proven to be very accurate in reproducing the altered neutron fluxes in a variety of cases in which global or localized perturbations are considered. Although the calculations were performed in two energy groups, this should not represent a limitation, the basic features of the method being preserved when more groups are assumed. Considering the speed with which the perturbation terms can be calculated, once the eigenvector basis has been formed, this method can be advantageously used in reactor system analysis. As far as future work is concerned, an effort is underway to write a 2-D code for dealing with more realistic cases. No major difficulties are anticipated, apart from the obvious fact that, due to the double indexing of the eigenfunctions, coincident eigenvalues may occur, necessitating some care in separating the corresponding eigenvectors.

Acknowledyements--The authors acknowledge the Laborat6rio Nacional de Computa~o Cientifica, for providing computer facilities. The authors would like to thank the Financiadora de Estudos e Projetos and the Conselho Nacional de Desenvolvimento Cientifico e Tecnol6gico for financial support.

106

F.C. DA SILVAel al. REFERENCES

Filippone W. L. and Conversano R. (1983) Nucl. Sci. Engn9 83, 75. Gandini A. (1978) Nucl. Sci. Engn9 67, 347. Gandini A. (1981) Nucl. Sci. Enyn9 79, 426.

Gomit J. M., Planchard J. and Sargeni A. (1982) Report HI 4345-07, D6partment Physique de R6acteurs, Electricit6 de France (1982). Henry A. (1975) Nuclear Reactor Analysis. MIT Press, Cambridge, Mass. Palmiotti G. (1983) Nucl. Sci. Engng 83, 281.