- Email: [email protected]
Bifurcation analysis of out-of-phase oscillations in boiling water reactors using multipoint neutron kinetics
Progress in Nuclear Energy 120 (2020) 103218
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: http://www.elsevier.com/locate/pnucene
Bifurcation analysis of out-of-phase oscillations in boiling water reactors using multipoint neutron kinetics Rohan Vora a, Abhishek Chakraborty b, Suneet Singh a, * a b
Department of Energy Science and Engineering, Indian Institute of Technology Bombay, Mumbai, 400076, Maharashtra, India Nuclear Power Corporation of India Limited, Anushaktinagar, Mumbai, 400094, Maharashtra, India
A R T I C L E I N F O
A B S T R A C T
Keywords: Stability analysis Multipoint reactor kinetics Hopf bifurcation Out-of-phase oscillations
Stability analysis is important for identifying the safe operating regimes in Boiling Water Reactors (BWRs). The power instabilities, due to the fuel and coolant void feedbacks, can be either in-phase or out-of-phase. The models used in earlier studies were quite complicated, hence, not suitable for bifurcation analysis. In this paper, a simpler model based on neutron multipoint kinetics coupled with coolant temperature feedbacks, with reason able accuracy is presented for bifurcation analysis of spatial oscillations in BWRs. The stability maps in different parametric planes are constructed using MATCONT. Both supercritical and subcritical Hopf bifurcations are observed and the Generalized Hopf points are also identified. Identification of subcritical Hopf is important as the linearly stable region may have unstable limit cycles near the stability boundary. The region in the stability map where in-phase oscillations are stable whereas the out of phase oscillations are unstable is identified and verified by numerical simulations.
(ii) Out-of-phase or regional power oscillations. 1. Introduction Boiling Water Reactor (BWR) instabilities such as periodic oscilla tions in power, density, flow rates and temperature have been experi mentally observed and studied. BWRs behave as linear systems under normal operating conditions (Waranpera and Anderson, 1981). How ever, operating at high power levels and low flow conditions, BWRs can be susceptible to limit cycle power oscillations. These oscillations are typically a nonlinear phenomenon (March-Leuba, 1989) and hence BWRs can behave non-linearly at high power levels and low flow con ditions. The instabilities in BWR is a complex phenomenon and can be classified into three types, namely plant instability, reactivity instability and channel thermal hydraulic instabilities (March-Leuba, 1984). Plant instability occurs from external disturbance; reactivity instability occurs when the feedback reactivity coefficients (namely coolant void and Doppler) are strong enough to cause power oscillation when a pertur bation is applied; and channel thermal hydraulic instabilities such as density wave oscillation which occurs due to the dynamics of two phase flow in a heated channel. Neutronics coupled with thermal hydraulics can result in power oscillations of two kinds: (i) In-phase or global power oscillations, and
In global power oscillations, the total power of the reactor oscillates in-phase whereas in out-of-phase oscillations, regional power oscilla tions take place. Thermal hydraulic instabilities lead to local power oscillations due to interactions with density reactivity coefficient feed back, flow rate and pressure drop in the channel whereas, the global oscillations are not much affected (Dokhane, 2004). The in-phase or global power oscillations can be studied using point reactor kinetics coupled with thermal hydraulics as given in (March- Leuba et al., 1986; Rizwan-Uddin, 2006; KarveRizwan-uddin and ~ oz-Cobo et al., 2000;Mun ~ oz-Cobo et al., 2016). In Dorning, 1997; Mun point reactor kinetics model, the reactor is assumed to be a point and spatial variation is not considered i.e., it is characterized by a time in dependent spatial shape function and a time dependent amplitude fac tor. When the flux shape undergoes significant change with time, this model is unable to capture the changes. However, for studying the out of phase oscillations, thermal hydraulics coupled with space time depen dent neutronics are used for modelling the system dynamics. Different methods for solution of space time dependent neutron kinetic equations have been used like Improved Quasi-Static (IQS) (Dutta and Doshi, 2009), modal kinetics (Dokhane, 2004; Mu~ noz-Cobo et al., 2016; Dokhane et al., 2007a, 2014; Lange et al., 2011), multi region neutron
* Corresponding author. E-mail address: [email protected] (S. Singh). https://doi.org/10.1016/j.pnucene.2019.103218 Received 1 June 2019; Received in revised form 23 October 2019; Accepted 29 November 2019 Available online 23 December 2019 0149-1970/© 2019 Elsevier Ltd. All rights reserved.
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Analytical Bifurcation method (Hassard et al., 1981) and numerical continuation method (Kuznetsov, 2004) are most commonly used for carrying out non-linear stability analysis. Semi analytical Bifurcation method is used in (Dokhane, 2004; Rizwan-Uddin, 2006; KarveR izwan-uddin and Dorning, 1997; Dokhane et al., 2014; Chakraborty et al., 2019) and continuation method is used in (Pandey and Singh, 2016; Chakraborty et al., 2018). The reduced order models as developed by (KarveRizwan-uddin and ~ oz-Cobo et al., 1996), based Dorning, 1997; Dokhane et al., 2007a; Mun on modal kinetics coupled with thermal hydraulics equations, involve use of system codes for approximating the cross coupling reactivities for each mode due to Doppler and coolant void feedbacks. The equations represented even with two thermal hydraulic channel model becomes very complicated. The thermal hydraulics coupled with multi region kinetics also lead to differential equations which are very cumbersome for analysis. The model developed by March-Leuba et al. (1986) to describe the behaviour of BWRs is very simple as compared to the other models involving explicit thermal hydraulics. The neutron kinetics in this model is based on point reactor kinetics and certain assumptions were made to obtain a simple thermal hydraulic model. The coolant is assumed to be at saturation temperature when it enters the core and the complete recirculation loop is treated as a one path where fluid travels with constant mass flow rate. This model was successful in the predic tion of the total power oscillations in BWRs. Detailed bifurcation anal ysis using this model using point kinetics was carried out using MATCONT in (Pandey and Singh, 2019). In this paper, multipoint kinetics are used for modelling the out of phase oscillations in a BWR and the equations representing the thermal hydraulic feedbacks are based on the model given in (March-Leuba, 1984). Linear and non-linear stability analysis has been performed using numerical continuation code MATCONT (Dhooge et al., 2004) and the results have been compared with one point model. Even though this model is very simple, it captures the dynamics and the results are in line with those predicted using more detailed thermal hydraulic codes.
Nomenclature nðtÞ N0 cðtÞ TðtÞ
ρ Λ β λ
ρα D
τ
H v0 G0
α12
Aij, Aji Dj dij
νi
Vi
Excess neutron population normalized to steady state neutron population Steady state neutron population Excess delayed precursors population normalized to steady state neutron population Excess fuel temperature Reactivity feedback Prompt neutron generation time Delayed neutron fraction Decay constant of delayed neutron precursors Coolant void feedback reactivity Doppler coefficient of reactivity Residence time of bubbles in the fluid channel Height of the reactor Propagation velocity Coolant mass flux Coupling coefficient between region 1 and 2 Interface area between two adjacent regions “i" and “j" One group diffusion coefficient for region ‘j’ Distance between the centroids of the two regions “i" and “j" Neutron speed in region “i" Volume of the region “i"
2. Mathematical model of a Boiling Water Reactor 2.1. Space independent model
Fig. 1. Block diagram for BWR model (March-Leuba, 1984).
The dynamic processes in a BWR involves the production of heat i.e., a neutronics loop and its removal from the core by light water in two phase (as coolant), represented by a thermal hydraulics loop. The main feedbacks are due to the rate of formation of coolant void and change in fuel temperature. The different processes can be represented in a block diagram form and is given in Fig. 1. The space independent dynamics of a BWR can be modelled using a fifth order model consisting of coupled nonlinear ordinary differential equations as given in March-Leuba et al. (1986) and are given below:
kinetics coupled with thermal hydraulics (UEHIRO et al., 1996; Lee and Pan, 2005). In order to identify stable and unstable regimes of operation, stability analysis is required. A straightforward approach for carrying out sta bility analysis is by directly solving the neutronics coupled with thermal hydraulics equations using system codes (Dutta and Doshi, 2009; Lange ~ oz-Cobo et al., 1996, 2000). However, et al., 2011; Wulf et al., 1984; Mun for identifying stable and unstable regimes, the simulations need to be carried out at different parametric values which is very computationally expensive and tedious since the entire range of parameters have to be spanned. In this regard, reduced order models are developed which are simpler than the explicit models but are very useful in capturing the phenomena. These reduced order models can be used for carrying out stability analysis for identification of stable and unstable regimes. Standard techniques such as frequency domain approach (March-Leuba, 1984) and eigenvalue approach (Rizwan-Uddin, 2006) can be used for arriving at the linearly stable regimes. Since BWRs behave as non-linear systems, linear stability analysis is not sufficient for identifying stable regimes. Non-linear stability analysis is also required because a “line arly” stable operating point can be “non-linearly” unstable. This is due to the fact that linear stability assesses the system dynamics for small perturbations and thereby for large perturbations, the system behaviour may be different. Hence, the identification and characterisation of such phenomena and the existence of Generalized Hopf bifurcation is very important in the identification of the stability of the system. Semi
dnðtÞ ρðtÞ β ρðtÞ ¼ nðtÞ þ λcðtÞ þ dt Λ Λ
(1)
dcðtÞ β ¼ nðtÞ dt Λ
λcðtÞ
(2)
dTðtÞ ¼ a1 nðtÞ dt
a2 TðtÞ
(3)
d2 ρα ðtÞ dρ ðtÞ þ a3 α þ a4 ρα ðtÞ ¼ kTðtÞ dt2 dt where nðtÞ ¼ ðNðtÞ NO Þ=NO and cðtÞ ¼ ðCðtÞ The reactivity feedback can be expressed as ðtÞ ¼ ρα ðtÞ
DTðtÞ
(4) CO Þ=NO . (5)
Eq. (4) is a second order differential equation and can be converted into two first order equations by the following substitution, 2
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Table 1 Model parameter values for Vermont Yankee test 7N (March-Leuba, 1984). Parameter
Value
Units
k
0.0037
K 1s
D
0.0000252
K
0.0056
–
Λ
0.00004
s
λ
0.08
s
a1
25.04
Ks
a2
0.23
s
a3
2.25
s
6.82
s
(6)
Substituting Eq. (6) in Eq. (4), we get
2
dρt ðtÞ þ a3 ρt ðtÞ þ a4 ρα ðtÞ ¼ kTðtÞ dt
1
β
a4
dρα ðtÞ ¼ ρt ðtÞ dt
(7)
The parameter k is related to fuel temperature and void coefficient of reactivity (March-Leuba, 1984). The parameters a3 and a4 scale as 1=τ and 1=τ2 , where τ is the residence time of bubbles in the fluid channel (March-Leuba, 1984). The residence time, τ ¼ H=v0 , H being the height of the reactor and v0 is the propagation velocity.v0 is directly propor tional to the coolant mass flux G0 as given in (March-Leuba, 1984) and hence τ is inversely proportional to G0 .The values of the parameters can be found in Table 1.
1 1 1
1 1 2
Fig. 2. Regions considered for analysis.
Fig. 3. Sub Critical and Super critical Hopf Bifurcation. 3
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
2.2. Space - time dependent model – multipoint kinetics
dT2 ðtÞ ¼ a1 n2 ðtÞ dt
In this paper, multipoint kinetics technique is used for solution of space-time dependent neutron diffusion equation. In this method, the core is divided into finite number of volumes with point kinetic equa tions for each region representing the dynamics and a coupling term between the adjacent regions which accounts for the neutron leakage. Since it has been observed that the oscillations in BWRs take place in the azimuthal direction (Dokhane, 2004; Dutta and Doshi, 2009; Dhooge et al., 2004) the reactor core is divided into two regions for the present analysis as shown in Fig. 2. The subscript 1 indicates the left region and the subscript 2 indicates the right region of the core with α12 repre senting the coupling between them. A fifth order model given in March-Leuba et al. (1986) is then adopted for both halves of the core separately. The coupling coefficients, αij ; representing the coupling between the ith and jth regions are calculated using the formula (Chakraborty et al., 2019)
αij ¼
(8)
ρ2 ðtÞ ¼ ρα2 ðtÞ þ DT2 ðtÞ
kc h
(9)
dhðtÞ ¼ kf ðn1 þ n2 Þ dt
ρ1 ðtÞ ¼ ρα1 ðtÞ þ DT1 ðtÞ
n1 ðtÞ þ
α12 Λ
n2 ðtÞ þ λc1 ðtÞ þ
ρ1 ðtÞ Λ
λc1 ðtÞ
(13)
dT1 ðtÞ ¼ a1 n1 ðtÞ dt
a2 T1 ðtÞ
(14)
dρα1 ðtÞ ¼ ρt1 ðtÞ dt
(15)
dρt1 ðtÞ þ a3 ρt1 ðtÞ þ a4 ρα1 ðtÞ ¼ kT1 ðtÞ dt
(16)
dn2 ðtÞ ρ ðtÞ β ¼ 2 dt Λ
α12
ρ2 ðtÞ ¼ ρα2 ðtÞ þ DT2 ðtÞ dc2 ðtÞ β ¼ n2 ðtÞ dt Λ
λc2 ðtÞ
n2 ðtÞ þ kc h
α12 Λ
n1 ðtÞ þ λc2 ðtÞ þ
ρ2 ðtÞ Λ
dhðtÞ ¼ kf ðn1 þ n2 Þ dt
(23)
(24)
where x is n dimensional state and r is the bifurcation parameter. The equilibrium solution is described as follows Fðx0 ; rÞ ¼ 0
(11)
dc1 ðtÞ β ¼ n1 ðtÞ dt Λ
(22)
x_ ¼ Fðx; rÞ
(25)
The characteristic roots of the Jacobian matrix, Jðx0 ; rÞ for the above system help determine the stability of local fixed points. Conditions required to be satisfied for Hopf bifurcation are as follows.
(12)
kc h
dρt2 ðtÞ þ a3 ρt2 ðtÞ þ a4 ρα2 ðtÞ ¼ kT2 ðtÞ dt
Hopf bifurcation is a codimension-1 bifurcation where there is only one free parameter for a bifurcation to occur. Changes in the parameter value of the system results in the change of its eigenvalues. These ei genvalues help determine the linear stability of a system in steady state. Positive real part of any one of the eigenvalues makes the system un stable whereas, when real parts of all the eigenvalues are negative, the system is stable. Hopf bifurcation occurs at the point at which a pair of complex eigenvalues cross the imaginary axis. Linear stability of the system is lost at this critical parameter value and a family of limit cycle bifurcates at that point, which is a non-linear phenomenon. Consider an autonomous system of ordinary differential equations
(10)
α12
(21)
3.1. Hopf Bifurcation
BWR system can now be modelled with multipoint kinetics as follows dn1 ðtÞ ρ ðtÞ β ¼ 1 dt Λ
dρα2 ðtÞ ¼ ρt2 ðtÞ dt
Hopf proved the existence of stable and unstable periodic oscillations for a non-linear dynamic system. Hopf theorem, sometimes also called as the Poincar�e–Andronov–Hopf (PAH) bifurcation theorem guarantees the existence of periodic solutions under certain conditions (Hassard et al., 1981; Guckenheimer and Holmes, 1983). The BWR system modelled by March – Leuba et al. (March-Leuba et al., 1986) meets all the conditions required by Hopf theorem. Since only Hopf bifurcation is observed in this study, a brief review of it and its mathematical background is given in this section. A detailed study on it can be found in (Kuznetsov, 2004) for more insights.
For the estimation of αij , the different constants are required which depend on the core geometry and fissile inventory. So for a two-region core, the coupling coefficient is inversely proportional to the distance between the two regions (dij ). This implies that a smaller core will have higher value of coupling coefficient since it is inversely proportional to the distance between the regions and vice versa. It is assumed that the total power oscillations are stable and the model only studies the spatial oscillations. In the study presented by ~ oz-Cobo et al. (Dhooge et al., 2004) constant pressure boundary Mun conditions were enforced to observe out-of-phase oscillations. Since in this study a simpler March-Leuba model (March-Leuba et al., 1986) is adopted, a feedback term is introduced to reactivity equation to observe out-of-phase oscillations. The feedback term introduced is kc h as shown in Eq. (8) and Eq. (9) where h is defined as shown in Eq. (10). kc h
(20)
3. Mathematical background
Di νi ΛAij dij Vj
ρ1 ðtÞ ¼ ρα1 ðtÞ þ DT1 ðtÞ
a2 T2 ðtÞ
1) Eigenvalues of Jðx0 ; rc Þ has at least one pair of imaginary eigen values, � iω, lying on the imaginary axis. Here rc represents the parameter’s critical value where Hopf bifurcation occurs. 2) The eigenvalues of Jðx0 ; rc Þ are of the form σðrÞ � iω ðrÞ and dσ ðr ¼ rc Þ=dλ 6¼ 0, i.e., the real parts of the eigenvalues of Jðx0 ; rc Þ cross the imaginary axis with non-zero speed as r crosses rc . 3.2. Supercritical and subcritical Hopf bifurcation Hopf theorem establishes the coexistence of limit cycles and steady state solutions. Limit cycles are isolated closed trajectories where neighbouring trajectories spiral towards or away from the limit cycles. Limit cycles are classified and as stable and unstable limit cycles. In stable limit cycles, neighbouring trajectories are attracted towards it, whereas in unstable limit cycles the neighbouring trajectories are repelled. Hopf bifurcation can again be classified as supercritical or subcritical Hopf bifurcation as shown in Fig. 3. In a subcritical bifur cation (Fig. 3a) unstable limit cycle bifurcates from a stable fixed point.
(17) (18) (19)
4
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
4. Results and discussions
Table 2 Additional parameter values for multipoint model. Parameter
Value
kc
0.01
kf
50
α12
0.0035
4.1. Model verification 4.1.1. Stability maps Bifurcation curves are generated using MATCONT (Dhooge et al., 2004) software. Stability boundaries for the BWR system described by multipoint model represented by Eq. (11) - Eq. (23) are identified. The parameters k and a3 are chosen as the bifurcation parameters, similar to the work presented in (Pandey and Singh, 2016). These two parameters are varied keeping the other parameters from Tables 1 and 2 constant. The bifurcation plot for k vs a3 given in Fig. 4 represents the region where the system changes its stability. As the parameter k is changed, the real part of eigenvalue changes its sign thereby implying Hopf bifurcation. Apart from identifying the Hopf points, first Lyapunov co efficients at each point are estimated using MATCONT which indicates the nature of the Hopf Bifurcation. Negative first Lyapunov coefficient indicates a supercritical Hopf bifurcation whereas subcritical Hopf bifurcation is indicated by positive first Lyapunov coefficient. Fig. 5 shows that the first Lyapunov coefficient is negative in the plotted region implying the nature of Hopf bifurcation is supercritical. The study presented in Lee and Pan (2005) concludes that for strong
Any perturbation with amplitude less than that of the limit cycle causes neighbouring trajectories to be attracted to the stable fixed point. Any perturbation with amplitude higher than that of limit cycle also repels neighbouring trajectories from the unstable limit cycle. Hence, in subcritical condition after bifurcation the trajectories must to jump to a distant attractor. In a supercritical Hopf (Fig. 3b) bifurcation system the physical system initially settles down to equilibrium through damped oscillations. As the parameter value changes, the decay decreases and slowly loses its stability resulting in small amplitude, sinusoidal limit cycles oscillation about the former steady state. Hence, in supercritical Hopf bifurcation a stable limit cycles bifurcates from an unstable fixed point. Any trajectory with lesser or higher amplitude than that of the limit cycle is attracted towards it.
Fig. 4. Bifurcation plot for k vs a3 at α12 ¼ 0.0035.
Fig. 5. First Lyapunov coefficient for k vs a3 bifurcation plot at α12 ¼ 0.0035. 5
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 6. Bifurcation plot for k vs α12
Fig. 7. First Lyapunov coefficient for k vs α12 bifurcation plot.
Fig. 8. Bifurcation plot of k vs a3 at α12 ¼ 0.0063.
6
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 9. First Lyapunov coefficient for k vs a3 bifurcation plot at α12 ¼ 0.0063.
Fig. 10. Bifurcation plot for k vs a3 for different values of α12
illustrated by values of first Lyapunov coefficient shown in Fig. 7. With a higher value of coupling coefficient α12 in the multipoint kinetics model, both supercritical and subcritical Hopf bifurcation are observed. Fig. 8 shows the occurrence of supercritical and subcritical Hopf bifurcation when α12 ¼ 0.0063. First Lyapunov coefficient for this bifurcation plot is shown in Fig. 9. Fig. 10 compares the bifurcation plot for k vs a3 for different values of coupling coefficient α12 . For lower values of α12 , only supercritical Hopf bifurcation is observed. By increasing the value of the coupling coeffi cient, both supercritical and subcritical Hopf bifurcation are observed, thus confirming the results presented in Lee and Pan (2005) The stability map can be divided into four regions I, II, III and IV. Since the coupling coefficient is inversely proportional to the size of the reactor core, the lowest value corresponds to the highest value of coupling coefficient and vice versa. So, it can be seen from the figure that for α12 ¼ 0.0010, the regions I, II and III are unstable. However for α12 ¼ 0.0063, which represents a smaller core size, only Region I is unstable. Hence, higher the value of the coupling coefficient α12 , larger is the parametric space where the system is stable.
Table 3 Parameter values at various points for numerical simulations. Point number
a3
a4
k
kc
kf
α12
1
3.1
6.82
1 � 10
2
2.1
6.82
8:1 � 10
3
0.01
1
0.0035
0.01
1
3
3.24
1.705
4:3 � 10
0.0035
3
0.01
1
4
3.24
1.705
4:3 � 10
0.0035
3
0.01
50
0.0035
2
Table 4 Updated values of a3 and a4 for multipoint model. Parameter
Value
Units
a3;m ¼ (a3;one )/2
1.125
s
a4;m ¼ (a4;one )/4
1.705
s
1 2
neutron interaction, which corresponds to a higher value of coupling coefficient α12 , both supercritical and subcritical Hopf bifurcations are observed. Similar results are obtained observed in this study by adopting the multipoint kinetics model. Bifurcation plot of k vs α12 in Fig. 6 shows the occurrence of Generalized Hopf (GH) point. GH points indicates the change in nature of Hopf bifurcation from supercritical to subcritical as
4.1.2. Numerical simulations The dynamics of n1 and n2 are obtained by solving Eq. (11) - Eq. (23). It can be seen that the fixed point comes out to be [n1 ðtÞ, c1 ðtÞ, T1 ðtÞ, ρα1 ðtÞ, ρt1 ðtÞ, n2 ðtÞ, c2 ðtÞ, T2 ðtÞ, ρα2 ðtÞ, ρt2 ðtÞ, h] ¼ ½0; 0; 0; 0; 0; 0; 0; 0; 0; 7
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 11. Trajectories of n1 and n2 operating at point number 1.
Fig. 12. Trajectories of n1 and n2 operating at point number 2.
Fig. 13. Bifurcation plot k vs a3 for multipoint kinetics.
8
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 14. First Lyapunov coefficient vs k for multipoint model.
Fig. 15. Bifurcation plot k vs a3 for one point kinetics.
Fig. 16. First Lyapunov coefficient vs k for one point model. 9
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 17. Comparison of stability boundaries for one and multipoint kinetics model.
Fig. 18. Trajectory of total n at point number 3 using one point model.
0; 0�. The system is then perturbed such that the initial point at t ¼ 0 is ½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0�. Numerical simulations are carried out at points whose parameter values are given in Tables 3 and 4. The parameters k and a3 are chosen close to the stability boundary shown in Fig. 4. Point number 1 from Table 3 is in the stable region of bifurcation plot. When the system is given a perturbation, trajectories of n1 and n2 decay to return to its steady state. This is illustrated in Fig. 11. Point number 2 from Table 3 is in unstable region with supercritical Hopf bifurcation. The perturbation applied to the system causes the neutron population to undergo a limit cycle with periodic oscillations as shown in Fig. 12. The excess neutron population n1 and n2 undergoes about 8–9 cycles of oscillations every 20 s. This is similar to the fre ~ oz-Cobo et al. (Dhooge et al., 2004) with 9–10 quencies reported in Mun cycles of oscillations every 20 s. Hence with appropriate parameter values, multipoint kinetics based approach can be used to model the BWR system.
model. With the mass flux being the half compared to that of one point model, the propagation velocity V0 becomes half (March-Leuba, 1984) and thereby the residence time of the bubbles in the channel τ doubles. Hence in the multipoint model, the parameter a3;m is half of a3;one and a4;m is one fourth to that of one point model, a4;one (as given in Table 4). The other parameters k, a1 and a2 which are defined by March-Leuba (1984) remain unchanged when the core is divided into two halves. 4.2.1. Stability maps for one point and multipoint kinetics model The stability boundaries for the BWR system modelled with multi point kinetics are again generated using MATCONT with new values of a3 and a4 to compare the boundaries with one point kinetics model. Parameters k and a3 are chosen as the bifurcation parameters. The bifurcation plot for k vs a3 is given in Fig. 13. The first Lyapunov coef ficient is illustrated in Fig. 14 which indicates that the nature of Hopf bifurcation is supercritical. Similar bifurcation curve can be generated for BWR system with one point kinetics model which is represented by Eq. (1)–(7). Parameter values for one point model are obtained from Table 1. The stability boundary is illustrated in Fig. 15 again by varying k and a3 and the first Lyapunov coefficient plot is presented in Fig. 16. Hence the one point kinetics model suggests the system undergoes subcritical Hopf bifurcation along with supercritical Hopf bifurcation.
4.2. Comparison of one point and multipoint kinetics model Since the core is divided into two equal halves, the coolant mass flux in multipoint model is half of one point model. Hence the parameters a3 and a4 will be different for multipoint model as compared to one point 10
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 19. Trajectories of n1 and n2 at point number 3 using multipoint model.
Fig. 20. Trajectory of n1 and n2 at point number 3 with low perturbation.
Fig. 21. Phase portrait of Limit cycles. 11
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 22. Trajectory of total n at point number 3.
Fig. 23. Bifurcation plot k vs a3 for multipoint kinetics at kf ¼ 1 and kf ¼ 50
Fig. 24. Bifurcation plot for k vs kf indicating Hopf point occurs at same value of k
12
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 25. Trajectory of total n at point number 4 where kf ¼ 50.
Fig. 26. Trajectories of n1 and n2 for different kf.
13
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 27. Comparison of n1 vs n2 limit cycle for different kf.
The multipoint model indicates that only supercritical Hopf bifurcation exists. While in previous studies, Dokhane et al. (2007a) pointed out existence of subcritical bifurcation alone, in more recent tests as pre sented in Dokhane et al. (2014), the nature of limit cycles observed is supercritical only. However, parameters such as kc or α12 in multipoint model can be adjusted to obtain GH points as well. Fig. 17 shows the comparison of stability boundaries for one and multipoint kinetics model. As seen in Fig. 17, the stability boundary for two point model is right to the stability boundary of one point kinetics model. In region I, both one point and multipoint models are unstable. In region II, one point model is stable whereas multipoint model in unstable and in region III, both one point and multipoint models are stable. Hence in region II, which is between the two boundaries, the total power oscillations are stable where as the spatial oscillations are unstable which is consistent with the assumption as given in Section 2.2.
low perturbation results in a stable point and high perturbation results in jumping to a distant limit cycle. Fig. 20 shows the limit cycles along which n1 and n2 oscillate in the time domain and the phase portrait of the limit cycle in n1 vs T1 and n1 vs n2 plane is given in Fig. 21. 4.2.3. Effect of kf on stability boundary Fig. 22 shows the dynamics total n obtained by adding n1 and n2 obtained from multipoint model as opposed to directly simulating one point model shown in Fig. 18. Fig. 22 indicates that though the total power is being supressed, it is still not constant. Hence values of kf are varied to get desired results. It is observed that the values kf does not change the stability boundary. This is illustrated in Fig. 23 and Fig. 24. Fig. 23 indicates that bifurcation plot for k vs a3 is same for different values of kf whereas, Fig. 24 shows that the bifurcation plot does not change with kf for a given value of k. The total n trajectory at kf ¼ 50 and other parameter values from point number 4 in Table 3 for the same initial perturbation t ¼ 0, ½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0� is shown in Fig. 25. This normalised total excess neutron density is oscillating at zero with negligible amplitude. Changing the value of kf does not change the out of phase behaviour of n1 and n2 as shown in Fig. 26 and Fig. 27. Increasing the value of kf , n1 and n2 become more anti symmetrical as seen in Fig. 27, thereby cancelling each other out. This results in regional out of phase oscillations where the total oscillations remains constant. Hence with a suitable values of the coupling coefficient, α12 and feedback coefficients, kc and kf in Eq. (13) - Eq. (25), regional out of phase oscillations can be observed.
4.2.2. Numerical simulations for one point and multipoint kinetics model The fixed point for multipoint model is [n1 ðtÞ, c1 ðtÞ, T1 ðtÞ, ρα1 ðtÞ, ρt1 ðtÞ, n2 ðtÞ, c2 ðtÞ, T2 ðtÞ, ρα2 ðtÞ, ρt2 ðtÞ, h] ¼ ½0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0�. The system is then perturbed such that the initial point at t ¼ 0 is ½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0�. Similarly for one point model, from Eqs. (1)– (7), it can be observed that the fixed point comes out to be [nðtÞ, cðtÞ, TðtÞ, ρα ðtÞ, ρt ðtÞ] ¼ ½0; 0; 0; 0; 0�. The dynamics of n is then observed by creating a perturbation such that the initial point at t ¼ 0 is ½1; 1 ;10 ;0; 0�. One point model represents the dynamics of total n and multipoint model represents the dynamics of n1 and n2 . To verify the stability map shown in Fig. 18, simulations are carried out at point in between the two stability boundaries which is presented as point number 3 in Table 3. Fig. 18 is obtained by perturbing the one point model. The trajectory of n eventually reaches to its steady state value. Fig. 19 illustrates the dynamics of n1 and n2 are obtained by solving the multipoint model when the perturbation is applied with parameter values given at point number 3 in Table 3. It can be observed that n1 and n2 are out-of-phase and oscillate in a limit cycle. Even with low perturbation such as the initial point at t ¼ 0 is 0:001*½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0�, the trajectories of n1 and n2 grow and reaches the same limit cycle. This supports the claim that the Hopf bifurcation is supercritical where high and low perturbation results in the same limit cycle as opposed to subcritical Hopf bifurcation where
5. Conclusions A simple model for studying out of phase oscillations in a BWR is developed based on multipoint kinetics. The core is divided into two parts which is coupled together by the coupling coefficient α12 . Two thermal hydraulic channels, one representing each region is modelled. The feedback reactivities due to coolant void and fuel temperature feedback are modelled using equations similar to March-Leuba (1989). Generally, constant pressure boundary conditions (Mu~ noz-Cobo et al., 1996) or constant total mass flow rate conditions (Lee and Pan, 2005) are enforced for controlling in-phase oscillations. Since the model pre sented in this paper does not explicitly model the thermal hydraulics, a 14
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
feedback reactivity is used to control the total power. ~ oz-Cobo et al. The results of the model are compared with Mun (Dhooge et al., 2004) and are found to be consistent. However, it is to be noted that the parameters used in the analysis are dependent on the system parameters and have to be derived for each BWR as given in March-Leuba (1989). But once these constants are derived, then this simple model becomes sufficient for studying these power oscillations. The developed model is further used to perform bifurcation analysis using MATCONT. The bifurcation analysis shows that the BWR system undergoes Hopf bifurcation as the parameters are varied. Both super critical and subcritical Hopf bifurcation were also observed in different parametric planes. Limit cycles have been observed by plotting the trajectory of n1 and n2 , the excess neutron populations normalized to the steady state neutron population. The stability boundaries of one point and two - point kinetics model are also compared. The stability boundary for one point model represents the stability characteristics of global oscillations whereas that for the two point model gives the same for out of phase oscillations. It is found that there exists a region where the global oscillations are stable but local oscillations are unstable.
Dokhane, A., Ferroukhi, H., Pautz, A., 2014. On out-of-phase higher mode oscillations with rotation and oscillation. Ann. Nucl. Energy 67, 21–30. Dutta, G., Doshi, J.B., 2009. Nonlinear analysis of nuclear coupled density wave instability in time domain for a boiling water reactor core undergoing core-wide and regional modes of oscillations. Prog. Nucl. Energy 51 (8), 769–787. Guckenheimer, J., Holmes, P., 1983. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer-Verlag. Hassard, B., Kazarinoff, N., Wan, Y., 1981. Theory and Applications of Hopf Bifurcation. Cambridge University Press, New York. Karve, A.A., Rizwan-uddin, Dorning, J., 1997. Stability analysis of BWR nuclear-coupled thermal-hydraulics using a simple model. Nucl. Eng. Des. 177, 155–177. Kuznetsov, Y., 2004. Elements of applied bifurcation theory. In: Applied Mathematical Sciences, vol. 112. Springer. Lange, C., Hennig, D., Chair, A.H., 2011. An advanced reduced order model for BWR stability analysis. Progress Nuclear Energy J. 53 (1), 139–160. Lee, J.D., Pan, C., 2005. Dynamic analysis of multiple nuclear-coupled boiling channels based on a multi-point reactor model. Nucl. Eng. Des. 235. March-Leuba, J., 1984. Dynamic Behaviour of Boiling Water Reactors. PhD Thesis. The University of Tennessee, Knoxville. March-Leuba, J., 1989. Non linear dynamic and chaos in boiling water reactors. In: Noise and Nonlinear Phenomena in Nuclear Systems, vol. 192. March-Leuba, J., Cacuci, D., Perez, R., 1986. Nonlinear dynamics and stability of boiling water reactors: part 1 –qualitative analysis. Nucl. Sci. Eng. 23, 93–111. Mu~ noz-Cobo, J.L., Pdrezi, R.B., Ginestar, D., Escriv, A., Verdu, G., 1996. Non linear analysis of out of phase oscillations in boiling water reactors. Ann. Nucl. Energy 23 (16), 1301–1335. Mu~ noz-Cobo, J.L., Rosello, O., Miro, R., Escriva, A., Ginestar, D., Verdu, G., 2000. Coupling of density wave oscillations in parallel channels with high order modal kinetics: application to BWR out of phase oscillations. Ann. Nucl. Energy 27, 1345–1371. Mu~ noz-Cobo, J.L., Escriv� a, A., Domingo, M.-D., 2016. In-phase instabilities in BWR with sub-cooled boiling, direct heating, and spacers effects. Ann. Nucl. Energy 87, 671–686. Pandey, V., Singh, S., 2016. A bifurcation analysis of boiling water reactor on large domain of parametric spaces. Commun. Nonlinear Sci. Numer. Simul. 38, 30–44. Pandey, V., Singh, S., 2019. Bifurcations emerging from a double Hopf bifurcation for a BWR. Prog. Nucl. Energy 117. Rizwan-Uddin, 2006. Turning points and sub- and supercritical bifurcations in a simple BWR model. Nucl. Eng. Des. 236 (3), 267–283. Uehiro, M., Rao, Y.F., Fukuda, K., 1996. Linear stability analysis on instabilities of inphase and out-of-phase modes in boiling water reactors. J. Nucl. Sci. Technol. 33, 628–635. Waranpera, Y., Anderson, S., 1981. BWR stability tests reaching limit cycle threshold at natural circulation. Trans. Am. Nucl. Soc. 39, 868. Wulf, W., Cheng, H., Diamond, D., Khatib-Rahbar, M., 1984. Description and Assessment of RAMONA-3B Mod. 0 Cycle 4: a Computer Code with Three-Dimensional Neutron Kinetics for BWR System Transients. NUREG/CR-3664, United States.
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.pnucene.2019.103218. References Chakraborty, A., Singh, S., Fernando, M.P.S., 2018. A novel approach for bifurcation analysis of out of phase xenon oscillations using multipoint reactor kinetics. Nucl. Eng. Des. 328, 333–344. Chakraborty, A., Singh, S., Fernando, M.P.S., 2019. An improved reduced order model for nonlinear stability analysis of spatial xenon oscillations. Prog. Nucl. Energy 116, 62–75. Dhooge, A., Govaerts, W., Kuznetsov, Y.A., 2004. MATCONT: a matlab package for numerical bifurcation analysis of ODEs. SIGSAM Bull. 38 (1), 21–22. Dokhane, A., 2004. BWR Stability and Bifurcation Analysis Using a Novel Reduced Order � Model and the System Code RAMONA. PhD Thesis. Ecole Polytechnique F�ed� erale de Lausanne, Switzerland. Dokhane, A., Hennig, D., Rizwan-uddin, Chawla, R., 2007. BWR stability and bifurcation analysis using reduced order models and system codes: identification of a subcritical Hopf bifurcation using RAMONA. Ann. Nucl. Energy 34, 792–802.
15
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: http://www.elsevier.com/locate/pnucene
Bifurcation analysis of out-of-phase oscillations in boiling water reactors using multipoint neutron kinetics Rohan Vora a, Abhishek Chakraborty b, Suneet Singh a, * a b
Department of Energy Science and Engineering, Indian Institute of Technology Bombay, Mumbai, 400076, Maharashtra, India Nuclear Power Corporation of India Limited, Anushaktinagar, Mumbai, 400094, Maharashtra, India
A R T I C L E I N F O
A B S T R A C T
Keywords: Stability analysis Multipoint reactor kinetics Hopf bifurcation Out-of-phase oscillations
Stability analysis is important for identifying the safe operating regimes in Boiling Water Reactors (BWRs). The power instabilities, due to the fuel and coolant void feedbacks, can be either in-phase or out-of-phase. The models used in earlier studies were quite complicated, hence, not suitable for bifurcation analysis. In this paper, a simpler model based on neutron multipoint kinetics coupled with coolant temperature feedbacks, with reason able accuracy is presented for bifurcation analysis of spatial oscillations in BWRs. The stability maps in different parametric planes are constructed using MATCONT. Both supercritical and subcritical Hopf bifurcations are observed and the Generalized Hopf points are also identified. Identification of subcritical Hopf is important as the linearly stable region may have unstable limit cycles near the stability boundary. The region in the stability map where in-phase oscillations are stable whereas the out of phase oscillations are unstable is identified and verified by numerical simulations.
(ii) Out-of-phase or regional power oscillations. 1. Introduction Boiling Water Reactor (BWR) instabilities such as periodic oscilla tions in power, density, flow rates and temperature have been experi mentally observed and studied. BWRs behave as linear systems under normal operating conditions (Waranpera and Anderson, 1981). How ever, operating at high power levels and low flow conditions, BWRs can be susceptible to limit cycle power oscillations. These oscillations are typically a nonlinear phenomenon (March-Leuba, 1989) and hence BWRs can behave non-linearly at high power levels and low flow con ditions. The instabilities in BWR is a complex phenomenon and can be classified into three types, namely plant instability, reactivity instability and channel thermal hydraulic instabilities (March-Leuba, 1984). Plant instability occurs from external disturbance; reactivity instability occurs when the feedback reactivity coefficients (namely coolant void and Doppler) are strong enough to cause power oscillation when a pertur bation is applied; and channel thermal hydraulic instabilities such as density wave oscillation which occurs due to the dynamics of two phase flow in a heated channel. Neutronics coupled with thermal hydraulics can result in power oscillations of two kinds: (i) In-phase or global power oscillations, and
In global power oscillations, the total power of the reactor oscillates in-phase whereas in out-of-phase oscillations, regional power oscilla tions take place. Thermal hydraulic instabilities lead to local power oscillations due to interactions with density reactivity coefficient feed back, flow rate and pressure drop in the channel whereas, the global oscillations are not much affected (Dokhane, 2004). The in-phase or global power oscillations can be studied using point reactor kinetics coupled with thermal hydraulics as given in (March- Leuba et al., 1986; Rizwan-Uddin, 2006; KarveRizwan-uddin and ~ oz-Cobo et al., 2000;Mun ~ oz-Cobo et al., 2016). In Dorning, 1997; Mun point reactor kinetics model, the reactor is assumed to be a point and spatial variation is not considered i.e., it is characterized by a time in dependent spatial shape function and a time dependent amplitude fac tor. When the flux shape undergoes significant change with time, this model is unable to capture the changes. However, for studying the out of phase oscillations, thermal hydraulics coupled with space time depen dent neutronics are used for modelling the system dynamics. Different methods for solution of space time dependent neutron kinetic equations have been used like Improved Quasi-Static (IQS) (Dutta and Doshi, 2009), modal kinetics (Dokhane, 2004; Mu~ noz-Cobo et al., 2016; Dokhane et al., 2007a, 2014; Lange et al., 2011), multi region neutron
* Corresponding author. E-mail address: [email protected] (S. Singh). https://doi.org/10.1016/j.pnucene.2019.103218 Received 1 June 2019; Received in revised form 23 October 2019; Accepted 29 November 2019 Available online 23 December 2019 0149-1970/© 2019 Elsevier Ltd. All rights reserved.
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Analytical Bifurcation method (Hassard et al., 1981) and numerical continuation method (Kuznetsov, 2004) are most commonly used for carrying out non-linear stability analysis. Semi analytical Bifurcation method is used in (Dokhane, 2004; Rizwan-Uddin, 2006; KarveR izwan-uddin and Dorning, 1997; Dokhane et al., 2014; Chakraborty et al., 2019) and continuation method is used in (Pandey and Singh, 2016; Chakraborty et al., 2018). The reduced order models as developed by (KarveRizwan-uddin and ~ oz-Cobo et al., 1996), based Dorning, 1997; Dokhane et al., 2007a; Mun on modal kinetics coupled with thermal hydraulics equations, involve use of system codes for approximating the cross coupling reactivities for each mode due to Doppler and coolant void feedbacks. The equations represented even with two thermal hydraulic channel model becomes very complicated. The thermal hydraulics coupled with multi region kinetics also lead to differential equations which are very cumbersome for analysis. The model developed by March-Leuba et al. (1986) to describe the behaviour of BWRs is very simple as compared to the other models involving explicit thermal hydraulics. The neutron kinetics in this model is based on point reactor kinetics and certain assumptions were made to obtain a simple thermal hydraulic model. The coolant is assumed to be at saturation temperature when it enters the core and the complete recirculation loop is treated as a one path where fluid travels with constant mass flow rate. This model was successful in the predic tion of the total power oscillations in BWRs. Detailed bifurcation anal ysis using this model using point kinetics was carried out using MATCONT in (Pandey and Singh, 2019). In this paper, multipoint kinetics are used for modelling the out of phase oscillations in a BWR and the equations representing the thermal hydraulic feedbacks are based on the model given in (March-Leuba, 1984). Linear and non-linear stability analysis has been performed using numerical continuation code MATCONT (Dhooge et al., 2004) and the results have been compared with one point model. Even though this model is very simple, it captures the dynamics and the results are in line with those predicted using more detailed thermal hydraulic codes.
Nomenclature nðtÞ N0 cðtÞ TðtÞ
ρ Λ β λ
ρα D
τ
H v0 G0
α12
Aij, Aji Dj dij
νi
Vi
Excess neutron population normalized to steady state neutron population Steady state neutron population Excess delayed precursors population normalized to steady state neutron population Excess fuel temperature Reactivity feedback Prompt neutron generation time Delayed neutron fraction Decay constant of delayed neutron precursors Coolant void feedback reactivity Doppler coefficient of reactivity Residence time of bubbles in the fluid channel Height of the reactor Propagation velocity Coolant mass flux Coupling coefficient between region 1 and 2 Interface area between two adjacent regions “i" and “j" One group diffusion coefficient for region ‘j’ Distance between the centroids of the two regions “i" and “j" Neutron speed in region “i" Volume of the region “i"
2. Mathematical model of a Boiling Water Reactor 2.1. Space independent model
Fig. 1. Block diagram for BWR model (March-Leuba, 1984).
The dynamic processes in a BWR involves the production of heat i.e., a neutronics loop and its removal from the core by light water in two phase (as coolant), represented by a thermal hydraulics loop. The main feedbacks are due to the rate of formation of coolant void and change in fuel temperature. The different processes can be represented in a block diagram form and is given in Fig. 1. The space independent dynamics of a BWR can be modelled using a fifth order model consisting of coupled nonlinear ordinary differential equations as given in March-Leuba et al. (1986) and are given below:
kinetics coupled with thermal hydraulics (UEHIRO et al., 1996; Lee and Pan, 2005). In order to identify stable and unstable regimes of operation, stability analysis is required. A straightforward approach for carrying out sta bility analysis is by directly solving the neutronics coupled with thermal hydraulics equations using system codes (Dutta and Doshi, 2009; Lange ~ oz-Cobo et al., 1996, 2000). However, et al., 2011; Wulf et al., 1984; Mun for identifying stable and unstable regimes, the simulations need to be carried out at different parametric values which is very computationally expensive and tedious since the entire range of parameters have to be spanned. In this regard, reduced order models are developed which are simpler than the explicit models but are very useful in capturing the phenomena. These reduced order models can be used for carrying out stability analysis for identification of stable and unstable regimes. Standard techniques such as frequency domain approach (March-Leuba, 1984) and eigenvalue approach (Rizwan-Uddin, 2006) can be used for arriving at the linearly stable regimes. Since BWRs behave as non-linear systems, linear stability analysis is not sufficient for identifying stable regimes. Non-linear stability analysis is also required because a “line arly” stable operating point can be “non-linearly” unstable. This is due to the fact that linear stability assesses the system dynamics for small perturbations and thereby for large perturbations, the system behaviour may be different. Hence, the identification and characterisation of such phenomena and the existence of Generalized Hopf bifurcation is very important in the identification of the stability of the system. Semi
dnðtÞ ρðtÞ β ρðtÞ ¼ nðtÞ þ λcðtÞ þ dt Λ Λ
(1)
dcðtÞ β ¼ nðtÞ dt Λ
λcðtÞ
(2)
dTðtÞ ¼ a1 nðtÞ dt
a2 TðtÞ
(3)
d2 ρα ðtÞ dρ ðtÞ þ a3 α þ a4 ρα ðtÞ ¼ kTðtÞ dt2 dt where nðtÞ ¼ ðNðtÞ NO Þ=NO and cðtÞ ¼ ðCðtÞ The reactivity feedback can be expressed as ðtÞ ¼ ρα ðtÞ
DTðtÞ
(4) CO Þ=NO . (5)
Eq. (4) is a second order differential equation and can be converted into two first order equations by the following substitution, 2
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Table 1 Model parameter values for Vermont Yankee test 7N (March-Leuba, 1984). Parameter
Value
Units
k
0.0037
K 1s
D
0.0000252
K
0.0056
–
Λ
0.00004
s
λ
0.08
s
a1
25.04
Ks
a2
0.23
s
a3
2.25
s
6.82
s
(6)
Substituting Eq. (6) in Eq. (4), we get
2
dρt ðtÞ þ a3 ρt ðtÞ þ a4 ρα ðtÞ ¼ kTðtÞ dt
1
β
a4
dρα ðtÞ ¼ ρt ðtÞ dt
(7)
The parameter k is related to fuel temperature and void coefficient of reactivity (March-Leuba, 1984). The parameters a3 and a4 scale as 1=τ and 1=τ2 , where τ is the residence time of bubbles in the fluid channel (March-Leuba, 1984). The residence time, τ ¼ H=v0 , H being the height of the reactor and v0 is the propagation velocity.v0 is directly propor tional to the coolant mass flux G0 as given in (March-Leuba, 1984) and hence τ is inversely proportional to G0 .The values of the parameters can be found in Table 1.
1 1 1
1 1 2
Fig. 2. Regions considered for analysis.
Fig. 3. Sub Critical and Super critical Hopf Bifurcation. 3
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
2.2. Space - time dependent model – multipoint kinetics
dT2 ðtÞ ¼ a1 n2 ðtÞ dt
In this paper, multipoint kinetics technique is used for solution of space-time dependent neutron diffusion equation. In this method, the core is divided into finite number of volumes with point kinetic equa tions for each region representing the dynamics and a coupling term between the adjacent regions which accounts for the neutron leakage. Since it has been observed that the oscillations in BWRs take place in the azimuthal direction (Dokhane, 2004; Dutta and Doshi, 2009; Dhooge et al., 2004) the reactor core is divided into two regions for the present analysis as shown in Fig. 2. The subscript 1 indicates the left region and the subscript 2 indicates the right region of the core with α12 repre senting the coupling between them. A fifth order model given in March-Leuba et al. (1986) is then adopted for both halves of the core separately. The coupling coefficients, αij ; representing the coupling between the ith and jth regions are calculated using the formula (Chakraborty et al., 2019)
αij ¼
(8)
ρ2 ðtÞ ¼ ρα2 ðtÞ þ DT2 ðtÞ
kc h
(9)
dhðtÞ ¼ kf ðn1 þ n2 Þ dt
ρ1 ðtÞ ¼ ρα1 ðtÞ þ DT1 ðtÞ
n1 ðtÞ þ
α12 Λ
n2 ðtÞ þ λc1 ðtÞ þ
ρ1 ðtÞ Λ
λc1 ðtÞ
(13)
dT1 ðtÞ ¼ a1 n1 ðtÞ dt
a2 T1 ðtÞ
(14)
dρα1 ðtÞ ¼ ρt1 ðtÞ dt
(15)
dρt1 ðtÞ þ a3 ρt1 ðtÞ þ a4 ρα1 ðtÞ ¼ kT1 ðtÞ dt
(16)
dn2 ðtÞ ρ ðtÞ β ¼ 2 dt Λ
α12
ρ2 ðtÞ ¼ ρα2 ðtÞ þ DT2 ðtÞ dc2 ðtÞ β ¼ n2 ðtÞ dt Λ
λc2 ðtÞ
n2 ðtÞ þ kc h
α12 Λ
n1 ðtÞ þ λc2 ðtÞ þ
ρ2 ðtÞ Λ
dhðtÞ ¼ kf ðn1 þ n2 Þ dt
(23)
(24)
where x is n dimensional state and r is the bifurcation parameter. The equilibrium solution is described as follows Fðx0 ; rÞ ¼ 0
(11)
dc1 ðtÞ β ¼ n1 ðtÞ dt Λ
(22)
x_ ¼ Fðx; rÞ
(25)
The characteristic roots of the Jacobian matrix, Jðx0 ; rÞ for the above system help determine the stability of local fixed points. Conditions required to be satisfied for Hopf bifurcation are as follows.
(12)
kc h
dρt2 ðtÞ þ a3 ρt2 ðtÞ þ a4 ρα2 ðtÞ ¼ kT2 ðtÞ dt
Hopf bifurcation is a codimension-1 bifurcation where there is only one free parameter for a bifurcation to occur. Changes in the parameter value of the system results in the change of its eigenvalues. These ei genvalues help determine the linear stability of a system in steady state. Positive real part of any one of the eigenvalues makes the system un stable whereas, when real parts of all the eigenvalues are negative, the system is stable. Hopf bifurcation occurs at the point at which a pair of complex eigenvalues cross the imaginary axis. Linear stability of the system is lost at this critical parameter value and a family of limit cycle bifurcates at that point, which is a non-linear phenomenon. Consider an autonomous system of ordinary differential equations
(10)
α12
(21)
3.1. Hopf Bifurcation
BWR system can now be modelled with multipoint kinetics as follows dn1 ðtÞ ρ ðtÞ β ¼ 1 dt Λ
dρα2 ðtÞ ¼ ρt2 ðtÞ dt
Hopf proved the existence of stable and unstable periodic oscillations for a non-linear dynamic system. Hopf theorem, sometimes also called as the Poincar�e–Andronov–Hopf (PAH) bifurcation theorem guarantees the existence of periodic solutions under certain conditions (Hassard et al., 1981; Guckenheimer and Holmes, 1983). The BWR system modelled by March – Leuba et al. (March-Leuba et al., 1986) meets all the conditions required by Hopf theorem. Since only Hopf bifurcation is observed in this study, a brief review of it and its mathematical background is given in this section. A detailed study on it can be found in (Kuznetsov, 2004) for more insights.
For the estimation of αij , the different constants are required which depend on the core geometry and fissile inventory. So for a two-region core, the coupling coefficient is inversely proportional to the distance between the two regions (dij ). This implies that a smaller core will have higher value of coupling coefficient since it is inversely proportional to the distance between the regions and vice versa. It is assumed that the total power oscillations are stable and the model only studies the spatial oscillations. In the study presented by ~ oz-Cobo et al. (Dhooge et al., 2004) constant pressure boundary Mun conditions were enforced to observe out-of-phase oscillations. Since in this study a simpler March-Leuba model (March-Leuba et al., 1986) is adopted, a feedback term is introduced to reactivity equation to observe out-of-phase oscillations. The feedback term introduced is kc h as shown in Eq. (8) and Eq. (9) where h is defined as shown in Eq. (10). kc h
(20)
3. Mathematical background
Di νi ΛAij dij Vj
ρ1 ðtÞ ¼ ρα1 ðtÞ þ DT1 ðtÞ
a2 T2 ðtÞ
1) Eigenvalues of Jðx0 ; rc Þ has at least one pair of imaginary eigen values, � iω, lying on the imaginary axis. Here rc represents the parameter’s critical value where Hopf bifurcation occurs. 2) The eigenvalues of Jðx0 ; rc Þ are of the form σðrÞ � iω ðrÞ and dσ ðr ¼ rc Þ=dλ 6¼ 0, i.e., the real parts of the eigenvalues of Jðx0 ; rc Þ cross the imaginary axis with non-zero speed as r crosses rc . 3.2. Supercritical and subcritical Hopf bifurcation Hopf theorem establishes the coexistence of limit cycles and steady state solutions. Limit cycles are isolated closed trajectories where neighbouring trajectories spiral towards or away from the limit cycles. Limit cycles are classified and as stable and unstable limit cycles. In stable limit cycles, neighbouring trajectories are attracted towards it, whereas in unstable limit cycles the neighbouring trajectories are repelled. Hopf bifurcation can again be classified as supercritical or subcritical Hopf bifurcation as shown in Fig. 3. In a subcritical bifur cation (Fig. 3a) unstable limit cycle bifurcates from a stable fixed point.
(17) (18) (19)
4
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
4. Results and discussions
Table 2 Additional parameter values for multipoint model. Parameter
Value
kc
0.01
kf
50
α12
0.0035
4.1. Model verification 4.1.1. Stability maps Bifurcation curves are generated using MATCONT (Dhooge et al., 2004) software. Stability boundaries for the BWR system described by multipoint model represented by Eq. (11) - Eq. (23) are identified. The parameters k and a3 are chosen as the bifurcation parameters, similar to the work presented in (Pandey and Singh, 2016). These two parameters are varied keeping the other parameters from Tables 1 and 2 constant. The bifurcation plot for k vs a3 given in Fig. 4 represents the region where the system changes its stability. As the parameter k is changed, the real part of eigenvalue changes its sign thereby implying Hopf bifurcation. Apart from identifying the Hopf points, first Lyapunov co efficients at each point are estimated using MATCONT which indicates the nature of the Hopf Bifurcation. Negative first Lyapunov coefficient indicates a supercritical Hopf bifurcation whereas subcritical Hopf bifurcation is indicated by positive first Lyapunov coefficient. Fig. 5 shows that the first Lyapunov coefficient is negative in the plotted region implying the nature of Hopf bifurcation is supercritical. The study presented in Lee and Pan (2005) concludes that for strong
Any perturbation with amplitude less than that of the limit cycle causes neighbouring trajectories to be attracted to the stable fixed point. Any perturbation with amplitude higher than that of limit cycle also repels neighbouring trajectories from the unstable limit cycle. Hence, in subcritical condition after bifurcation the trajectories must to jump to a distant attractor. In a supercritical Hopf (Fig. 3b) bifurcation system the physical system initially settles down to equilibrium through damped oscillations. As the parameter value changes, the decay decreases and slowly loses its stability resulting in small amplitude, sinusoidal limit cycles oscillation about the former steady state. Hence, in supercritical Hopf bifurcation a stable limit cycles bifurcates from an unstable fixed point. Any trajectory with lesser or higher amplitude than that of the limit cycle is attracted towards it.
Fig. 4. Bifurcation plot for k vs a3 at α12 ¼ 0.0035.
Fig. 5. First Lyapunov coefficient for k vs a3 bifurcation plot at α12 ¼ 0.0035. 5
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 6. Bifurcation plot for k vs α12
Fig. 7. First Lyapunov coefficient for k vs α12 bifurcation plot.
Fig. 8. Bifurcation plot of k vs a3 at α12 ¼ 0.0063.
6
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 9. First Lyapunov coefficient for k vs a3 bifurcation plot at α12 ¼ 0.0063.
Fig. 10. Bifurcation plot for k vs a3 for different values of α12
illustrated by values of first Lyapunov coefficient shown in Fig. 7. With a higher value of coupling coefficient α12 in the multipoint kinetics model, both supercritical and subcritical Hopf bifurcation are observed. Fig. 8 shows the occurrence of supercritical and subcritical Hopf bifurcation when α12 ¼ 0.0063. First Lyapunov coefficient for this bifurcation plot is shown in Fig. 9. Fig. 10 compares the bifurcation plot for k vs a3 for different values of coupling coefficient α12 . For lower values of α12 , only supercritical Hopf bifurcation is observed. By increasing the value of the coupling coeffi cient, both supercritical and subcritical Hopf bifurcation are observed, thus confirming the results presented in Lee and Pan (2005) The stability map can be divided into four regions I, II, III and IV. Since the coupling coefficient is inversely proportional to the size of the reactor core, the lowest value corresponds to the highest value of coupling coefficient and vice versa. So, it can be seen from the figure that for α12 ¼ 0.0010, the regions I, II and III are unstable. However for α12 ¼ 0.0063, which represents a smaller core size, only Region I is unstable. Hence, higher the value of the coupling coefficient α12 , larger is the parametric space where the system is stable.
Table 3 Parameter values at various points for numerical simulations. Point number
a3
a4
k
kc
kf
α12
1
3.1
6.82
1 � 10
2
2.1
6.82
8:1 � 10
3
0.01
1
0.0035
0.01
1
3
3.24
1.705
4:3 � 10
0.0035
3
0.01
1
4
3.24
1.705
4:3 � 10
0.0035
3
0.01
50
0.0035
2
Table 4 Updated values of a3 and a4 for multipoint model. Parameter
Value
Units
a3;m ¼ (a3;one )/2
1.125
s
a4;m ¼ (a4;one )/4
1.705
s
1 2
neutron interaction, which corresponds to a higher value of coupling coefficient α12 , both supercritical and subcritical Hopf bifurcations are observed. Similar results are obtained observed in this study by adopting the multipoint kinetics model. Bifurcation plot of k vs α12 in Fig. 6 shows the occurrence of Generalized Hopf (GH) point. GH points indicates the change in nature of Hopf bifurcation from supercritical to subcritical as
4.1.2. Numerical simulations The dynamics of n1 and n2 are obtained by solving Eq. (11) - Eq. (23). It can be seen that the fixed point comes out to be [n1 ðtÞ, c1 ðtÞ, T1 ðtÞ, ρα1 ðtÞ, ρt1 ðtÞ, n2 ðtÞ, c2 ðtÞ, T2 ðtÞ, ρα2 ðtÞ, ρt2 ðtÞ, h] ¼ ½0; 0; 0; 0; 0; 0; 0; 0; 0; 7
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 11. Trajectories of n1 and n2 operating at point number 1.
Fig. 12. Trajectories of n1 and n2 operating at point number 2.
Fig. 13. Bifurcation plot k vs a3 for multipoint kinetics.
8
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 14. First Lyapunov coefficient vs k for multipoint model.
Fig. 15. Bifurcation plot k vs a3 for one point kinetics.
Fig. 16. First Lyapunov coefficient vs k for one point model. 9
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 17. Comparison of stability boundaries for one and multipoint kinetics model.
Fig. 18. Trajectory of total n at point number 3 using one point model.
0; 0�. The system is then perturbed such that the initial point at t ¼ 0 is ½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0�. Numerical simulations are carried out at points whose parameter values are given in Tables 3 and 4. The parameters k and a3 are chosen close to the stability boundary shown in Fig. 4. Point number 1 from Table 3 is in the stable region of bifurcation plot. When the system is given a perturbation, trajectories of n1 and n2 decay to return to its steady state. This is illustrated in Fig. 11. Point number 2 from Table 3 is in unstable region with supercritical Hopf bifurcation. The perturbation applied to the system causes the neutron population to undergo a limit cycle with periodic oscillations as shown in Fig. 12. The excess neutron population n1 and n2 undergoes about 8–9 cycles of oscillations every 20 s. This is similar to the fre ~ oz-Cobo et al. (Dhooge et al., 2004) with 9–10 quencies reported in Mun cycles of oscillations every 20 s. Hence with appropriate parameter values, multipoint kinetics based approach can be used to model the BWR system.
model. With the mass flux being the half compared to that of one point model, the propagation velocity V0 becomes half (March-Leuba, 1984) and thereby the residence time of the bubbles in the channel τ doubles. Hence in the multipoint model, the parameter a3;m is half of a3;one and a4;m is one fourth to that of one point model, a4;one (as given in Table 4). The other parameters k, a1 and a2 which are defined by March-Leuba (1984) remain unchanged when the core is divided into two halves. 4.2.1. Stability maps for one point and multipoint kinetics model The stability boundaries for the BWR system modelled with multi point kinetics are again generated using MATCONT with new values of a3 and a4 to compare the boundaries with one point kinetics model. Parameters k and a3 are chosen as the bifurcation parameters. The bifurcation plot for k vs a3 is given in Fig. 13. The first Lyapunov coef ficient is illustrated in Fig. 14 which indicates that the nature of Hopf bifurcation is supercritical. Similar bifurcation curve can be generated for BWR system with one point kinetics model which is represented by Eq. (1)–(7). Parameter values for one point model are obtained from Table 1. The stability boundary is illustrated in Fig. 15 again by varying k and a3 and the first Lyapunov coefficient plot is presented in Fig. 16. Hence the one point kinetics model suggests the system undergoes subcritical Hopf bifurcation along with supercritical Hopf bifurcation.
4.2. Comparison of one point and multipoint kinetics model Since the core is divided into two equal halves, the coolant mass flux in multipoint model is half of one point model. Hence the parameters a3 and a4 will be different for multipoint model as compared to one point 10
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 19. Trajectories of n1 and n2 at point number 3 using multipoint model.
Fig. 20. Trajectory of n1 and n2 at point number 3 with low perturbation.
Fig. 21. Phase portrait of Limit cycles. 11
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 22. Trajectory of total n at point number 3.
Fig. 23. Bifurcation plot k vs a3 for multipoint kinetics at kf ¼ 1 and kf ¼ 50
Fig. 24. Bifurcation plot for k vs kf indicating Hopf point occurs at same value of k
12
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 25. Trajectory of total n at point number 4 where kf ¼ 50.
Fig. 26. Trajectories of n1 and n2 for different kf.
13
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
Fig. 27. Comparison of n1 vs n2 limit cycle for different kf.
The multipoint model indicates that only supercritical Hopf bifurcation exists. While in previous studies, Dokhane et al. (2007a) pointed out existence of subcritical bifurcation alone, in more recent tests as pre sented in Dokhane et al. (2014), the nature of limit cycles observed is supercritical only. However, parameters such as kc or α12 in multipoint model can be adjusted to obtain GH points as well. Fig. 17 shows the comparison of stability boundaries for one and multipoint kinetics model. As seen in Fig. 17, the stability boundary for two point model is right to the stability boundary of one point kinetics model. In region I, both one point and multipoint models are unstable. In region II, one point model is stable whereas multipoint model in unstable and in region III, both one point and multipoint models are stable. Hence in region II, which is between the two boundaries, the total power oscillations are stable where as the spatial oscillations are unstable which is consistent with the assumption as given in Section 2.2.
low perturbation results in a stable point and high perturbation results in jumping to a distant limit cycle. Fig. 20 shows the limit cycles along which n1 and n2 oscillate in the time domain and the phase portrait of the limit cycle in n1 vs T1 and n1 vs n2 plane is given in Fig. 21. 4.2.3. Effect of kf on stability boundary Fig. 22 shows the dynamics total n obtained by adding n1 and n2 obtained from multipoint model as opposed to directly simulating one point model shown in Fig. 18. Fig. 22 indicates that though the total power is being supressed, it is still not constant. Hence values of kf are varied to get desired results. It is observed that the values kf does not change the stability boundary. This is illustrated in Fig. 23 and Fig. 24. Fig. 23 indicates that bifurcation plot for k vs a3 is same for different values of kf whereas, Fig. 24 shows that the bifurcation plot does not change with kf for a given value of k. The total n trajectory at kf ¼ 50 and other parameter values from point number 4 in Table 3 for the same initial perturbation t ¼ 0, ½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0� is shown in Fig. 25. This normalised total excess neutron density is oscillating at zero with negligible amplitude. Changing the value of kf does not change the out of phase behaviour of n1 and n2 as shown in Fig. 26 and Fig. 27. Increasing the value of kf , n1 and n2 become more anti symmetrical as seen in Fig. 27, thereby cancelling each other out. This results in regional out of phase oscillations where the total oscillations remains constant. Hence with a suitable values of the coupling coefficient, α12 and feedback coefficients, kc and kf in Eq. (13) - Eq. (25), regional out of phase oscillations can be observed.
4.2.2. Numerical simulations for one point and multipoint kinetics model The fixed point for multipoint model is [n1 ðtÞ, c1 ðtÞ, T1 ðtÞ, ρα1 ðtÞ, ρt1 ðtÞ, n2 ðtÞ, c2 ðtÞ, T2 ðtÞ, ρα2 ðtÞ, ρt2 ðtÞ, h] ¼ ½0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0�. The system is then perturbed such that the initial point at t ¼ 0 is ½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0�. Similarly for one point model, from Eqs. (1)– (7), it can be observed that the fixed point comes out to be [nðtÞ, cðtÞ, TðtÞ, ρα ðtÞ, ρt ðtÞ] ¼ ½0; 0; 0; 0; 0�. The dynamics of n is then observed by creating a perturbation such that the initial point at t ¼ 0 is ½1; 1 ;10 ;0; 0�. One point model represents the dynamics of total n and multipoint model represents the dynamics of n1 and n2 . To verify the stability map shown in Fig. 18, simulations are carried out at point in between the two stability boundaries which is presented as point number 3 in Table 3. Fig. 18 is obtained by perturbing the one point model. The trajectory of n eventually reaches to its steady state value. Fig. 19 illustrates the dynamics of n1 and n2 are obtained by solving the multipoint model when the perturbation is applied with parameter values given at point number 3 in Table 3. It can be observed that n1 and n2 are out-of-phase and oscillate in a limit cycle. Even with low perturbation such as the initial point at t ¼ 0 is 0:001*½1; 1; 10; 0; 0; 0; 0; 0; 0; 0; 0�, the trajectories of n1 and n2 grow and reaches the same limit cycle. This supports the claim that the Hopf bifurcation is supercritical where high and low perturbation results in the same limit cycle as opposed to subcritical Hopf bifurcation where
5. Conclusions A simple model for studying out of phase oscillations in a BWR is developed based on multipoint kinetics. The core is divided into two parts which is coupled together by the coupling coefficient α12 . Two thermal hydraulic channels, one representing each region is modelled. The feedback reactivities due to coolant void and fuel temperature feedback are modelled using equations similar to March-Leuba (1989). Generally, constant pressure boundary conditions (Mu~ noz-Cobo et al., 1996) or constant total mass flow rate conditions (Lee and Pan, 2005) are enforced for controlling in-phase oscillations. Since the model pre sented in this paper does not explicitly model the thermal hydraulics, a 14
R. Vora et al.
Progress in Nuclear Energy 120 (2020) 103218
feedback reactivity is used to control the total power. ~ oz-Cobo et al. The results of the model are compared with Mun (Dhooge et al., 2004) and are found to be consistent. However, it is to be noted that the parameters used in the analysis are dependent on the system parameters and have to be derived for each BWR as given in March-Leuba (1989). But once these constants are derived, then this simple model becomes sufficient for studying these power oscillations. The developed model is further used to perform bifurcation analysis using MATCONT. The bifurcation analysis shows that the BWR system undergoes Hopf bifurcation as the parameters are varied. Both super critical and subcritical Hopf bifurcation were also observed in different parametric planes. Limit cycles have been observed by plotting the trajectory of n1 and n2 , the excess neutron populations normalized to the steady state neutron population. The stability boundaries of one point and two - point kinetics model are also compared. The stability boundary for one point model represents the stability characteristics of global oscillations whereas that for the two point model gives the same for out of phase oscillations. It is found that there exists a region where the global oscillations are stable but local oscillations are unstable.
Dokhane, A., Ferroukhi, H., Pautz, A., 2014. On out-of-phase higher mode oscillations with rotation and oscillation. Ann. Nucl. Energy 67, 21–30. Dutta, G., Doshi, J.B., 2009. Nonlinear analysis of nuclear coupled density wave instability in time domain for a boiling water reactor core undergoing core-wide and regional modes of oscillations. Prog. Nucl. Energy 51 (8), 769–787. Guckenheimer, J., Holmes, P., 1983. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer-Verlag. Hassard, B., Kazarinoff, N., Wan, Y., 1981. Theory and Applications of Hopf Bifurcation. Cambridge University Press, New York. Karve, A.A., Rizwan-uddin, Dorning, J., 1997. Stability analysis of BWR nuclear-coupled thermal-hydraulics using a simple model. Nucl. Eng. Des. 177, 155–177. Kuznetsov, Y., 2004. Elements of applied bifurcation theory. In: Applied Mathematical Sciences, vol. 112. Springer. Lange, C., Hennig, D., Chair, A.H., 2011. An advanced reduced order model for BWR stability analysis. Progress Nuclear Energy J. 53 (1), 139–160. Lee, J.D., Pan, C., 2005. Dynamic analysis of multiple nuclear-coupled boiling channels based on a multi-point reactor model. Nucl. Eng. Des. 235. March-Leuba, J., 1984. Dynamic Behaviour of Boiling Water Reactors. PhD Thesis. The University of Tennessee, Knoxville. March-Leuba, J., 1989. Non linear dynamic and chaos in boiling water reactors. In: Noise and Nonlinear Phenomena in Nuclear Systems, vol. 192. March-Leuba, J., Cacuci, D., Perez, R., 1986. Nonlinear dynamics and stability of boiling water reactors: part 1 –qualitative analysis. Nucl. Sci. Eng. 23, 93–111. Mu~ noz-Cobo, J.L., Pdrezi, R.B., Ginestar, D., Escriv, A., Verdu, G., 1996. Non linear analysis of out of phase oscillations in boiling water reactors. Ann. Nucl. Energy 23 (16), 1301–1335. Mu~ noz-Cobo, J.L., Rosello, O., Miro, R., Escriva, A., Ginestar, D., Verdu, G., 2000. Coupling of density wave oscillations in parallel channels with high order modal kinetics: application to BWR out of phase oscillations. Ann. Nucl. Energy 27, 1345–1371. Mu~ noz-Cobo, J.L., Escriv� a, A., Domingo, M.-D., 2016. In-phase instabilities in BWR with sub-cooled boiling, direct heating, and spacers effects. Ann. Nucl. Energy 87, 671–686. Pandey, V., Singh, S., 2016. A bifurcation analysis of boiling water reactor on large domain of parametric spaces. Commun. Nonlinear Sci. Numer. Simul. 38, 30–44. Pandey, V., Singh, S., 2019. Bifurcations emerging from a double Hopf bifurcation for a BWR. Prog. Nucl. Energy 117. Rizwan-Uddin, 2006. Turning points and sub- and supercritical bifurcations in a simple BWR model. Nucl. Eng. Des. 236 (3), 267–283. Uehiro, M., Rao, Y.F., Fukuda, K., 1996. Linear stability analysis on instabilities of inphase and out-of-phase modes in boiling water reactors. J. Nucl. Sci. Technol. 33, 628–635. Waranpera, Y., Anderson, S., 1981. BWR stability tests reaching limit cycle threshold at natural circulation. Trans. Am. Nucl. Soc. 39, 868. Wulf, W., Cheng, H., Diamond, D., Khatib-Rahbar, M., 1984. Description and Assessment of RAMONA-3B Mod. 0 Cycle 4: a Computer Code with Three-Dimensional Neutron Kinetics for BWR System Transients. NUREG/CR-3664, United States.
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.pnucene.2019.103218. References Chakraborty, A., Singh, S., Fernando, M.P.S., 2018. A novel approach for bifurcation analysis of out of phase xenon oscillations using multipoint reactor kinetics. Nucl. Eng. Des. 328, 333–344. Chakraborty, A., Singh, S., Fernando, M.P.S., 2019. An improved reduced order model for nonlinear stability analysis of spatial xenon oscillations. Prog. Nucl. Energy 116, 62–75. Dhooge, A., Govaerts, W., Kuznetsov, Y.A., 2004. MATCONT: a matlab package for numerical bifurcation analysis of ODEs. SIGSAM Bull. 38 (1), 21–22. Dokhane, A., 2004. BWR Stability and Bifurcation Analysis Using a Novel Reduced Order � Model and the System Code RAMONA. PhD Thesis. Ecole Polytechnique F�ed� erale de Lausanne, Switzerland. Dokhane, A., Hennig, D., Rizwan-uddin, Chawla, R., 2007. BWR stability and bifurcation analysis using reduced order models and system codes: identification of a subcritical Hopf bifurcation using RAMONA. Ann. Nucl. Energy 34, 792–802.
15