On butterfly effect in the dynamics of Oskarshamn-2 instability event of 1999

The current research is a further investigation on the OECD/NEA benchmark on the Oskarshamn-2 stability event of cycle 24, occurred on 25th February 1999, in order to study the effect of minor core design changes on both steady state and dynamic behavior of the reactor core, using SIMULATE-3K. The analysis of steady state and transient have been performed along with a systematic comparison between results of the reference and modified core. Results are presented in terms of different parameters such as power/flow ratio, 2D core maps for power and flow for steady state; and channels ’ power time series, 2D core maps for maximum power, decay ratio, and resonance frequency, for the transient. The results show that, the introduction of minor changes in the core through the reshuffling of the two hottest channels, and even less limiting channels, has small effects on the steady state core behavior. However, a significant change in the core dynamic behavior is observed through a drastic reduction of the power oscillation amplitude and the simulation demonstrates that, if such minor core design changes had been introduced in the beginning of cycle 24, the stability event of February 25, 1999 might have been potentially avoided. This drastic change in dynamical behavior is explained by carrying out a detailed stability and bifurcation analysis, which concluded that the reason behind such significant reduction in the power oscillation amplitudes is due to a shift in the stability boundary, making the core more stable and also due to the fact that the event took place very close to the stability boundary, region in which the system, unexpectedly, found very sensitive to small changes in the unstable region where the supercritical Hopf bifurcation is expected.


Introduction
Boiling Water Reactors (BWRs) are highly nonlinear complex systems that can become vulnerable to instabilities, under specific conditions, due to the strong coupling between the neutronic and thermal hydraulic components, via the void and Doppler reactivity feedbacks.Even though extensive amounts of studies have been carried out in the field of stability analysis of BWRs, many open questions of fundamental nature are still open and need further efforts using more advanced mathematical and simulations tools.
In recent years, an OECD/NEA Benchmark has been launched in order to analyze a real plant transient (stability event) undergone in the Swedish plant, Oskarshamn-2, on 25th February 1999, during Cycle 24 (Kozlowski, et al., 1999).Among other institutions and universities (Gajev et al., 2013;Kozlowski, 2014;Balestra et al., 2013), the Paul Scherrer Institute (PSI) was an active participant and in this context efforts have been dedicated to study this interesting and challenging stability event (Dokhane et al., 2016;Dokhane et al., 2017).In that context, a stability methodology has been developed and validated, at PSI, based on the three-dimensional reactor kinetics state-of-the-art code SIMULATE-3K (S3K) (Dokhane et al., 2016).Results of these studies demonstrated the capability of S3K not only to simulate such complex behavior of O2 stability event but also to investigate possible scenarios that could have been occurred whether with or without the intervention of the operator.As an example, the S3K results indicated that if the reactor scram had not been introduced, the observed high oscillation amplitudes would have decayed and hence the reactor would have gone back to a stable state.
Within a recent study performed, using TRACE/PARCS, it has been reported that when the O2 core design is slightly modified, the original high amplitude oscillations of the power damp drastically although the resonance frequency value remains identical (Walser, 2016).However, no explanation was given to the observed phenomenon, which could be unphysical and simply a result of numerical effects.The findings of these previous studies suggest that, this specific stability event is rich of complex and interesting behaviours of high interest that needs further and deeper analysis, which should provide deeper insights into the underlying physical mechanisms.
Therefore, the main goal of the research is twofold.First, to carry out an independent study, using SIMULATE-3K being one of the most reliable codes for stability analysis, in order to verify that the behavior observed in (Walser, 2016), i.e. high sensitivity of the reactor dynamics to minor core design changes, is a physical phenomenon and not simply due to numerical diffusion.Second, more importantly, to carry out a detailed and an in-depth stability and bifurcation analysis in order to explain the reasons behind such drastic change of the core dynamics.
The outcome of such study could help core designers to better optimize the core loading pattern, which could avoid or mitigate such instabilities and therefore could improve the safety of BWRs.

Overview on Oskarshamn-2 25 February 1999 event
The Oskarshamn-2 reactor, located in Oskarshamn, Sweden, and designed by ABB ATOM, is an external-loop type BWR comprising four recirculation loops and a core with 444 fuel assemblies, belonging to four different bundle designs (Kozlowski, et al., 1999).It began operation in 1975 with a 1700 MW thermal power, then was uprated to 1800 MW in the early 1980 (Kozlowski, et al., 1999).
During a maintenance work on the switchyard outside of the O-2 reactor, On 25th of February 1999, a load rejection signal was transmitted to the turbine and triggered a turbine trip, which was not transferred to the reactor and therefore kept on running at nominal power.As a result, the feedwater preheaters were no longer operational and the feedwater temperature started decreasing, which leaded to the entrance of colder water in the reactor vessel and created a positive reactivity feedback and hence an increase in the core power level.This leaded the automatic system to reduce the main recirculation flow, which consequently reduced the power.But due to the continuous decrease in the feedwater temperature, this sequence was repeated twice before the operator partially inserted control rods (see Fig. 1).However, the flow of the colder feedwater continued, causing the reactor power to increase rapidly in an oscillatory manner and then the reactor scrammed due to high power, at about 132 % (see Fig. 1).More details about the event can be found in (Kozlowski, et al., 1999).

PSI stability analysis methodology
At PSI, a stability analysis methodology has been developed based on the S3K code (Dokhane et al., 2013), which is a best-estimate LWR core dynamics code consisting of a full core coupling between neutronics and thermal hydraulics.The neutronics is based on three-dimensional two-Fig.1.Power Oscillation of Oskarshamn-2 Event (Kozlowski, et al., 1999).group neutron diffusion solver using a transverse-integrated nodal method for the spatial integration and the frequency transform method for the time integration during the transient.The thermal hydraulic model is based on a six-equation solver: mass, energy and momentum conservation equations for both liquid and vapour.It should be noted that, instead of solving the separate phasic momentum conservation equations, a mixture momentum equation and a weighted difference of the liquid and vapour momentum equations.For more details about S3K models, see (Grandi and Belblidia, 2011).
The PSI stability analysis methodology has the following particular aspects.First, the S3K core model is directly fetched from upstream validated reference models of the lattice code CASMO and the static solver SIMULATE-3 of the Swiss reactors, allowing a full consistency in terms of neutronic methods (Ferroukhi et al., 2008).Second, the core model is complemented by a hydraulic model for the vessel, including all recirculation loop components.As a result, the main principle of the PSI methodology is that the core and vessel models along with the associated physical models and code options are employed in a systematic and consistent manner for all operating conditions and all cycles, covering all types of oscillation modes (e.g.(Dokhane et al., 2013;Dokhane et al., 2014).
For the current study, a similar methodology is applied for the Oskarshamn-2 core, as shown on the Fig. 2.However, because the few group homogenized nuclear data were provided as part of the specifications, the only difference is that no CASMO models were developed.On that basis, the CMSLINK and SIMULATE-3 models were developed, using for the latter, the nodal history distributions, also provided as part of the specifications.

Analyses and results
This section presents, first, an overview about the modification introduced to the original core design in order to study the associated impact on the core behavior.After that, since most of the results of the reference core design have been presented in previous works (Dokhane et al., 2016), the focus here will be mainly on a systematic and a detailed comparison of steady state and transient results between the modified and the original (reference) core models.
It should be noted that, in the framework of the current study, a Matlab platform has been developed in order to carry out a detailed, systematic and automated analysis.More specifically, after extraction of relevant data from the S3K output file of the corresponding case, the Matlab programs produce the steady state and the transient results.The steady-state results are presented in terms of 2D core maps of core power, void, and flow distributions, while transient results are illustrated in terms of time series of core power and power of any selected FA, the 2D core maps of the decay ratio (DR), resonance frequency (RF), maximum oscillation amplitude and oscillation phase-shift between assemblies.

Core design modification
The analysis of the original core design using S3K showed that, the fuel assemblies 146 and 299, located in (8,12) and (15,11) in the core map, see Fig. 3, are the hottest channels in the core at nominal core conditions, before the transient.In addition, the same two channels have also the highest void fraction and the lowest mass flow, which results in  highest power/flow ratio (See Table 1 and left side of Fig. 5).This means that these two channels should be the most vulnerable to instabilities and consequently have more impact on the global stability of the core compared to the remaining channels, as will be shown in the transient analysis, in Section 4.3.On this basis, the core design modification is performed in such a way to exchange the two hottest fuel assemblies with two other FAs close to the core periphery, while keeping the core symmetry, as illustrated in Fig. 3.This is done by exchanging the hottest FAs 146 and 299 with those located in 44 and 401, respectively, in a similar way done in the study performed in (Walser, 2016), allowing the verification of the results of this latter.It should be noted that, in Table 1, the ranking of each channel, involved in the swap, is based on the power/flow ratio value.For instance, the hottest channels 146 and 299 are ranked 2 and 1, respectively, among the 444 core channels, while channels 44 and 401 are ranked 291 and 292, respectively.In addition, the table shows also the burnup values of the two involved pairs of fuel assemblies, where, as expected, it is observed that channels 146 and 299 are less burned than channels 44 and 401.
The modified core design is modeled by adapting the SIMULATE-3 core model, according to the above fuel assemblies swap, while the SIMULATE-3K model does not need to be modified since, as explained in Sec. 3, the S3K core model is directly transferred from SIMULATE-3 through the restart file.
In order to carry out the systematic comparison of the O-2 core, for both steady state and transient, between the original and modified (called also reshuffled or Swap 1) core, the simulation of the complete O-2 stability event for the modified core is carried out under identical operating and boundary conditions as for the reference model.

Steady-state analysis
In this section, the S3K steady state results of the modified core model are presented along with those of the original core model in Table 2 and Figs. 4 and 5, in order to perform a systematic comparison  A. Dokhane et al. and to evaluate the effect of the small core design change due to the reshuffling of only 4 out of 444 fuel assemblies.As can be seen in Table 2, summarizing the main steady state parameters, the introduction of the small changes in the core design has negligible effects on most parameters, with a very slight reduction in the keff of aboub 60 pcm.As expected, the swap of the two hottest channels from the central to the  peripheral part of the core has an effect on the radial peaking factor with a slight reduction, i.e. from 1.71 to 1.63.In Fig. 4, representing the 2D core maps of the radial normalized power and active mass flow, it is shown that the effect of the core design change has again a minor effect on the radial power and flow values and this minor impact is limited to the vicinity of the locations of the swapped 4 fuel assemblies, where the originally hottest fuel assemblies with the lowest mass flow lose these properties once they are moved close to the periphery.
Fig. 5 presents the power/flow ratio distribution of both core design, where it is shown that for the original core the smooth distribution is broken by the two hottest channels, which are distinguished by the largest power/flow ratio among the other channels, while for the modified core the distribution is clearly smooth.

Transient analysis
The S3K simulation of the O2 1999 stability event is presented in this section for both reference and reshuffled cores with a systematic comparison in Figs. 6, 7, 8 and 9.As stated above, the main focus here will be on the comparison between the results of the two cores since the reference core has been extensively studied in an earlier work (Dokhane et al., 2016).
Fig. 6 presents the power time-series of selected fuel assemblies, including the four swapped ones ((146,299) vs (44.401)), during the O2 event for both cores.As can be seen, the dynamics of each of the swapped channels changes significantly despite the minor change introduced in the core.For instance, for the reference core, the power oscillations in the hottest channels (146 and 299) grow rapidly and reach a maximum amplitude of almost 180 % of the respective steady state value.Then, due to the event BCs (continuous decrease of feedwater flow and temperature) these oscillations decay and reach a stable state.However, the same channels when swapped show a drastic reduction in the oscillation amplitudes with a very slow growing rate with a power level that remains below the steady state value.More interesting, even the channels that were not involved in the reshuffling, e.g. 1 and 444, have undergone the same drastic change in their dynamics, i.e. significant reduction in oscillation amplitude.
This drastic change in the dynamics, observed for the selected channels above, is found to be a global core trend, when scanning the status of each channel, as illustrated in Fig. 7.The figure represents the 2D core map of maximum power amplitude and shows that, most of the channels in the original core reach very large oscillation amplitudes, with a maximum value about 180 % of its respective steady-state value, corresponding to the hottest channels, as expected.This is a clear indication of a strong instability at most channels in the core.However, for the slightly modified core, the complete picture changes drastically with very small oscillation amplitudes and maximum power levels during the instability that do not reach even the reference steady state value.In addition, as can be seen, the high amplitude region shifts from the central zone (original core) to the four outer corner zones of the core (reshuffled core).The results in terms of the stability parameters, i.e. the decay ratio and resonance frequency, are presented in Figs. 8 and 9.It should be noted that, the evaluation of DR and RF is performed on the power time series of each assembly in the interval of the oscillation amplitude growing, i.e. 235-255 s, using the PSI signal analysis methodology based on ARMA modeling approach (Dokhane et al., 2006).The current results confirm, again, the strong impact of the minor changes, introduced in the core.As can be seen, while the DR at every FA is in the range 1.35-1.38 for the original core, it is significantly reduced, about 1.16, for the FAs of reshuffled core.However, the resonance frequency seems robust to this minor change.
Finally, the overall effect during the stability event on the total reactor power is presented in Fig. 10, along with the measured power.As can be observed, while during the first transient period, i.e. until about 230 s, the behavior of both cores is exactly the same, beyond 230 s the dynamics of the two cores start to diverge in an exponential manner despite the minor core design differences, exactly similar to the butterfly effect, known in Chaos theory, referring to the sensitive dependence on initial conditions in a nonlinear dynamical system (Gleick, 2011).More specifically, for the reference core design, the total power oscillations grow in the first transient period and reach a maximum amplitude of about 170 % then damp until reaching a stable state.As explained in (Dokhane et al., 2016), this is due to a double crossing of the stability boundary, i.e. from stable to unstable then again to a stable region.However, for the reshuffled core, the total power oscillations grow also during the first transient period but very slowly and while the amplitude is still very small, the oscillations start to damp and reach finally a stable state.Here again, similar to the reference core, it is clear that the core state has crossed the stability boundary twice, i.e. stable-unstable-stable regions, however in the current case, the amount of time spent inside the unstable region seems to be shorter, which does not allow to the oscillation amplitude to reach very large values.
The main conclusion here is that although the minor changes introduced in the O2 core have very small effect on the steady state behavior, their effect on the dynamic behavior is very significant.This means that, if such minor core design changes had been introduced in the beginning of Cycle 24, the stability event of February 25, 1999 might potentially have been avoided.Of course, this statement is valid assuming that the reshuffled core would not deteriorate the steady-state core performance along the cycle depletion including, e.g.Minimum Critical Power Ratio (MCPR), something that can evidently not be verified here.
It should be emphasized that, the complete analysis was focused on the swap of FAs No. 146 and 299 with FAs No. 44 and 401, respectively.As stated in section 4.1, the choice of this swap is based on the swap performed in (Walser, 2016), which allows a systematic verification and comparison of the results.However, this is not the only swap that show the high sensitivity of the core stability behavior to small core design.Within an additional investigated, it has been found that there are many other swaps resulting in the same effects.In addition, it is also found that it is not necessary that the swap of only hottest channels is the reason behind the drastic change in the dynamics, as illustrated for Swaps 2 and 3, in Fig. 10.Where in Swap 2, FAs 147 and 298, ranked 26 and 25, respectively, based on power/flow ratio (see Table 4, in Appendix 1), are exchanged with FAs 64 and 381, respectively; while in Swap 3, FAs 188 and 257, ranked 8 and 7, respectively, based on power/flow ratio (see Table 4, Appendix 1), are interchanged with FAs 57 and 388, respectively.It should be noted that, Appendix 1 comprises additional results related to Swap 2 and 3. Furthermore, many other swaps found to have small effect on the core stability behavior.This may suggest that individual fuel assemblies contribute differently to the core stability, where some act towards destabilizing the core while others towards stabilizing the core.The challenge then is to quantify the contribution of each FA to the stability response of the core through a parameter that can help in predicting the effect of any FA swap on the core dynamic response.However, this needs additional efforts, which are beyond the scope of the current research but planned to be carried out in a future work.
In the following section, further analysis is carried out in order to explain the reason behind such high sensitivity of the core dynamics to such minor change in core design.

Stability and bifurcation analysis
The main goal in this section is to get more insights on the reasons behind the observed high sensitivity of the core dynamics to such minor core design change.To this aim, the stability behavior is analyzed for both cores, i.e. original and reshuffled, in the neighborhood of the operating point at which the O2 instability was triggered.This operating point is located between points 10 and 11 in Fig. 1 and characterized by 65 % total power and 33 % total flow conditions.In this context, a detailed local stability and bifurcation investigation is carried out for both cores in order to study how the solution manifold of the system varies as a function of the mass flow, selected here as the bifurcation parameter.
The results of this detailed stability investigation are collected in Fig. 11 and show the different solutions obtained for both cores at different mass flow conditions, while the rated power is kept constant at 65 %.
With regards to the original core, as can be seen in the left-hand side of the figure, while the oscillation amplitudes are decaying with 33.5 % flow condition, indicating a stable core (stable fixed point solution), a reduction of the flow by 0.3 % (at 33.2 %) the solution type changes completely, indicating the occurrence of a bifurcation, and becomes a stable limit cycle with a very small amplitude.
The alteration of the solution type, from a stable fixed point (at 33.5 %) to a stable limit cycle (at 33.2 %), is comprehensible and it is a signature that a stability boundary, separating stable solutions from unstable ones, is located between these two values, exactly at about 33.4 % flow, as illustrated in Fig. 11 (left-hand side).In addition, the appearance of the stable limit cycle is an indication that the system loses its stability though a supercritical Hopf bifurcation.This deduction is based on the experience accumulated along the years using both reduced order models (ROMs) and system codes where the so-called correspondence hypothesis was proposed (Dokhane et al., 2007).The hypothesis A. Dokhane et al. defines the correspondence between the obtained limit cycle (stable or unstable) using system codes and the underlying Hopf bifurcation, being the main type of bifurcation encountered so far in BWR instabilities.
More interesting, when the mass flow is slightly reduced to about 32.8 %, a drastic change in the behavior of the core dynamics is observed.More specifically, the power that was oscillating around the

Table 3
Burnup, Power, Flow, Void and Power/Flow values for Swap 1, 2 and 3. steady state value, i.e. 65 %, with a very small amplitude (see case with 33.2 % flow, Fig. 11 (left)), it is now oscillating with very large amplitudes and reach a maximum value almost twice as much as that of the steady state, i.e. about 120 %.This illustrates that the high sensitivity of the dynamics of O2 core is not limited to minor changes of core design but also to slight changes of other parameters like the mass flow.It should be noted however that, this drastic change in the unstable region is uncommon, based on the past experience on bifurcation using system codes, in the sense that, for a supercritical Hopf bifurcation, when a small change in the bifurcation parameter is introduced, the limit cycle amplitude change would be rather small (Dokhane, 2004).Hence, it is suggested that a high-resolution bifurcation investigation should be carried out to understand how this drastic change in the limit cycle amplitude arises and check if more complicated bifurcation types exist in this parameter region, which is beyond the current scope.
With regards to the reshuffled core, based on Fig. 11 (right-hand side), as can be observed, the solution manifold landscape and type of bifurcation is qualitatively similar to that of the original core, except a noticeable shift in the stability boundary towards lower values of mass flow, from 33.4 % (for original core) to about 32.6 %.Hence, the effect of the change introduced in the core design has a stabilizing effect because the original design core that was unstable in the interval 32.6 %-33.4 %, it is now stable when reshuffled.
More interestingly, if the behavior of the two cores is analyzed at the same operating point, e.g. with 32.8 % flow condition (Fig. 11), a significant difference between the two behaviors is clearly seen, i.e. from a highly unstable core with maximum oscillation amplitude reaching 120 % of rated power, for original core, to a stable state for the reshuffled core.This is due to the fact that this operating point is located in the stable side for the reshuffled core but situated in the unstable region for the original core and due also to the high sensitivity of the system within small parameter intervals in the unstable region where the supercritical Hopf bifurcation occurs.
Hence, based on the above detailed stability and bifurcation analysis, it can be concluded that, the reason behind the drastic reduction in the power oscillation amplitudes in the O2 event, following a minor core design change, is mainly due to the shift in the stability boundary making the core more stable and also this reduction is also due to the fact that the event is taking place very close to stability boundary, region in which the system, unexpectedly, found very sensitive to small changes in the unstable region where supercritical Hopf bifurcation is expected.

Summary and conclusions
The current research is an additional study carried out on the OECD/ NEA benchmark on the Oskarhsamn-2 Stability event, occurred on 25 February 1999, in order to explore hypothetical scenarios that could have helped in controlling, or mitigating, or even preventing the instability and hence the reactor scram.In this context, the current study has investigated the effect of minor core design changes on the overall core behavior, where the two hottest channels of the reference core have been swapped with two channels located closer to the core periphery.Both steady state and transient analysis have been performed and results were presented in terms of many parameters such as power/flow ratio, 2D core maps for power and flow in steady state; and for the transient, channels' power time series, 2D core maps for decay ratio and resonance frequency, maximum oscillation amplitude and oscillation phase-shift.
The results show that while the introduction of such minor changes in the core design has negligible effects on the steady state core behavior, a drastic change in the dynamic behavior has been observed with a huge reduction of the power oscillation amplitude.As a result, it has been demonstrated that, if such minor core design changes had been introduced in the beginning of cycle 24, the stability event of February  A. Dokhane et al. 25, 1999 and associated reactor SCRAM might have been avoided.It should be emphasized that, similar results have been obtained even when the swap involves less limiting channels as the case with Swap 2 and 3, suggesting that the core dynamics is highly sensitive to even minor changes.
Such drastic change in the core dynamic behavior can only be explained based on a detailed stability and bifurcation analysis, which concluded that such high sensitivity is not limited to the minor core design changes but also to some parameter changes, like the mass flow.This analysis has revealed the reason behind the drastic reduction in the    In the current Appendix, additional results are presented for Swap 2 and 3.In Fig. 12, the location of the two pairs of channels involved in the different swaps are shown.
Table 3 presents the burnup, power, flow, void and power/flow ratio of two pairs of assemblies involved in each of the studied swaps.In addition, a parameter, called Rank, is showed representing a ranking of the fuel assemblies based on the value of the power/flow ratio, where fuel assembly with max value of power/flow ratio is ranked 1, while the assembly with minimum value is ranked 444 (number of FAs in the core).
Table 4 summarizes the main steady state parameters for the Original core and the three different reshuffled cores based on Swap 1, 2 and 3.
Figs. 13 and 14 represents the core FA power/flow ratio value distribution before and after the performance of Swap 2 and 3, respectively, pointing out the two pairs of assemblies involved in the swapping.It should be noted here that, in the two swaps, even the hottest channels of the original core (channels #146 and 299) are not involved in the swapping, their power/flow ratio values are affected with a clear decrease.
Figs. 15 and 16 represents the stability parameters results in terms of 2D distribution of DR and RF for Swap 2 and 3, respectively.

Fig. 15 .
Fig. 15.2D O2 Core Map of Decay Ratio and Resonance Frequency for Swap 2.

Fig. 16 .
Fig. 16.2D O2 Core Map of Decay Ratio and Resonance Frequency for Swap 3.
A.Dokhane et al.Appendix 1. Results summary related to Swap 2 and 3

Table 1
Limiting Core Values of Burnup, Power, Flow, Void and Power/Flow Ratio.

Table 2
Main Steady State Parameters.

Table 4
Main Steady State Parameters for Original core and different swaps.