Complex BWR dynamics from the bifurcation theory point of view

Annals of Nuclear Energy 67 (2014) 91–108

Contents lists available at ScienceDirect

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

Complex BWR dynamics from the bifurcation theory point of view Carsten Lange ⇑, Dieter Hennig, Manuel Schulze, Antonio Hurtado Hydrogen and Nuclear Energy, Institute of Power Engineering, Faculty of Mechanical Engineering, Technische Universität Dresden, 01069 Dresden, Germany

a r t i c l e

i n f o

Article history: Received 17 May 2013 Accepted 22 July 2013 Available online 16 August 2013 Keywords: Hopf bifurcation Asymptotic decay ratio In-phase and out-of-phase power oscillations Stable and unstable limit cycles

a b s t r a c t In a recently published (motivation) paper Lange et al., 2013 we analysed the BWR mode oscillation phenomenon from the bifurcation theory point of view. We demonstrated that some stability properties can be explained only in terms of nonlinear dynamics. In the present paper we shall continue the application of elements of the nonlinear stability analysis and analyse stability states which are characterized by the existence of simultaneous Hopf bifurcations (so-called Hopf–Hopf-bifurcations). From the bifurcation theory it is well-known that the solution manifold of the DEs describing the BWR dynamics can be surprisingly extensive in the vicinity of such bifurcation points. In the present paper we analyse the dynamics of a real BWR and search for a Hopf–Hopf point by appropriate parameter variation. The objective of this paper is to investigate and demonstrate the system dynamics which we have to expect if a Hopf– Hopf bifurcation exists. In a future paper we will analyse these different stability states reflecting the solution manifold near a Hopf–Hopf point concerning their technical relevance. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The dynamical behavior of boiling water reactors (BWRs) which is described by a system of coupled nonlinear partial differential equations is governed by strong feedback mechanisms. From the nonlinear dynamics it is well known that dynamical systems can show a very complex behavior under specific parameter configurations which is reflected in its solution manifold (Neyfeh and Balachandran, 1995; Guckenheimer and Holmes, 1984; Seydel, 2010; Strogatz, 1994; Kuznetsov, 1998). In particular, such systems can generate a multitude of single and multi-stability states which cannot be explained by linear stability analysis methods and is difficult to interpret merely by nonlinear system code analyses results. In order to understand these nonlinear phenomena in detail, the solution manifold of the equation set describing the BWR stability behavior must be examined. In a previous paper (Lange et al., 2013) we discussed an especial phenomenon of the BWR stability behavior, the existence of socalled in-phase and out-of-phase power oscillations (March-Leuba and Blakeman, 1991; Miro et al., 2000; Zhou and Uddin, 2002) from the nonlinear analysis point of view. According to the interpretation of March-Leuba and Blakeman (March-Leuba and Blakeman, 1991) in the first case the fundamental neutron flux (spatial) mode (Fig. 1, left picture) is excited to oscillate in a stable or unstable manner and, in the second case, the first and second azimuthal spatial modes (Fig. 1, right picture) are excited (there are two pairs of azimuthal modes with a phase shift of p/2). The ⇑ Corresponding author. Tel.: +49 351 463 32083; fax: +49 351 463 37161. E-mail address: [email protected] (C. Lange). 0306-4549/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2013.07.034

mode excitation picture (Henry, 1975) bases on the spatial mode expansion method applied in many fields of the physics. A physically reasonable function can be expanded in a series of a complete set of orthogonal functions (Henry, 1975; Weinberg and Schweinler, 1948; Weinberg and Wigner, 1958; Verdù et al., 1994; Miro et al., 2002). Hence the time- and space-dependent neutron flux or the reactor power can be expanded as

Uð~ r; tÞ ¼

X nk ðtÞWk ð~ rÞ;

k ¼ 0; 1; 2; . . .

ð1Þ

k

rÞ the spatial modes, e.g. the k-modes (Verdù et al., 1994; with Wk ð~ Miro et al., 2002), and nk ðtÞ the amplitude functions of the mode k. Because of the subcriticality of the higher modes (in a critical reactor) the higher modes are decaying in time. But in BWRs we encounter the phenomenon that by a mode coupling mechanism the fundamental mode (k = 0) and the first/s azimuthal modes (k = 1, 2) could become linear unstable and oscillate with constant or increasing amplitudes (Fig. 3). The question that arises is what means a stability boundary in the framework of the mode expansion picture? The system stability boundary is given by the first (linear) unstable flux mode (in an in-phase state the fundamental mode and in an out-of phase state the first and/or second azimuthal modes constitute the dominant limit cycle oscillation pattern). However, it should be noted that because of the solution structure of systems of linear differential equations all modes are excited if the stability boundary is reached (Zhou and Uddin, 2002; Dokhane et al., 2004). The mode excitation explanation is given for the first time by the authors in March-Leuba and Blakeman (1991) but they treat e.g. an unstable in-phase state as a ‘‘pure’’ fundamental mode oscillation and overlooked

92

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

Fig. 1. In-phase and out-of-phase modes.

Fig. 2. Simplest visualization of the appearance of a Hopf bifurcation (see Section 3 for mathematical explanation). If the nonlinear dynamical system under control parameter variation encounters a Hopf bifurcation point (reaches the stability boundary in part b)), an periodic orbit (limit cycle) is born and the bifurcation analysis must be used to estimate the (nonlinear) limit cycle stability (the stability indicator is given by b2 Kuznetsov, 1998; Hassard et al., 1981, see Section 3 of this paper).

that due to the different time-constants and amplitudes of the excited mode oscillations we see, e.g. in an in-phase state, first of all, an increasing fundamental mode oscillation (higher mode oscillations are decaying) but the azimuthal modes are also growing up after sufficient time (see the example in Fig. 3, and expression (4) in Section 1.1) depending on the system’s time constants. In the framework of the paper (Lange et al., 2013) we investigated in detail the mathematical background of the in-phase and out-of-phase power oscillation state. In particular, it was shown that both the in-phase state and the out-of-phase state can consecutively experience Hopf bifurcations under parameter variation which generate co-existing in-phase and out-of-phase limit cycle oscillation states. In the scope of the current paper, in particular, we will study the nonlinear dynamics of the system when both, fundamental mode and the azimuthal modes, experience a Hopf bifurcation simultaneously.

1.1. Stable and unstable mode oscillation (spatial mode expansion approach) More in-depth investigations reveal that the eigenvalues and the magnitude of the components of the eigenvectors of the system matrix (the Jacobian matrix) play the essential role in both scenarios. If, for a selected operating point (OP) of the dynamical system, the dominant pair (strictly defined below, see (Hassard et al., 1981) of the complex conjugated eigenvalues reaches, under parameter variation, the imaginary axis and the Hopf conditions (Kuznetsov, 1998; Hassard et al., 1981) are fulfilled, a stable or unstable limit cycle occurs (Fig. 2b) Neyfeh and Balachandran, 1995; Guckenheimer and Holmes, 1984; Seydel, 2010; Strogatz, 1994; Kuznetsov, 1998. Strictly speaking, the dynamical system moves discontinuously to a multi-stability state (stable limit cycle and unstable fixed points, Fig. 2c or unstable limit cycle and stable fixed points, Fig. 2a). If this pair of eigenvalues has crossed the imaginary axis

93

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

n0(t)-n0,steady

Time evolutions of the fundamental and first azimuthal mode and the channel inlet velocities 0.015

fundamental mode

0.000 -0.015

1

n (t)

0.15

first azimuthal mode

0.00

vinlet, i(t) /v0, inlet

-0.15

v inlet,1 (t)

1.025

v inlet,2 (t)

1.000 0

100

200

300

400

500

600

s

700

800

Time Fig. 3. The OP NPP Ringhals, cycle14, record9, in an out-of-phase mode (2-channel BWR-ROM results) (Lange et al., 2009, 2011, 2012a,b,c). The time evolutions of the fundamental and first azimuthal mode n0 and n1 (different ordinate scales!) show clearly that the system is in a linearly unstable state (stable limit cycle), this means the stability boundary is given by the first azimuthal mode stability, thus the system is in an out-of-phase state (below the channel inlet velocities v1;inlet ðtÞ and v2;inlet ðtÞ are added).

the system will be linearly unstable and all components of the solution will increase with time or a periodic orbit will appear with a constant amplitude oscillation. Thereby the amplitude of the limit cycle can change under parameter variation (Fig. 2c, supercritical bifurcation). However, we repeat: Although we assume that all eigenvalues different from the dominant pair are negative (or have negative real parts), the time evolution of all components of the solution vector are, first of all, decreasing but later increasing (Fig. 3). We have to note that in the framework of the spatial neutron kinetic mode description (mode expansion model) the flux (or power) modes reflect a mathematical picture (the modes are used as a mathematical vehicle) to understand and visualize an experimental observed (stable or unstable) in-phase and out-of-phase oscillatory state. If one of the modes becomes linear unstable then the whole system is linear unstable which means that all components of the solution vector contain an unstable term (Kuznetsov, 1998; Dokhane et al., 2004) and, in case there is at least one pair of complex conjugated eigenvalues of the Jacobi matrix, all spatial modes oscillate with different amplitudes indicating the different mode dominance (Lange et al., 2013). Thereby the dominance of an in-phase or out-of-phase stability state, however, is determined by the eigenvalues with the largest real part (Lange et al., 2013; Zhou and Uddin, 2002; Dokhane et al., 2004) (see also definitions below). The linear stability boundary is given by the fixed points with zero real parts of the complex eigenvalues (we find it by variation of a control parameter (Lange et al., 2013). In accordance with the theorem of Hartman and Grobman (Kuznetsov, 1998) (roughly speaking) the time behavior of the nonlinear system in the close neighborhood of a hyperbolic fixed point can be described by the linearized system (topological equivalence of the nonlinear and the linearized system (Kuznetsov, 1998). In (Lange et al., 2013) we demonstrated the spontaneous change in the stability properties and dominance of in-phase and out-of-phase power oscillations under continuous parameter variation. To this end, we introduced some important definitions which will be summarized as follows. The BWR stability can be described by the autonomous dynamical system

!

@ X ðtÞ ! ! ¼ F ðX ; ~ cÞ: @t !

ð2Þ

!

!

!

X ðtÞ with X 2 Rn is a state vector, F (with F : Rn  Rm ! Rn is C1) is a vector field, ~ c 2 Rm is a parameter vector (also called control parameter vector!with m components) and t 2 R is the time. The general solution X ðtÞ of the linearized system around the steady ! state solution X 0 which is satisfying equation !

! ! dX 0 ¼ 0 ¼ F ðX 0 ; ~ cÞ dt

ð3Þ

can be written as !

!

X ðtÞ ¼ X 0 þ

n X ci~ pi eki t i¼1

0

1 0 1 0 1 x1 ðtÞ x01 ðtÞ pi1 n B. C B . C X B . C kt i B. C¼B . Cþ B ci @ .. C @. A @ . A Ae ; i¼1 xn ðtÞ x0n ðtÞ pin

ð4Þ

pi are the eigenvectors of the Jacobian matrix J with where ~

0 @F

1

@x1

B B J ¼ B ... @

@F n @x1

 ..

.



@F 1 @xn

1

C .. C ; . C A

ð5Þ

@F n @xn

ki are the corresponding eigenvalues and ci are constants which can be calculated from the initial conditions. All vector components of ! the solution vector X ðtÞ are formed by a sum of terms containing all eigenvalues (superposition principle, see Eq. (4)). We define (Lange et al., 2013; Zhou and Uddin, 2002; Dokhane et al., 2004):  Eigenvalue or pair of complex conjugated eigenvalues (k1 ) of the Jacobian matrix of the system (2) with the largest real part: Reðk1 Þ > Reðki Þ 8ki with i–1.  Eigenvalue or pair of complex conjugated eigenvalues (k2 ) with the second largest real part: Reðk1 Þ > Reðk2 Þ > Reðki Þ 8ki with i–1; 2.

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

pin and ~ pout of the linearized BWR system  The eigenvectors ~ (ROM) correspond to specific reactor eigenstates which are referred to as in-phase state and out-of-phase state, respectively. The eigenvalues corresponding to ~ pin and ~ pout are denoted by kin and kout . In case of a two dimensional system, where the state variables are the amplitude functions (n0 and n1 ) of the fundamental spatial mode W0 and first azimuthal spatial mode W1 , the general solution reads. !

!

X ðtÞ ¼ X 0 þ cin~ pin ekin t þ cout~ pout ekout t

ð6Þ

or in components (n00(t)and n01 ðtÞ are the steady states)



n0 ðtÞ



n1 ðtÞ

 ¼

n00 ðtÞ



n01 ðtÞ

þ cin

pPin0 pPin1

!

ekin t þ cout

0 pPout P1 pout

!

Signals of LPRM 9 and LPRM 31 of KKLc7_rec4 70 66 64 62 60 58 56 54 52 50

ekout t ;

48

ð7Þ

LPRM 9 LPRM 31 NF = 0.58 Hz

68

Relative Amplitude

94

46

70 65 60 55 50 730

732

400

 The system is linear unstable in the in-phase state, if Reðkin Þ > 0; 8 ki –kin ; Reðki Þ < 0. The system is linear unstable in the out-of-phase state, if Reðkout Þ > 0; 8 ki –kout ; Reðki Þ < 0.  If Reðkin Þ > Reðki Þ, 8 ki –kin , the asymptotic solution (see Eqs. (4) and (6)) of the BWR can be approximated by !

!

pin ekin t ; X ðtÞ ¼ X 0 þ cin~

ð8Þ

this means the fundamental power mode oscillation is dominant.  If Reðkout Þ > Reðki Þ, 8 ki –kout , the asymptotic solution (see Eqs. (4) and (6)) of the BWR can be approximated by !

!

X ðtÞ ¼ X 0 þ cout~ pout ekout t ;

ð9Þ

this means the first/s azimuthal power mode oscillation is dominant.  The dominance of the eigenvector ~ pin is reflected in the following properties: (1) The element corresponding to the fundamental mode pPin0 is larger (at least one order of magnitude) than the element corresponding to the first azimuthal mode pPin1 and (2) the elements corresponding to the thermal–hydraulic and heat conduction variables in the first core1 half have the same sign and the same absolute values as the corresponding elements of the second core half (this was discussed in Zhou and Uddin (2002), Dokhane et al. (2004).  The dominance of the eigenvector ~ pout is reflected in the following properties: (1) The element corresponding to the first azi1 muthal mode pPout is larger (at least one order of magnitude) 0 than the element corresponding to the fundamental mode pPout and (2) the elements corresponding to the thermal–hydraulic and heat conduction variables in the first core half have the opposite sign but have the same absolute values as the corresponding elements of the second core half (this was also discussed in Zhou and Uddin (2002), Dokhane et al. (2004)). These definitions allow a clear distinction of an in-phase oscillation state and an out-of-phase oscillation state. Due to the solution structure (4) and (7), where each solution component depends on all eigenvalues of the Jacobian matrix, if one of the eigenvalues has a real part larger than zero, all solution components will asymptotically grow with time. This means, for instance, if the BWR encountered a supercritical Hopf bifurcation at ~ cc where a stable out-of-phase power oscillation is generated, the solution (Hassard et al., 1981)

h i2pt i ! pout þ 0ðe2 Þ X ðt; ~ cÞ ¼ X 0 ð~ cc Þ þ eRe eTðeÞ~ !

1

ð10Þ

The symmetry line of the first azimuthal mode subdivides the core into two core halves. In our ROM, the thermal–hydraulic behavior of these two core halves is described by two representative flow channels which include also a fuel heat conduction model, respectively (Lange et al., 2009, 2011).

734

500

736

738

600

740

700

s

800

Time Fig. 4. Time evolution of the measured LPRM signals (LPRM 9 and 31) KKL Cycle7, stability test, record #4 (23:39).

predicts also an oscillation of all solution components in the neighborhood of ~ cc . In particular, the fundamental mode n0 ðtÞ oscillates even though the BWR is in an out-of-phase power oscillation state (Lange et al., 2013; Zhou and Uddin, 2002; Dokhane et al., 2004). This is shown in Fig. 3. In solution (10) ~ cc is the critical parameter at which the supercritical Hopf bifurcation appears, e is the small amplitude and TðeÞ the period of the oscillation (Hassard et al., 1981). 2. ROM analyses of the KKL, cycle7, stability test During the stability test KKLcycle7 an out-of-phase limit cycle was observed and comprehensively analysed (Lange et al., 2013; Aguirre and Hennig, 2002; Blomstrand et al., 1992; Wiktor et al., 1997). Hence, these stability test results are well suitable to be used for extended nonlinear stability analysis. From the previous analyses we know that at least one Hopf bifurcation point exists (Lange et al., 2013). In the framework of the current investigations we are searching for the existence of Hopf–Hopf bifurcation points. From the nonlinear literature we know that the system dynamics in an environment of a Hopf–Hopf point can be very complex. Hence it is worth from the author’s point of view to look for the physical conditions of the existence of double Hopf dynamics in real dynamical systems. To this end, we used the reduced order model (Dokhane et al., 2004; Lange et al., 2009) and selected suitable system parameters for parameter variations. In the thesis of Lange (Lange et al., 2009) the ROM version developed at PSI (Dokhane et al., 2003a,b) was qualified for BWR stability analysis among others by the model extension to an external loop (Appendix A). Because of the simplicity of the loop model many parametric investigations were necessary to demonstrate the equivalence of this model to the real loop dynamics. In this framework flow cross section variations were performed which prove also successful for model-specific options of extended nonlinear studies. Hence, we are searching for double Hopf states on the basis of external loop parameter variations. The investigations are embedded into the KKLc7, record4, stability analysis. The ROM is published in Lange et al. (2013, 2009, 2011, 2012a), Dokhane et al. (2004), Dokhane et al. (2003a,b) and the recirculation loop model is summarized in Appendix A. In the next section the KKL stability test c7 is briefly presented. 2.1. KKL stability test The stability test was performed when the plant was still in its power ascension phase (Aguirre and Hennig, 2002; Blomstrand

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

95

Fig. 5. Numerical integration of the ROM equations for Aol ¼ 12:68.

Power Flow Map for Aol = 12.68

The test was initiated by pulling control rods. The onset of undamped regional power oscillations occurred at approximately 57% power. Further increase of power to 59% caused the oscillation to grow to approximately up to 30% peak-to-average on several LPRMs. Notice, a continued power increase using control rod positions in this case (Aguirre and Hennig, 2002; Blomstrand et al., 1992; Wiktor et al., 1997) made the radial power distribution bowl-shaped with low-powered central region. This is by experience a typical power or neutron flux pattern to have a disposition to regional oscillation if the system becomes linear unstable. Since the oscillations were growing up it was decided to suppress them by reducing the reactor power (see example in Fig. 4).

80

Thermal Power

70

60

l in

e

u ar

n

bl e sta lin

ea

reg r

ion

bl e sta

reg

ion

50

KKLc7 record4 OP: Thermal Power (59.5%) Core Flow (36.5%)

40

SB (subcritical PAH-B) 112% Rod Line exclusion region

30 0.2

0.3

0.4

0.5

2.2. ROM analyses: Influence of the recirculation loop (see Appendix A: Recirculation loop model (Lange et al., 2013, 2009, 2011) 0.6

Core Flow Fig. 6. Stability boundary (SB) transformed into the power flow map.

et al., 1992; Wiktor et al., 1997). The goal was to gather stability data for Cycle7 and compare them with the stability performance demonstrated during initial core stability tests (Cycle1). The KKL core in Cycle7 contained only 8  8 fuel assembly types supplied by GE (General Electric) apart from four (4) SVEA-64 fuel assemblies supplied by ABB Atom and four (4) KWU-9  9 fuel assemblies from Siemens. On September 6th 1990 (see Aguirre and Hennig, 2002), during the reactor start-up stage, after nine hours holding thermal power at 70% of rated power,2 the core flow was reduced to about 37% core rated flow.3 Initially the flow control valves (FCVs) were closed to their minimum position and the recirculation pumps were switched to the lower rotational speed. The desired core flow was then managed by opening the FCVs back to their maximum position.

2 3

100% Power = 3138 MW. 100% Flow = 11151 kg/s

In the following investigations the NPP Leibstadt cycle 7 record #4 (stability test, KKLc7_rec4 (Aguirre and Hennig, 2002; Blomstrand et al., 1992; Wiktor et al., 1997) data were analyzed. To this end, semi-analytical bifurcation analyss analyses with the ROM were performed for different values of Aol ¼ Ainlet =Adoc which is a very significant external loop design parameter (see Appendix A, channel flow cross section/downcomer flow cross section). The selection of the reference OP and the procedure to calculate the ROM input is not presented here (see Lange et al., 2009, 2011). In the scope of the bifurcation analysis, N sub and DP ext are defined to be the iteration and bifurcation parameters, respectively, and the downcomer friction is neglected. In order to study the effect of the recirculation loop on the BWR stability behavior, the ratio Aol was varied in a physically and technically reasonable domain. Thereby, a change of Aol ¼ Ainlet =Adoc corresponds to a change of the downcomer flow cross section Adoc . For the calculation of the ratio Aol from the RAMONA5 model, the flow cross section of downcomer 2 (DC2, see Studsvik/Sc, 2003; Wulff et al., 1984) was taken into account. The result is Aol ¼ 12:68 and is considered to be the reference value for this analysis. At first numerical integration was carried out at the reference OP for Aol ¼ 12:68. The transient was initiated by introducing per-

96

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

Fig. 7. Numerical integration of the ROM equations for Aol ¼ 8:00.

Fig. 8. Numerical integration of the ROM equations for Aol ¼ 7:00.

turbations in the channel inlet velocities with opposite sign (an inphase oscillation is triggered). The result is shown in Fig. 5. It can be seen that the fundamental mode decays and the first azimuthal mode oscillates. While the nonlinear stability of this OP cannot be deduced from the picture it can be said that the fundamental mode is (linearly) stable and the first azimuthal mode is not, i.e. it is either divergent or reaches a stable limit cycle. Fig. 6 shows the additionally computed stability boundary (SB) for the first azimuthal mode in the power flow map where the 112% rod-line and the exclusion region for cycle 7 are included.

The analysis has shown that the pair of complex conjugated eigenvalues with the largest real part corresponds to the out-ofphase oscillation state (k1 ¼ kout with Reðk1 Þ > Reðki Þ 8 ki ; i–1). For this parameter configuration, all the other eigenvalues have strictly negative real parts. Then the ratio Aol was changed to Aol ¼ 8:00 and numerical integration has been carried out again. The analysis result is shown in Fig. 7. Again, at the reference OP only the out-of-phase mode is excited, because Reðkout Þ > 0 and Reðkin Þ < 0 with k2 ¼ kin , but the fundamental mode has become less stable.

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

97

Fig. 9. Numerical integration of the ROM equations for Aol ¼ 6:75.

Fig. 10. Numerical integration of the ROM equations for Aol ¼ 6:0.

Furthermore, numerical integration for the ratios Aol ¼ 7:00, Aol ¼ 6:75 and Aol ¼ 6:00 have been performed. The time series for these cases are presented in Figs. 8–10, respectively. Now we can see that the fundamental mode becomes (linearly) unstable as well. This means that an additional stability boundary for the fundamental mode was crossed during the variation of Aol . This is confirmed by bifurcation analysis for Aol ¼ 6:00 depicted in Fig. 11. The regions with different configurations of the real parts of kout and kin are labeled. For Aol ¼ 8:00, for instance, the reference OP is not located in the region where the in-phase state

is linearly unstable. For Aol ¼ 12:68 a region where Reðkin Þ > 0 does not exist. Aol-variation causes a change of the stability characteristics of the in-phase state while the out-of-phase state is not affected. The larger Aol the smaller the region where Reðkin Þ > 0. The analysis has shown that the downcomer flow cross section variation does not have an effect on the out-of-phase oscillation state because the Aol variation does not change the location of operational points where Reðkout Þ ¼ 0 in the N sub  DP ext -parameter map. On the other hand, the Aol (Adoc ) variation has a significant ef-

98

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

Fig. 14. Bifurcation diagram in the (DPext , N sub ) plane for the ROM (Aol ¼ 6:00).

Fig. 11. Stability boundary for Aol ¼ 6:00.

fect on the in-phase oscillation state. If Aol is decreased (Adoc increased) the area of the region for which the in-phase state is excited will grow. I.e. if the downcomer flow cross section is increased the region where Reðkin Þ P 0 (in the N sub  DP ext -parameter space) is growing. The stability characteristics of the generated limit cycles are shown in form of the first Lyapunov coefficients (Neyfeh and Balachandran, 1995; Kuznetsov, 1998)4 for the first azimuthal mode and fundamental mode in Figs. 12 and 13 respectively. The starting point is the left intersection (HHB1) between the Hopf bifurcation curves of the both modes. It can be seen that stable limit cycles in the linear unstable region are always created for the first azimuthal mode and also for the fundamental mode with the exception of a small region close to the intersections between the two modes’ stability boundaries where two so-called generalized Hopf bifurcations (Lange et al., 2013; Kuznetsov, 1998) occur (see discussion in Section 3). The two intersections labeled by HHB1 and HHB2 in Figs. 11 and 14 will be the subject of the following discussion. 3. ROM dynamics near Hopf–Hopf bifurcation points (for a short textbook introduction of this phenomenon see Appendix B))

Fig. 12. First Lyapunov coefficient along the Hopf bifurcation curve of the first azimuthal mode.

Fig. 13. First Lyapunov coefficient along the Hopf bifurcation curve of the fundamental mode.

In the following we demonstrate the surprising solution manifold in the vicinity of double Hopf bifurcations. It should motivate to search for parameter regions where these kinds of bifurcations exist and find them by system codes too. Since each stability boundary in Fig. 11 represents a Hopf bifurcation curve, an intersection between two of them means that a Hopf–Hopf bifurcation occurs (in Appendix B we demonstrate that in this case the center manifold reduction of the ROM equation set results in a 4-dimensional center manifold and, accordingly, we find a 4-dimensional normal form). In order to study the dynamics of the ROM in the vicinity of a Hopf–Hopf bifurcation the ROM was coupled with the bifurcation analysis code MatCont (Dhooge et al., 2007, 2006) (for a short summary see Appendix C). MatCont has been selected for this task because it is, to the author’s knowledge, the only available code that can compute Hopf–Hopf bifurcation normal form coefficients. Because the ROM is implemented in MATLAB, MatCont’s implementation in the same language is also advantageous. Using a ROM version coupled with MatCont the two stability boundaries were computed for Aol ¼ 6:00. When selecting a sufficiently small continuation step size MatCont automatically recog4 The first Lyapunov coefficient is directly connected with the second Floquet exponent b2 (Guckenheimer and Holmes, 1984; Seydel, 2010), this means if the first Ljapunov coefficient is zero, b2 = 0 too and a generalized Hopf bifurcation occurs (see Lange et al., 2013 for details).

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

99

Fig. 15. Region (ii), torus.

Fig. 16. Region (ii), fixed point.

Fig. 17. Region (ii), torus.

nizes the occurrence of the Hopf–Hopf bifurcations at the two points of intersection and prints the required normal form coefficients. These can be used for further studies in an analytical representation of the Hopf–Hopf bifurcation’s normal form. We can

conclude from the normal form that a Hopf–Hopf bifurcation induces additional bifurcations (Neyfeh and Balachandran, 1995; Guckenheimer and Holmes, 1984; Seydel, 2010; Strogatz, 1994; Kuznetsov, 1998) whose curves can also be computed by MatCont.

100

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

Fig. 18. Region (ii), fixed point.

Fig. 19. Region (iii).

Fig. 20. Region (iv).

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

Fig. 21. Region (v).

Fig. 22. Region (ix).

Fig. 23. Region (ix) in the (n0 , n1 ) plane.

101

102

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

Fig. 24. Hopf–Hopf-bifurcation close to the origin (linear approximation).

Fig. 25. Hopf–Hopf-bifurcation, numerical solution.

Namely it generates two Neimark-Sacker bifurcation curves T 1;2 (Appendix B). The Hopf–Hopf bifurcation may also create bifurcation curves C and J whose meaning is discussed in Appendix B. Contrary to the existence of the curves T 1;2 the existence of C and J is not guaranteed by the normal form. Although MatCont is not able to compute their route in the bifurcation diagram, their existence was proven indirectly as explained below. The full bifurcation diagram for the ROM (Aol ¼ 6:00) is depicted in Fig. 14. For the clear visualization of these bifurcation curves it is convenient to compute phase space diagrams at parameter configurations in particular regions separated by the bifurcation curves. These can then be compared with results gained from the normal form. Selected (n0 , n1 , V inlet;1 ) phase space diagrams depicted in Figs. 15–23 will be discussed in the following. The labeling of the regions is consistent with the labeling used in Fig. 14 (see Figs. 24 and 25). Region (ii): Starting sufficiently near the Hopf–Hopf bifurcation point in this region leads to the system evolving onto a densely filled torus with incommensurable frequencies. As one departs further from the bifurcation point the torus gets deformed and eventually destroyed and the system evolves to a fixed point. This means that tori exist even in the linear stable region, although only near Hopf–Hopf bifurcations. Regions (iii–v): The system is constrained to a torus as predicted by the normal form (see Appendix B). The tori’s amplitudes are small for the selected parameters. Region (ix): In this region, which does not appear in the normal form, the system moves on a strange attractor with relatively large amplitudes that is not described by the normal form in Appendix B. It should be noted that this attractor covers all state variables and is not restricted to neutron kinetic or thermo-hydraulic state variables alone. The existence of this attractor can be explained as follows. Suppose curve C exists and has been passed in region (ix) (see Appendix B). This means an unstable three-dimensional torus has been created at C. According to the Ruelle–Takens–Newhouse theorem (Guckenheimer and Holmes, 1984; Kuznetsov, 1998) the

torus is unstable under certain parameter perturbations and decays to a strange attractor. This explains the existence of the strange attractor but it also means that the existence of C can be reasonably assumed. This shows the existence of J as well, since the three-dimensional torus created at C has to disappear (at J). Additionally the attractor’s Lyapunov dimension has been calculated with the AnT 4.669 software package (AnT 4.669, 2003). The Lyapunov dimension is 4.05 which means the attractor’s fractality is not very pronounced. To prove the attractor’s strangeness the Lyapunov exponents have been calculated (also with AnT 4.669). There’s one positive exponent (0.041) which proves the strangeness of the attractor. Because its close to zero the exponential departure on the torus is not very pronounced as well which means that trajectories starting close to each other remain so for a relatively long time. However, the negative exponents span six orders of magnitude which shows that the attractor attracts weakly in some and strongly in other directions. At the Neimark–Sacker bifurcation curves two-dimensional tori are created or annihilated. There are points along these curves where tori with rational winding numbers are created. By departing from the Neimark–Sacker curves through parameter variation regions are found that still exhibit this so called frequency-locking (Arnold tongues). This explains tori with commensurable frequencies in regions (iii and v). March-Leuba already reported on chaotic dynamics in BWRs in March-Leuba et al. (1986), March-Leu (1984). However, he showed a route to chaos via period doubling and the chaotic states found lie deep inside the linear unstable region. To the authors’ knowledge there have not been any reports yet regarding a Ruelle–Takens–Newhouse scenario, describing chaotic states close to the stability boundary, induced by a Hopf–Hopf bifurcation. 4. Summary and conclusions In the present work the BWR dynamics was examined in the environment of a so-called double Hopf bifurcation point. The Hopf–Hopf points were searched with the aid of the bifurcation code MatCont by variation of an external loop design parameter. We found a surprising complex dynamics represented by the solution manifold of the BWR ROM equations appearing in the vicinity of double Hopf points. The objective of this paper is to demonstrate this solution manifold and motivate an in-depth analysis of the BWR dynamics in the sense of the RAM–ROM method (Lange et al., 2013, 2009, 2011). Acknowledgements We want to thank VGB Power Tech Service GmbH for supporting the Chair of Hydrogen and Nuclear Energy of the Technische Universität Dresden. Appendix A First a short description of the recirculation loop model, which is a sub-model of the ROM, is given. A complete description of the ROM is given in Lange et al. (2009). In the notation used here an asterisk on a variable or parameter indicates the original dimensional quantity while any quantity without an asterisk is dimensionless. The mass balance of the downcomer can be written as

@  _ ðtÞ ¼ 0 m @z tot _ tot ðtÞ and the total mass flow m

ð11Þ

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

" # X X      _ tot ðtÞ ¼ _ n ðtÞ ¼ qf An;inlet vn;inlet ðtÞ m m n

ð12Þ

n

can accordingly be expressed by the sum of the core channel mass _ n ðtÞ, because the coolant in all hydraulic regions is considflows m ered to be incompressible. Here, n is the channel number and An;inlet is the flow cross section of the nth heated channel. The energy balance of the downcomer is reduced to a boundary   condition hinlet ¼ hdoc inlet ¼ const: because energy gain and loss are neglected. The momentum balance of the downcomer can be written as



      _ tot ðtÞ _ tot ðtÞ 2 @pdoc m @ m f1U ; þ ¼ þ g  qf ;     @t Adoc 2qf Ddoc Adoc @z

ð13Þ

where the term on the left hand side describes the pressure drop in the downcomer, the first term on the right hand side describes the pressure drop due to inertial effects of the coolant, the second term states the downcomer friction and the last term is the gravity term. Adoc is the downcomer flow cross section, Ddoc is the hydraulic diam_ tot is the total mass flow, qf is the liquid eter of the downcomer, m density and f1U is the single phase friction factor. Substituting and transforming into dimensionless form leads to



! " #2 X@ X @pdoc ¼ Aol vn;inlet ðtÞ  Nf 1U A2ol Dol vn;inlet ðtÞ @t @z n n þ Fr 1 ;

ð14Þ

wherein Aol and Dol are defined as

 Aol ¼ Ainlet Adoc Dh



ð15Þ

Ddoc

(Dh

and Dol ¼ hydraulic diameter of the heated channel). The ODEs for the channel inlet velocities vn;inlet ðtÞ are determined by

I

@p dz ¼ 0 ¼ @z

Z

1

0

@ph ch dz þ @z

Z

0

1

@precirc dz @z

ð16Þ

in which the pressure drop over the recirculation loop is given by

Z

0

1

@precirc dz ¼ @z

Z 1

o

@pdoc dz  DPhead : @z

ð17Þ

The evaluation of Eq. (16) with expression (17) was performed by using the symbolic toolbox of MATLAB. The final ODE for the nth heated channel can be written as

Xd d vn;inlet ðtÞ ¼ An ðtÞ  Bn ðtÞ  vn;inlet ðtÞ dt dt n  Bn ðtÞNf 1U Aol Dol

!

X vn;inlet ðtÞ

103

because the contribution of the downcomer friction is very small in (18). It can be shown that the effect of the downcomer friction on the thermal–hydraulic stability behavior is very small. Hence, it can be neglected for BWR stability analyses. Eq. (20) relates the steady state external pressure drop to the pump head. In addition to that it can be seen that the steady state core inlet velocity vn;inlet;0 does not depend on the downcomer flow cross section. This means, vn;inlet;0 depends on DPhead or DP ext only. As expected, each ODE for vn;inlet ðtÞ is hydraulically coupled because of (12) with all the other heated channels. Eq. (18) is the ODE for the channel inlet velocity vn;inlet ðtÞ of the nth heated channel. The ODEs for vn;inlet ðtÞ can be written separately for one and the two heated channel(s). The result for one heated channel is

d 1 vinlet ¼ ½A þ BNf 1U Aol Dol ðvinlet Þ2  dt 1þB

ð21Þ

and the ODEs for vn;inlet ðtÞ for two heated channels (n ¼ 1; 2) can be written as

d A2  B1 þ A1  B2 þ A1 B1 Nf ;1U Aol Dol v1;inlet ¼ þ ½v1;inlet þ v2;inlet 2 ; dt 1 þ B1 þ B2 1 þ B1 þ B2 d A1  B2 þ A2  B1 þ A2 B2 Nf ;1U Aol Dol v2;inlet ¼ þ ½v1;inlet þ v2;inlet 2 : dt 1 þ B1 þ B2 1 þ B1 þ B2 ð22Þ According to Eq. (18) (second and third term on the right hand side) v_ n;inlet ðtÞ depends on all heated channels. The inertial term contributes to the mass flow changes of all heated channels. Thus, the inertial term of the downcomer momentum balance describes the impact of all heated channels on the nth heated channel. From the physical point of view, if the downcomer flow cross section is increased (Aol decreases), the inertial effects of the downcomer mass flow decrease. For Adoc ! 1 the ratio Aol ! 0 which corresponds to the constant external pressure drop boundary condition. In this case (Aol ¼ 0), the inertial term in the downcomer momentum balance vanishes and the nth heated channel is independent of all other heated channels. This means the change of the mass flow in the kth heated channel does not affect the nth heated channel. Consequently, if the ratio Aol is zero (Aol ¼ 0) the inertial effects of the downcomer vanish and (18) is reduced to v_ n;inlet ¼ An which is the result of

0

Z 0

1

@ph ch dz  DPext ; @z

ð23Þ

evaluated in Dokhane et al. (2004) (final ODEs for the heated channel inlet velocities presented by Dokhane (Dokhane et al., 2004) with DP ext ¼ Fr 1 þ DPhead ).

!2 ;

ð18Þ

n

wherein An ðtÞ is defined as

 dN n;pch ðtÞ dl ðtÞ An ðtÞ ¼ ffn;11 ðtÞ n þ ffn;12 ðtÞ ffn;14 ðtÞ dt dt  þffn;13 ðtÞ þ 1 þ Fr  DPhead FrDPext |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}

Appendix B B.1. Highly complex bifurcation type: Hopf–Hopf bifurcation

1

ð19Þ

and Bn ðtÞ is defined as Bn ðtÞ ¼ Fr Aol =ffn;14 ðtÞ. The time dependent intermediate terms ffn;11 ðtÞ, ffn;12 ðtÞ, ffn;13 ðtÞ and ffn;14 ðtÞ are calculated in Dokhane et al. (2004), Lange et al. (2009). DP ext is the steady state external pressure drop with DPext ¼ Fr1 þ DP head , where downcomer friction is neglected. In the steady state case Eq. (16) can be approximated as

Z 1

0

@precirc dz  Fr 1  DPhead ¼ DP ext ; @z

ð20Þ

From the theory of nonlinear dynamical systems it is wellknown that we can transform Eq. (2) into the so-called normal form. The transformation can be carried out in two different ways: (1) by using a multi-scale perturbation theory and (2) by using near identity transformations. The normal form of the bifurcation type is of central meaning in the framework of bifurcation theory. From the Hopf bifurcation phenomenon we know that we can reduce an n-dimensional (nonlinear) equation system to a 2-dimensional one by center manifold reduction. The Hopf bifurcation emerges if one pair of complex eigenvalues (of the Jacobian matrix) under control parameter variation reaches the imaginary axis (and the Hopf conditions (Guckenheimer and Holmes, 1984) are fulfilled). Now we shall look what happens if more than one pair of

104

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

(

Introducing ri ¼ .2i simplifies the planar system as follows:

r_ 1 r_ 2

 

n1 ¼ p11 r 1 ; complex conjugated eigenvalues reach (and pass) the imaginary axis (called Hopf–Hopf or double Hopf bifurcation). Calculating the Hopf–Hopf-bifurcation normal form coefficients for the ROM via MatCont (Appendix C) allows an analytical representation of the normal form in order to study the dynamical system’s solution in a neighborhood of the bifurcation. It is possible to deduce important features of the system’s solution by looking at the normal form. The following discussion follows Kuznetzov (Kuznetsov, 1998). The Hopf–Hopf-bifurcation’s normal form is a four dimensional dynamical system with two parameters and its representation in polar coordinates ð.1 ; .2 ; u1 ; u2 Þ reads as follows:

¼ .1 ðl1 þ p11 ðlÞ.21 þ p12 ðlÞ.22 þ s1 ðlÞ.42 Þ þU1 ð.1 ; .2 ; u1 ; u2 ; lÞ; ¼ .2 ðl2 þ p21 ðlÞ.21 þ p22 ðlÞ.22 þ s2 ðlÞ.41 Þ þU2 ð.1 ; .2 ; u1 ; u2 ; lÞ;

ð24Þ

s ¼ 2t:

n2 ¼ p22 r2 ;

The signs are chosen such that n1;2 are non-negative variables. Additionally, introduce the following new parameters (l1;2 remain the system’s bifurcation parameters):

h¼

p12 ; p22

d¼

p21 ; p11

H¼

s1 ; p222

D¼

s2 : p211

These new parameters have been computed by MatCont for the left Hopf–Hopf-bifurcation point of the ROM and read as follows

h ¼ 0:6696;

d ¼ 6:0379;

H ¼ 39624;

D ¼ 239961:

The coefficients for the cubic terms H and D have been set to +1 and 1 respectively, retaining their correct signs. The system now reads as follows:

(

¼ x1 ðlÞ þ W1 ð.1 ; .2 ; u1 ; u2 ; lÞ;

n_ 1 n_ 2





l1  n1 þ hn2 þ Hn22 ;  ¼ n2 l2  dn1 þ n2 þ Dn21 : ¼ n1

ð26Þ

Looking at the system’s fixed points

¼ x2 ðlÞ þ W2 ð.1 ; .2 ; u1 ; u2 ; lÞ;

(

wherein the system’s parameters l ¼ ðl1 ; l2 Þ are the real and x1;2 are the imaginary parts of the eigenvalues

k1;4 ¼ l1 ix1 ;

ð25Þ

The criterion p11 p22 70 defines two significantly different classes of bifurcation diagrams for the system. For all following calculations the ROM normal form coefficients of the left Hopf–Hopfbifurcation point as computed by MatCont have been used. This revealed that only the case p11 < 0, p22 > 0, i.e. p11 p22 < 0, is relevant here. Some additional substitutions have to be made to treat the system effectively. Introduce the following non-negative variables and perform a time scaling:

Fig. 26. Region (i), the system exhibits a stable spiral and a limit cycle.

8 ._ 1 > > > > > > > > < ._ 2 > > > > > > _1 >u > : u_ 2



l1 þ p11 r1 þ p12 r2 þ s1 r22 ; ¼ 2r 1 l2 þ p21 r 1 þ p22 r 2 þ s2 r21 : ¼ 2r 1

n_ 1 n_ 2

¼ 0; ¼0

;

reveals a trivial fixed point

k2;3 ¼ l2 ix2 3

of the system. Furthermore U1;2 ¼ Oðð.21 þ .22 Þ Þ and W1;2 ¼ Oð1Þ are 2p – periodic functions in /1;2 and pij , si are certain partial derivatives of the original system’s right-hand side. In the following, the normal form is to be discussed without the higher order terms U1;2 , W1;2 . Removing these terms causes that the first pair of equations (the planar equations) becomes independent of the second pair, or angular equations, which describes two independent rotations with angular velocities x1;2 in the planes .2;1 ¼ 0. This means that the system’s dynamics are essentially given by the planar equations. Removing these terms also means that the normal form will not be locally topologically equivalent anymore to any (generic) system exhibiting a Hopf–Hopf bifurcation. Still, important features can be learned from the normal form. The following relations between the planar equations and the full system hold:  Fixed point .1 ¼ .2 ¼ 0 of the planar equations ? fixed point in the origin of the full system.  Fixed point with either .1 ¼ 0 or .2 ¼ 0 of the planar equations ? limit cycle of the full system with a frequency of either x2 or x1 .  Fixed point with .1;2 > 0 of the planar equations ? two-dimensional torus of the full system with frequencies x1;2 .  Limit cycle of the planar equations ? three-dimensional torus of the full system with frequencies x1;2 and the frequency of the limit cycle.

E0 ¼ ð0; 0Þ;

ð27Þ

as well as fixed points on the coordinate axes:

E1 ¼ ðl1 ; 0Þ;

l1 > 0

ð28Þ

and

E2 ¼ ð0; l2 Þ;

l2 < 0:

ð29Þ

These fixed points cannot exist in the whole parameter plane because n1;2 must not become negative. Instead they are created at the curves

H1 ¼ fðl1 ; l2 Þ : l1 ¼ 0g

ð30Þ

and

H2 ¼ fðl1 ; l2 Þ : l2 ¼ 0g:

ð31Þ

Looking at the relations between the planar and full system above shows that these curves create limit cycles of the full system, i.e. they are two independent Hopf-bifurcation-curves of the full system. While E0;1;2 are fixed points one can find by simply looking at the planar equations there is an additional non-trivial fixed point whose linear approximation is revealed by linearizing the bracket terms of the planar system (26) and finding their root:



0 ¼ l1  n1 þ hn2 ; 0 ¼ l2  dn1 þ n2 :

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

105

Fig. 31. Region (vi), the system exhibits a stable spiral, the unstable fixed point E2 and the saddle E3 .

Fig. 27. Region (ii), the system exhibits a stable spiral.

Fig. 28. Region (iii), the system exhibits a stable spiral and the saddle E1 on n1 . Fig. 32. Region (vii), the system exhibits a stable spiral, the unstable fixed point E2 and an invisible limit cycle.

Fig. 29. Region (iv), the system exhibits a stable spiral, the saddle E1 and the unstable fixed point E2 on n2 . Fig. 33. Region (viii), the system exhibits a stable spiral point, the unstable fixed point E2 and a limit cycle.

relevant case dh < 1, E3 exists for hl2 < l1 < l2 =d. This results in two curves in between which E3 exists:

  T 1 ¼ ðl1 ; l2 Þ : l1 ¼ hl2 þ O l22 ; l2 < 0

ð33Þ

and

  T 2 ¼ ðl1 ; l2 Þ : l2 ¼ dl1 þ O l21 ; l1 > 0 : Fig. 30. Region (v), the system exhibits a stable spiral, the stable fixed point E1 , the unstable fixed point E2 and the saddle E3 .

Solving this system yields the fourth fixed point

  l  hl2 dl  l2 E3 ¼  1 þ Oðklk2 Þ;  1 þ Oðklk2 Þ : dh  1 dh  1

ð32Þ

This fixed point corresponds to a two-dimensional torus of the full system. Again, the coordinates must not be negative. For the

ð34Þ

Along these curves E3 vanishes in E2;1 . In the full system this corresponds to a creation of a two-dimensional torus from a limit cycle and vice versa, i.e., T 1;2 are Neimark–Sacker-bifurcationcurves of the full system. Furthermore, E3 can encounter a Hopf-bifurcation. By calculating the eigenvalues of the Jacobian of the planar system with respect to E3 (in linear approximation) it is shown below that their real parts vanish along the curve

C¼



d1 l ðl1 ; l2 Þ : l2 ¼  l1 þ Oðl21 Þ; hl2 < l1 < 2 : h1 d

ð35Þ

106

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

Fig. 34. Adjacency graph of bifurcation types covered by MatCont (Dhooge et al., 2006). The meaning of all abbreviations is explained in next figure.

The real parts of its eigenvalues are

Rek1 ðE3 Þ ¼ Rek2 ðE3 Þ ¼

Fig. 35. Table of all bifurcation types considered by MatCont (Dhooge et al., 2007, 2006).

Regarding the full system C describes the creation of a threedimensional torus from a two-dimensional one. The Jacobian is

df ¼ dn



l1  2n1 þ hn2

hn1

dn2

l2  dn1 þ 2n2

With respect to E3 it is

 df  ¼ dnn¼E3

l1 hl2 dh1 dðdl1 l2 Þ dh1

1 hl2 Þ  hðldh1 1 l2  dldh1

! :

 :

1 1 ðð1  dÞl1 þ ð1  hÞl2 Þ: 2 dh  1

Either one of the real parts can be set equal to zero (necessary condition for a Hopf bifurcation) and brought in the above form l2 ðl1 Þ for C. In our case C creates a heteroclinic orbit that connects E1 with E0 . This scenario can also be thought of as an invisible unstable limit cycle partly lying on the n2 -coordinate axis. There is an additional curve J at which this limit cycle becomes visible by detaching itself from the n2 -axis, i.e. the heteroclinic orbit disappears. This phenomenon is also referred to as a cycle blow-up. In the full system J describes the blow-up of a three-dimensional torus. Numerical studies on the Hopf–Hopf-bifurcation with MatCont showed that the above-mentioned linearisations are valid only in a small region. Upon leaving this region higher order terms gain importance and the curves start to bend and may leave the quadrant they started in. This leads to a complex bifurcation diagram for the Hopf–Hopf-bifurcation, which is to be discussed in the following. When starting in region (i) the limit cycle already exists and disappears when crossing C towards region (ii). The quadrant consisting of regions (i and ii) equals the linear stable region since there are no fixed points on the coordinate axes. Crossing H1 towards region (iii) creates the saddle E1 on n1 . This corresponds to the excitation of the fundamental mode in the ROM. Similarly crossing H2 towards region (iv) creates the unstable fixed point E2 on n2 , which corresponds to the excitation of the first azimuthal mode. Passing T 2 towards region (v) creates the non-trivial saddle E3 and arriving in region (vi) annihilates E1 at H1 . Crossing C results in the creation of the invisible limit cycle around E3 and in a change of stability of this fixed point. When traversing T 1 towards region (vii) E3 disappears. The limit cycle’s blow-up happens when crossing J and the cycle becomes visible in region (viii). The tour ends in region (i) with the annihilation of E2 at H2 . It should be noted that for a generic system like the ROM the normal form only guarantees the existence of H1;2 and T 1;2 , but not C and J, which is due to the influence of higher order terms. Yet their existence in the ROM was shown indirectly (see above). In the following the Hopf–Hopf-bifurcation dynamics’ practical consequences for the full system will be discussed separately for

107

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

every region. Here, practical consequence means the consideration of asymptotic dynamics. The corresponding phase spaces are shown in Figs. 26–33. These diagrams are sketches of phase spaces that have been computed by setting l1 and l2 to their appropriate values and performing manual numerical integration for several initial conditions for n1 and n2 . To learn about the unstable fixed point and saddles MatConts backward integration capabilities, i. e. integration with reversed (negative) time, have been used.  Region (i): When starting inside of the unstable three-dimensional torus the stable two-dimensional torus is reached. When starting outside of it the system is divergent in the neighborhood of the Hopf–Hopf-bifurcation but may converge to a stable invariant set that’s not described by the normal form.  Regions (ii–iv): The system reaches a stable two-dimensional torus not described by the normal form, but determining the system’s asymptotic behavior. This torus corresponds to an additional non-trivial fixed point of the planar system that eludes a linear approximation.  Region (v): The system evolves to either the stable torus or the stable limit cycle on the n1 -axis.  Region (vi): The system lands on the stable torus or the trivial stable fixed point.  Regions (vii–viii): Starting inside of the unstable three-dimensional torus the system moves towards the stable two-dimensional torus. Starting outside of it means reaching the trivial fixed point. This discussion shows that the asymptotic dynamics of the normal form are governed by an additional non-trivial fixed point in the regions (ii–vi) not described by the normal form itself.

Appendix C C.1. MatCont summary (Dhooge et al., 2007, 2006) In the framework of in depth analysis of dynamical systems, in particular stability analysis, bifurcation analysis with appropriate software is common practice. In the current study, the continuation toolbox MatCont embedded in the MATLAB environment is used. The objective of MatCont is that  it covers as many bifurcations in ODEs as possible (see Fig. 34 and 35), in particular all bifurcations with two control parameters can be revealed (Dhooge et al., 2006) and  it implements robust and efficient numerical methods for all computations (Dhooge et al., 2006).  Furthermore It computes normal form coefficients of the most bifurcation types like of the Hopf–Hopf bifurcation In the framework of continuation we consider the mathematical ! ! problem 0 ¼ F ðX ; ~ cÞ. The objective of numerical continuation is to estimate branches in the neighborhood of bifurcations. These branches are approximated by a consecutive sequence of points which are estimated by the continuation solver. MatCont employs a predictor–corrector method to estimate the points on the branch. ! All points X i , i ¼ 1; 2; . . . forming the branch have to satisfy a se! ! lected tolerance criterion k F ðX i!; ~ cÞk 6 e for some e > 0 and an additional accuracy condition kdX i k 6 e0 where e0 > 0. The crucial point of application of the continuation technique is to find the!first point on the specific bifurcation branch. If the initial point X i on the curve have been found the normalized tangent

!

!

yi at X i will be computed, yi ¼ 0 and vector ~ where F x ðX i ; ~ cÞ~ ! hyi ; yi i ¼ 1. The estimation of X iþ1 requires to steps: (1) the prediction of the new point and (2) the correction of the predicted point. Appendix D D.1. Model reduction (ROM) A BWR loop constitutes a nonlinear dynamical system described by a set of partial differential equations !

! @ X ðtÞ ! n ! ¼ F ð@ X ðtÞ; . . . ; X ðtÞ; ~ cÞ @t

ð36Þ

in the case of distributed parameter systems like system codes5 or ordinary differential equations !

@ X ðtÞ ! ! ¼ F ðX ; ~ cÞ @t

ð37Þ

in the case of lumped parameter systems (reduced order model equation sets6).

r; tÞ @xi ð~ r; tÞg; f@ n xj ð~ r; tÞg; cÞ; ¼ F i ðfxj ð~ @t ¼ 1; . . . ; N:

Or in components :

i; j

We assume a complet set of orthonormal function fum ð~ rÞg exist, the so - called spatial modes. Expansion xi ð~ r; tÞ in this basis:

xi ð~ r; tÞ ¼

M X xi;m ðtÞum ð~ rÞ yields m¼1

dxi;m ðtÞ ¼ dt

*

uym ; F i

d~ xðtÞ ! ¼ F ð~ x; cÞ; dt

( X m

) ( ) !+ X xj;m ðtÞ@ n um ð~ rÞ ; c

xj;m ðtÞum ð~ rÞ ;

m

autonomous system hf ð~ rÞi ¼

Z

f ð~ rÞd~ r

ð38Þ

Vol

Model reduction in the sense above (elimination of the spatial variable) leads to a set of N  M coupled nonlinear ordinary differential equations (ODE’s). The coupling reflects that the excitation of the spatial modes um occurs not independently. The modes should be chosen on the basis of physical criteria. The Reduced Order Model (ROM) must be adequately to the origin model regarding his stability behavior. If this is true it is easier, in a first step, to analyze the ROM equations concerning their stability behavior in order to get an overview about the stability properties of the origin dynamical system (Zhou and Uddin, 2002; Dokhane et al., 2004; Lange et al., 2009, 2011). The last group of formulas demonstrate the mode expansion approach and the conversion of the set of coupled, nonlinear PDEs solved in integral system codes to the ROMs, a set of coupled, nonlinear ODEs. ROMs are, as mentioned, in essential reduced regarding the space dependence of the dynamical system. The spacedependence is taken into account approximately by the mode expansion terms. Hence ROMs are more convenient and of highly physical transparency in comparison to the PDF sets in system codes. In addition ROMs could be couplet straightforward with special methods of the nonlinear dynamics like bifurcation analysis codes (Hassard et al., 1981; Dhooge et al., 2006). It was proved in 5 System code: integrated BWR loop model, the coupled neutron kinetic and thermal-hydraulic PDE‘s are solved. 6 Reduced order model: spatial averaged PDE but the stability properties of the resulting ordinary DE (ODE) are very similar to the PDE of the system code equations.

108

C. Lange et al. / Annals of Nuclear Energy 67 (2014) 91–108

the past it is helpful to use system codes and ROMs complementary for BWR stability analysis (we call it RAM-ROM method (Lange et al., 2009, 2011). References Aguirre, C., Hennig, D., February 2002. Stability Measurement at NPP Leibstadt Cycle 7, September 6th 1990, NACUSP, FIKS-CT-2000-00041. AnT 4.669-Description, 2003. University of Stuttgart. Blomstrand, J., 1992. The KKL Core Stability Test, Conducted on September 6th, 1990, ABB, Report BR-254. Dhooge, A., Govaerts, W., Kuznetsov, Yu.A., Mestrom, W., Riet, A.M., 2006. MATCONT and CL-MATCONT: Continuation Toolboxes in MATLAB. Utrecht University, The Netherlands. Dhooge, A., Govaerts, W., Kuznetsov, YU.A., Meijer, H.G.E., Sautois, B., 2007. New features of the software MatCont for bifurcation analysis of dynamical systems. Mathematical and Computer Modeling of Dynamical Systems. Dokhane, A., Hennig, D., Chawla, R., Rizwan-uddin, 2003. Nuclear-coupled thermalhydraulic nonlinear stability analysis using a novel BWR reduced order model: Part 1 – The effects of using drift flux versus homogeneous equilibrium models. In: Proc. Of 11th International Conference on Nuclear Engineering, Tokyo, Japan. Dokhane, A., Hennig, D., Chawla, R., Rizwan-uddin, 2003b. Nonlinear stability analysis of boiling water reactor on the basis of system codes and reduced order models. Proceedings of Annual Meeting on Nuclear Technology 126, 15–19. Dokhane, A., 2004. BWR Stability and Bifurcation Analysis using a Novel Reduced Order Model and the System Code RAMONA. Doctoral Thesis, EPFL, Switzerland. Guckenheimer, Holmes, P., 1984. Nonlinear oscillations, dynamical systems and bifurcation of vector fields. In: Applied Mathematical Sciences. Springer Verlag, 42. Hassard, B.D., Kazarinoff, N.D., Wan, Y.H., 1981. Theory and Application of Hopf Bifurcation. Cambridge University Press, Cambridge, London, New York, Melbourne, Sydney. Henry, A.F., 1975. Nuclear Reactor Analysis. MIT Press, Cambridge; Massachusetts; USA. Kuznetsov, Y.A., 1998. Elements of applied bifurcation theory, . second ed. Applied Mathematical Science second ed., vol. 112 Springer-Verlag, New York. Lange, C., 2009. Advanced Nonlinear Stability Analysis of Boiling Water Nuclear Reactors. PhD Thesis, TU Dresden, Germany. Lange, C., Hennig, D., Hurtado, A., 2011. An advanced reduced order model for BWR stability analysis. Progress in Nuclear Energy 53, 139–160. Lange, C., Hennig, D., Hurtado, A., 2012a. A saddle-node bifurcation of cycles found in the parameter space of a BWR. International Journal of Bifurcation and Chaos 22 (2).

Lange, C., Hennig, D., Hurtado, A., Schuster, R., Lukas, B., Aguirre, C., 2012b. Remarks on boiling water reactor stability analysis – Part 1: Theory. Kerntechnik 77, 5. Lange, C., Hennig, D., Hurtado, A., Schuster, A., Lukas, B., Aguirre, C., 2012c. Remarks on boiling water reactor stability analysis – Part 2: Stability monitoring. Kerntechnik 77, 6. Lange, C., Hennig, D., Schulze, Manuel, Hurtado, A., 2013. Comments on the application of bifurcation analysis in BWR stability analysis. Progress in Nuclear Energy 68, 1–15. March-Leuba, J., 1984. Dynamic Behaviour of Boiling Water Reactors. Ph.D. Dissertation; University of Tennessee; USA. March-Leuba, José, Blakeman, E.D., 1991. A mechanism for out-of-phase power instabilities in boiling water reactors. Nuclear Science and Engineering 107 (2), 173–179. March-Leuba, J., Cacuci, Dan G., Perez, R.B., 1986. Nonlinear dynamics and stability of boiling water reactors: qualitative analysis. Nuclear Science and Engineering 93, Part1: pp. 111–123; Part2: pp. 124–136. Miro, M., Ginestar, D., Hennig, D., Verdu, G., 2000. On the regional oscillation phenomenon in BWR’s. Progress in Nuclear Energy 36 (2), 189–229. Miro, R., Ginestar, D., Verdù, G., Hennig, D., 2002. A nodal modal method for the neutron diffusion equation. Application to BWR instabilies analysis. Annals of Nuclear Energy 29, 1171–1194. Neyfeh, A.H., Balachandran, B., 1995. Applied Nonlinear Dynamics. John Wiley & Sons Inc., New York, Singapore. Seydel, R., 2010. Practical bifurcation and stability analysis, . third ed. Interdisciplinary Applied Mathematics third ed., vol. 5 Springer, New York. Strogatz, S.H, 1994. Nonlinear Dynamics and Chaos. Addison-Wesley. Studsvik/Scandpower, 2003. RAMONA5.5 Manual. Verdù, G., Ginestar, D., Vidal, V., Munoz-Cobo, J.L., 1994. 3D k-modes of the neutron diffusion equation. Ann. of Nuclear Energy 21 (7), 221–405. Weinberg, A.M., Schweinler, H.C., 1948. Theory of oscillating absorber in a chain reactor. Physical Review, Second Series 74, 8. Weinberg, A.M., Wigner, E.P., 1958. Physical Theory of Neutron Chain Reactors. University of Chicago Press, Chicago, Illinois. Wiktor, C.G., 1997. KKL Cycle7 Stability Test-Core Simulation Input Data, KKL, Report BET/97/111. Wulff, W., Cheng, H.S., Diamond, D.J., Khatib-Rhabar, M., 1984. A Description and Assessment of RAMONA-3B Mod.0 Cycle4: A Computer Code with ThreeDimensional Neutron Kinetic for BWR System Transients. NUREC/CR-3664, BNL-NUREC-51746, Brookhaven National Laboratory. Zhou, Q., Uddin, R., 2002. Bifurcation analyses of in-phase and out-of-phase oscillations in BWRs. In: Proceedings of the International Conference of the New Frontiers of Nuclear Technology Reactor Physics, Safety and High Performance Computing (Physor2002), Seoul, Korea, October 7–10, CD-ROM, 2002.