Influence of subcooled boiling on out-of-phase oscillations in boiling water reactors

Nuclear Engineering and Design 235 (2005) 1267–1282

Influence of subcooled boiling on out-of-phase oscillations in boiling water reactors J.L. Mu˜noz-Cobo a,∗ , S. Chiva b , A. Escriv´a a a

Department of Chemical and Nuclear Engineering, Polytechnic University of Valencia, P.O. Box 22 012, 46071 Valencia, Spain b Department of Technology, Fluid Mechanics Area, Jaume I University, Campus del Riu Sec, 12080 Castell´ on, Spain Received 13 January 2005; received in revised form 14 January 2005; accepted 31 January 2005

Abstract In this paper, we develop a reduced order model with modal kinetics for the study of the dynamic behavior of boiling water reactors. This model includes the subcooled boiling in the lower part of the reactor channels. New additional equations have been obtained for the following dynamics magnitudes: the effective inception length for subcooled boiling, the average void fraction in the subcooled boiling region, the average void fraction in the bulk-boiling region, the mass fluxes at the boiling boundary and the channel exit, respectively, and so on. Each channel has three nodes, one of liquid, one with subcooled boiling, and one with bulk boiling. The reduced order model includes also a modal kinetics with the fundamental mode and the first subcritical one, and two channels representing both halves of the reactor core. Also, in this paper, we perform a detailed study of the way to calculate the feedback reactivity parameters. The model displays out-of-phase oscillations when enough feedback gain is provided. The feedback gain that is necessary to self-sustain these oscillations is approximately one-half the gain that is needed when the subcooled boiling node is not included. © 2005 Elsevier B.V. All rights reserved.

1. Introduction The most advanced thermal–hydraulic codes for boiling water reactors, like RAMONA, TRAC-B, and TRAC-M, use a two fluid model, and solve the conservation equations of mass, energy, and momentum for each phase of the two-phase mixture (Wulff et al., 1984). These equations must be solved at each channel ∗ Corresponding author. Tel.: +34 963 877 631; fax: +34 963 877 639. E-mail address: [email protected] (J.L. Mu˜noz-Cobo).

0029-5493/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2005.01.018

of the reactor, and these reactor channels are coupled to the reactor vessel by proper boundary conditions. The channel thermal–hydraulics is modelled in 1D and neglecting transversal effects, however, the vessel is usually modelled in 3D. As a consequence, these codes are extremely time consuming, even with a small number of thermal–hydraulic nodes, and the application of these tools is very limited. A complementary effort has focused on the development of reduced order models, often consisting of a limited number of ordinary differential equations that represent the most important dynamical processes of a

1268

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

Nomenclature A cpl cpf Dh (r) fPm

channel area specific heat of the liquid specific heat of the liquid at saturation hydraulic diameter power distribution factor for the mth mode and the rth core region Fs fraction of the energy transferred to the channel that is invested in steam production Gin,r (t) mass flux at the inlet of channel r Gbb,r (t) mass flux at the bulk-boiling boundary of channel r Gex,r (t) mass flux at the exit of channel r G2φ  average mass flux in the two-phase region hin specific enthalpy of the liquid at the entrance of the channel hl,z1 specific enthalpy of the liquid at the boiling inception point hl,sb average specific enthalpy in the subcooled boiling region hf specific enthalpy of the liquid at saturation conditions (hA)ch,r fuel to coolant heat transfer coefficient times the heat transfer area at channel r H1φ single-phase heat transfer coefficient Hc length of the channels in the reactor core kl heat conductivity of the liquid K1,r form loss coefficient at the channel inlet plus form loss coefficients of the fuel rod spacers in the mono-phasic region K2,r form loss coefficient of the fuel rod spacers in the bi-phasic region Kex,r form loss coefficients at the channel exit nm (t) oscillating normalized components of the neutron flux expansion in lambda modes p pressure P power P0 steady-state power q˙ f,r heat generation rate fluctuation at core region r Q heat flux ˙ f,r Q heat generation rate at fuel core region r

˙ ch,r Q Tlinc

heat transfer rate to the fluid of channel r liquid temperature at the inception point for boiling Twinc wall temperature at the inception point for boiling V volume of the core V(r) volume of the region r of the reactor core (r) Wnm reactivity weighting distribution factors WKP square power weighting factors x dynamic quality of the two-phase flow mixture xbb (t) dynamic quality at the boiling boundary xex,r (t) dynamic quality at the exit of channel r Z1,r (t) effective boiling inception length at channel r Zbb,r (t) bulk-boiling boundary length at channel r Greek letters αbb (t) void fraction at the boiling boundary αbbR average void fraction at the bulk-boiling region αsb average void fraction in the subcooled boiling region β fraction of delayed neutron precursors p pressure drop Φ2 two-phase pressure drop multiplier λs bubble decay constant due to bubbles collapsing in the subcooled core of the channel λ disintegration constant of delayed neutron precursors Λ neutron generation time θ f,r lumped temperature fluctuations in the fuel of channel r ρg steam density ρfg difference between gas and liquid density at saturation conditions (ρg − ρf ) difference between gas and liquid denρlg sity at subcooled conditions (ρg − ρl ) ρl liquid density ρf liquid density at saturation conditions F ρnm feedback reactivity for the nth and mth modes ρl,sb average liquid density in the subcooled boiling region

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

BWR. These fast running codes give, in general, worse predictions than the 3D codes, but allow performing sensitivity analysis and give qualitatively correct predictions for a wide range of operating and design parameters. These reduced order models use generally (Van Bragt and Van der Hagen, 1998; Van Braght, 1998; Mu˜noz-Cobo et al., 2000; Lee and Onyemaechi, 1989) a two nodes model for each channel of the reactor core. The first node is a single-phase node, where the liquid is uniformly heated from the channel entrance up to the boiling boundary where saturation is reached. The following node is a two-phase node, where the two phases are in thermodynamic equilibrium and all the heat transferred to the channel is invested in steam production. However, in real cases there is a region called the subcooled boiling region, where there is not thermodynamic equilibrium between the phases, and part of the power transferred to the channel is invested to heat the liquid, while the rest is invested in steam production. In this subcooled boiling region, the bubbles are formed at the walls (clad surface), and detach from the walls to travel to the subcooled region of the channel, where they collapse. Because reactivity feedback is very sensitive to void variations, and is axially weighted by a distribution factor that depends on the square of the power distribution, then, the void fraction variations in the subcooled boiling region will have a big contribution to the void feedback reactivity. In this paper, we have developed a reduced order model where each channel is represented by three nodes: one subcooled liquid node from the channel entrance to the boiling inception point, where subcooled boiling starts; the second node starts at the inception of subcooled boiling and ends at the boiling boundary, where the average enthalpy of the two-phase mixture attains saturation conditions; finally, the last node, or bulk-boiling region, extends from the bulk-boiling boundary to the channel exit. The goal of this paper is to study the influence of the subcooled boiling on the feedback mechanisms that lead to the development of out-of-phase instabilities in boiling water reactors, and to make a consistent lumped model that includes the subcooled boiling in a realistic way. The studies performed on the void feedback reactivity (Wulff et al., 1984; Mu˜noz-Cobo et al., 1994) have concluded that the typical functionalization of the void reactivity as a second-degree polynomial in the void

1269

fraction, must be multiplied by a weighting factor that depends on the square of the power distribution. Therefore, the void feedback reactivity is enhanced in the reactor regions where the power is higher. The new fuel designs tend to make the power distribution more peaked at the reactor bottom, or lower part of the core. Therefore, the reactivity feedback at the subcooled boiling region will become more important with these new designs because is enhanced by the bottom peaked axial power distribution. This is the main reason to add a subcooling boiling node to the classical reduced order models with only two nodes, and where the reactivity feedback in the subcooled boiling region is not considered. As it is well known, the out-of-phase instabilities appear because one or two subcritical modes are excited. The mechanism to excite these modes is to provide enough reactivity feedback to overcome the eigenvalue separation between the fundamental mode and the subcritical one (Mu˜noz-Cobo et al., 2000). For out-ofphase oscillations, the inlet mass flow rate to the reactor core remains constant, and the two oscillating core regions adjust their flows to maintain approximately constant the pressure drop across the core (March-Leuba and Rey, 1993). The paper has been organized as follows: Section 2 is devoted to the development of the model equations. Section 3 is devoted to the results of the model and the discussion about the influence of subcooling boiling on the onset of the instabilities. Finally, Section 4 is devoted to analyze the main conclusions of this paper.

2. The reduced order model of a BWR with subcooled boiling In this section, we explain the main characteristics of the DWOS M SU (density wave oscillation with modal kinetics and subcooled boiling) lumped parameter model. One of the defects of the previous lumped models is that they do not include the subcooling void fraction. Therefore, the void fraction along the channel is smaller in previous lumped models than the real one, mainly at the beginning of the channel. This effect has a strong influence on the void reactivity feedback, which depends on the void fraction and the power distribution through the square power distribution factor. Now, in

1270

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

Fig. 1. DWOS lumped parameter model. This model is a two channel model, with three regions at each channel, the liquid region from Z = 0 to the boiling inception height Z1 the subcooled boiling region from Z1 , to the boiling boundary length Zbb , and the bulk-boiling region from Zbb to the channel exit.

the lower part of the channel, the power is higher due to the fact that the density of thermal neutrons is higher. This is a direct consequence of the bigger moderation of fast neutrons produced by the higher density of the water (less voids) in this region. Therefore, the square power distribution factor is higher in the lower part of the channel and, therefore, the consideration of the subcooled boiling will have a strong influence on the void reactivity feedback. The core fuel assemblies are grouped into two core regions that are represented by two averaged channels. The DWOS lumped parameter model of each channel has three regions or nodes, as it is displayed in Fig. 1. 2.1. Calculation of the effective inception temperature Tlinc , and the dynamic equation of the inception length Z1 for subcooled boiling The most important thing, in any effective subcooling boiling model, is the ability of the model to be able to predict accurately where significant void fraction appears, this location of the void departure will be denoted by Z1 . Several criteria can be used to predict the inception point (Lahey and Moody, 1993). However, we have checked several criteria for the effective inception point and the differences were very small, so we have chosen a criterion that is obtained equating the

single-phase forced convection heat flux, at Z = Z1 , to a Jens–Lottes type heat flux. This criterion has been successfully used by the LAPUR code (Otaduy, 1979), for many years. This criterion lead to the following expression for the liquid temperature at the inception point, Tlinc , in terms of the heat flux, Q , at the fuel walls: Tlinc = C2−0.25 Q

0.25

+ Tsat −

Q , H1φ

(1)

where H1φ is the single-phase forced convection heat transfer coefficient, given by the Dittus–Boelter formula. Tsat is the liquid saturation temperature, and the constant C2 is given by:
 4p C2 = 2.5454 exp . (2) 6, 207, 385 Therefore, the liquid enthalpy at Z1 is given by: hlz1 = hf − cpl (Tsat − Tlinc ).

(3)

Therefore, at steady-state conditions, the effective inception boundary length can be computed by means of the following obvious formula, deduced assuming the uniform heating of the channel: Z1,0 =

(hlz1 − hin )Gin,0 A , ˙ ch,0 /Hc Q

(4)

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

where hin is the enthalpy of the subcooled liquid at the entrance of the channel, A the channel flow area, Gin,0 the mass flux at steady-state conditions, Qch,0 the power transferred to the channel fluid at steady-state conditions, and finally, Hc is the total channel length. To get the dynamical equation governing the evolution of Z1,r , with time at the rth channel, we integrate the energy conservation equation in the single-phase region. This calculation yields: ˙ ch,r dZ1,r 2Q Z1,r = dt AHc ρl,in hin − ρl,z1 hl,z1 + 2Gin,r

hin − hl,z1 ρl,in hin − ρl,z1 hl,z1

(5)

where Qch,r is the heat transferred per unit time to the channel fluid, ρl,z1 the liquid density at the inception point Z1 , Gin,r the mass flux at the entrance of the rth channel, and finally, hl,z1 is the liquid enthalpy at the boundary of the effective inception point for subcooled boiling. 2.2. Dynamic equation governing the average void fraction in the subcooled boiling region To obtain the dynamic equation for the average void fraction in the subcooled boiling region, we start from the mass conservation equation of the steam in this region, given by: ˙ Fs Q ∂ ∂ ch,r (ρg αr ) = −λs ρg αr − (xr Gr ) + ∂t ∂z Ahfg

(6)

where αr is the void fraction in the rth channel, Gr the mass flux in the rth channel, λs the bubble decay constant due to the bubble collapse in the subcooled region of the channel, and finally, Fs is the fraction of the energy transferred to the channel that is invested in steam formation. Integration of Eq. (6), between the subcooled boiling inception length Z1 and the bulk-boiling boundary length Zbb of the rth channel, yields:  Zbb,r ∂ dz (αρg ) ∂t Z1,r = −λs ρg (Zbb,r − Z1,r )αsb,r − xbb,r Gbb,r (t) ˙ ch,r F¯ s Q + (Zbb,r − Z1,r ), Ahfg Hc

where αsb,r is the average void fraction in the subcooled boiling region of the rth channel, xbb,r the dynamic quality at the boiling boundary of the rth channel, and finally, Fs  is the average fraction of the energy transferred to the channel that is invested in steam formation in the subcooled boiling region. Now, we apply the Leibnitz rule to the left-hand side of Eq. (7), and on account of the definition of the average void fraction in the subcooled-boiling region, given by:  Zbb,r 1 αsb,r = dz α. (8) Zbb,r − Z1,r Z1,r It is obtained, from Eq. (7), the following equation for the evolution of the average void fraction, αsb,r , in the subcooled boiling region of the rth channel: αsb,r d d Z1,r αsb,r = Zbb,r − Z1,r dt dt +

αbb,r − αsb,r d Zbb,r − λs αsb,r Zbb,r − Z1,r dt

−

˙ ch,r xbb,r Gbb,r Fs Q + . (Zbb,r − Z1,r )ρg Ahfg Hc ρg

(9)

We observe that the evolution of the average void fraction in the subcooled region depends positively on the average of the transferred energy invested in steam formation, and negatively of the collapsing of bubbles term that is proportional to the average void fraction in this region. 2.3. Bulk-boiling boundary dynamics The equation for the dynamic behavior of the boiling boundary length, Zbb,r (t), at the rth representative channel is given by (Van Bragt and Van der Hagen, 1998; Lee and Onyemaechi, 1989): ˙ ch,r Zbb,r (t) dZbb,r (t) 2Gin,r (t) 2Q = − dt ρf (hf − hin )AHc ρf

(10)

where Qch,r is the fuel to coolant heat transfer rate per fuel assembly of the rth core region, which can be written as: ˙ ch,r = P0,chr (t) + (hA)ch,r θf,r , Q

(7)

1271

(11)

where P0,chr is the steady-state power per fuel assembly in the rth channel, (hA)ch,r the fuel to coolant heat trans-

1272

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

fer coefficient times the heat transfer area in a typical or average channel of the rth region, and θ f,r denotes the lumped temperature fluctuation in the fuel of the rth channel. 2.4. The equation for the mass flux at the boiling boundary The continuity equation in the subcooled boiling region is given by:

Gbb,r (t) =

[(ρg hg − ρl,sb hl,sb ) − (ρlg )sb hl,z1 ]r [ρg hg − ρl,sb hl,sb − ρlg,sb (hl,bb + xbb hlg,bb )]r −

∂ ∂ (αr ρg + (1 − αr )ρl ) = − Gr . ∂t ∂z

where Gz1,r is the liquid mass flow rate at the inception subcooling boiling boundary, hl,z1 the liquid enthalpy at the inception point for subcooled boiling, hl,bb the enthalpy of the liquid at the boiling boundary. We have assumed that ρl,sb hl,sb is approximately time independent. Next, we eliminate the integral in Eq. (15), with the help of Eq. (13), and we make the approximation Gin,r ∼ = Gz1,r . From the resulting equation, and after some algebra, it is obtained the following expression for the mass flux at the boiling boundary: Gin,r (t)

˙ ch (Zbb − Z1 )] [ρlg,sb Q r AHc [ρg hg − ρl,sb hl,sb − ρlg,sb (hl,bb + xbb hlg,bb )]r

(16)

where we have introduced the following definitions: (12)

ρlg,sb = ρg − ρl,sb ,

(17)

The integration of Eq. (12), between the subcooled boiling inception boundary length, Z1,r , and the boiling boundary length, Zbb,r , yields:  Zbb,r ∂α Gin,r − Gbb,r dz , (13) = ∂t ρg − ρl,sb Z1,r

ρlg,bb = ρg − ρl,bb ,

(18)

hlg,bb = hg − hl,bb .

(19)

where ρl,sb is the average density of the liquid in the subcooling boiling region, and we have made the approximation Gz1,r ∼ = Gin,r . Finally, Gbb,r is the total mass flux at the boiling boundary length of the rth channel. The next step is to integrate the energy equation from Z1,r to Zbb,r . The two-phase energy equation is given by: ∂ ((1 − α)ρl hl + αρg hg ) ∂t ˙ ch ∂ Q = − ((1 − x)Ghl + xGhg ) + (14) ∂z AHc Assuming that that the channel is uniformly heated and integrating Eq. (14) from Z1,r to Zbb,r it is obtained:  Zbb,r ∂α (ρg hg − ρl,sb hl,sb ) dz ∂t Z1,r = (Gz1 hl,z1 )r − {(1 − xbb )hl,bb + xbb hg }r Gbb,r (t) +

˙ ch,r Q (Zbb,r − Z1,r ) Ar H c

(15)

Eq. (16) relates the mass flux at the boiling boundary with the mass flux at the channel inlet and the rate of heat transfer to the subcooled boiling region. 2.5. The equation of the mass flux at the channel exit First, we integrate the mass conservation Eq. (12) from the boiling boundary length Zbb,r to the total height Hc of the channel, this calculation yields:  Hc ∂α Gbb,r − Gex,r dz = . (20) ∂t ρg − ρ f Zbb,r Then, we integrate the energy conservation equation between the same limits; this calculation yields:  Zbb,r ∂α dz (ρg hg − ρf hf ) ∂t Z1,r = {(1 − xbb )hl,bb + xbb hg }r Gbb,r (t) − {(1 − x)ex hf + xex hg }r Gex,r +

˙ ch,r Q (Hc − Zbb,r ). Ar H c

(21)

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

Now, we eliminate the integral in Eq. (21), with the help of Eq. (20). Then, after some algebra, it was obtained the following expression for the mass flux at the channel exit: Gex,r (t) =

1273

where p = pin − pex is the difference between the inlet and outlet channel pressures, that is channel independent, because all the channels have commons lower and upper plena.

[((1 − xbb )hl,bb + xbb hg )ρfg − (ρg hg − ρf hf )]r Gbb,r (t) {[(1 − x)ex hf + xex hg ]ρfg − (ρg hg − ρf hf )}r +

˙ ch (Hc − Zbb )] [ρfg Q r AHc {[(1 − x)ex hf + xex hg ]ρfg − (ρg hg − ρf hf )}r

Eq. (22) gives the mass flux at the exit of the channel in terms of the mass flux at the boiling boundary and the heat transfer rate to the bulk-boiling region. 2.6. Equation for the evolution of the average void fraction in the bulk-boiling region (bbR) The equation for the evolution of the average void fraction, in the bulk-boiling region, is easily obtained by integrating the mass conservation equation between Zbb,r and Hc . This calculation yields after some arrangements:

pacc,r represents the pressure drop due to the fluid acceleration in the channel:   2 G2 xex (1 − xex )2 G2ex Gin,r ex pacc,r = + − αex ρg (1 − αex )ρf ρl,in r (25) pg,r represents the gravitational pressure drop, which is given by the following expression: pg,r = gρl Z1,r (t) + [g(1 − αsb,r )ρl,sb + gαsb,r ρg ](Zbb,r (t) − Z1 (t)) + [g(1 − αbbR,r )ρf + gαbbR,r ρg ]

d Gbb,r (t) − Gex,r (t) αbbR,r = dt ρfg (Hc − Zbb,r ) (αbb − αbbR )r dZbb,r − , Hc − Zbb,r dt

(22)

× (Hc − ZbbR,r (t)), (23)

where the mass fluxes Gbb,r and Gex,r , at the region boundaries, are easily calculated by means of Eqs. (16) and (22), respectively. 2.7. Momentum conservation equation A key assumption in many thermal-hydraulic models of out-of-phase oscillations is that the pressure drop across the channels remains approximately constant during the out-of-phase oscillation transient (Van Bragt and Van der Hagen, 1998; Van Braght, 1998; Mu˜nozCobo et al., 2000; Lee and Onyemaechi, 1989; MarchLeuba and Rey, 1993). Assuming that the channel area is constant, we get, integrating the momentum equation along the channel, the following result:  H ∂Gr (z, t) dz = p − pacc,r − pg,r − pf,r , ∂t 0 (24)

(26)

where αsb,r and αbbR,r are the average void fractions, in the subcooled boiling region and in the bulkboiling region of the rth channel, respectively. Finally, pf,r gives the friction pressure drop, which is calculated by means of the following expression:
 2 Z1,r (t) Gin,r (t) pf,r = K1,r + fr 2ρl DH + K2,r Φ2r ΩJ r + fr

G2φ 2r 2ρf

G2φ 2r Hc − Z1,r (t) Φ2r ΩJ r DH 2ρf

+ Kex,r Φ2ex,r ΩJ r

G2ex,r (t) , 2ρf

(27)

where the first term of expression (27) gives the pressure drop in the single-phase region of the rth channel, with K1,r being the form loss coefficient due to the losses at the channel inlet and the fuel rod spacers, located in the single-phase region. The second term gives

1274

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

the channel pressure drop due to the rod spacers located in the biphasic region. The third term gives the pressure drop due to friction with the walls. Finally, the last term gives the pressure drop at the channel exit. In these terms, Φ2 is the two-phase pressure drop multiplier, and Ωr is the Jones and Dight (1963) correction factor. To compute the pressure drop in the channels, we neglect the time derivative term (Lee and Onyemaechi, 1989), and we assume that the average flux in the biphasic region at a given time can be approximated by: Gr = 0.5(Gin,r (t) + Gex,r (t)).

(28)

Therefore, we write p0 = pacc,r + pg,r + pf,r .

2.8.1. Fraction of power invested in steam production in the subcooled boiling region According to Lahey and Moody (1993) and Otaduy (1979), the energy transferred from the fuel to the subcooled boiling fluid can be split into three components: (i) formation of steam bubbles near the heating surface which may detach into the main flow stream, (ii) pumping of the liquid mass out of the control volume by the expanding action of the steam bubble formation, (iii) single-phase convective heating through the parts of the heating surface not generating bubbles. The investigations performed have concluded that the steam formation and the pumping process are predominant over the normal convective process. Therefore, we can write:  qevap  qevap

 + qpump

=

1 1+

 qpump  qevap

ε=

 qpump  qevap

=

ρl (hf − hl ) . ρg hfg

(31)

On account of Eqs. (30) and (31), and the LAPUR code (Otaduy, 1979), we use for Fs the following expression: Fs =

1

,

(32)

hf −hl hfg

gives the degree of subcooling of

Hl 1 + fp 1−η

where Hl =

s the liquid, and η = ρf ρ−ρ . Finally, fp is an adjustable f factor, equal to 1.3, that will allow for better to correlate predictions with measurements.

(29)

2.8. Constitutive equations and calculation procedures

Fs =

and therefore we have:

,

(30)

 is the part of the energy flux at the surwhere qevap  is face associated to the steam formation, while qpump the part of the energy flux associated to the pumping process. According to Rouhani and Axelsson (1970), the ratio of the heat fluxes, due to pumping and to steam formation, can be calculated with the approximation that the liquid that leaves the control volume is at saturation

2.8.2. Steam decay constant in the subcooled boiling region The model used is based in the model of Jones and Dight (1963) and LAPUR (Otaduy, 1979), the results of this model have been compared with the expression used by the RELAP5 code that gives similar results. The Jones model uses the following expression for the bubble decay ratio: λs = cλ0 φHl2 ,

(33)

where Hl is the degree of subcooling of the liquid, φ and λ0 are given by the following expressions:
2 hfg Hw2 φ= and λ0 = , (34) cpf (Tc0 − Tsat ) kf ρl cpf where Tc0 is the clad temperature at the inception of subcooled boiling, Hw the single-phase heat transfer coefficient, and cpf is the specific heat of water at saturation. Finally, we have checked that assuming the parameter c equal to 0.005 gives values of the subcooled boiling void fraction that are close to the experimental ones, when only one subcooled node is used. The problem is that when we have many nodes in the subcooled boiling region, for instance, 50, then in the three or four nodes near the inception point, λs is very big and then becomes much more smaller in the rest of nodes. We have checked that when using c equal to 0.005, we get a value for the average subcooling void fraction that is very close to the average value obtained using 30 nodes in this region.

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

2.8.3. Relation between the quality and the void fraction To obtain the value of the quality at the boiling boundary, xbb (t), and at the channel exit, xex (t), we start from the formula that relates the void fraction with the quality (Thom, 1964): α=

γx 1 + (γ − 1)x

(35)

where in the HEM model (Todreas and Kazimi, 1990) γ = vg /vl . The average void fraction in the subcooled boiling region of the rth channel can be obtained as follows:  xbb,r (t) 1 α(t)sb,r = α dx. (36) xbb,r (t) 0 Direct substitution of Eq. (35) in (36) yields:  γ 1 α(t)sb,r = 1− γ −1 (γ − 1)xbb,r (t)  × ln(1 + (γ−1)xbb,r (t)) . (37) The average void fraction in the bulk-boiling region can be easily obtained as follows:  xex,r (t) 1 α(t)bbR,r = α dx. (38) xex,r (t) − xbb,r (t) xbb,r (t) Direct substitution of the expression (35), which relates the void fraction with the quality, followed by integration yields:  γ 1 α(t)bbR,r = 1− γ−1 (γ − 1)(xexit,r (t) − xbb,r (t))  1 + (γ − 1)xex,r × ln . (39) 1 + (γ − 1)xbb,r Eq. (39) gives the average void fraction in the bulkboiling region, in terms of the quality at the channel exit and the quality at the boiling boundary. When solving the set of dynamics equations, we get at each time step the average void fractions in the subcooled boiling region, and in the bulk-boiling region. To get the mass fluxes Gbb,r (t) and Gexit,r (t), we need to know the qualities at the boiling boundary, xbb (t), and at the channel exit, xexit (t). We get these qualities by iteration in Eqs. (37) and (39).

1275

2.9. Neutronic model and feedback reactivity The normalized components nm (t) of the neutron flux expansion in λ modes obey, for the case of only two modes, the following set of coupled differential equations (Mu˜noz-Cobo et al., 2000; Hashimoto, 1993): ρF (t) ρF (t) − β dn0 n0 (t) + 00 = 00 dt Λ0 Λ0 F (t) ρ01 + n1 (t) + λc0 (t) Λ0 F (t) − β ρF (t) ρs + ρ11 dn1 n1 (t) + 10 = 1 dt Λ1 Λ1 , F ρ (t) + 10 n0 (t) + λc1 (t) Λ1 dc0 β = n0 (t) − λc0 (t) dt Λ0 dc1 β n1 (t) − λc1 (t) = dt Λ1

(40)

where cm (t) are the oscillating normalized components of the expansion of the delayed neutron precursor concentrations in terms of harmonic λ-modes, Λi the neutron generation time for the ith mode, ρ1s the subcritical F are reactivity of the first subcritical mode. Finally, ρmn the feedback reactivities for a given configuration of the reactor core, and a given position of the reactor control rods. The main feedback reactivity contributions are the void and Doppler feedback reactivities, so we write: F V D ρmn (t) = ρmn (t) + ρmn (t).

(41)

V (t) are The modal void feedback reactivities ρmn computed adding up the void reactivity contributions of the various reactor core regions: V ρmn (t) =

 r

V,r ρmn (t),

(42)

V,r (t) from the where the void reactivity contribution ρmn rth core region can be obtained from the reactivity (r) weighting distribution factors Wm,n (Mu˜noz-Cobo et al., 2000) as follows: V,r r (r) ρmn (t) = RV 00 (α )Wm,n ,

(43)

1276

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

r where RV 00 (α )is the scaled void reactivity for the rth subcore region defined by:

r RV 00 (α )

= =

Φ∗0 ,

 K

∂M

(r)  ∂L ∂α − ∂α δαΦ0 ∗ (r) Φ0 , M0 Φ0

(C1 + C2 αrk

+ C3 (αrk )2 )WKP δαrK ,

Cofrentes case 1 instability event 1991Ringhals test

(44)

where M and L are the production and removal operators, and Φ0 and Φ∗0 are the direct and adjoint fluxes, · means integration over a given region. We have assumed that the ratio of Eq. (44) for the rth core region scales with the void fraction αrk as in the full core, i.e. with the same polynomial fit constants C1 , C2 , and C3 , the index k denotes the axial nodes in the fitting. WKP are the typical square power weighting factors: WKP

P2 = K2 K PK

(45)

The modal reactivity weighting factors can be approximated by (Mu˜noz-Cobo et al., 2000):

(r) Wmn

(r) Φ∗m , Φn = Φ∗m , Φm

(46)

The Doppler reactivity contribution of the rth core region can be obtained by a similar method (Mu˜nozCobo et al., 2000). The reactivity coefficients C1 , C2 , C3 , displayed in Table 1, have been obtained by means of a consistent multidimensional methodology (Mu˜noz-Cobo et al., 1994). The modal reactivity weighting factors, defined at expression (46) were computed with the code LAMBDA (Mir´o, 2002), and the results are displayed in Table 2. The LAMBDA code permits us to solve the 3Deigenvalue equation: LΦn = λn MΦn ,

(47)

Table 1 Void feedback reactivity coefficients for Cofrentes event January 1991 C1 C2 C3

−0.23256 0.64433 −0.92516

Table 2 Reactivity weighting factors for a two regions model with two modes the fundamental one and the first harmonic

1 = 0.4995 W 2 = 0.5005 W00 00

1 = 0.4985 W 2 = 0.5014 W00 00

1 = −0.4529W 2 = 0.4529 W01 01

1 = −0.4299W 2 = 0.4312 W01 01

1 W10 1 W11

= 0.4541

1 = −0.4321W 2 = 0.4298 W10 10

= 0.4901

1 = 0.4748 W 2 = 0.4783 W11 11

2 = −0.4541W10 2 = 0.4918 W11

where L is the differential operator for the absorption and leakage of neutrons given by:  

· (D1 ∇)

+ Σa1 + Σ12 −∇ 0 L= ;

· (D1 ∇)

+ Σa2 −Σ12 −∇ (48) D1 and D2 are the diffusion coefficients for the fast and thermal groups, respectively; Σ a1 and Σ a2 are the absorption cross-sections of the fast and thermal group, respectively; Σ 12 = Σ r is the removal or downscattering cross-section from the fast to thermal group; M is the two-group production operator   νΣf1 νΣf2 M= ; (49) 0 0 Φ = column[Φ1 , Φ2 ] and λn = 1/kn are the nth eigenvector and the nth Lambda eigenvalue, respectively. The method used by the LAMBDA code is to transform the eigenvalue problem, associated to the differential operator L, into an algebraic eigenvalue problem. This step is performed by discretizing the space in parallelepiped cells, followed by a Legendre expansion of the fluxes in the cells. This method was originally developed by Herbert (1987), and applied successfully by Ginestar et al. (2002), and belongs to the class of nodal collocation methods. The previous method allows us to compute the eigenvalues and eigenvectors of the two group diffusion equation in three dimensions. In the Case of Cofrentes NPP, the application of this method for the determination of the eigenvalues of the out-of-phase modes yields the results displayed at Table 3. Finally, in Fig. 2, we display the first subcritical harmonic mode of Cofrentes Reactor computed for the conditions of the instability event of 1991, or case 1 of the report D16f (Escriv´a et al., 2003).

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282 Table 3 Eigenvalues of the Cofrentes nuclear reactor and residual errors for the case 1 of report D16f (Escriv´a et al., 2003), with 30 × 30 × 27 nodes and a Legendre polynomial expansion of second order

k0 k1 k2

Eigenvalues

Residual error

0.99724 0.99090 0.98850

2.23 × 10−7 1.43 × 10−4

This result shows us that if the first harmonic mode is excited one-half of the reactor core will increase its power while the other half will decrease its power. In Table 2, we display the reactivity weighting factors of Cofrentes NPP, when the reactor core is divided in two regions and only the first harmonic mode is excited. We display also in Table 2, the values of the reactivity weighting factors computed for Ringhals test (Mu˜noz-Cobo et al., 2000). We can see that in both nuclear plants the reactivity weighting factors have similar values. Let us analyze first the physical meaning of this reactivity weighting factors. The physical meaning of the reactivity weighting factors is obvious from Table 2. For case 1, the core is split in two regions, approximately of the same size. In region (1), the first harmonic mode is negative, while in region (2), is positive. Because both region have the same size and are practically symmetrical, the contribution of each region to the total reactivity of the fundamental mode is approx-

1277

imately one-half. Therefore, the reactivity weighting factors, for the fundamental mode should be approximately equal to 0.5. This is the value obtained for Cofrentes and Ringhals, using Eq. (46), and the fundamental eigen-modes for both reactors computed with the LAMBDA code. For a cubic reactor core with two regions, an elemental calculation shows that the reactivity weighting (1) (2) factors W01 and W01 are equal to −4/(3π) = −0.424 and 4/(3π) = 0.424. Obviously, the reactors are not cubic but these values are very close to the values displayed in Table 2, with the true geometry and the true eigen-modes. Because for a given rth region the model only computes the average void fraction in that particular region, then we make the approximation: αrk (t) = αk0 + δαr (t),

(50)

i.e. we have assumed that the void fraction perturbations are the same in all the axial nodes belonging to the same region. αk0 is the axial steady-state void fraction distribution, computed with a 1D average channel model for LAPURX (Otaduy, 1979). We have computed this distribution with an average channel model with 25 axial nodes. Finally, the total reactivity ρmn is written in the form:  F r D r r r = Kgmn (RV ρmn (51) 00 (α ) + R00 (α , Tf ))Wmn , r

Fig. 2. First subcritical harmonic mode for Cofrentes NPP. Case 1 of the Nacusp Project (Escriv´a et al., 2003).

1278

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

where Kgmn is a bifurcation parameter to increase the r r feedback reactivity of the system, and RD 00 (α , Tf ) is the scaled Doppler reactivity (Mu˜noz-Cobo et al., 2000). 2.10. Power generation and fuel dynamics The power generated in a given region can be obtained integrating the reactor heat generation rate per unit volume over the volume of this region. This calculation yields (Mu˜noz-Cobo et al., 2000): ˙ f,r (t) = P0 f (r) + Q P0

∞ 

(r)

P0 fPm nm (t)

m=0

˙ (r) (t) + q˙ f,r (t), =Q f,0

(52)

Cobo et al., 2000): cf,r Mf,r

d θf,r (t) = q˙ fch,r (t) − (hA)ch,r θf,r (t), dt

(54)

where Mf,r is the fuel mass contained in a fuel channel assembly of the rth region, q˙ fch,r the fluctuation in the heat generation rate in the fuel of one fuel channel assembly of the rth region, cf,r the specific heat of the fuel, and (hA)ch,r is the effective heat transfer coefficient times the heat transfer area. The fluctuation q˙ fch,r (t) in the heat generation rate for one fuel channel assembly can be obtained from expression (52), dividing this expression by the number of fuel channel assemblies Nch (r), of a given region r. This calculation yields:

∞ (r) P0 fPm nm (t) q˙ fch,r (t) = m=0 (55) Nchan (r)

(r)

where P0 is the steady-state power, and fPm are the power distribution factors for the mth mode, which are defined as follows:  V (r) (Σf1 Φm1 + Σf2 Φm2 ) dV (r) fPm =  , (53) V (Σf1 Φ01 + Σf2 Φ02 ) dV where V(r) is the volume of the rth region of the reactor core and V is the total volume of the core, Σ f1 and Σ f2 are the fission cross-sections for the fast and thermal group, respectively. Φm1 and Φm2 are the fast and thermal components of the mth lambda eigenfunction. These coefficients have been computed with the code LAMDA (Mir´o, 2002), and the result of the calculation is displayed in Table 4. Dividing the power generated in this region by the number of fuel assemblies of this region we get the power generated per fuel channel assembly of the rth ˙ fch,r . Then, core region, this magnitude is denoted by Q the governing equation for the average fuel temperature fluctuation θ f,r , in a fuel channel assembly of the rth core region, is given by the following equation (Mu˜noz-

Table 4 Power distribution factors of Cofrentes NPP for case 1, when the core is divided in two region and we consider a model with only two modes Region 1

Region 2

(1)

fP0 = 0.4995

(1)

fP1 = 0.4960

Fundamental mode

fP0 = 0.50042

First harmonic mode

fP1 = −0.5039

(2) (2)

3. Results and discussion 3.1. Steady-state results of DWOS model for Cofrentes NPP The steady-state and the dynamic equations described in this paper were implemented in a FORTRAN code denoted: DWOS M SU. The input data for the Cofrentes model were obtained from (Escriv´a et al., 2003). In Table 5, we display the main parameter values of this model. Then, we solve the steady-state equations with the parameter values given in Table 5, the main results of the DWOS steady-state calculations are displayed in Table 6. Also, in this table, we display some results of the SIMULATE code for the same case. Table 5 Parameter values used for Cofrentes instability event Parameter

Value

Parameter

Value

p (Pa) hin (J/kg) Hc (m) A (m2 ) P0 (W) Gin0 (kg/m2 s) Nch Tl,in K1 K2

6.6478 × 106

Kex f Dh ρ1s Λ (s) λ β τ (s) cf M f

0.9 0.0266 0.01306 −0.006398 3.591 × 10−4 6.636 × 10−2 0.556 × 10−2 4.11 82270

0.111 × 107 3.81 0.009783 1.121 × 109 515.97 624 529 57.43 2.81

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282 Table 6 Steady-state effective inception length Z1 , bulk-boiling boundary length Zbb , average void fraction in the subcooled boiling region αsb , average void fraction in the bulk-boiling region αbbR , average void fraction in the biphasic region α2φ , average void fraction in the channel αc , quality at the exit of the channel, xex , computed with DWOS, quality in the upper plenum computed with SIMULATE, xUP(SIMULATE) Z1,0 Zbb,0 αsb,0 αbbR,0 αc α2φ xex xUP(SIMULATE)

0.345 1.479 0.216 0.67 0.47 0.522 0.169 0.146

We observe that the quality at the exit of the channel computed with DWOS is 0.169. We have displayed the quality at the upper plenum computed by the SIMULATE code for this same case; however, this quality is smaller due to the fact that, in the upper plenum, the two-phase mixture coming from the fuel channel assemblies mixes with the water coming from the bypass, this fact reduces the channel exit quality about a 10%, so the DWOS result and the SIMULATE result agrees. Then, in Table 7, we display the pressure drop in the channel computed with DWOS and SIMULATE at steady-state conditions. We observe that the pressure drops computed with both codes are very similar in spite of the simplifications of the DWOS code. Also, we observe that the effective inception point for subcooled boiling starts at 0.345 m from the beginning of the channel, with an average void fraction of 0.216. This means that the subcooled boiling region has an average void fraction that cannot be neglected in reduced order models. Table 7 Total pressure drop in the channel computed with DWOS, p, total pressure drop in the channel computed with SIMULATE, pSIMULATE , acceleration pressure drop, pacc , gravity pressure drop, pg , friction pressure drop, pf p (Pa) pSIMULATE (Pa) pacc (Pa) pg (Pa) pf (Pa)

0.380 × 105 0.392 × 105 0.127 × 104 0.155 × 105 0.211 × 105

1279

3.2. Transient results To simplify the model, we have reduced the number of feedback gain parameters to only two (Mu˜noz-Cobo et al., 2000), the feedback gain for the fundamental mode denoted by Kg0 : Kg0 = Kg00 = Kg01

(56)

And the feedback gain for the first harmonic reactivity denoted by Kg1 : Kg1 = Kg11 = Kg10

(57)

The feedback reactivity coefficients of Table 4 were calculated for a model with 25 axial nodes, following a consistent methodology (Mu˜noz-Cobo et al., 1994). Because in this particular model, the number of axial nodes is three (one for the subcooled liquid region, one for the subcooled boiling region, and one for the bulk-boiling region), it is expected that the feedback gain necessary to achieve limit cycle oscillations should be bigger than one. These results have been recently proved by Ginestar et al. (1999); these authors have proved that with the reduction of the number of nodes, the feedback gain must be increased to get the critical value. To get self-sustained out-of-phase nuclear-coupled density wave oscillations, we fix the feedback gain of the fundamental mode at a value of one, and we increase the reactivity feedback gain Kg1 of the first harmonic mode. The critical value is attained at Kg1 = 3.1. In this case, the model displays out-of-phase oscillations with a period of 2.95 s, and a frequency of 0.34 Hz. At this point, we must remark that the feedback gain necessary to achieve out-of-phase oscillations, when the subcooling boiling is not included in the DWOS model, is more than two times the gain that it is necessary when the subcooling boiling is included. Let us analyze now some results of the DWOS code. Fig. 3 displays the bulk-boiling boundary lengths Zbb,1 and Zbb,2 , when out-of-phase oscillations are excited with a feedback gain of 3.2. We observe that when in one-half of the reactor core the bulk-boiling boundary length Zbb,1 attains its maximum value, in the other half of the rector core the bulk-boiling boundary length Zbb,2 reaches its minimum value. Other important parameters are the mass fluxes (kg/m2 s) at channel 1, we have observed that the mass

1280

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

Fig. 3. Bulk-boiling boundary lengths, Zbb,1 and Zbb,2 , vs. time when out-of-phase oscillation are excited with a feedback gain of 3.2.

fluxes Gbb,1 (t) and Gex,1 (t), at the bulk-boiling boundary of channel 1 and at the exit of channel 1, are delayed with respect to the mass flux at the entrance of channel 1. Also, we have noticed that the mass flux Gex,1 , at the exit of channel 1, oscillates with a delay of 180◦ with respect to the mass flux at the entrance of this same channel. This behavior is typical of density wave oscillations in unstable channels. This behavior is displayed in Fig. 4, which shows the mass flux at the entrance and exit of channel 1 versus time.

Next, we must remark that the average void fraction in the bulk-boiling region of channel 1 oscillates outof-phase of the average void fraction in the bulk-boiling region of channel 2. This behavior is typical of out-ofphase instabilities. Finally, the powers transferred to channels 1 and 2 are out-of-phase, i.e. while in one-half of the reactor core the power attains its maximum value, at the other half of the reactor core the power attains its minimum value. This is a direct consequence of the

Fig. 4. Mass fluxes at the inlet of channel 1, Gin,1 (t), and the exit of channel 1, Gex,1 (t), vs. time.

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

oscillation of the normalized amplitude of the first subcritical mode n1 (t) with time and a period of 2.95 s, equivalent to a frequency of 0.34 Hz. Because for out-of-phase oscillations n1 (t) is bigger than n0 (t), then the power in both halves of the reactor core oscillates out-of-phase. However, the fundamental mode still has oscillations of small amplitude due to the driven term F ρ01 Λ n1 (t).

4. Conclusions The fact that parallel channels of a BWR have common lower and upper plena impose the same pressure drop p to all the reactor fuel channel assemblies, but this pressure drop can change with time if the number of channels is low (Mu˜noz-Cobo et al., 2002). However, when the number of fuel channel assemblies is very large, as in BWR, the pressure drop is practically constant and the oscillations are negligible (Hennig, 2001). The second boundary condition for out-of-phase oscillations is that the mass flow rate coming from the downcomer, and entering into the reactor core, is practically constant but it oscillates at the inlet of each half of reactor core. In this paper, we have improved the DWOS model, adding one additional node to each channel, this extra node contains the subcooling boiling model. Therefore, the new DWOS model is a two channel thermal–hydraulic model of the reactor core coupled to high order modal kinetics, with three nodes at each channel. These nodes are: the liquid node, where all the heat is invested in the heating of the liquid, the subcooling boiling node, where a fraction of the heat is invested in steam production, and the bulk-boiling region where all the heat is invested in steam production. The coupling, between the neutronic and the thermal–hydraulic parts of the model, is performed through the modal feedback reactivities and the power distribution factors (Mu˜noz-Cobo et al., 2000). The new DWOS model contains four additional differential equations in comparison with the old one (Mu˜noz-Cobo et al., 2000), two for the dynamics of the inception boundary lengths Z1,1 and Z1,2 , and two for the dynamics of the average void fraction in the subcooled boiling region of each channel.

1281

The steady-state results of the DWOS code and the SIMULATE code are very similar for the pressure drop along the channels and the exit qualities, as we display in Section 3. Increasing the feedback gain of the first harmonic mode beyond the critical value, the model displays outof-phase oscillations. We have found that the frequency of these oscillations is 0.34 Hz, which is lower than the experimental value found from the signal analysis of 0.47 Hz, and the value found with the LAPUR 5.2 code of 0.41 Hz, this is due to use a homogeneous model in this paper (Aguirre et al., 2005). Also, it is interesting to point out that the mass flux oscillations at the bulkboiling boundary are delayed with respect to the mass flux oscillations at the channel inlet. Also, the mass flux oscillations at the channel exit are delayed with respect to the mass flux oscillations at the bulk-boiling boundary and with respect to the oscillations at the channel inlet. The delay between the mass flux oscillation at the channel exit and the channel inlet is half a period or 180◦ , this fact is typical of density wave oscillations, and suggests that the density wave mechanism is also at the root of out-of-phase oscillations. The oscillations of the power transferred to the channels are out-of-phase. Nevertheless, when we have out-of-phase instabilities, the normalized oscillating component n1 (t) displays a behavior practically symmetric with time, and therefore as a consequence the power at each half of the reactor core also shows this behavior. This fact has been confirmed experimentally and numerically with 3D calculations performed with the codes RAMONA, MODKIN, and LAMBDA REACT. Therefore, some of the characteristics of the out-of-phase oscillations can be studied with the DWOS code. Acknowledgments The authors of this paper are indebted with the members of the European project NACUSP, and with the MCYT by their support under the Contract BMF20012690.

References Aguirre, C., Caruge, D., Castrillo, F., Dominicus, G., Geutjes, A.J., Saldo, V., van der Hagen, T.H.J.J., Hennig, D., Huggenberger, M.,

1282

J.L. Mu˜noz-Cobo et al. / Nuclear Engineering and Design 235 (2005) 1267–1282

Ketelaar, K.C.J., Manera, A., Munoz-Cobo, J.L., Prasser, H.-M., Rohde, U., Royer, E., Yadigaroglu, G., 2005. Natural circulation and stability performance of BWRs (NACUSP). Nucl. Eng. Des. 235, 401–409. Escriv´a, A., Mu˜noz-Cobo, J.L., Ballester, D., Sancho, J., Chiva, S., July 2003. UPV Results Cofrentes NPP. Report EVOLNACUSP-D016f, Project NACUSP EC FIKS-CT-2000-00041. Ginestar, D., Verd´u, G., March-Leuba, J., 1999. Thermal–hydraulics oscillations and numerical integration. In: Proceedings of the 9th International Topical Meeting on Nuclear Reactor Thermal–Hydraulics (NURETH-9), American Nuclear Society, IL, USA. Ginestar, D., Mir´o, R., Verd´u, G., Henning, D., 2002. A transient modal analysis of a BWR instability event. J. Nucl. Sci. Technol. 39 (4), 554–563. Hashimoto, K., 1993. Linear model analysis of out-of-phase instability in boiling water reactors. Ann. Nucl. Energy 20 (12), 789–797. Hennig, D., 2001. Personal communication. Herbert, A., 1987. Development of the nodal collocation method for solving the neutron diffusion equation. Ann. Nucl. Energy 14 (10), 527–541. Jones, A.B., Dight, D.G., 1963. Hydrodynamic Stability of a Boiling Channel. KAPL-2290. Lahey, R.T., Moody, F.J., 1993. The Thermal–Hydraulics of Boiling Water Reactors, second ed. American Nuclear Society, Illinois. Lee, J.C., Onyemaechi, A., 1989. Phase plane analysis of nuclear coupled density wave oscillations. In: Noise and Non Linear Phenomena in Nuclear Systems, Series B, vol. 192. Editorial Plenum, New York. March-Leuba, J., Rey, J.M., 1993. Coupled thermal–hydraulic neutronic instabilities in boiling water reactors: a review of the state of the art. Nucl. Eng. Des. 145, 97–111. Mir´o, R., 2002. Modal methods to study instabilities in boiling water reactors. Ph.D. Thesis. Polytechnic University of Valencia, Spain.

Mu˜noz-Cobo, J.L., Verd´u, G., Pereira, C., Escriv´a, A., R´odenas, G., Castrillo, F., Serra, J., 1994. Consistent generation and functionalization of one dimensional cross sections for TRAC-BF1. Nucl. Technol. 107, 125–137. Mu˜noz-Cobo, J.L., Rosell´o, O., Mir´o, R., Escriv´a, A., Ginestar, D., Verd´u, G., 2000. Coupling of density wave oscillations in parallel channels with high order modal kinetics: application to out-ofphase oscillations. Ann. Nucl. Energy 27, 1345–1371. Mu˜noz-Cobo, J.L., Podowski, M.Z., Chiva, S., 2002. Parallel channel instabilities in boiling water reactor systems: boundary conditions for out-of-phase oscillations. Ann. Nucl. Energy 29, 1891–1917. Otaduy, P.J., 1979. Modeling of the dynamic behaviour of large boiling water reactor cores. Ph.D. Thesis. University of Florida. Rouhani, S.Z., Axelsson, E., 1970. Calculation of void volume fraction in sub-cooled and quality boiling regions. Int. J. Heat Mass Transfer 13, 383–393. Thom, J.R.S., 1964. Prediction of pressure drop during forced circulation of boiling water. Int. J. Heat Mass Transfer 7, 709– 724. Todreas, N.E., Kazimi, M.S., 1990. Nuclear Systems I—Thermal– Hydraulics Fundamentals. Hemisphere Publishing Corporation, New York. Van Bragt, D.D.B., Van der Hagen, T.H.J.J., 1998. Stability of natural circulation boiling water reactors; Part I: description stability model and theoretical analysis in terms of dimensionless groups. Nucl. Technol. 121, 40–51. Van Braght, D.D.B., 1998. Analytical modeling of boiling water reactor dynamics. Ph.D. Thesis. Delft University, Delft University Press, ISBN-90-407-1719-2/CIP. Wulff, W., Cheng, H.S., Diamond, D.J., Khatib-Rahbar, M., 1984. A Description and Assessment of Ramona 3B Mod 0 Cycle 4: A Computer Code With 3D Neutron Kinetics for BWR System Transient. NUREG/CR-3664.