A finite element formulation for perturbation theory calculations

A FINITE ELEMENT

FORMULATION FOR

PERTURBATION THEORY

CALCULATIONS

Bilge Özgener, Sibel Kaluç

Energy Institute

1. UNDERLYING THEORY

 When the configuration of a nuclear system is to be modified, it is desirable to know the ensuing change in reactivity.

 If the effective multiplication factor of the original system ( ) is predetermined by a fission source iteration, the determination of the change in

reactivity () requires the calculation of the effective multiplication factor of the modified

system ( ) by another fission source iteration

eff

k

eff



(1)

 During fission source iteration, some

numerical model (i.e.finite differences, finite elements, boundary elements etc.) must be employed for the solution of the multigroup neutron diffusion equations

.

eff eff k 1 k 1     

 Assuming multigroup neutron diffusion theory constitutes a sufficiently accurate model of the nuclear system, this method for calculating the change in reactivity, which we'll call,"the diffusion theory method" or DTM is capable of calculating  almost exactly provided that a sufficiently fine spatial mesh is used in the numerical model.

 When the number of modifications in the system

configuration we must consider is large, the

computational load of DTM can become excessive since each new modification requires a seperate determination by a new fission source iteration.

eff

k

 On the other hand, if the modification in the

system configuration is neutronically not overly significant, the approximate multigroup

perturbation theory expression can be used for the approximate calculation of . The use of the multigroup perturbation theory expression requires only the knowledge of the effective

multiplication factor( ) , the group flux vector

 (2) eff k

 

r



₁

 

r ₂

 

r _g

 

r _G

 

r



T        _{    } 

 The group adjoint flux vector

(3)

 The modifications in the matricial multigroup

operators and : (4)

 

r



₁

 

r ₂

 

r _g

 

r _G

 

r



T             _{    } 

 

_{ }

 

_{ }

 

_{ }                                                 r G , r 1 g G , s 1 G , s r g , r g 1 g , s r 1 , r 1 r D r D r D                      M  _F

(5)

 Now we define the inner product of two

functions f and g as:

 (6)

 







 





 



 







 





 



 







 





 



                                                   r r r r r r r r r G , f G G g , f g G 1 , f 1 G G , f G g g , f g g 1 , f 1 g G , f G 1 g , f g 1 1 , f 1 1                    

 

f ,g f

   

r g r dV V  





 It can be written as:

(7)



In contrast to (1), (7) requires no

knowledge about the effective multiplication

factor of the modified system ( ).







 





 



          _   F , M , F , k / 1 _eff eff k

 Once and the group flux - adjoint flux vectors

of the original system is determined by fission source iteration, the approximate change in

reactivity for any modification of the original

system can be calculated simply by carrying out the indicated integrations in (7) without any

further fission source iteration. The use of (7) for the approximate calculation of  will be called "the perturbation theory method" or PTM.

 The use of both DTM and PTM requires the

solution of the generalized eigenvalue problem for the original system:

eff

k

(8)

 The use of PTM also requires the

determination of the group adjoint flux vector of the original system by solving the adjoint

generalized eigenvalue problem whose

eigenvalues coincide with those of (8), namely:

 (9)    F k 1 M eff   _{ }  + eff + F k 1 M



_{and are the matricial adjoints of}

and respectively. The use of DTM also requires the solution of the generalized

eigenvalue problem of the modified system for the determination of namely:

(10) (11) (12)  M F eff k M F         _F k 1 M eff M M M    F F F   



_{being the group flux vector of the modified}

system.

 The solution of (8), (9) or (10) by fission

source iteration requires the solution of the

within group diffusion equations consecutively. If we designate the iteration number in fission source iteration by n, the within group diffusion equation for the g'th energy group can be

written generically both for the forward and adjoint calculations as:

(13)





 

r

 

r

   

r r q

 

r D_g   _g  _r,g  _g  _g           



₍₁₃₎ where (14) (15)

 

r

 

r

   

r r q

 

r D_g   _g  _r,g  _g  _g           

 

 _{ }         _ ns calculatio adjoint for ns calculatio forward for r n g n g g 

 

 

 

 

 

 

 

 

 

 

 

 

   

 

                     









                         G 1 g g G 1 g 1 n g g 1 n g , f g n g g g , s 1 g 1 g G 1 g 1 n g g , f g 1 n g n g g g , s g r k r r r r r k r r r q         

 In this study, we tried to assess the merit of the approximate PTM relative to the rigorous DTM. For ease of application, the system geometry is taken as one dimensional cylindrical geometry. In case of 1D cylindrical geometry (13) simplifies to:  (16)

 where R is the system radius. We also assume

that the zero-flux vacuum boundary condition is valid.

 

   

r r q

 

r , 0 r R dr d r rD dr d r 1 g g g , r g g             

 We solve (16) numerically by the linear finite

element method. We first divide the cylindrical system into N homogeneous annular rings

(finite elements) with outer radii ,n=1,…N. Using a Galerkin type weighted residual finite element formulation [1], we demand that:

(17) where are the linear finite element basis

functions. is the basis function of node n.

n

R

 

rdr

     

r r w r rdr q

   

r w r rdr 0 dr dw dr d r D R 0 g R 0 g g , r R 0 g g      







 

r h_n

 

r h_n

 It assumes the value of 1 at the n'th node and

vanishes at all other nodes. It varies linearly

 (18)  (19)  radially over the (n-1)'th and n'th elements

(annular rings for the 1D cylindrical geometry) whom it connects. It is zero elsewhere.

 

_

 

  N 0 n n n r w h r w

 

_

 

    N 0 n gn n g r h r

 The requiring (17) to be valid for all possible

weight functions of form (18) results in the N dimensional linear system:

 (19)

 _{where the matrix elements and right hand}

side (RHS) vector are of the form:

 (20) g g g q A  

 



_

R _

 

 

     

_ 0 m n g , r m n g nm g r h r h r rdr dr dh dr dh r D a



₍₂₁₎

 The form of the RHS vector in (22) is obtained

by employing the "lumped source approximation" .

 



_

R

   

0 n g n g q r h r rdr q

2. APPLICATIONS

 A FORTRAN program called PERTURB has

been developed to calculate the change in reactivity for modifications in system

configuration. The program is capable of

calculating and the group flux distributions within the context of multigroup diffusion theory using linear finite element discretization of

Galerkin type. eff

 Using PERTURB, the reactivity change can be

calculated either by the rigorous diffusion theory approach (DTM) or the approximate perturbation theory approach (PTM).The program has been developed in the MS

FORTRAN environment under the WINDOWS XP operating system.

 The program has been validated via

comparison with analytical solutions.The relative merit of the PTM has also been assessed by means of this comparison.

eff

k_eff k

 Then, PERTURB with both DTM and PTM

options has been run for the determination of the reactivity change ensuing by the

replacement of a fuel rod by water in a

reactor.Thus assessment in realistic situations has also been rendered possible.

 The first problem we consider is a

bare,homogeneous ,critical,cylindrical reactor with radius R=36.07238336536454 cm.The one group constants for this reactor are:

 D=0.9 cm, ,

For cases when the absorption cross section is increased by

( ) and

( ) in the central cylindrical section of radius =R/20 (small) and

=R/10(large),analytical absolute values of the change in reactivity ( ) have been

obtained

and given below in table below.

1 a 0.066cm    1 f 0.07cm     1 a 0.0007574cm     % 15 . 1 / _a a     1 a 0.007574cm     % 5 . 11 / _a a     0

R

0

R

. ana  

 The absolute values of the changes in

reactivity calculated by PERTURB using DTM and PTM have also been listed as

and respectively in the same table. The per cent errors associated with DTM and PTM with respect to the analytical result are given underneath in paranthesis.

 The PERTURB runs for both DTM and PTM

have been made with an extraorinarily large number of finite elements (320 elements of

constant radial thickness) to annul the effect of discretization errors. DTM   PTM  

 The table shows that the reactivity values

calculated by DTM are almost identical to the analytical reactivities(with errors on the order one in ten thousand), indicating the

convergence of the finite element method.

 When both the size and magnitude of the

perturbation is small (column one in the table), perturbation theory is capable of calculating

reactivities which are quite accurate(with error on the order of one in thousand).

 When both the size and magnitude of the

 When both the size and magnitude of the

perturbation is large , the error in the reactivity computed by PTM is too large (almost 8 %) rendering that approach useless.



The next application we'll consider

involves the TRIGA Mark II Reactor of the

Institute of Energy at Istanbul Technical

University .

 The reactor consists of the A,B,C,D,E,F rings

and the graphite reflector, extending in that order radially outward from the center of the cylindrical reactor.

 All rings consist of a number of fuel cells,water

gaps and graphite element cells. A ring is

simply the central thimble.B ring consists of 6 fuel cells while 11 fuel cells and a water gap constitute the C ring. D ring consists of 17 fuel cells and a water gap while 23 fuel cells and a water gap constitute the E ring. F ring consists of 12 fuel cells, 2 water gaps and 16 graphite element cells. The graphite reflector surrounds the outermost F ring.



Spectrum calculations have been carried

out seperately for the fuel cell, water gap,

graphite element cell and the graphite

reflector using the WIMS-D/4 program

and two group homogenized cross

sections for these have been obtained.



Then, by the use of these cross sections in

simple geometric volume homogenization,

the homogenized cross sections of the

 The modification in this application is the

replacement of one fuel cell by water gap in the B,C and D rings respectively.The

homogenized cross sections for the modified ring is then recalculated.



Finally, PERTURB runs with both DTM

and PTM have been made.To minimize

discretization errors, a large number (560)

of finite elements has been used.



The calculated changes in reactivity is

 The calculated changes in reactivity is given in

the table .

 Fourth column gives the per cent difference

between the second and the third columns.

 The change in reactivity increases as we go to

the outer rings as expected.

 The use of PTM seems unwarranted in this

case. This is quite expected since the

replacement of a fuel cell by water is not a minor modification.