Application of the modified neutron source multiplication method to the prototype FBR Monju

Application of the modified neutron source multiplication method to the prototype FBR Monju

Annals of Nuclear Energy 51 (2013) 94–106 Contents lists available at SciVerse ScienceDirect Annals of Nuclear Energy journal homepage: www.elsevier...

3MB Sizes 0 Downloads 52 Views

Annals of Nuclear Energy 51 (2013) 94–106

Contents lists available at SciVerse ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Application of the modified neutron source multiplication method to the prototype FBR Monju Guillaume Truchet a, W.F.G. van Rooijen b,⇑, Y. Shimazu b, K. Yamaguchi c a

Institut National des Sciences et Techniques Nucléaires, Centre CEA de Saclay, F-91191 Gif-sur-Yvette Cedex, France Research Institute of Nuclear Engineering, University of Fukui, Kanawa-cho 1-2-4, Tsuruga-shi, T914-0055 Fukui-ken, Japan c Japan Atomic Energy Agency, FBR Plant Engineering Center, Shiraki 1, Tsuruga-shi, T919-1279 Fukui-ken, Japan b

a r t i c l e

i n f o

Article history: Received 6 April 2012 Received in revised form 25 July 2012 Accepted 27 July 2012 Available online 9 October 2012 Keywords: Monju ERANOS Sodium cooled fast reactor Mode extraction Reactivity of sub-critical systems Modified neutron source method

a b s t r a c t The Modified Neutron Source Method (MNSM) is applied to the prototype FBR Monju. This static method to estimate the reactivity in sub-critical conditions has already given good results on commercial pressurized water reactors. The MNSM consists both in the extraction of the fundamental mode seen by a detector to avoid the effect of higher modes near sources, and the correction of flux distortion effects due to control rod movement between different reactor states. Among Monju’s particularities that have a big influence on the MNSM factors are: the presence of two californium sources near the core and the position of the detector, which is located far from the core outside of the reactor vessel. The importance of spontaneous fission and (a, n) reactions, which have increased during the shutdown period of 15 years, will also be discussed. In order to evaluate the detector count rate, a propagation calculation has been conducted from the reactor vessel to the external detector. For two sub-critical states, an estimation of the reactivity has been made and compared to experimental data obtained in the restart experiments at Monju (2010). Results indicate a good agreement between the MNSM reactivity and the reactivity measured with other methods. The reactivity dependence of the correction to apply to the point kinetic equation is discussed. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction A fundamental problem exists when trying to estimate the reactivity of a nuclear reactor if one uses a static method, with a constant neutron source in the reactor. Indeed, the well-known relation to calculate the reactivities of different states of the reactor, associating the reactivity ratio to the inverse of a measured count rate ratio, does not take into account the effects due to the change of the flux shape between the different reactor states. The Modified Neutron Source Multiplication Method (MNSM) can be used to reconstruct the link between the reactivity defined with the eigenvalue of the critical system and the count rate observed by a detector placed in a source driven system. It is shown that the correction comes from three physical effects: 1. a spatial effect due to the distortion of the critical flux shape between different states of the reactor (for instance due to control rod movement),

⇑ Corresponding author. Tel.: +81 776 27 9871; fax: +81 776 27 9892. E-mail addresses: [email protected] (G. Truchet), [email protected] (W.F.G. van Rooijen), [email protected] (Y. Shimazu), yamaguchi.katsuhisa@ jaea.go.jp (K. Yamaguchi). 0306-4549/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2012.07.040

2. an extraction effect due to the difference between the critical flux shape and the actual source driven flux shape, 3. an importance effect due to the distortion of the adjoint flux: a neutron released by a source may not have the same importance in the fission chain depending on the reactivity of the reactor. These three correction factors are mathematically introduced in Section 2 of this paper. In Section 3, the method is applied to the prototype FBR Monju. The MNSM correction factors mainly depend on the source and detector positions, and thus the sources and detectors must be characterized with proper accuracy. In order to precisely estimate the intensity of the independent sources coming from californium, spontaneous fission in the fuel and (a, n) reactions, a calculation has been done taking into account the roughly 15 years of aging of the fuel, corresponding to the shutdown period of Monju between December 1995 and May 2010. The necessity to calculate the flux both in the core and at the detector position causes a computer memory problem if one desires to make a model which takes into account the whole geometry. In Monju, the detector is located far away from the core, outside of the vessel, in a nitrogen atmosphere. It has been chosen to treat this problem in two steps: first, a model of the core limited to the primary vessel has been created, and second, a propagation

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

of the flux from the boundary of the vessel to the detector position is done. Two ways of propagating the flux are explored: 1. an analytical way considering nitrogen as void, 2. a propagation based on the adjoint source-detector problem, from the detector to the vessel. The last section presents the results obtained for a detector placed anywhere in the primary vessel, as well as for the actual position of the detector after propagation of the flux. These results are compared with results of dynamic measurements done by JAEA. In general, the MNSM results are different but within the error margin of the measured reactivities. The reactivity dependence of the MNSM correction is also analyzed in order to provide a direct expression of reactivity. 2. The modified neutron source multiplication method Before we start our discussion, a note about definitions. A nuclear reactor is deemed to be critical if, in the absence of an independent source of neutrons, there is a steady state neutron population in the reactor. This implies a balance between the number of neutrons produced by fissions and the number of neutrons lost by leakage and absorption. In calculations one usually does not deal with a critical reactor, and an artificial parameter is introduced into the neutron balance equation so that an artificial steady state is obtained. Mathematically the neutron balance equation is transformed into an eigenvalue equation. The most common choice is to apply the correction factor (the eigenvalue) to the fission term in the neutron balance equation. In such a case, one often speaks of the k-eigenvalue and the corresponding eigenmodes are called k-modes. Related to the k-eigenvalue is the multiplication factor k:

k

1 k

ð1Þ

In the framework of the Point Kinetics Equations (PKEs), there is yet another parameter in use: the reactivity q, which is related as follows:

q

k1 1k k

ð2Þ

Note here that technically one can only speak of a reactivity in the framework of the point kinetics equations. However, in the present paper we will find expressions similar to Eq. (2) and we will loosely use the term reactivity for these. This section presents the three correction factors of the modified neutron source multiplication method using the approach of Tsuji et al. (2003). When a reactor is sub-critical and if there is a constant source of neutrons, the neutron population in the reactor is constant. Such a steady state of the reactor, with a constant reactivity, leads to a constant neutron population in the reactor and this yields a constant count rate in a detector. This paper refers to such a situation as a reactivity state. Let us start with the classic static point kinetic relation in Eq. (3) where q is the reactivity, and M is the count rate. Suppose there are two distinct reactivity states for the reactor (e.g. two different control rod patterns are present in the reactor, to create two distinct reactivity states). The detector count rates are observed in the two reactivity states. Let one of the reactivity states be defined as the reference state (indexed as ref). The other state is an arbitrary reactivity state of the reactor different from the reference state. rd is the detector position. If the reactivity of the reference state is known (for instance by performing a measurement), then the reactivity q of any other state can be calculated from the observed count rate in that other state, as follows:



Mref ðrd Þ q Mðrd Þ ref

95

ð3Þ

Such a reactivity is called the estimated reactivity. The spatial distortion of the flux as a function of the reactivity is not taken into account in Eq. (3), which is defined for a point reactor model. The MNSM Method will introduce three coefficients to correct this equation. The MNSM equation is given in Eq. (4), with Cim the importance correction, Csp the spatial correction and Cext the extraction correction.

q ¼ C im C sp C ext

Mref ðrd Þ q Mðrd Þ ref

ð4Þ

 The importance correction Cim reflects the difference of the importance between a source neutron in state ref and the state presently under consideration. For example, in a deeply sub-critical state the source neutrons are relatively more important than in a nearly critical state.  Csp reflects the change of the flux distribution due to the physical change in the reactor between the reference state and the state under observation, e.g. the insertion of one or more control rods.  Cext is necessary to extract the fundamental mode component from the actual (source driven) flux distribution in the reactor. The fundamental mode depends on the state of the reactor, and thus the fundamental mode will change (slightly) between the reference state and the state under observation. Note that the reactivity in the MNSM approach is not the same concept as the reactivity in the PKE, as will be demonstrated later. The task at hand is now to estimate, given some arbitrary sub-critical state of the reactor, the multiplication factor, which is defined using the eigenvalue k = 1/k of the criticality equation (Eq. (5)), given that the reactor is source-driven, i.e. the flux corresponds to the solution of the source-driven equation (Eq. (7)). Introduce the forward criticality equation:

L/c ðrÞ ¼

1 c F/ ðrÞ k

ð5Þ

It is adjoint equation:

Ly /cy ðrÞ ¼

 y 1 Fy /cy ðrÞ k

ð6Þ

and the source-driven equation:

L/s ðrÞ ¼ F/s ðrÞ þ sðrÞ

ð7Þ

In these equations, L is defined as the destruction and scatter operator, F as the production (fission) operator and s as the independent source. We use the multi-group formalism, i.e. the operators L and F are matrices, and the source s and flux /c,s are vectors. The superscript c indicates critical, i.e. Eq. (5) and the superscript s indicates source-driven, i.e. Eq. (7). It is known that Eq. (5) has a spectrum of eigenvalues ki ¼ k1 and corresponding eigenfunctions i /ci . The smallest (in magnitude) eigenvalue k1 is the fundamental eigenvalue, the corresponding 1/k1 is the multiplication factor, and the corresponding eigenvector /c1 is the fundamental mode. All discussions of the multiplication factor of the (sub-critical) reactor deal with k1 and /c1 . We now introduce the assumption that the eigenfunctions form a complete, orthogonal basis (there are some non-trivial technicalities involved with this assumption, see Section 2.1 for more details). Then, we can write:

/s ðrÞ ¼ A1 /c1 ðrÞ þ

1 X Ai /ci ðrÞ i¼2

ð8Þ

96

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

If we can determine A1, then we know how much the fundamental mode contributes to the source-driven flux /s. We can operate with F on both sides of this equation, and take the scalar product with the adjoint fundamental mode /cy 1 . We find:



s /cy 1 ; F/



¼



c A1 /cy 1 ; F/1



1 X   c Ai /cy þ 1 ; F/i

ð9Þ

i¼2

  c As discussed in Section 2.1, /cy 1 ; F/i ¼ 0 for i – 1, we find for A1:



 s /cy 1 ; F/ A1 ¼  cy  /1 ; F/c1

ð10Þ

We can establish a relation between A1 and the multiplication factor k1 as follows. Take the scalar product of the source-driven equation (Eq. (7)) with the fundamental mode adjoint, and use the properties of adjoint operators:



  s y cy   cy  /s ; Ly /cy ¼ / ; F /1 þ /1 ; s 1 Since k1 



y k1 ,

ð11Þ

we can write:

    cy  1 s  1 /cy ¼ /1 ; s 1 ; F/ k1

ð12Þ

 cy    cy  / ;s / ;s k1 1 A1 ¼  cy1 c  ¼  s  cy1 c  1  k1 /1 ; F/1 q /1 ; F/1

ð13Þ

where we can loosely define qs as the MNSM reactivity, corresponding to the multiplication factor of the sub-critical reactor in the source-driven state. Given our knowledge of A1, we can rewrite Eq. (8) as follows: 1 X Ai /ci ðrÞ

ð14Þ

i¼2

with /s1 ðrÞ the fundamental mode component of the source-driven flux, defined by:

   /c1 ðrÞ 1   s /cy  1 ; s  cy q /1 ; F/c1

ð15Þ

Suppose one has a detector at location rd typified by a cross section Rd. If the reactor is sub-critical with a constant source, the count rate will be:

Mðrd Þ ¼ Rd /s ðrd Þ

ð16Þ

If only the fundamental mode would exist in the reactor, the count rate would be:

M 1 ðrd Þ ¼ Rd /c1 ðrd Þ

ð17Þ s

/c1 ðrd Þ.

One can calculate / (rd) and Thus, given a count rate observed in an experiment M, the contribution from the fundamental mode can be calculated from:

M1 ¼

ð20Þ

Thus, the MNSM reactivity of an arbitrary state l of the reactor is found:

3 2 R /c  d 1;l 2  cy 3 7 6 /cy ;F/c 7 1;l 1;l qsl 4 /1;l ; s 5 6 7 6 ¼  c s 6 Rd /1;ref 7     cy qref 4 5 /1;ref ; s cy |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} /1;ref ;F/c1;ref |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} C im 2

cy /1;ref ;F/sref

C sp

3

c

R / 6 d 1;ref 7 6 cy Mref 7 7 M 6 /1;ref ;F/c1;ref 7 6 ref 6  7 7 Ml 6 cy s c / ;F/ 7 6 R / l 4  1;l d 1;l 5 /1;l ;F/c1;l

ð21Þ

Ml

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} C ext



/s1 ðrÞ ¼ A1 /c1 ðrÞ ¼

  cy c  s  R /cy ; s d CM /c1;l /1;ref ; F/1;ref 1;l qref Ml ðrd Þ 1;ref   ¼ M Mref ðrd Þ qsl Rd /cy ; s /c1;ref /cy ; F/c C 1;l 1;l 1;ref 1;l

cy

from which we find that:

/s ðrÞ ¼ /s1 ðrÞ þ

Suppose we observe the count rate in some other, arbitrary subcritical state of the reactor (call it state l), and express the ratio:

/c1 ðrd Þ M ¼ CM 1 ðrd ÞM /s ðrd Þ

ð18Þ

Now suppose that the reactor is in a given state, which we will label the reference state, and suppose that the reactivity of this state is known (from some experiment or some other reliable source). Then the fundamental mode contribution to the count rate is found as: s M1;ref ðrd Þ ¼ C M 1;ref ðrd ÞM ref ðrd Þ ¼ Rd /1;ref   /c1;ref ðrÞ 1  /cy ¼ Rd  s 1;ref ; s c qref /cy 1;ref ; F/1;ref

ð19Þ

The three MNSM factors indicated in Eq. (21) can be calculated with a computer based on an adequate reactor model. Note that several calculations need to be done to obtain the fluxes appearing in the scalar products, namely forward and adjoint eigenvalue calculations in the reference state and the perturbed state, as well as a source-driven calculation of the flux in the reference and the perturbed state. Note also that if the reference and perturbed states are due to control rod movement, a code needs to be used which is capable of explicitly representing such control rod movement in the calculational mesh. One may wonder if the estimated reactivity with the MNSM method is really different from the reactivity as introduced in the PKE. One may argue that one needs the reference reactivity qsref , and if this reference reactivity is determined from measurements assuming the validity of the PKE, then the MNSM reactivity is nothing but an improved PKE-reactivity. There are two points to take into account, though. First, the reference reactivity may be measured in a configuration for which the PKE are valid, but the MNSM-reactivity remains valid even for reactor states where the PKE are not valid. Secondly, the reference reactivity may be derived from other sources than the PKE. For example, one could imagine high-fidelity Monte-Carlo calculations to determine the multiplication factor of some reference state of the reactor, and from that infer a reactivity. The MNSM-reactivity is thus a more general concept than the PKE-reactivity, even though the two will coincide in certain cases. That being said, one may in fact write the entire MNSM-approach in terms of k1 rather than q, and not refer to PKE at all. 2.1. Pesky mathematical details In the derivation of the MNSM factors, we have used some properties of the flux, most notably the fact that the eigenmodes are a complete set. This point is briefly discussed by Bell and Glasstone (1985, p. 45, symbols as in original text) for transport theory: It should be noted that the k eigenfunctions are not a complete set of functions for expansion of solutions of the transport equation. In some cases of one-speed problems, however, it has been found that when the k eigenfunctions are integrated over X, they do form a complete set for expanding functions of r only. and in the discussion of the multi-group diffusion equation (Bell and Glasstone, 1985, p. 189, symbols as in original text):

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

There are presumably other eigenvalues, smaller in magnitude than k0, for which the corresponding eigenfunctions are sometimes negative or even complex, but no single one of these higher modes can be realized physically. Although these higher modes can be found explicitly for simple cases, such as the onespeed approximation in simple geometry, little is known about such modes in general. In (Henry, 1975, p. 347, symbols as in original text) we find the following: Little is known about the nature of the eigenvalues kn. In particular it is not known whether the eigenfunctions belonging to these eigenvalues form a complete set (in the sense that any function obeying the boundary and continuity conditions of W can be expanded in terms of them). It would be an advantage to be able to prove completeness, since formal proofs of various properties of solutions to the transport equations could then be obtained using eigenfunction expansions. It is known that the eigenfunctions belonging to the eigenvalues kn are not complete unless the eigenvalue k = 0 is included. Moreover there are an infinite number of linearly independent eigenfunctions Wð0Þ n corresponding to that eigenvalue. They are solutions of

M Wnð0Þ ¼ 0

ð22Þ

where the subscript n for this case refers to independent eigenfunctions of the particular eigenvalue zero. Despite the lack of positive evidence, it is current practice to assume that, with the Wð0Þ included, the set of eigenfunctions n are complete. Thus it is assumed that any arbitrary function F(r, X, E) obeying the continuity and boundary conditions imposed on the directional flux can be expanded as

Fðr; X; EÞ ¼

1 X

an WnðkÞ ðr; X; EÞ þ

n¼0

1 X a0n Wð0Þ n ðr; X; EÞ

ð23Þ

n¼0

where the parts of the expansion involving eigenfunctions belonging to nonzero and zero eigenvalues have been written separately. In (Henry, 1975), M is the fission operator. Henry discusses the issue in terms of transport theory. From his comments, it is clear that the completeness of the sets of eigenfunctions is not proven mathematically, but it is assumed. Therefore we feel somewhat confident that the series expansion is allowed. As for the equivalence of the fundamental forward and adjoint eigenvalues, we refer to Williams (1986) and Henry (1975). Note that for the higher order modes it is in general not possible to show that the forward and adjoint eigenvalues are the same. One final point needs to be made, and that is that in general:



 c /cy 1 ; F/i ðrÞ ¼ 0;

i–1

ð24Þ

Again, we refer to (Henry, 1975, p. 348, symbols as in original text), where it is stated that:

D

E

D

E

ðkÞ ðkÞ  ðkÞ Wm jMjW0 ¼ 0 ¼ WðkÞ ; n jM jW0

m; n – 0

ð25Þ

due to the physical properties of the fundamental forward and adjoint fluxes (i.e. the fundamental forward and adjoint fluxes are non-negative everywhere in the calculational domain). 3. Monju application On 8 May 2010, criticality was reached of the prototype FBR Monju. Between May and August 2010, various reactor physical parameters were measured. Details of this measurement campaign

97

and results of the experiments are detailed in various publicly available sources (JAEA, 2011; Jo et al., 2011) These measurements included the determination of the reactivity of the reactor in various sub-critical configurations. The control rod worth curve of the central control rod was measured using a period method. This method assumes the validity of the familiar point-kinetics equations to establish a link between the reactivity and the measured period of the flux increase or decrease. All other control rod worth curves were measured with the compensation method, using the central control rod as a calibrated rod (JAEA, 2011). In our work, the control rod worth of the central control rod is investigated. The locations of the detector and of the sources are essential to determine MNSM factors. The Monju model has been prepared based on publicly available sources (Takashita et al., 2000; NISA, 2007). The model takes into account the actual 3D hexagonal geometry and uses a loading pattern of five types of fuel subassemblies, whose enrichments and Pu-vectors have been estimated from publicly available documents (Taro et al., 1994; Takashita et al., 2000; NISA, 2007). Actual composition data are classified. It should be noted that the fuel was irradiated to approximately 30 Equivalent Full Power Days (EFPDs) during the first set of commissioning tests of Monju in 1994–1995. This depletion history is not taken into account in the present work for two reasons: first, the fuel composition of the original fuel is not known very accurately anyway (rather, it is estimated based on available sources), and secondly, many fuel assemblies were reloaded or refueled prior to the criticality of May 2010. Modeling is done with the well-established ERANOS 2.0 code (Rimpault et al., 2002). Since the compositions of the fuel have been determined for each type at a certain production date, a dedicated module in ERANOS was used to simulate the aging of the fuel until May 2010. See Fig. 1 for an overview of the May 2010 core. The VARIANT module (Palmiotti et al., 1995) of ERANOS has been used for the analysis at the core level. The full core calculations are done with diffusion theory and six neutron energy groups. In order to take into account the deep absorption in the center of the control rods, the control rod sub-assemblies have been homogenized following the reactivity equivalence method described in the ERANOS manual (Rimpault et al., 1997). The aim of this method is to get the same reactivity in the homogeneous cell as the original control rod cell surrounded by fuel. Such a treatment is necessary to obtain a good correlation between the control rod position and the reactivity variations. For a given insertion depth of the central control rod (all other rods being completely withdrawn from the core), the reactivity calculated by ERANOS is, unsurprisingly, slightly different from the experimental reactivity (difference of several hundred pcm, attributed to modeling errors and imprecise knowledge of the fuel composition). Therefore, slight adjustments to the cross-section data have been made to obtain the same reactivity in the calculations as experimentally observed, i.e. q = 0 is calculated with the control rod in the simulation at the measured critical position. With the reactivity equivalence method, it was possible to reproduce with ERANOS the measured control rod worth curve of the central control rod with good accuracy (JAEA, 2011). Of course, there are several sources of errors that should be taken into account, as well as their potential impact on the results of the analysis. As stated before, the actual fuel composition of the Monju 2010 core is classified. More detailed information about the fuel composition and how it was estimated is given in Tamagno and Rooijen (2012). Results from that work certainly give credibility to the core model. The depletion of the fuel during the trial operation of Monju in 1994/1995 is neglected, but the reactivity decrease due to the decay of Pu-241 during the 15 years of shutdown is much more important, and this decay is taken into account. The six energy groups of the core calculations use the

98

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

Fig. 1. Monju loading 2010. The core loading has a

2 3

p symmetry, but the presence of the two californium sources breaks this symmetry; symmetry over p remains.

recommended energy structure for Super-Phénix, according to the ERANOS manual. Diffusion theory is used because the original paper of Tsuji et al. (2003) uses diffusion theory for a much smaller reactor, and in that paper good results are obtained. Of course, to obtain k = 1.000 in our model with the control rod positions from JAEA (2011), some adjustments were needed. However, for the present research the most important point is that the calculation of the multiplication factor as a function of the central control rod position is correct. The fact that the control rod curve from JAEA (2011) could be calculated satisfactorily with our model gives credibility to our work. Besides, even if a bias is present in the model, it is not really a problem because we are dealing with reactivity ratios and reactivity differences, which are relatively unaffected by modeling bias. We believe that our modeling is reasonable given the objectives of the study. We could have used more ‘‘high-fidelity’’ modeling options, but they may give a false sense of accuracy. The results are not only dependent on the core modeling options, but also on the modeling of the detectors, the sources, and the propagation of the flux, all of which have some associated approximations and errors. Besides, the correction factors are calculated as scalar products over the entire reactor space, so that ‘‘local’’ errors may be averaged out.

using ORIGEN-S by considering that a-decay occurs with the MOX-fuel as surrounding matrix. The majority of (a, n) reactions happen on O-18. Elements that contribute the most to this production of neutrons are presented in Fig. 2b. From Fig. 2, one can deduce that the Am-241 (produced from the decay of Pu-241 over a period of 15 years) presently represents 23% of the total distributed neutron source. The distribution of all neutron sources in the reactor is plotted in Fig. 3. It must be noted that the total neutron source in the fuel is equal to 1.1  109 n/s; more than five times the californium source strength! This high ratio is due to the decay of the californium sources which have lost about 80% of their strength between 1993 and 2010.

3.2. Detector position and energy response cross-section Now that the independent sources have been calculated, it is required to analyze the detectors in the Monju reactor and to

3.1. Source evaluation There are three types of sources in Monju. Firstly, there are two sources, see Fig. 1, each one made of 2.81 mg of pure Cf-252 on 12 October 1993, and located in the blanket region. Secondly there is the source due to spontaneous fissions of the fuel, and finally the source due to (a, n) reactions, mainly on oxygen. The two last kinds of source are present in the fuel sub-assemblies and are called distributed sources in this paper. Let us first evaluate the Cf-252 source. Cf-252 decays into Cm248 (TCf252 = 2.645 years). Given the production date, 12 October 1993, one finds an activity of 727 MBq in May 2010. A Cf-252 decay can be either alpha decay (96.91%) or spontaneous fission (3.09%). The average number of neutrons released upon spontaneous fission is 3.756 (JEFF2.2). Thus one finds a source strength of 84.6  106 n/s due to spontaneous fissions in the Cf-sources in May 2010. This value has been confirmed using ORIGEN-S (Gauld et al., 2009) and the energy spectrum of the neutron released from the californium source has been condensed into the six energy groups used in the core calculations. Spontaneous fission occurs in the fuel as well and was evaluated with ORIGEN-S. Fig. 2a indicates the relative contribution of isotopes to the total neutron source due to spontaneous fissions. Thirdly, a distributed source of neutrons is due to (a, n) reactions in the fuel. The strength of this source has been determined

Fig. 2. Neutron sources in Monju’s fuel, type Inner Fuel 1; the most important contributions from individual isotopes are indicated.

99

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

M ¼ Rd /d ¼

G X

Rdi /di

R3

92

calculate their energy response. BF3-detectors are located outside of the vessel, in a nitrogen atmosphere. Fig. 4 illustrates the two kinds of detectors in use: the Source Range Monitor (SRM) and the Experimental Nuclear Instrumentation System (E-NIS). Note: The count rate measured with the E-NIS detector is less accurate than the count rate measured with the SRM detector. Therefore, this paper focuses only on data measured with the source range monitor. It is required to evaluate the detector cross-section Rdg, i.e., the count rate of the detector if a neutron of energy group g hits the surface of the detector. The detector geometry is presented in Fig. 5. To determine Rdg, a point source located very far from the detector in a void medium is considered. The number of source neutrons is distributed uniformly over an energy group g and the angular distribution is set so that all neutrons will propagate towards the detector. Then the number of (n, a) reactions that occurs due to all source neutrons of energy group g is calculated. This gives the detector response cross-section Rdg. It is assumed that the detector count rate is linearly proportional with the number of (n, a) reactions. Using the Monte Carlo code TRIPOLI (Petit et al., 2008), Rdg has been determined for each energy group g. Then if the group-wise flux /g is known at the detector position, one can calculate the total count rate as follows:

ð26Þ

i¼1

Fig. 6 provides a global view of the collision points in the detector. All outer boundaries are vacuum boundaries (neutrons are lost if they cross the outer boundary). The fact that only a point source has been used and not the actual flux corresponding to each configuration of the reactor certainly does not constitute a problematic hypothesis as long as the flux /d is correctly evaluated (taking for example into account the fact that a neutron may scatter somewhere outside of the detector, come back to the detector and make a contribution to the count rate). The SRM detector energy response cross-section is given in Table 1. The SRM detector gives a high count rate for fast neutrons: the graphite moderates the neutrons in order for them to have a high probability to cause an (n, a) reaction. 3.3. Propagation methods From Eq. (21), it is clear that one needs to know the flux at the detector position. Modeling the reactor and the surrounding volume up to the detector position is unfeasible with VARIANT in ERANOS. For this reason, the reactor and surrounding areas up to

Fig. 4. Overview of Monju detectors for flux measurements. Indicated are detectors used in the present work, i.e. one of the Source Range Monitors (SRMs) and the ENIS detector. Also indicated is the position and orientation of the core relative to the detectors.

the vessel have been modeled in ERANOS. Then some calculation is necessary to propagate the neutrons from the vessel surface to the detector position. Two ways of propagating the flux are described in this paper: (1) an analytical propagation and (2) an adjoint calculation from the detector to the core vessel. 3.3.1. Analytical propagation The analytical propagation assumes that the concrete walls (see Fig. 4) do not reflect neutrons, and that the neutrons can only propagate directly from the vessel to the detector. In Fig. 7, the total flux calculated by VARIANT and the energy spectrum in six groups is plotted, from the center of the core to the vessel boundary. The energy spectrum shows that the energy spectrum at the vessel surface is dominated by the lowest energy group. Now let us move on to the equations used to propagate the scalar flux, given by the ERANOS calculation, from the vessel boundary to the detector location. An illustration of the analytical propagation is given in Fig. 8. Considering the transport equation in the nitrogen (nitrogen as void):

~  rW þ 0 ¼ 0 X

ð27Þ

From this equation and for any point ~ r in the nitrogen:

~Þ ¼ Wð~ ~; X ~Þ Wð~ r; X r  RX

Fig. 3. Neutron sources in the core of Monju. Clearly visible are the two californium sources, as well as a rather homogeneously distributed source in the fuel assemblies.

ð28Þ

R is the distance along the flight path of the neutrons. Considering that neutrons at the vessel boundary are low energy neutrons, they have undergone many scattering reactions and have therefore a uniform probability to have any direction. Of course once the neutrons propagate into the nitrogen there is no possibility to scatter back to the vessel. From this interpretation, it has been assumed that the outgoing angular flux is constant as a function of X in each

100

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

Fig. 5. Detector geometry for the determination of Rdg (SRM detector).

Table 1 Response of SRM detector using TRIPOLI. Group

Rdi

Uncertainties (%)

1 2 3 4 5 6

2.94E003 2.99E003 1.33E003 1.31E003 1.32E003 1.28E003

0.13 0.16 0.16 0.15 0.16 0.20

point r 2 S, S being the ERANOS boundary and n the normal at this point to the surface:

/s ð~ rÞ ¼ ¼

Z 4p

ws ðr; XÞdX

Z

ws ðr; XÞdX þ

XnP0

Z

ws ðr; XÞdX ¼ 2pws ðrÞ þ 0

ð29Þ

Xn<0

where /s(r) is the scalar flux taken from ERANOS. In order to obtain the contribution of all nodes on the surface of the VARIANT model to the detector response, one has to sum all the contributions:

/d ¼

Fig. 6. Collision points in the SRM detector. The source of neutron is placed on the right of the drawings.

point on the boundary of the ERANOS model. A mathematical interpretation of this half isotropic phenomenon is given in Fig. 8. For a

1 Vd

Z r0 2V d

Z

/s ðrÞ r2S

2pkr0  rk2

dS dV

ð30Þ

3.3.2. Adjoint source-detector propagation This time the effect of the concrete shielding on the detector count rate will be taken into account. One can understand that if there is reflection on the concrete walls, scattered neutrons may also contribute to the detector count rate. First, we show how the reflection can affect the final count rate of the detector as a function of the reactivity of the reactor. Refer to Fig. 9. In a deep sub-critical state the flux is highly peaked around the sources. The consequence is a rather non-uniform flux on the vessel boundary. On the other hand, in a critical state, the flux distribution on the vessel is more or less uniform. The difference of the neutron distribution on the vessel implies that the contribution of scattered neutrons to the detector count rate is a function of the

101

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

Fig. 7. Radial evolution of the energy spectrum (left), and the total flux (right), from the center of the core to the vessel boundary. Note that on the vessel boundary almost all neutrons are in group 6, i.e. at a relatively low energy because of scattering collisions.

Ω.n < 0 ⇒ Φ(Ω) = 0 Ω.n .n > 0 ⇒ Φ(Ω) = c2

0 0 200 200 Fig. 8. Propagation scheme from the vessel boundary to the detector position. Illustrated is the geometrical model used in the VARIANT simulations, i.e. a hexagonal geometry. The actual vessel is cylindrical.

reactivity. Therefore, if reflections are possible on the concrete walls, it is important to know from which point on the vessel surface a neutron originates before it interacts in the detector. Using TRIPOLI, this effect has been analyzed in 2D. To find which parts of the vessel contributes most to the count rate of the detector, the basic idea is to use adjoint source-detector problem, i.e.: start neutrons in the detector and propagate them to the vessel. To solve this kind of propagation, TRIPOLI uses a pseudo-adjoint calculation. By calculating a direct propagation from the vessel to the detector and saving the history of the neutrons, one can reconstruct the adjoint flux (because initial positions of the neutrons are saved, it is possible to determine from which point of the vessel a detected neutron had been emitted). The vessel has been angularly divided in 16 emission regions in TRIPOLI and neutrons reaching the blue surface of Fig. 9 are counted. The importance of each of the 16 zones is given in Fig. 10. These results indicate that a neutron emitted anywhere on the vessel may reach the detector. In the case of the analytical propagation discussed before, the emitting part of the vessel has been constrained between an angle of 210° and 270°. Therefore only 64% of neutrons that actually interact with the detector are taken into account. By weighting the flux coming from ERANOS with the importance function, one will obtain a value of /d.

Fig. 9. TRIPOLI model for the adjoint propagation.

102

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

4. Results 4.1. Measurements Three reactivity states have been chosen corresponding to three control rod positions, see Fig. 13 and Table 2. The Coarse Control Rod 1, located at the center of the core, is withdrawn progressively. Reactivities indicated in Table 2 are deduced by JAEA using the total control rod worth, measured using a kinetic method (Jo et al., 2011). Associated static count rates have also been provided by JAEA and are also given in Table 2 (JAEA, 2011). To express the MNSM reactivity, the MNSM correction factors must be determined. This determination is done with the VARIANT module of ERANOS. This code calculates the necessary flux distribution to calculate the MNSM factors, using diffusion theory and six energy groups. If the correction factors are needed for a detector outside of the vessel, two ways to propagate the flux to the desired position have been described in Section 3.3. The actual Monju detector is located outside of the vessel and propagation is necessary. After presenting briefly some flux distributions obtained with VARIANT, this section will analyze the different effects observed on the correction factors inside of the Monju vessel (for a hypothetical detector inside the vessel). Then using the propagation methods described in Section 3.3, MNSM coefficients and estimated reactivities are compared with the measured reactivities provided by JAEA. The influence of the source strength is also discussed. To conclude this section, the evolution of the correction factors is analyzed as a function of reactivity in order to directly express the estimated reactivity as a function of the count rate ratio. 4.2. Flux solution As illustrated in Fig. 11, the source driven flux in energy group 1 (high energy) peaks locally around the sources, and decreases at the CCR1 position (the exact center of the core). In Fig. 12, the critical flux in group 6, i.e. the thermal group, is represented. It is seen that the flux is locally thermalized by control rod channels (filled with sodium) and shield regions (mainly constituted of steel). 4.3. Correction factors inside the vessel c s Because /cy 1 ; /1 and / are known everywhere inside the core, it is possible to imagine a hypothetical detector placed anywhere. Using these fluxes, the relevant correction factors in a horizontal or in a vertical cut of the Monju reactor can be plotted. Fig. 13 shows the spatial dependence of the spatial correction factor Csp. It is seen that the spatial correction due to physical changes inside the reactor is only relevant close to the location of the perturbation (i.e. the control rod in this case). The movement of the absorbent

Fig. 10. Normalized contribution of flux on the vessel surface to the flux at the detector position (importance function) as a function of the angular region on the vessel surface.

Table 2 Count rates data, for three reactivity states of Monju. The different states corresponds to different insertion depths of CCR1.

CCR1 (mm) Reactivity ( ) SRM⁄ detector (cpm) Error (cpm) E-NIS⁄ detector (cpm) Error (cpm)

State 2

State 1

Reference

488 3.0 ± 0.4 1367.37 2.61 168.38 0.92

455 19.9 ± 0.5 226.3 1.06 27.33 0.37

396 50.8 ± 0.9 91.56 0.68 10.74 0.23

part of the control rod distorts axially the Csp, as seen in Fig. 13. It basically means that the correction for a detector placed very close to the control rod would be important. The correction Csp to apply elsewhere is nearly a constant. It means that the spatial distribution of /c1 between two states is the same far away from the perturbation caused by the control rod. Even if the spatial distribution of /c1 are the same between two reactivity states, it is not the case of their amplitudes, which are inversely proportional to the reactivity. The magnitude of Csp far away from the perturbation is therefore never equal to unity (except when the two states considered are the same). The importance coefficient Cim is constant everywhere, and in our case, calculation shows that it is very close to unity. It reflects the fact that /cy 1 does not change much with the chosen reactivity. The Cext coefficient reflects the dissimilarity between /c1 and /s. Since this dissimilarity is maximal at the source position, it explains the two dips observed in the spatial representation of Cext in Fig. 14. The SRM detector is located in a region affected by Cext because it is behind a source (looking from the core). In Fig. 15, the contour graph of the Cext has been superposed with the geometry of the core. The correction factor for the E-NIS detector can be calculated directly using the value of ERANOS whereas the flux needs to be propagated to reach the SRM detector. 4.4. Correction factor for the SRM detector In order to compare the measurements presented in Section 4.1 with the correction factors, it has been chosen to calculate fMNSM and its error DfMNSM:

fMNSM ¼ C im C sp C ext ¼ qq MM ref ref M DfMNSM ¼ Dq q Mref þ Dqref q2qM M ref ref

þ DM q þ DMref qM2 M q qref Mref ref ref ref ð31Þ

The global factor fMNSM calculated from the JAEA measurements M , the one obtained with the analytical propagation is denoted fMNSM Prop is denoted fMNSM ; the one obtained with the adjoint propagation is Adjoint denoted fMNSM . Finally, a last fMNSM corresponds to the estimation of the global factor at the boundary of the ERANOS model: the average of fMNSM on the surface C of Fig. 15. This third fMNSM is deERANOS noted fMNSM . All results for the SRM detector are gathered in Table 3. A difference is observed between the measured and calculated coefficients, but this difference does not exceed the experimental uncertainties. Table 4 shows the correction factors for the hypothetical case of a stronger californium source. The influence of the source is as expected: a stronger source will cause a larger distortion between the source driven flux and the critical flux (for the same sub-criticality) and thus the correction factors are expected to be larger, as attested by the results in the table. The corrections calculated using different kinds of propagation are very similar. In fact, even if no propagation is taken into account, results stay close to each other.

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

103

Fig. 11. Source driven flux, group 1 (high energy), horizontal cut through the midplane of the core. Note the presence of the two californium sources emitting highly energetic neutrons. Also notice the flux depression in the core center due to the inserted central control rod.

Fig. 12. Critical flux, thermal group (group 6), horizontal cut through the midplane of the core. Note the effect of the sodium-filled control rod channels. The thermal flux in the outer core and blanket region is reduced due to absorption.

4.5. Evolution of the correction factor as a function of reactivity If one knows the evolution of the total correction factor fMNSM as a function of the reactivity, one can express the correction for any reactivity state. By plotting the correction factors as a function of the reactivity, let us see if an evolution law can be extracted. In Fig. 16, the total correction factor as a function of reactivity has been plotted. Note that a further point has been added which corresponds to the reactivity state with CCR1 fully inserted: q = 1.83$ (calculated value, corroborated by the results published in Jo et al. (2011)). Whatever propagation considered, the total coefficient fMNSM evolves more or less linearly with reactivity. Please note that the reactivities (expressed in pcm) in Fig. 16 are calculated by ERANOS and correspond to the control rod positions given in Table 3.

Let us compare the slight differences observed between the propagation methods. The best straight line is obtained for the evoProp lution of fMNSM : Prop fMNSM ðqpcm-unit Þ ¼ 0:982  0:0001qpcm-unit

r2 ¼ 0:99903

ð32Þ

Adjoint The worst straight line is obtained for the evolution of fMNSM . Even if differences are small and might not be taken into account between the propagation methods, the nonlinear effect of the reflections of the flux on the concrete walls of the containment is observed in Fig. 16. If Eq. (32) is satisfied, one can also rewrite the expression of the reactivity:



0:982 MrefMqref 1 þ 104 M

M ref qref

ð33Þ

104

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

Fig. 13. Spatial coefficient Csp: vertical cut along line B–B indicated in Fig. 12. Indicated also is the position of the central control in the three reactivity states (reference states, state 1 and state 2). The surface on the left represents Csp for state 1. Note that this factor has a large magnitude very near the location of the perturbation. Between the reference state and state 1, the control rod is withdrawn. Thus the flux below the control rod will increase and the flux at the top of the control will decrease, as indeed reflected in this graphical representation of Csp.

Fig. 14. Extraction coefficient Cext: horizontal cut through the core, going through the midplane of a source. Note that the magnitude of this coefficient is influenced by the presence of the source.

It is seen that the correction goes from 0.98 for CCR1 fully extracted, to 1.04 for CCR1 fully inserted. In that case there will be 6% difference between a control rod worth estimated without correction and a control rod worth estimated with MNSM. It should be noted that in the present measurements, the difference between the classic estimation of the reactivity and the MNSM-estimation is around 10% in state 2, and thus the 6% correction on the ‘‘full stroke reactivity’’ is negligible to the measurement (10% over a much smaller stroke). However, given that the MNSM is a general method, it is always the preferred method to measure the reactivity of a sub-critical reactor.

Fig. 15. Extraction coefficient Cext: horizontal cut going through the midplane of the source. Superimposed is the overall geometry. It is seen that the E-NIS detector is in the ‘‘wake’’ of the source, but the SRM detector should have less of an influence due to the source. Also note that if the detectors would be moved clockwise, the influence of Cext would be negligible.

105

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106 Table 3 Correction factors for the SRM detector (channel 1). SRM detector

State 2

State 1

Reference

CCR1 (mm) Reactivity ( ) M fMNSM

488 3.0 ± 0.4 0.882 ± 0.124

455 19.9 ± 0.5 0.968 ± 0.051

396 50.8 ± 0.9 1.000 ± 0.018

Prop fMNSM

0.982

0.989

1.000

Adjoint fMNSM

0.983

0.989

1.000

ERANOS fMNSM

0.987

0.992

1.000

ERANOS

1.002 ± 0.002

1.001 ± 0.001

1.000

spERANOS

0.989 ± 0.002

0.993 ± 0.001

1.000

0.996 ± 0.002

0.997 ± 0.001

1.000

C im C

ERANOS

C ext

Table 4 Correction factors for the SRM detector (ch1) and a Cf-source strength of 1E9 n/s instead of 84.6E6 n/s. SRM detector

State 2

State 1

Reference

CCR1 (mm) Reactivity ( ) M fMNSM

488 3.0 ± 0.4 0.882 ± 0.124

455 19.9 ± 0.5 0.968 ± 0.051

396 50.8 ± 0.9 1.000 ± 0.018

Prop fMNSM

0.954

0.975

1.000

Adjoint fMNSM ERANOS fMNSM

0.961

0.976

1.000

0.963

0.981

1.000

ERANOS

1.000 ± 0.003

1.000 ± 0.001

1.000

spERANOS

0.989 ± 0.003

0.993 ± 0.001

1.000

0.977 ± 0.003

0.988 ± 0.001

1.000

C im C

ERANOS

C ext

Fig. 16. Plot of the core reactivity and the (calculated) magnitude of the corresponding MNSM correction factor. For the three reactivity states considered in this paper (to the right part of the figure), the MNSM factor is nearly unity. If the central control rod is completely inserted, the correction becomes larger.

5. Conclusion In order to improve the evaluation of the reactivity in a sub-critical state, the MNSM method was successfully applied for the first time to the Monju reactor. A model in two steps has been built: the ERANOS code has been used to model the core up to the primary vessel and a propagation of the flux is then realized to determine the flux at the detector location. TRIPOLI has been used to calculate the energy response of the SRM detector of Monju in order to express the count rate. During the shutdown period of 15 years, the Cf-sources have decreased in strength while the source due Am-241 has increased (Am241 decays by a-emission which in turn causes neutron production

from the (a, n)-reaction on O-18). In May 2010, the distributed neutron source was five times stronger than the Cf-sources. The spatial evolution of the MNSM factors has been analyzed everywhere in the prototype FBR Monju. The shape of these correction factors has been plotted and may be used to improve source and detector placement in a reactor. The reactivity determined with MNSM is somewhat larger in magnitude compared to measured data, but within the error margins of the experiment. Because of the source-detector configuration of Monju with the detector far away from the core and the low contribution to the count rate of neutrons of the Cf-sources, the correction factors are close to unity. The measured reactivity was obtained with a kinetic method, i.e. a method assuming that

106

G. Truchet et al. / Annals of Nuclear Energy 51 (2013) 94–106

the point kinetics equations apply. As stated before, the MNSM takes into account the fundamental mode extraction, i.e. does not necessarily rely on the validity of point kinetics. The MNSM-reactivity is a more general concept than the PKE-reactivity, and it is therefore our opinion that the MNSM method yields the best estimate of the reactivity. For the present work, the reactivities estimated with the MNSM are therefore interpreted as being the correct reactivities for the sub-critical states of the Monju reactor. Acknowledgments The work presented in this paper was performed in the scope of an internship for G. Truchet’s M.Sc. degree through an agreement between INSTN, CEA, JAEA, and the Research Institute of Nuclear Engineering, University of Fukui. The authors would like to acknowledge the financial support from the European Nuclear Education Network (ENEN), Electricité de France (EdF), as well as the material support from JAEA to enable the stay of G. Truchet in Japan. G. Truchet thanks T. Hazama of JAEA for his help and assistance. References Bell, G.I., Glasstone, S., 1985. Nuclear Reactor Theory. Robert E. Krieger Publishing Company. Gauld, I., Hermann, O., Westfall, R., 2009. ORIGEN-S: SCALE System Module to Calculate Fuel Depletion, Actinide Transmutation, Fission Product Buildup and Decay, and Associated Radiation Source Terms. Oak Ridge National Laboratory. Henry, A.F., 1975. Nuclear-Reactor Analysis. The MIT Press. JAEA, 2011. External Evaluation on Monju Core Confirmation Tests in FY2010 (the Technical Committee on Monju Research Utilization). Tech. Rep. JAEAEvaluation 2011-002. (in Japanese).

Jo, T., Goto, T., Yabuki, K., Ikegami, K., Miyagawa, T., Mori, T., Kubo, A., Kitano, A., Nakagawa, H., Kawamura, Y., Ida, T., Muranaka, M., Morohashi, Y., Kato, Y., 2011. Core Confirmation Test in System Startup Test of the Fast Breeder Reactor Monju. Tech. Rep. JAEA-Technology 2010-052. (in Japanese). NISA, 2007. Application for the Change of a Nuclear Reactor Facility – JAEA Fast Reactor Development Center – About the Nuclear Design. Tech. Rep. 111A-1-2, Nuclear and Industrial Safety Agency. (in Japanese). Palmiotti, G., Lewis, E.E., Carrico, C.B., 1995. VARIANT: VARIational Anisotropic Nodal Transport. Tech. Rep. Argonne National Laboratory. Petit, O., Hugot, F.-X., Lee, Y.-K., Jouanne, C., Mazzolo, A., 2008. TRIPOLI-4 Version 4, Manuel de l’Utilisateur. CEA. Rimpault, G., Plisson, D., Tommasi, J., Jacqmin, R., Rieunier, J.-M., Verrier, D., Biron, D., 2002. The ERANOS code and data system for fast reactor neutronic analyses. In: PHYSOR 2002. ANS, Seoul, South Korea. Rimpault, G., Smith, P., Jacqmin, R., Malvagi, F., Rieunier, J.-M., Honde, D., Buzzi, G., Finck, P., 1997. Note Technique – Schéma de calcul de réfÉrence du Formulaire ERANOS et Orientations Pour le Schéma de Calcul de Projet. CEA/NT/SPRC/LEPH. Takashita, H., Higuchi, M., Togashi, N., Hayashi, T., 2000. Report on Neutronic Design Calculational Methods (Technical Report). Tech. Rep. JNC TN8410 2000-011, Japan Nuclear Cycle Development Institute (JNC). (in Japanese). Tamagno, P., Van Rooijen, W.F.G., 2012. Uncertainty analysis of the prototype FBR Monju with the JENDL-4.0 nuclear data set. Annals of Nuclear Energy 51, 257– 273 (in Japanese). Taro, K., Takeaki, O., Nobuhiro, O., Jin, Y., Koichiro, K., 1994. Optimization of the fuel design of the prototype FBR Monju based on the startup test results. Tech. Rep. PNC-TJ1678-95-006, Advanced Reactor Technology Development Corporation (PNC). (in Japanese). Tsuji, M., Suzuki, N., Shimazu, Y., 2003. Subcriticality measurement by neutron source multiplication method with a fundamental mode extraction. Journal of Nuclear Science and Technology 40 (3), 158–169. Williams, M.L., 1986. Chapter ‘‘perturbation theory for nuclear reactor analysis. In: Ronen, Y. (Ed.), Handbook of Nuclear Reactor Calculations. CRC Press, pp. 63– 188.