Reactivity calculation in the frequency domain for BWR stability codes

Available online at www.sciencedirect.com annals of

NUCLEAR ENERGY Annals of Nuclear Energy 35 (2008) 1185–1198 www.elsevier.com/locate/anucene

Reactivity calculation in the frequency domain for BWR stability codes J.L. Munoz-Cobo a,*, C. Garcı´a a, A. Escriva´ a, J. Melara b a

Department of Chemical and Nuclear Engineering and Institute of Energy Engineering, Universidad Polite´cnica, Camino de Vera 14, Valencia 46022, Spain b IBERDROLA Ingenierı´a y Construccio´n, C/ Jose´ Bardasano Baos 9, 28016 Madrid, Spain Received 24 July 2007; received in revised form 27 December 2007; accepted 27 December 2007 Available online 4 March 2008

Abstract This paper deals with the problem of computing the feedback reactivity in the frequency domain codes as the LAPUR code. First, we explain how to calculate the feedback reactivity in the frequency domain using slab-geometry (1D) kinetics, also we show how to perform the coupling of the 1D kinetics with the thermal–hydraulic part of the LAPUR code in order to obtain the density to reactivity feedback coefficients, the power to reactivity feedback coefficients and the inlet temperature to reactivity feedback coefficients for the different channels. Also we explain how to solve the lambda eigenvalue equation in the frequency domain by the spectral nodal method and the comparison between the reactivity obtained by solving the eigenvalue equation in the frequency domain and the reactivity obtained from first order perturbation theory. Because LAPUR works in the linear regime, it is assumed that in general the perturbations are small. There is also a section devoted to the reactivity weighting factors used to couple the reactivity contribution from the different channels to the reactivity of the entire reactor core in point kinetics and 1D kinetics. Also we display a comparison of the reactivity obtained using point kinetics and 1D kinetics for the different kind of perturbations and we check the consistency of the results. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction The calculation of stability in nuclear reactors is an important issue from the very beginning of nuclear energy. There are two main ways to compute the reactor stability: in the time domain and in the frequency domain. The first way is used by the MATSTAB (Hangii, 2001) and RAMONA codes, while the second way is used by LAPUR 5 (March-Leuba and Otaduy, 1983; Escriva´ and MarchLeuba, 2001), NUFREQ (Peng et al., 1985), and STAIF. MATSTAB linearizes the equations in the time domain and computes the eigenvalues and eigenvectors of the Jacobian matrix, and from the dominant eigenvalue of this matrix the decay ratio is computed. Frequency domain stability codes proceed in two steps: first they solve the steady state thermal–hydraulic equa-

*

Corresponding author. E-mail address: [email protected] (J.L. Munoz-Cobo).

0306-4549/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2007.12.016

tions, to obtain the steady state values of the different magnitudes that define the reactor state at a given operating point of the power-flow map. Then, in a second step these codes solve in the frequency domain the linearized equations of mass, energy, and momentum, of the liquid and gas phases, the heat conduction equations, and the neutron kinetics equations, to obtain the reactor response to small perturbations performed around the steady state. Usually these codes use point kinetics to compute the reactivity to power transfer function. These codes can perform predictions of the reactor stability parameters like the global decay ratio and the oscillation frequency. Some of them as LAPUR incorporate also the reactivity to power transfer function of the first subcritical mode and are able to predict the decay ratio of the reactor for out-of-phase oscillations. Otaduy (1979) developed the LAPUR code to obtain the stability parameters as decay ratio of a BWR reactor. March-Leuba (1984) added the transfer function of outof-phase modes and coupled this transfer function to the

1186

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

Nomenclature Latin symbols B2r geometrical radial buckling concentration of delayed neutron precursors for Ck the kth group Dg diffusion coefficient for the gth group eix flow to reactivity transfer function for fuel channel ix fix power to reactivity transfer function for fuel channel ix hix inlet subcooling temperature to reactivity transfer function for fuel channel ix H specific liquid enthalpy in the subcooled region J neutron current k multiplication constant ~kðxÞ multiplication constant in the frequency domain L removal operator M production operator Nax number of axial nodes Nch(ix) number of fuel channels of type ix NCH total number of core fuel channels Sf fission source d Te in inlet temperature perturbation V volume of the reactor core Wi,ix reactivity weighting factor for the axial node i, and the channel ix ðaxialÞ Wi axial weighting factor ðradialÞ W ix radial weighting factor

thermal–hydraulics. The new version of LAPUR was able to get good estimations of the DR for in-phase and outphase instabilities. This paper deals with the problem of implementing 1D kinetics in the frequency domain stability code LAPUR and how to couple the 1D kinetics with the thermal– hydraulic part of this code in order to get the different transfer functions that make possible the calculation of the stability parameters of the reactor. The paper has been organized as follows: Section 2.1 deals with the calculation of the reactivity in the time domain and in the frequency domain. Section 2.2 deals with the reactivity weighting factors in point kinetics and 1D kinetics when the reactor core is partitioned in different regions, with different thermal– hydraulic characteristics as different axial power profile, power per channel, frictional resistance in the spacers and so on. Section 2.3 explains how is performed the calculation of the following transfer function: (i) power to reactivity, (ii) flow to reactivity and (iii) inlet subcooling temperature to reactivity, in point kinetics and 1D kinetics.

d~y in z

perturbation in the inlet normalized flow axial coordinate

Greek symbols a void fraction b fraction of delayed neutron precursor per fission e perturbation parameter / neutron flux /+ adjoint neutron flux k lambda eigenvalue of the neutron diffusion equation kk disintegration constant of the delayed neutron precursors for group k qfb feedback reactivity ~fb ðxÞ feedback reactivity in the frequency domain q q0s steam density q0f saturated liquid density R macroscopic cross section R12 fast to thermal scattering cross section Subscripts ax axial dens density fb feedback g neutron group in inlet 0 steady state or unperturbed magnitudes 1D one dimensional kinetics

Section 3 explains the numerical solution of the neutron diffusion equations with two-energy groups for complex cross sections. Section 4 gives the reactivity calculation in the frequency domain, for different channels, different frequencies and different perturbations, when using point kinetics and 1D kinetics. To do this comparison the 1D kinetics reactivity subroutines were implemented in the LAPUR code, and the change in reactivity was computed by both methods. Finally Section 5 discusses the main conclusions of the paper and the future work in this area.

2. Reactivity calculation in the time domain and in the frequency domain 2.1. Basic concepts The equation that models the neutron flux evolution in slab geometry with K groups of delayed neutron precursors and two-energy groups is given in the diffusion approximation by

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

1 o /ðz; tÞ ¼ Lðz; tÞ/ðz; tÞ þ ð1  bÞMðz; tÞ/ðz; tÞ v ot K X þ kk C k ðz; tÞ

ð2:1Þ

1

oC k ¼ bk Mðz; tÞ/ðz; tÞ  kk C k ðz; tÞ ð2:2Þ ot where L(z, t) and M(z, t) are the standard 1D neutron removal and production operators given in standard matrix notation by: " Lðz; tÞ ¼

This equation describes a critical reactor that possesses identical nuclear properties to those of the actual reactor at the instant t, and is made critical by a fictitious value of the effective multiplication constant k(t). Clearly as pointed out by Akcasu, both k(t), and /k(t), will depend on time parametrically because the eigenvalue equation will be different at different times. The reactivity change between time 0 and time t, will be given by #

 ozo D1 ðz; tÞ ozo þ Ra1 ðz; tÞ þ DB2r;1 þ R12 ðz; tÞ

0

R12

 ozo D2 ðz; tÞ ozo þ Ra2 ðz; tÞ þ DB2r;2

 Mðz; tÞ ¼

mRf 1 ðz; tÞ mRf 2 ðz; tÞ 0



 ð2:4Þ

0

Finally /(z, t) = (/1(z, t), /2(z, t))T is the neutron flux column vector. If we assume that the production and removal operators at time t can be expressed as the sum of a time independent part plus a time dependent one, then: Lðz; tÞ ¼ L0 ðzÞ þ e dLðz; tÞ

ð2:5Þ

Mðz; tÞ ¼ M 0 ðzÞ þ e dMðz; tÞ

ð2:6Þ

Where we assume that a perturbation has been performed in the reactor system at time t = 0, and the subscript zero denotes the production and removal operators before the perturbation; e is an order parameter that we set equal to one. At time t = 0 before the perturbation takes place the eigenvalue diffusion problem is easily obtained from Eqs. (2.1) and (2.2) by multiplying the production operator by the fundamental lambda eigenvalue k0 = 1/k0 and setting the time derivatives to zero. Thus we obtain:  L0 ðzÞ/0 ðzÞ þ ð1  bÞk0 M 0 /0 ðzÞ þ

K X

kk C k;0 ðzÞ ¼ 0

ð2:7Þ

k¼1

bk k0 M 0 ðzÞ/0 ðzÞ ¼ kk C k;0 ðzÞ

ð2:8Þ

After some trivial simplifications Eqs. (2.7) and (2.8) reduce to: L0 ðzÞ/0 ðzÞ ¼ k0 M 0 ðzÞ/0 ðzÞ

ð2:9Þ

If we introduce perturbations at time t = 0 in such a way that the production and removal operators change with time as displayed in Eqs. (2.5) and (2.6). For instance this will be the case during an in-phase oscillation where the reactor power and hence the void fractions are changing at the same time (in-phase) in the full reactor core, producing as a result, a variation in the cross sections and in the production and removal operators. Then, as pointed by Akcasu et al. (1971), we can write at time t the following equation: 1 Lðz; tÞ/kðtÞ ðzÞ ¼ Mðz; tÞ/kðtÞ ðzÞ kðtÞ

ð2:10Þ

1187

dqðtÞ ¼ dk ¼ 

1 1  kðtÞ k 0

ð2:3Þ

 ¼

kðtÞ  k 0 kðtÞk 0

ð2:11Þ

If the change in reactivity is small, then Eq. (2.11) can be approximated (Mun˜oz-Cobo et al., 2006) by the following expression that is additive and commutative:   kðtÞ dqðtÞ ¼ log ð2:12Þ k0 So solving the eigenvalue Eq. (2.10), with the nuclear properties at the instant t, allows computing the reactivity change, by means of Eqs. (2.11) or (2.12). Also it is frequently used the following exact expression to compute the reactivity (Lewins, 1978):     dMðtÞ /þ ;  dLðtÞ /0 0 k0 dq ¼ dk ¼ ð2:13Þ 0 0 ð/þ 0 ;M / Þ where the prime means that the neutron flux or the operators are in the perturbed state. For small perturbations, it is possible to compute the reactivity change by using first order perturbation theory, although this method can be easily extended to higher order in the perturbation parameters, (Stacey, 2001). In this case the flux solution and the eigenvalue are written as a series expansion in the perturbation parameter as follow: / ¼ /0 þ e d/ð1Þ þ e2 d/ð2Þ þ oðe2 Þ

ð2:14Þ

k ¼ k0 þ e dkð1Þ þ e2 dkð2Þ þ oðe2 Þ

ð2:15Þ

Substituting Eqs. (2.5), (2.6), (2.14) and (2.15), into Eq. (2.10), retaining only up to first order terms in the perturbation parameter and using Eq. (2.9) to cancel the zero order terms, we obtain the following result: dL/0 þ L0 d/kðtÞ ¼ k0 dM/0 þ k0 M 0 d/kðtÞ þ dkð1Þ M 0 /0

ð2:16Þ

Furthermore we take the inner product by the adjoint flux /þ 0 of Eq. (2.16), this calculation yields, after rearranging some terms the following result:

1188

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

þ ð/þ 0 ; ðL0  k0 M 0 Þd/kðtÞ Þ  ð/0 ; ðk0 dM  dLÞ/0 Þ

¼ dkð1Þ ð/þ 0 ; M 0 /0 Þ " e Lðz; xÞ ¼

e a1 ðz; xÞ þ R e 12 ðz; xÞ e 1 ðz; xÞ o þ R  ozo D oz e 12 ðz; xÞ R

ð2:17Þ

0

0

where the superscript T means transpose. The adjoint flux /þ 0 ðzÞ is the solution of the adjoint eigenvalue equation: Lþ 0

¼

k0 M þ 0 ðzÞ/0 ðzÞ

# ð2:21Þ

e a2 ðz; xÞ e 2 ðz; tÞ o þ R  ozo D oz

where the inner product of the vectors a and b has been defined in the usual way as: Z H ðaðz; tÞ; bðz; tÞÞ ¼ dzaT ðz; tÞbðz; tÞ ð2:18Þ

Lþ 0 ðzÞ/0 ðzÞ

e x ðz; xÞ, cross sections. These cross sections, denoted as R will depend in general on the position and frequency. The removal operator e Lðz; xÞ appears as:

ð2:19Þ

Mþ 0

where and are the adjoints of the removal and production operators, respectively. The first term of the left hand side of Eq. (2.17) cancels by virtue of the commutation relation of the adjoint operators, and on account of Eq. (2.19). Therefore, from Eq. (2.17), one can obtain the reactivity change in first order perturbation theory that is given by:     /þ ; dMðtÞ  dLðtÞ /0 0 k  0þ dqð1Þ ¼ dkð1Þ ¼ ð2:20Þ /0 ; M 0 /0 The integrations that appear in Eq. (2.20) are performed only over the spatial variable. In the frequency domain stability codes, like LAPUR, one performs three different types of perturbations in the following magnitudes: (i) mass flow rate entering the channels, (ii) inlet temperature subcooling and (iii) power. The first one, known as the flow perturbation, is a unitary flow perturbation performed in the mass flow rate entering the reactor channel; this particular perturbation permits to compute the flow to reactivity transfer function. These frequency domain codes solve the linearized mass, energy and momentum conservation equations for each perturbation and each frequency, in the reactor channel nodes. The solution of these equations is a set of perturbations in the fluid temperature and the void fraction for each axial node, of each particular channel. Because we are working in the frequency domain, these quantities are complex numbers, therefore the 1D cross sections that depend in general on the fluid temperature and the void fraction will be also complex numbers. The same thing happens with the reactivity that depends on the production and removal operators that can be expressed in terms of cross sections. Therefore, in the frequency domain codes the cross sections and the reactivities will be, in general, complex numbers. If we have the BWR core cross section data set parameterized in terms of the void fraction, the moderator temperature and the fuel temperature, then, each particular perturbation will change the void fractions and hence the

where the absorption cross sections include the leakage terms. The frequency dependent production operator will be given by: " # e f 1 ðz; xÞ m R e f 2 ðz; xÞ m R e ðz; xÞ ¼ M ð2:22Þ 0 0 Next we look for the complex eigenvalue ~kðxÞ ¼ 1=~kðxÞ of the following equation: e Lðz; xÞ/~k ðz; xÞ ¼

1 e M ðz; xÞ/~k ðz; xÞ ~kðxÞ

ð2:23Þ

Because the production and removal operators depend on the frequency, then the eigenvector /~k ðz; xÞ and the eigenvalue ~kðxÞ will depend on the frequency parametrically. The method to solve this equation in the complex domain is described in Section 3. The reactivity change in the frequency domain for small perturbations will be given by:   ~kðxÞ  k 0 d~ qðxÞ ¼  ~kðxÞ  k0 ¼ ð2:24Þ ~kðxÞk 0 For small perturbations, one can use the following definition of reactivity, which is additive and commutative: d~ qðxÞ ¼ logð~kðxÞ=k 0 Þ

ð2:25Þ

However, it is interesting to show that the eigenvalue of Eq. (2.23) is closely related to the reactivity change in the frequency domain by means of Eq. (2.24). To do this we perform the following perturbation expansions in the frequency domain of the production and removal operators, the eigenfunction and the eigenvalue: ~kðxÞ ¼ k0 þ e d~kð1Þ ðxÞ þ Oðe2 Þ

ð2:26Þ

e Lðz; xÞ ¼ L0 ðzÞ þ e d e Lðz; xÞ þ Oðe Þ

ð2:27Þ

e ðz; xÞ ¼ L0 ðzÞ þ e d M e ðz; xÞ þ Oðe2 Þ M

ð2:28Þ

/~k ðz; xÞ ¼ /0 ðzÞ þ d/...k ðz; xÞ þ Oðe2 Þ

ð2:29Þ

2

then we substitute Eqs. (2.26)–(2.29) in the frequency domain eigenvalue Eq. (2.23), and we neglect second and higher order terms. After some elementary algebra on account of Eq. (2.9) we obtain the following result: de Lðz; xÞ/0 ðzÞ þ L0 ðzÞd/~k ðz; xÞ e /0 ðzÞ þ k0 M 0 ðzÞd/~k ðz; xÞ ¼ dk~ð1Þ ðxÞM 0 ðzÞ/0 ðzÞ þ k0 d M ð2:30Þ

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

Performing the inner product of Eq. (2.30) by the unperturbed adjoint flux /þ 0 ðzÞ yields on account of Eq. (2.19): dkð1Þ ðxÞ ¼ ð~ kðxÞ  k0 Þ þ Oðe2 Þ     e ðxÞ  d e /þ ; k d M LðxÞ /0 0 0  þ ¼ ð2:31Þ /0 ; M 0 /0 If we apply the Laplace transform to Eq. (2.20), we obtain the following expression for the reactivity change d~ qð1Þ ðxÞ in the frequency domain and in first order perturbation theory:  ~    dMðxÞ e /þ  d LðxÞ /0 0; k  0þ d~ qð1Þ ðxÞ ¼ d~ kð1Þ ðxÞ ¼ ð2:32Þ /0 ; M 0 /0 where the variation of the production and removal operators in the frequency domain are given by: e ðxÞ ¼ M e ðz; xÞ  M 0 ðzÞ dM de LðxÞ ¼ e Lðz; xÞ  L0 ðzÞ

ð2:33Þ ð2:34Þ

Comparison of Eqs. (2.31) and (2.32), show that for small reactivity perturbations, the reactivity change in the frequency domain computed using first order perturbation theory is, in first order approximation, equal to the reactivity computed by solving the eigenvalue equation in the frequency domain. 2.2. Reactivity weighting factors in 1D kinetics In this section we describe how reactivity contributions of the different channel regions are weighted to obtain the total reactivity of the entire core in the LAPUR code, also we explain how we have generalized this procedure to obtain the total reactivity of the core using 1D kinetics. In the LAPUR code, the entire reactor core is divided in different regions denoted by the index ix, each one of these regions has several reactor fuel channels that are characterized by the same power axial profile and the same type of fuel elements. Physically this corresponds to fuel elements of the same type located within the same radial region. We denote by Nch(ix), the number of reactor channels belonging to the same reactor fuel region ix and by NCH, the total number of channels of the entire core. A fuel channel of type ix is characterized as representative of a group of Nch(ix) channels with the same axial power profile and the same type of fuel element. We begin with the well known exact reactivity expression in the time domain in 3D:  R þT dM /  dL /0 dV 0 V k0 R þT 0 0 dq ¼ ð2:35Þ /0 M / dV V where V denotes the core volume, and the primes denote the perturbed values. Now we assume that the neutron flux can be factorized in a radial part and an axial part for each jth channel

1189

/  fp;j ðr; hÞ/ðjÞ ðz; tÞ for r; h 2 V j

ð2:36Þ

To compute the weighting factors in 1D kinetics we assume one-energy group diffusion theory where the direct and adjoint fluxes are the same and therefore we can express Eq. (2.35) using consistent cross section for the removal operator, i.e, adding the radial buckling terms to the absorption cross sections (Munoz-Cobo et al., 1993; Munoz-Cobo et al., 2006). Because the void fraction and the coolant temperatures change from the hot to the peripheral channels, the cross sections and the neutron fluxes vary from one channel to the others. If the index j runs over the different fuel S channels and we express the total core volume as: V = jVj, then Eq. (2.35) can be written as follows (Munoz-Cobo et al., 1993, 1996): R ðjÞ   PR H þ 0 2 dM dAf ðr;hÞ dz/ ðzÞ  dL ðz; tÞ / consistent p;j 0 j Aj 0 k0 dq  R ðjÞ PR H 2 0 0 dz/þ 0 ðzÞM / ðz; tÞ j Aj dAf p;j ðr; hÞ 0 ð2:37Þ

Eq. (2.37) can be expressed in the form: R ðjÞ   H 0 dM N CH dz/þ X 0 ðzÞ k 0  dLconsistent / ðz; tÞ 0 dq  F 2p;j R ðjÞ H þT 0 0 j¼1 dz/ M / 0 0 R ðjÞ H 0 0 dz/þ 0M / 0  ð2:38Þ ðkÞ PN CH 2 R H 0 0 F dz/ ðzÞM / ðz; tÞ 0 k¼1 p;k 0 where Fp,j is approximately the fraction of power that goes to the jth fuel channel and the first term is the reactivity variation of the jth channel, so we can recast this expression as follows: dq 

N CH X

F 2p;j

j¼1

R H 

0

PN CH

2 k¼1 F p;k

R

0 0 dz/þT 0 M /

ðjÞ

ðjÞ ðkÞ dq dz/0 ðzÞM / ðz; tÞ

H 0

0

ð2:39Þ

0

The axial power profile and the thermal–hydraulic properties are the same for all the channels belonging to one specific region, denoted by ix in the LAPUR code. Then if we denote by Nch(ix), the number of fuel channels belonging to the ixth region, then because the reactivity variation is the same for all the channels belonging to the same region, we can denote this reactivity by dq(i) = dq(ix) for i 2 ix, and Eq. (2.39) can be expressed in the form: X dq  N ch ðixÞF 2p;ix ix

R

ðixÞ 0 0 dz/þT 0 M / ðixÞ  R ðkxÞ dq P H þ 2 0 0 dz/0 ðzÞM / ðz; tÞ kx N ch ðkxÞF p;kx 0 H 0

ð2:40Þ

1190

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

Finally, by transforming the integral in Eq. (2.40), in a sum over the axial nodes we can write the following approximate expression for the reactivity variation when using 1D kinetics to calculate the reactivity of the channels: dq 

X

dqN ch ðixÞ

N ax X

ix

W i;ix

ð2:41Þ

i¼1

where Wi,ix are the weighting factors that according to Eq. (2.40), and neglecting perturbation terms are given by: ðixÞ

W i;ix ¼ P

F 2p;ix ðP i Þ2 PN ax ðkxÞ 2 2 kx N ch ðkxÞF p;kx j¼1 ðP j Þ

ð2:42Þ

ðixÞ

where P i is the axial power distribution in the axial nodes denoted by the subscript i of the core region ix. Finally, expression (2.42) can be expressed in terms of radial and axial weighting factors. To do this, let us define the weighting factor for a given region ix, as: ðradÞ

W ix

¼P

F 2p;ix 2 ix N ch ðixÞF p;ix

ð2:43Þ

Also we define the axial weighting factor for the axial node i, belonging to the ix region, in the form: ðixÞ

ðaxialÞ

Wi

ðP Þ2 ðixÞ ¼ PN ax i ðixÞ 2 j¼1 ðP j Þ

ð2:44Þ

In terms of the radial and axial weighting factors, the reactivity weighting factors Wi,ix given by expression (2.42), can be expressed in the following form: PN ax ðixÞ 2 ðradÞ ðaxialÞ W ix W i ðixÞ k¼1 ðP k Þ W i;ix ¼ P ð2:45Þ ðradÞ PN ax ðkxÞ 2 kx N ch ðkxÞW kx j¼1 ðP j Þ At this point we remark that the weighting factors, defined by Eq. (2.42), satisfy the following normalization condition: N ax XX ix

N ch ðixÞW i;ix ¼ 1

ð2:46Þ

i¼1

2.3. Feedback reactivity with point kinetics and 1D kinetics To understand the way in which the 1D kinetics has been implemented into the LAPUR code, we will first describe the way in which the point kinetics is implemented in LAPUR 5.2, then we will explain the development of the 1D kinetics and its coupling with the LAPUR code. In LAPUR 5.2 the calculation of the flow to density reactivity feedback coefficient eix at a specified channel type (ix) is computed in the FREQ subroutine performing a unit perturbation of normalized flow at the channel inlet, while the power and the inlet subcooling temperature perturbations are set to zero: eix ¼ bd~ qixfb ðxÞc d~y in ¼1 d~ q¼0

de T in ¼0

ð2:47Þ

qixfb ðxÞ that where eix gives the feedback reactivity change d~ is produced in one channel of type ix, when a unit perturbation of normalized flow d~y in (ix) is performed in the flow entering this channel. The LAPUR 5 code solves the linearized and discretized mass, energy and momentum equations in the frequency domain and computes the void fraction perturbations d~aix ðiÞ and also the liquid enthalpy perturbations in the nodes of the channel ix. It also solves the heat conduction equations in the frequency domain to compute the perturbations in the fuel temperature. Then the transfer function eix is computed in the FREQ subroutine by means of the following expression:  N ax  X oqfb eix ¼ W i;ix ½d~ q0dens;i;ix  d~y in ¼1 ð2:48Þ 0 oq d~ q¼0 dens i;ix i¼1 de T in ¼0   oqfb is the density reactivity coefficient for the where oq0 dens

i;ix

ith node, of the ixth channel. The value of this coefficient for the ith node in LAPUR 5.2 is obtained from a set of tabulated density reactivity coefficients in terms of the average relative water density ðq0dens =q0f Þa;i;ix , the way to obtain these coefficients have been explained by Munoz-Cobo et al. (2006). Wi, ix is the reactivity weighting factor that depends on the square of the power distribution, and is given by Eqs. (2.42) or (2.45). Finally ½d~ q0dens;i;ix  d~y in ¼1 is the density variation in the ith d~ q¼0

de T in ¼0

node of the ix channel produced by the unit flow perturbation d~y ix . For the nodes belonging to the subcooled-boiling region, the average density variation in the ith node is given by:  0 oql 0 e a;i þ ðq0  q0 Þ d~aa;i d~ qdens;ai ¼ ð1  d~aa;i Þd H ð2:49Þ s f ai oH i where q0s and q0f are the densities of the steam and saturated liquid, respectively, the subscript a refers to the average vae a;i is the liquid lue of a physical magnitude in the ith node; d H enthalpy variation produced in the ith node by the unit flow e a;i ¼ 0, in perturbation. In the bulk boiling region we set d H Eq. (2.49), while in the nodes of the non-boiling region we have only variations in the liquid enthalpy. We have dropped the subscript ix in expression (2.49). The calculation of the power to density reactivity feedback coefficient fix at a specified channel type (ix), is computed also in the TRANS subroutine performing a unit power perturbation d~q: fix ¼ ½d~ qixfb  d~y in ¼0

ð2:50Þ

d~ q¼1

de T in ¼0

qixfb that is where fix gives the feedback reactivity change d~ produced when a unit perturbation of normalized power d~q is performed in the fuel power of the ix-channel. The LAPUR 5 code solves the linearized mass, energy and momentum equations of that channel in the frequency domain and computes the void fraction perturbations d~ aix ðiÞ, e a;i in the nodes of and the liquid enthalpy perturbations d H

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

the channel ix, produced by the power perturbation. It also solves the heat conduction equations in the frequency domain to compute the perturbations in the fuel temperatures. Then, the transfer function is computed in the FREQ subroutine by means of the following expression:  N ax  X oqfb fix ¼ W i;ix ½d~ q0dens;i;ix  d~y in ¼0 ð2:51Þ oq~0 dens i;ix d~ q¼1 i¼1 T in ¼0 de The calculation of the inlet temperature to density reactivity feedback coefficient hix at a specified channel type ix is also performed in the FREQ subroutine. To obtain hix, we perform a unitary perturbation in the inlet temperature of the subcooled liquid entering to one channel of the ixtype. Then the program solves the linearized mass, energy and momentum conservation equations in the frequency domain and obtains for each channel node the perturbations in the flow, void fraction and liquid enthalpy, i.e. e a;i , produced by this particular perturbation. d~y a;i ; d~ aa;i ; d H Finally, this feedback coefficient is also computed in the FREQ subroutine by means of the following expression:  N ax  X oqfb hix ¼ W i;ix ½dq~0 dens;i;ix  d~y in ¼0 ð2:52Þ o~ q0dens i;ix d~ q¼0 i¼1 T in ¼1 de The density reactivity change d~ qdens;fb in the entire reactor core produced by the set of perturbations d~y in ðixÞ, d Te in and d~ q, performed in the reactor channels, can be recasted in terms of the density reactivity feedback coefficients as follows: X X d~ q dens;fb ¼ eix N ch ðixÞd~y in ðixÞ þ fix N ch ðixÞd~ q ix

þ

X

ix

hix N ch ðixÞd Te in

ð2:53Þ

ix

where Nch(ix) is the number of channels of ix-type, and we have assumed perfect mixing of the fluid in the lower plenum i.e. d Te in ðixÞ ¼ d Te in . Now we describe the way to implement the 1D kinetics, that is as follows, for each perturbation, for instance in the flow d~y ix entering to the channels of region ix, the code solves the mass, energy and momentum equations, and computes the perturbations in the flow, void fraction and liquid enthalpy. The changes in the void fraction produce changes in the cross sections. These new cross sections are complex numbers because we are working in the frequency domain, then we solve the 1D diffusion eigenvalue problem given by Eq. (2.23), by the spectral nodal method (Barros et al., 2003, 2004), with the new perturbed cross section values, and we compute the complex eigenvalue ðixÞ ~kðxÞ, and the reactivity variation d~ q1D ðxÞ ¼ ð1=~kðxÞ ðixÞ 1=k 0 Þ . Also to check the consistency of the results we compute this same reactivity using first order perturbation theory and Eq. (2.31). Then the flow to density reactivity feedback coefficient eix,for a specified channel type ix, is computed by means of the expression:

ðixÞ

q1D ðxÞ d~y ix ¼1 eix ¼ ½d~ de T in ¼0

N ax X

W i;ix

1191

ð2:54Þ

i¼1

d~ q¼0

The power to density reactivity feedback coefficient fix at a specified channel type (ix) is computed also in the TRANS subroutine performing a unit power perturbation dq in channel ix, solving the mass, energy and momentum conservation equations in the channel, computing the variation in the cross section and solving the 1D diffusion eigenvalue problem given by Eq. (2.23), with the new cross section set, or computing the reactivity change by first order perturbation theory using Eq. (2.31). ðixÞ

fix ¼ ½d~ q1D ðxÞ d~y ix ¼0 de T in ¼0

N ax X

W i;ix

ð2:55Þ

i¼1

d~ q¼1

A similar method is used to compute the inlet temperature to density reactivity feedback coefficient hix. ðixÞ

hix ¼ ½d~ q1D ðxÞ d~y ix ¼0 de T in ¼1

N ax X

W i;ix

ð2:56Þ

i¼1

d~ q¼0

3. Numerical solution of two-energy group slab-geometry neutron diffusion equations by the spectral nodal method for complex cross sections In this section we describe a spectral nodal method for two-energy group slab-geometry diffusion eigenvalue problems. The essence of the spectral nodal method for multigroup slab-geometry diffusion eigenvalue problems is carefully described in the papers by Barros et al. (2003, 2004). We perform a spectral analysis of the two-group diffusion equations to obtain the local general solution in each spatial node, for a given estimate of the dominant eigenvalue in the power iterations. The numerical solution generated by the present spectral nodal diffusion (MEND) method for multigroup eigenvalue problems in slab geometry is completely free from spatial truncation errors. Let us consider the two-energy group slab-geometry neutron diffusion equation: dJ 1 ðzÞ 1 þ RT 1 ðzÞ/1 ðzÞ ¼ ½mRf 1 ðzÞ/1 ðzÞ þ mRf 2 ðzÞ/2 ðzÞ dz k d J 1 ðzÞ ¼ D1 ðzÞ /1 ðzÞ ð3:1Þ dx dJ 2 ðzÞ þ RT 2 ðzÞ/2 ðzÞ ¼ Rs12 ðzÞ/1 ðzÞ dz d ð3:2Þ J 2 ðzÞ ¼ D2 ðzÞ /2 ðzÞ; 0 6 z 6 H : dz Here the subscripts 1 and 2 represent the fast and the thermal energy groups, respectively, Dg is the group diffusion coefficient, k is the dominant eigenvalue, mRfg is the group g fission macroscopic cross section times the fission neutron yield, RT1 is the total macroscopic cross section for the fast energy group and RT2 is the total macroscopic cross section

1192

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

Fig. 1. One-dimensional spatial grid for the application of the spectral-nodal method.

for the thermal energy group, Rs12 is the group transfer cross section. We remark that in Eqs. (3.1) and (3.2) we have neglected the upscattering events and we have considered all the parameters as complex numbers. Considering Fig. 1, that illustrates an arbitrary spatial grid. To perform a spectral analysis of Eqs. (3.1) and (3.2) and determine the local general solution in node Xi, we seek elementary solutions of the forms: /g;n ðzÞ ¼ ag ðnÞ expðnzÞ;

ð3:3Þ

J g;n ðzÞ ¼ bg ðnÞ expðnzÞ;

z 2 Xi ;

g ¼ 1; 2

ð3:4Þ

whose first terms ag(n) and bg(n) carry the energy dependence and the exponential terms carry the spatial dependence. Substituting Eqs. (3.3) and (3.4) into Eqs. (3.1) and (3.2) we obtain a homogeneous system of four algebraic linear equations in four unknowns ag(n) and bg(n), g = 1:2. In order to obtain non-trivial solutions of the previous homogeneous system, the determinant of the coefficient matrix must be equal to zero. As a result, we obtain a polynomial equation of degree four, whose roots are the four values of n. It is simple to verify that this polynomial equation has only even powers of n. Therefore, the four values of n appear in ±pairs. Moreover, for practical applications, the material parameters generally keep algebraic relationships that yield simple roots to the polynomial equation.   mRf 1 D2 4 D1 D2 n þ  RT 1 D2  RT 2 D1 n2 þ RT 1 RT 2 k 1  ðmRf 1 RT 2 þ mRf 2 Rs1!2 Þ ¼ 0 ð3:5Þ k Once we have determined the four roots ±n1 and ±n2 of Eq. (3.5), we choose a1 ðn‘ Þ ¼ n‘ ;

‘ ¼ 1 : 2;

ð3:6Þ

therefore we obtain a2 ðnl Þ ¼

nl RS1;2;i ; ðRT 2;i  n2l D2;i ; Þ

‘ ¼ 1 : 2:

ð3:7aÞ

Moreover, we obtain bg ðnÞ ¼ nDg;i ag ðnÞ;

g ¼ 1 : 2;

ð3:7bÞ

Therefore, for each estimate of the dominant eigenvalue k in the outer iterations, the expressions for the local general solutions of Eqs. (3.1) and (3.2) in node Xi appear as /g ðzÞ ¼

4 X ‘¼1

a‘ /g;n‘ ðzÞ ¼

4 X ‘¼1

a‘ ag ðn‘ Þ expðn‘ zÞ

ð3:8aÞ

and J g ðzÞ ¼

4 X

b‘ J g;n‘ ðzÞ ¼

‘¼1

x 2 Xi ;

4 X

b‘ bg ðn‘ Þ expðn‘ zÞ;

‘¼1

g ¼ 1 : 2;

ð3:8bÞ

i ¼ 1 : I;

where a‘ and b‘ are arbitrary constants. At this point we describe a convergent numerical method that generates dominant numerical solution for Eqs. (3.1 and 3.2) that is completely free from spatial truncation errors, because it preserves the analytical general solutions (3.8a and 3.8b) and assure continuity at the node interfaces and satisfy the appropriate boundary conditions. Therefore, we apply the operator Z 1 ziþ1=2 dz ð3:9Þ hi zi1=2 to Eqs. (3.1) and (3.2) to obtain the two-group discretized spatial balance diffusion equations in the form J 1;iþ1=2  J 1;i1=2  1i ¼ 1 ½mRf 1i /  1i þ mRf 2i /  2i  þ RT 1i / k hi D1i J1i ¼  ð/1;iþ1=2  /1;i1=2 Þ hi J 2;iþ1=2  J 2;i1=2  1i ¼ Rs12 /  1i þ RT 2i / hi D2i J2i ¼  ð/2;iþ1=2  /2;i1=2 Þ: hi

ð3:10Þ

ð3:11Þ

Here hi is the width of node Xi, i = 1:I and we have used the following definition for the average of a magnitude over a node: Z 1 ziþ1=2 sg;i ¼ sg ðzÞdz; s ¼ / or J : ð3:12Þ hi zi1=2 In the system represented by Eqs. (3.10) and (3.11), there are four linear algebraic equations and eight unknowns  1i ; J 2;iþ1=2 ; for each spatial node, i.e., J 1;iþ1=2 ; /1;iþ1=2 ; J1i ; /    /2;iþ1=2 ; J 2i ; /2i . Therefore, we need four auxiliary equations that we write as  1i ¼ c11 ð/1;iþ1=2 þ /1;i1=2 Þ þ c21 / i i 2 ð/2;iþ1=2 þ /2;i1=2 Þ  2

ð3:13Þ

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

ð/2;iþ1=2 þ /2;i1=2 Þ þ c12 i 2 ð/1;iþ1=2 þ /1;i1=2 Þ  2 ðJ þ J 1;i1=2 Þ D1i 21 1;iþ1=2 J1i ¼ c11 þ c i 2 D2i i ðJ 2;iþ1=2 þ J 2;i1=2 Þ  2 ðJ þ J 2;i1=2 Þ D2i 12 2;iþ1=2 J2;i ¼ c22 þ c i 2 D1i i ðJ 1;iþ1=2 þ J 1;i1=2 Þ :  2

 2i ¼ c22 / i

ð3:14Þ

ð3:15Þ

ð1Þ

where the superscript L stands for the left end or bottom of the slab. Similarly, for z = zI+1/2 we write !    R  /1;Iþ1=2 J 1;Iþ1=2 0 a11 ¼ ; ð3:18Þ /2;Iþ1=2 J 2;Iþ1=2 aR aR 21 22 where the superscript R indicates that these are the righthand side or upper boundary conditions. Following the last procedure, we obtain two three-diagonal systems of algebraic linear equations for the group node-edge scalar fluxes. For the fast group (g = 1), we obtain 1 A1 ~ /1 ¼ ðF1~ /1 þ F2 ~ /2 Þ: k

ð3:19Þ

For the thermal group (g = 2), we obtain /2 ¼ B1~ /1 : A2 ~

iterative power scheme (Tyrtysshnikov, 1997) in the outer iterations. The solution procedure is initiated by estimating ð0Þ the fission source distribution S f , and an effective multiplication constant, k(0), and solving the group 1 equation ð1Þ (3.19) for the first iterate flux, /1 , here we used the Gaussian elimination with backward substitution to solve each tridiagonal system. Next, the group 2 equation (3.20) is ð1Þ solved for the first iterative flux, /2 , using the first iterative ð1Þ flux /1 . Once we determine all group fluxes we calculate the average fluxes using the Eqs. (3.13 and 3.14). The average fluxes are then used to compute a first iterate fission source:

ð3:16Þ

21 22 In the auxiliary Eqs. (3.13)–(3.16) the constants c11 i ,ci , ci 12 and ci are calculated in such a way that the general analytical solution (3.8a) and (3.8b) is preserved (Barros et al., 2003, 2004). We solve the system (3.10) and (3.11) by substituting the node average scalar fluxes and currents given by the auxiliary Eqs. (3.13)–(3.16) together with appropriate continuity and boundary conditions. Here we remark that in this particular case the resulting constants are complex numbers, because we have used cross sections with complex values. In addition, without upscattering, we write the boundary conditions for z = z1/2 in the form !    L  /1;1=2 J 1;1=2 0 a11 ¼ ; ð3:17Þ /2;1=2 J 2;1=2 aL21 aL22

ð3:20Þ

In Eqs. (3.19) and (3.20)Ag, Fg for g = 1, 2 and B1 are complex symmetric three-diagonal matrices of order I + 1, and ~ /1 , ~ /2 are the group node-edge scalar fluxes. We remark that, at each estimate of k in the outer iterations, we need to calculate the four values of n, and with these, the values of cg0 ;g;i , in order to calculate the entries of the complex tridiagonal matrices A1 and A2, which are not constant in the numerical scheme. To solve the resulting eigenvalue problem (3.19) and (3.20) for complex matrixes we have used the conventional

1193

S f ðrÞ ¼

2 X

mRgf ðrÞ/ð1Þ g ðrÞ:

ð3:21Þ

g¼1

And a first iterate effective multiplication constant is determined by means of the expression: R ð1Þ k ð0Þ drS f ðrÞ k ð1Þ ¼ R : ð3:22Þ ð0Þ drS f ðrÞ The iterative process is continued until the relative deviations of the estimates of the eigenvalues and eigenvectors generated in two successive iterations differ by less than a given convergence number. That is,

ðnÞ

k  k ðn1Þ



ð3:23Þ

k ðn1Þ < ek : and

ðnÞ

/  /ðn1Þ



i;g i;g

< ej :

ðn1Þ

/i;g

ð3:24Þ

Because the matrixes A1 and A2 are complex, we worked in complex arithmetic throughout and we obtained complex values for eigenvalues and eigenvectors. We used the well known library LAPACK (Anderson et al., 1999) to compute the algebraic eigenvalue problem in a complex field. 4. Comparison of reactivity calculations by point kinetics and 1D kinetics in the LAPUR code In this section we compare the reactivity calculations using point kinetics and 1D kinetics in the LAPUR code. We have performed these calculations for several cases of the BWR Ringhals stability benchmark of the OCDE (1996). Because LAPUR code works in the frequency domain, we have performed the calculations of the reactivity variations using point kinetics and 1D kinetics, for a set of frequencies ranging from 0.01 Hz, up to 1 Hz. Also these calculations have been performed for different types of perturbations, and for different regions of the reactor core. The reactivity calculations using 1D kinetics have been computed by two different methods. The first one is solving the eigenvalue equation in the frequency domain by the spectral nodal method explained in Section 3 and then computing the reactivity variation d~ qðixÞ ðxÞ ¼ ð1=~kðxÞ

1194

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

1=k 0 ÞðixÞ . Using the second method we computed the reactivity variation by first order perturbation theory. In this case the forward and adjoint fluxes of the unperturbed case were computed by the analytical nodal method (Henry, 1980). For comparison purposes of 1D kinetics results and point kinetics results we compare the reactivity variation associated to one given channel belonging to a particular region ix of the reactor core. Because the LAPUR code works in the linear regime then all the equations in the frequency domain have been linearized in this code. This can be easily checked by dividing the inlet flow perturbation by a factor Facflow, for instance 10, computing the flow to reactivity transfer function for this particular perturbation and multiplying the result by 10, we obtain the same result as performing a unit flow perturbation. However when calculating the reactivity

change by solving first the eigenvalue problem given by Eq. (2.23), and then computing the reactivity using Eq. (2.24), it happens that if the perturbations are not small then the reactivity variations are not the same as those generated by first order perturbation theory. Normally the temperature perturbations and the power perturbations are small within the range of linear theory and the reactivities computed using the MEND (Spectral Nodal) method and first order perturbation theory match exactly as displayed in Figs. 2 and 3, where we can see that these reactivities are practically the same for all the frequencies but are different from point kinetics. However both 1D and point kinetic reactivities have the same order of magnitude and the shape of both curves is similar. We remark that we have computed the feedback reactivity variation per channel of type ix, due to a unit perturba-

Channel 5 Inlet temperature perturbation-P Real 0.5 Lapur

Mend

Ktrac

0.4 0.3

Reactivity in %

0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 -0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 2. Comparison of the real parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for a unit perturbation in the inlet subcooling temperature to channel 5, case R15b1 of Ringhals Stability Benchmark.

Channel 5 Temp-P Imag 0.5 Lapur

Mend

Ktrac

0.4 0.3

Reactivity %

0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 -0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 3. Comparison of the imaginary parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using the MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for a unit perturbation in the inlet subcooling temperature to channel 5, case R15b1 of Ringhals Stability Benchmark.

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

tion in the inlet temperature i.e. hix ¼ ½d~ qixfb ðxÞ d~y in ¼0 , using d~ q¼0

T in ¼1 de Eq. (2.56) for 1D kinetics and Eq. (2.52) for point kinetics. For the flow variations because the mass flow rate entering one channel is normalized with the steady state mass flow rate entering that channel, then one normalized mass flow perturbation unit produces big changes in the thermal–hydraulic variables. This does not produce any effect in point kinetics theory, but when we use 1D kinetics and we solve the 1D eigenvalue equation we must check that the perturbations are within linear theory. The way to solve this problem is to divide the mass flow rate perturbation entering the channel by a factor, normally 10, seeking to work in the linear region, where the perturbations are small. Then once we have computed the reactivity we mul-

1195

tiply by this factor i.e. we multiply by 10 to obtain the feedback reactivity produced by a perturbation unit. We have computed the reactivity variations at different frequencies ranging from 0.01 Hz to 1 Hz, that are produced by a normalized inlet flow perturbation equal to 1/ 10 in the channel number 5. Figs. 4 and 5 display the real and imaginary parts of the reactivity variations produced by this perturbation and computed by using point kinetics (LAPUR). In addition these same figures display for the same perturbation and the same channel, the real and imaginary parts of the reactivity variations calculated using 1D kinetics (MEND method), and 1D kinetics (first order perturbation theory). In this case the three methods agree very well, and we note that for this particular case the point kinetics and the MEND methods match exactly for frequencies below 0.3 Hz and that the MEND and the first Comparison of 1D kinetics and point kinetics

Channel 5 Flux perturbation-P Real 8.0

Lapur

Mend

Ktrac

7.0 6.0

Reactivity in %

5.0 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 4. Comparison of the real parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using the MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for an 1/10 perturbation of the unit inlet mass flow entering channel 5, case R15b1 of Ringhals Stability Benchmark.

Channel 5 Flow-P Imag 3.0 Lapur

Mend

Ktrac

2.0

Reactivity %

1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 5. Comparison of the imaginary parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using the MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for an 1/10 perturbation of the unit inlet mass flow entering channel 5, case R15b1 of Ringhals Stability Benchmark.

1196

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198 Channel 4 Power-P Real case R15.c2 2.0 Lapur

Mend

Ktrac

1.0

Reactivity 1D %

0.0 -1.0 -2.0 -3.0 -4.0 -5.0 -6.0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 6. Comparison of the real parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using the MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for a unit perturbation in power in channel 4, case R15c2 of Ringhals Stability Benchmark.

Canal 4 Power P Imag Case R15.c2 5.0 Lapur

Mend

Ktrac

4.0

Reactivity 1D %

3.0 2.0 1.0 0.0 -1.0 -2.0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 7. Comparison of the imaginary parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using the MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for a unit perturbation in power in channel 4, case R15c2 of Ringhals Stability Benchmark.

Canal1 Temp-P Real Case R15.C2 0.5 Lapur

Mend

Ktrac

0.4 0.3

Reactivity 1D %

0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 -0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 8. Comparison of the real parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using the MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for a unit perturbation in the inlet subcooling temperature to channel region 1, case R15c2 of Ringhals Stability Benchmark.

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

1197

Canal1 Temp-P Imag 0.5 Lapur

Mend

Ktrac

0.4 0.3

Reactivity 1D %

0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 -0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Frequency [Hz]

Fig. 9. Comparison of the imaginary parts of the reactivity versus the frequency for point kinetics (LAPUR), 1D kinetics using the MEND method, and 1D kinetics using first order perturbation theory (KTRAC), for a unit perturbation in the inlet subcooling temperature to channel region 1, case R15c2 of Ringhals Stability Benchmark.

order perturbation theory methods match exactly above 0.3 Hz. Next we have computed the reactivity change versus the frequency produced by a unit perturbation in power performed at channel region 4 for case R15.c2 of Ringhals OCDE stability benchmark. We observe in Figs. 6 and 7 that the reactivities computed above 0.2 Hz agree for the three methods. However the real parts of the reactivities below 0.2 Hz display a better agreement between point kinetics and the MEND method, however the differences with the first order perturbation theory method are very small except at very low frequency i.e. 0.01 Hz. The imaginary parts of the reactivity for both 1D methods, i.e first order perturbation theory and the MEND method, agree very well at all frequencies for this case R15.c2 and channel 4. However for frequencies below 0.1 Hz some differences do exist for point kinetics and 1D kinetics, as displayed in Fig. 7. It is noteworthy to remark that both methods MEND and first order perturbation theory display in general a very consistent behaviour for the reactivity calculations in all type of channels for instance for the hot Channel (region 1), we have computed the reactivity variation for a perturbation in the inlet subcooling temperature of 1 degree Celsius. The real and imaginary parts of the reactivity variations produced by this perturbation, as calculated by the three methods in the channels of region 1, are displayed in Figs. 8 and 9. 5. Conclusions and final remarks Stability calculations in the frequency domain mainly when we have an axial variation produced by the presence of xenon can give important differences between point kinetics and 1D kinetics. For this reason we have developed the calculation of 1D kinetics in the complex domain. We have performed this calculation by two different meth-

ods. The first one is the spectral nodal method based on previous development by Barros et al. (2003, 2004) and that has been extended in this paper to the complex domain. The second one is first order perturbation theory using the analytical nodal method to compute the steady state neutron fluxes. These 1D methods in the frequency domain have been implemented in the LAPUR code. In order to couple the 1D kinetics in the complex domain to the LAPUR code, we have obtained the expressions of the weighting factors that need to be used in order to get the reactivity contribution of each region to the total reactivity of the reactor core. In general we conclude that for small perturbations within the linear regime, the results of the reactivity variation in the frequency domain solving the 1D eigenvalue equation by the MEND method or computing the reactivity variation by first order perturbation theory are very similar at all frequencies. Also these reactivity changes computed by using 1D kinetics in the reactor channels are in general very close to the reactivity changes computed by using point kinetics. This fact explains why the LAPUR predictions with point kinetics are in general very good, assuming that a good methodology to generate reactivity coefficient tables is used. The convergence of the MEND method in the complex domain to obtain the ~kðxÞ eigenvalue is very fast and the number of iterations that required for convergence is small. Only for large perturbations in the mass flow which are outside of the linear range we need more than 100 iterations. The method for coupling the 1D reactivity calculations by means of the reactivity weighting factor in order to get the different transfer functions eix, fix, hix seems to be well implemented since for the inlet temperature perturbations of several cases as the ones displayed in Figs. 8 and 9, the reactivity variations computed using 1D and point kinetics are practically the same. In other cases small differ-

1198

J.L. Munoz-Cobo et al. / Annals of Nuclear Energy 35 (2008) 1185–1198

ences between point kinetics and 1D kinetics do appear, but this is due to the limitations of point kinetics to predict correctly the reactivity when we have strong variations in the axial properties. Comparison of the decay ratio calculated using 1D complex kinetics in the program LAPUR for several reference cases supplied by the vendor for Cofrentes NPP show a better agreement than using point kinetics. Acknowledgements This research has been performed under the sponsorship of IBERDROLA GENERATION inside the framework of a R&D project called DROP (Decay Ratio On-line Predictor) to develop a predictive stability tool to couple with the Core Monitoring System of a BWR. The authors are particularly grateful to M. Albendea from IBERDROLA GENERATION, for the project coordination and continuous interest and support. References Akcasu, Z., Lellouche, G.S., Shotkin, L.M., 1971. Mathematical Methods in Nuclear Reactor Dynamics. In: Nuclear Science and Technology, vol. 7. Academic Press. Anderson, E., Bai, Z., Bischof, C. et al., 1999. LAPACK Users Guide. SIAM. Barros, R. et al., 2003. The application of spectral nodal methods to discrete ordinates and diffusion problems in Cartesian geometry for

neutron multiplying systems. Progress in Nuclear Energy 42 (4), 385– 426. Barros, R. et al., 2004. Recent advances in spectral nodal methods for numerically solving neutron diffusion Eigenvalue problem. In: Transport Theory and Statistical Physics, vol. 33. Taylor and Francis, pp. 331–346, ISSN: 0041-145. Escriva´, A., March-Leuba, J., 2001. LAPUR 5.2, Verification and User’s Manual. NUREG/CR-6696, ORNL/TM-2000/340. Hangii, P., 2001. Investigating BWR stability with a new linear frequency domain method and detailed 3D neutronics. PhD dissertation, ETH14416 Zurich. Henry, A.F., 1980. Nuclear Reactor Analysis. MIT Press (second printing). Lefvert, T., 1996. Ringhals 1 Stability Benchmark–Final Report NEA/ NSC/DOC(96)22 (November 1996). Lewins, J., 1978. Nuclear Reactor Kinetics and Control. Pergamon Press. March-Leuba, J., 1984. Dynamic behaviour of boiling water reactors. PhD dissertation, University of Tennessee, Knoxville. March–Leuba, J., Otaduy, J.P., 1983. A Comparison of BWR Stability Mesurements with Calculations Using the Code LAPUR-IV, NUREG/CR-2998, ORNL/TM-8546. Oak Ridge National Laboratory (January). Munoz-Cobo, J.L., Escriva´, A., Garcı´a, C., 2006. Consistent generation and point feedback reactivity parameters for stability and thermal– hydraulic codes. Annals of Nuclear Energy. Otaduy, P.J., 1979. Modelling of the dynamic behaviour of large boiling water reactor cores. PhD Thesis, University of Florida. Peng, S.J., Podowski, M.Z., Lahey, R.T., Becker, M., 1985. NUFREQNP: A Compute Code for the Linear Stability Analysis of Boiling Water Nuclear Reactors NUREG/CR-4116. Stacey, W., 2001. Nuclear Reactor Physics. John Wiley and Sons. Tyrtysshnikov, E.E., 1997. A Brief Introduction to Numerical Analysis. Birkhauser (Chapter 9).