Stability in boiling water reactors: Models and digital signal processing

Stability in boiling water reactors: Models and digital signal processing

1

All human knowledge begins with intuitions, proceeds from thence to concepts, and ends with ideas. Immanuel Kant, Critique of Pure Reason (1781), B 730.

1.1

Nuclear power plants and their impact in our world

Energy represents one of the most remarkable aspects of the contemporary world, being the basis of development and sustainability for each of the countries that make it up. Currently, population growth requires an increasing consumption of energy, which must be generated from an existing source. Added to this is the problem of global warming, which paints a complicated image for the future. Both problems, sufficient and clean energy, require a unified response from all the countries of the world because, obviously, this is no longer a local or regional problem. The solutions point to alternative energy sources, such as nuclear, and, in general, to renewable energy sources, as well as to radical changes in our lifestyles that allow us to minimize our energy consumption. Global warming has devastating consequences, manifested in an increase in the intensity of hurricanes, floods, glacial melting, and forest fires, as well as ecological structural changes that affect the migration of animal species and their extinction, to mention some important aspects of increasing the Earth’s temperature. Many countries are pursuing greenhouse gases (GHG) mitigation policies, which result in the increased use of renewable sources in the electricity sector to mitigate these CO2 emissions. Wind power is the clean source preferred option to mitigate these emissions. However, due to its intermittence, backup power is needed, which in most of the cases must be provided by combined cycle thermal plants using natural gas. An alternative strategy based on nuclear power can definitely be used as part of that policy for CO2 mitigation. Indeed, nuclear power plants have an essential characteristic: they can generate high energy density without emitting GHG (a nonemitting CO2 source). Also, they do not suffer from the problem of intermittency. However, their main drawback is the high investment required for deployment. Nuclear energy is a nearly carbon-free technology that has progressed through several generations of development and that can compete favorably with alternatives for base load electricity generation and enhance the security of energy supply. Periodic surveys give confidence that available natural resources and technological progress will ensure that nuclear Linear and Non-Linear Stability Analysis in Boiling Water Reactors. https://doi.org/10.1016/B978-0-08-102445-4.00001-1 Copyright © 2019 Elsevier Ltd. All rights reserved.

2

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

fuel requirements will be met in all scenarios of the development of nuclear energy (Alonso & del Valle, 2013). According to the economic analysis of Alonso and del Valle (2013) for an alternative strategy for CO2 mitigation, a strategy based on nuclear deployment is less expensive than a strategy that uses wind power with combined cycles as a backup. At the beginning of the 1950s, the commercial development of nuclear reactors began, highlighted by thermal reactors cooled by water, whose technology is based on boiling water reactors (BWRs) and pressurized water reactors (PWRs). Currently the BWR is of direct cycle and the PWR is of indirect cycle to produce electrical energy. Nearly 70 years ago, commercial nuclear energy appeared and approximately 500 nuclear reactors are now in operation in the world, distributed in approximately 30 countries, and more than 70 are under construction with advanced technology. Currently, there are about 80 BWRs operating all over the world that provide an important contribution to the overall electricity supply from the nuclear source. In this book, we are going to focus on the linear and nonlinear behavior of BWRs. BWRs have characteristics that distinguish them from the other types of reactors: refrigerant undergoes a phase change in the reactor core and is a direct cycle. The evolution of BWRs, which includes technological improvements and innovations, shows that these unique characteristics have been maintained into the present. Thus, the studies and advances in nuclear-thermalhydraulic instabilities in BWR cores are current and are still of scientific and technological interest. However, in the beginning, the first commercial BWR was of dual cycle and natural convection. The ABWR (advanced BWR) model is the most advanced of the BWRs in operation, whose main characteristic is the reactor internal pumps. The last model to be designed (that still has not come into operation) is the ESBWR (economic simplified BWR), whose significant characteristic is natural convection, and which is equipped with passive systems for safety. In Chapter 2 we will introduce and explain in detail the technological evolution of BWRs. The problems of radioactive waste and storage are latent risks, however, they can be reduced in important ways through the development of new technologies, both in the storage and recycling of fuel in nuclear reactors. These developments can include the design of nuclear fuels with burn cycles much greater than the current ones, including generating its own fuel while operating (Nun˜ez-Carrera, Lacouture, del Campo, & Espinosa-Paredes, 2008). The fourth-generation reactors contribute to this idea of sustainable nuclear reactors. The generation IV (G-IV) reactor must meet the following design criteria: economic competitiveness, inherent safety, minimization of waste, nonproliferation enhancement, and social acceptance. The Generation IV International Forum (GIF) seeks to develop a new generation of nuclear energy systems for commercial deployment by 2030–2040. These systems include both the reactors and their fuel-cycle facilities. The aim is to provide significant improvements in economics, safety, sustainability, and proliferation resistance. The systems selected for development are the very high-temperature gas-cooled reactor (VHTR), the sodium-cooled fast reactor (SFR), the gas-cooled fast reactor (GFR), the lead-cooled fast reactor (LFR), the molten salt reactor (MSR), and the super-critical water-cooled reactor (SCWR). We can observe that among the proposed G-IV technologies, only a

Stability in boiling water reactors: Models and digital signal processing 650

SCW CANDU

500

Temperature (°C)

3

SCWR (US)

SCW

R

25

M

Pa

Pa

Critical point 15

M

PWR

350

10

Pa Pa

CANDU-6

320 310

M

7M

BWR

285 280 265 3.0

3.5

4.0

4.5

5.0

5.5

6.0

6.5

Entropy (kJ/kg K)

Fig. 1.1 Operating conditions of current nuclear reactors and SCWRs. Based on Mokry, S., Lukomski, A., Pioro, I., Gabriel, K., & Naterer, G. (2012). Thermalhydraulic analysis and heat transfer correlation for an intermediate heat exchanger linking a SuperCritical Water-cooled Reactor and a Copper-Chlorine cycle for hydrogen co-generation. International Journal of Hydrogen Energy 37(21), 16542–16556.

water-cooled reactor was considered, which can be an important option when thinking about the possibility of a supercritical direct-cycle nuclear power plant as an option instead of a dual cycle. The SCWR is one of the most promising and innovative designs. This is a very high-pressure water-cooled reactor that will operate under conditions above the thermodynamic critical point. Water enters the reactor core and then exits without change of phase, that is, no water/steam separation is necessary. The thermal efficiency is on the order of 30%–35% for current nuclear power plants, and approximately 45%–50% for SCWRs. In Fig. 1.1 we show the difference in the operating conditions of current generation reactor systems in comparison to SCWRs. In SCWRs, the target is to increase the coolant pressure from 10–16 MPa to about 25 MPa, the inlet temperature to about 350°C, and the outlet temperature to about 625°C.

1.2

BWR and the stability issue

In December 2011, everyone in the world knew the word Fukushima. An accident occurred in the Fukushima Daiichi nuclear reactor located in this city of Japan due to the occurrence of a huge tsunami after a massive earthquake (Amano, 2011).

4

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

The earthquake situation affected the electrical power supply lines at the site and the tsunami caused massive damage to the site infrastructure, which resulted in the loss of the cooling system. Every reactor in the world is susceptible to a severe accident. Fukushima was not the first. In 1979 (NRC, 1979), a severe accident occurred at the Three Mile Island Unit 2 (TMI-2). A failure in the nonnuclear part of the plant triggered a series of automated responses in the reactor coolant system (RCS) and the relief valve at the top of the pressurizer failed to close when the pressure returned to a proper value. After that, the operators were unaware that cooling water was pouring out of the stuck open valve. The lack of proper water flow allowed the reactor core to become partially uncovered and severely damaged. As in the case of the TMI-2 and the Chernobyl accident (which occurred on April 26, 1986), public concern around the world over nuclear safety has been aroused by the Fukushima Daiichi nuclear accident. Because of its severity and enormous impact on public opinion, it is important to study the accident and its implications. Both the Chernobyl and Fukushima accidents were categorized as major accidents (Level 7) on the International Nuclear and Radiological Event Scales (INES), while TMI-2 was classified as Level 5 (Murray & Holbert, 2015). These severe accidents significantly destroyed the infrastructure and were potentially harmful to humans and the environment. Consequently, it is necessary to study and understand the progression of severe accidents to prevent their occurrence or mitigate their consequences in case they are unavoidable (Espinosa-Paredes, Batet, Nun˜ez-Carrera, & Sugimoto, 2012; Espinosa-Paredes, Camargo-Camargo, & Nun˜ez-Carrera, 2012). After the Chernobyl accident, power uprate represented one of the best options to obtain the major benefits of nuclear power plants (NPPs) because of the over regulation for the construction of new NPPs and because of public opinion against nuclear energy. However, plant owners have been performing power uprates since 1970. The existence of thermalhydraulic-neutronic instabilities in boiling water reactors have been known since the beginning of BWR research using prototype reactors. Indeed, the physical mechanism, in which a recursive flow change between an upper and a lower limit, whose behavior is of a periodic nature, was observed at the end of the 1930s in two-phase flow components and systems of industrial plants and boilers. In many BWRs, different power oscillation phenomena have been observed, especially the ones produced at low-flow and high-power operational conditions that have been classified as in phase or out of phase according to whether the neutron flux oscillates: in phase over the whole reactor core or out of phase between individual half core regions. Besides, in this book the term of partial out-of-phase oscillation is applied to establish modes of oscillation or stability that can be considered in the transition region between stable and unstable behavior. Another possibility is a mode oscillation that is not strictly defined as an in-phase or out-of-phase oscillation, but it is also possible that there is a combined mode of in-phase and out-of-phase oscillation with transitions between them. The simultaneous oscillations were interpreted as dual oscillations and this was considered as a novel instability phenomenon in which it was found that and in-phase and an out-phase mode were existing concurrently, with the same resonance frequency, but different stability properties and signal amplitudes.

Stability in boiling water reactors: Models and digital signal processing

5

There are different causes of instability in two-phase flow systems and some of them can be present in a BWR, for example, flow pattern transitions in two-phase, acoustic instabilities, and the density wave mechanism. Other sources of instabilities are the steam generation systems subjected to flow instabilities due to parametric fluctuations and inlet conditions, which may result in mechanical vibrations of components and system control problems. The BWR is a complex system in which the sources of the instabilities are diverse and whose nature is also diverse. For example, there is a particular anomalous behavior in the recirculation loop flow in the jet pump of the BWR known as bistable flow or bistable vortex, and, according to preliminary studies, it is induced by turbulence. Bistable flow has been attributed to random variations in the recirculation drive flow during steady state operation and is highly sensitive to the variation in the recirculation flow rate due to changes in the position of the control valve. The bistable flow induced by turbulence has been observed in the header cross of the BWR/3 through the BWR/6. These random variations produce variations in the core flow and in the reactor power from 0.1% up to 5%. As mentioned above, a two-phase flow system under certain conditions can present thermalhydraulic instabilities that can be manifested by periodic variations of flow and pressure. In this sense, we are talking strictly about pure thermalhydraulic instabilities, that is, in systems known as conventional. The explanations for these phenomena have been studied mainly at application scales (at system level), notwithstanding that the interfacial forces between the gas-phase and the liquid-phase and their hydrodynamic and thermal conditions largely govern the phenomena of instability. The model of two-phase determines the degree of precision of the representation of the phenomenon, and for this it is important to understand the scales involved with the hydrodynamic phenomena of interaction and feedback of the concerned system. The instabilities can be linear and nonlinear and require experiments or numerical models of simulation to understand. Along these lines, de Bertodano, Fullmer, Clausse, and Ransom (2017) recently presented a monograph on two-fluid model stability, simulation, and chaos, and this is an example of the scientific and technological importance on the state of the art about this subject. BWRs can develop pure thermalhydraulic instabilities, also in operating conditions in which neutron power is an important variable. The instability can be identified by periodic variations of the power, in addition to the variations of the core flow. There are important differences in the instabilities of the pure thermalhydraulic type and with neutron feedback. The effects of feedback can generate instabilities whose intensity is normally much greater than the pure instability in which the amplitude of the oscillation is much greater. In this book, the feedback effects are approached with theoretical depth with two-fluid models, with numerical experiments and linear analysis in the frequency domain and with the nonlinear analysis in the time domain, as well as neutron power data processing. Chapter 3 is dedicated to describing in detail the stability phenomena that can occur in BWRs. The neutron monitoring system, based on fission chamber-type detectors known as local power range monitors (LPRMs) and physically placed in the core, permits the measurement of the neutron flux across the core. In general, the data from these detectors are averaged outside the core, forming some signals called APRMs, which are

6

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

observed (and analyzed) in the control room, thus allowing the detection of a possible instability. However, an out-of-phase instability may not be observed correctly based on these APRM signals, because the LPRM signals could present a significant phase to each other, destroying the averaged signal. This scheme, through the APRMs, is more effective when instabilities are observed in phase. The data from LPRMs (or APRMs) have been analyzed for decades with signal processing techniques to understand the phenomena, especially in unusual situations in the operation of BWR reactors. In particular, we study the instability phenomenon due to a density wave, a typical characteristic present in BWRs by the single fact of operating with two-phase flows. It is still a challenge to predict this phenomenon with sufficient time to operate the reactor to safety conditions, it being essential to propose mathematical models capable of including the description of the propagation phenomena (governing the instability by wave density) covering both linear and nonlinear analysis. The only models for singleand two-phase flow proposed today still consider constitutive laws of the Fick type for mass transfer and for neutron diffusion theory, as well as the Fourier law for heat transfer, or Newton’s law of viscosity for the momentum transfer, including consideration, in a nonrigorous manner, of the sonic propagation velocities in the gas and liquid phases of the two-phase flows of an incompressible flow. Other simplifications, such as homogeneous flow (identical phase velocities) or thermodynamic equilibrium between the phases, are also included in these models. All these widely accepted hypotheses have an impact on the analysis of instability. The problem with these constitutive laws is that they are limited by not including a precise description of the propagation phenomena, still considering these phenomena as diffusive when the velocities of propagation of the waves are infinite. Based on the issues mentioned above, we think that it is necessary to change the paradigm of reduced order models (ROM) of type March-Leuba (1986) for dynamic analysis in a BWR. Indeed, the feedback phenomena between the thermalhydraulics in two-phase flow, heat transfer in fuel, neutronic process, and the reactivity by voids fraction and Doppler are perfectly defined for the dynamic analysis in the core of the reactor. The ROM typically considers neutron point kinetics equations, obtained from the neutron diffusion equation; nevertheless, it does not consider the parameters related to propagation phenomena, such as the relaxation time, whose physical interpretation is related to the response time to a disturbance. The fuel heat transfer dynamics based on the Fourier law considers that the speed of propagation of the heat transport is infinite because again the Fourier law is of diffusive type. The relaxation time in fuel heat transfer is crucial in the instability analysis. The nuclear fuel in current nuclear reactors is different than it was decades ago because the fuel cycles are larger and its dimensions (diameter of the fuel rod) smaller in order to improve its heat transfer characteristics, especially in extended power uprate in nuclear reactors. One of the most important factors in instabilities in nuclear reactors, and the one of highest incidence, is the fact that the geometrical configuration and nuclear characteristics of the nuclear fuel changes, modifying the speed of propagation by changing the time of relaxation of the fuel. If the current models do not consider the propagation velocities in the heat transfer, the results on the stability of the reactor will not be precise. This issue can be avoided (or improved) by including in the analysis the physical propagation of the heat transport wave.

Stability in boiling water reactors: Models and digital signal processing

7

In this book, we present our own mathematical models, established through fundamental principles and universal laws, with the idea of making new proposals based on a different paradigm in mathematical modeling that integrates effects that are not usually included as propagation phenomena for the modeling of the neutronic processes and processes of heat transfer in the fuel. Propagation phenomena related to feedback mechanisms in nuclear reactors are crucial, especially in BWRs, due to the phenomenon of phase change in the reactor core under the two-phase flow regime. We dedicate Chapter 4 to a presentation of the phenomena of propagation in BWR reactors. Ideas, such as neutron wave propagation velocities, are introduced, as well as propagation velocities of the wave of heat transport in the fuel and a special treatment on propagation velocities of the sonic waves in two-phase flow is carried out, including the void propagation speed waves in two-phase flow systems. Some of these ideas have already been tried in previous work, leading to a mathematical model that yields better results, but those proposals are very different from the ones presented in this book, in which the propagation phenomenon is widely addressed. Regarding the nonlinear analysis, the previous models can be simplified in order to obtain reduced order models as the proposed by March-Leuba (1986). These models allow us to explore nonlinear phenomena in the time domain in which the attractors can be studied at different conditions of instability, as well as to generate the bifurcations diagram and to apply techniques, such as the Lyapunov exponent, to understand the complex behavior of a system, such as a BWR. Fig. 1.2 depicts the cycle limit formed between neutron flux density n(t) and fuel temperature T(t) in a scenario of instability in a BWR. Following with these ideas, considering an approximation of the neutron wave propagation velocities, we can derive the neutron point kinetic equation directly with the P1 approximation of the transport equation that includes the time derivative of the 25

20

n(t)

15

10

5

0 −10

0

10

20

30

40 T(t)

50

60

70

Fig. 1.2 Example of a limit cycle in a scenario of instability in a BWR.

80

8

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

current vector, thus generating a velocity of propagation as well as a neutron relaxation time. In the case of the fuel heat transfer, the wave equation can be approximated using a constitutive law of type non-Fourier. We explore these ideas in Chapter 5.

1.3

Dynamical analysis in BWR: Introducing codes

During decades of research on the behavior of BWR reactors, numerical codes have been developed describing the integral behavior of the reactor, which includes internal devices, such as the jet pumps of the recirculation system, separators and steam dryers, and also the core of the reactor where the thermal and hydrodynamic processes of the refrigerant are carried out in one- and two-phase flow due to the fact that the lower part of the reactor core operates in a regime of single liquid phase while the rest is twophase flow with a coexisting boundary region between both regions. To describe these processes, mathematical models are based on the conservation equations of mass, energy, and quantity of momentum that depend on time and space in one, two, or three dimensions. The distribution of pressures, as well as the superficial or intrinsic velocities and temperatures, are essential in the analysis, but the variable that is crucial to consider to explain the feedback phenomena is the fraction of voids that represents the volume fraction of the gas phase in the reactor core. A typical approach is onedimensional in the axial direction, although there are fundamental works that describe the distribution of voids in the radial direction. The fraction of voids maintains a very close relationship with the neutrons, which are manifested by a decrease or increase in their power in the reactor under any perturbation because the coolant in two-phase flow constitutes a feedback void fraction as a unique aspect of the BWR. Then, a description of a mathematical model must accurately describe the void fraction, which is truly complicated. Current models are far from describing, with very good approximation, the flow topologies in two phases that come to be presented in different patterns of flow, because they depend on the interfacial forces to describe the interfacial speeds. At the moment, this is complicated and methods are used to approximate a two-phase flow model whose more advanced description in applications for BWR is the two fluid model, that is, average field equations for each of the phases that consider the phenomena of interaction between the phases. The model of heat transfer in nuclear fuel is another essential element in the description of a nuclear reactor. The approaches are varied, but a model of distributed parameters dependent on space and time is widely accepted, with three considered regions, such as the fuel region, the gap region, and the clad region, in which the thermal resistance is very high in the region of the gap because it is made of an inert gas. The model of heat transfer in the fuel describes the distribution of temperature in the radial direction with very good precision because in the axial direction the transfer of heat is usually considered negligible on account of the length of the fuel being much greater than its diameter. The effects of nonuniform heat flow in the axial direction are determined by the enrichment of the fuel rod. The transient temperature distribution obtained with this approach is crucial to consider the Doppler feedback due to the average temperature of the fuel, thus also allowing us to determine the wall

Stability in boiling water reactors: Models and digital signal processing

9

temperature of the clad to calculate the convective heat flux, that is, the heat flux transferred from the fuel to the refrigerant. The neutron processes in a nuclear reactor can be described with a point model of neutron kinetics, with a model based on neutron diffusion theory (with one or several energy groups), and even with transport theory that depends on time and space in one, two, or three dimensions. The neutron model calculates the source of heat of nuclear origin for the model of heat transfer in the fuel and the volumetric heat of gamma sources that is a source of direct heat to the refrigerant, and in turn the model of the neutron receives the volume fraction of the gas phase, coolant density, coolant temperature, and fuel temperature as fundamental interaction variables. Accordingly, the dynamic analysis of a BWR-type nuclear reactor requires three main models: l

l

l

Thermalhydraulic (single- and two-phase flow) Fuel heat transfer Neutronic process

With this model, the approximation for the dynamic analysis can be classified in: l

l

Frequency-domain analysis Time-domain analysis

The frequency-domain approach is for the purpose of linear stability analysis of BWRs. Most of the field equations describing the system dynamics are complicated nonlinear partial differential equations. These nonlinear equations can be turned into linear conservative equations in the frequency domain. Therefore, it is also called linear frequency-domain stability analysis. The system governing equations are linearized by a small perturbation about steady state, and transfer functions can be obtained between perturbed variables. Due to simplifications of the linearization, the frequency-domain analysis can be used to obtain the stability boundary, that is, the onset of flow instability, for the density wave oscillations based on certain stability criteria. The time domain approach is mainly used for transient analysis and safety analysis but can be applied to perform nonlinear analysis. Some codes in the frequency domain are presented by different authors (e.g., D’Auria, 2008; March-Leuba & Rey, 1993; Prasad, Pandey, & Kalra, 2007). The codes, such as FABLE/BYPSS (General Electric, USA), HIBLE (Hitachi, japan), K2 (Toshiba, Japan) LAPUR-5 (NRC and ORNL, USA), NUFREQ (General Electric), STAIF-PK (Nuclear Fuel Industries Ltd., Japan), and DITENU (CNSNS, Mexico), use the point neutron kinetics equations. On the other hand, NUFRQ-NP (RPI, USA), ODYSY (General Electric), STAIF (Siemens-KWU, Germany), and STAIF-PK (Nuclear Fuel Industries Ltd., Japan), use one-dimensional neutron kinetics, inclusive 2D and 3D neutron kinetics as is the case of NUFRQNP. One-dimensional parallel channel thermalhydraulic of the core is the main characteristic of these codes, except NUFREQ. LAPUR was developed at the Oak Ridge National Laboratory (ORNL) for the NRC and is currently used by NRC, ORNL, and others. LAPUR’s capabilities include both point kinetics and the first subcritical mode of the neutronics for out-of-phase

10

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

oscillations. The thermalhydraulic part is modeled to consider up to seven flow channels with inlet flows coupled dynamically at the upper and lower plena to satisfy the pressure drop boundary condition imposed by the recirculation loop. LAPUR’s main result is the open- and closed-loop reactivity-to-power transfer function, from which a decay ratio is estimated. Its current version is LAPUR-5. NUFREQ calculates reactor transfer functions for the fundamental oscillation mode. The main difference between them is their ability to model pressure as an independent variable (NUFREQ-NP) so that the pressure perturbation tests can be reproduced. NUFREQ-NPW is a proprietary version currently used by Asea Brown Boveri (ABB); its main feature is an improved fuel model that allows modeling of mixed cores. FABLE can model up to 24 radial thermalhydraulic regions that are coupled to point kinetics to estimate the reactor transfer function for the fundamental mode of oscillation. Time-domain codes are more widely used and include RAMONA-3B, TRAC-BF1, TRAC-G, RETRAN, EPA, SABRE, TRAB, TOSDYN-2, STANDY, and SPDA. Typically, when using time-domain codes, the thermalhydraulic solution requires orders of magnitude more computational time than the neutronics codes. Because of the large expense associated with the computational time, thermalhydraulic channels are often averaged into regions to reduce computational time. For example, the RAMONA-3B version was developed by BNL and has a full three-dimensional (3D) neutron kinetics model that is capable of coupling to the channel thermalhydraulics on a one-to-one basis. This code uses an integral momentum solution that significantly reduces the computational time and allows for the use of as many computational channels as necessary to accurately represent the core. TRAC has two versions currently used in BWR stability analysis. TRAC-BF1 is the open version used mostly by Idaho National Engineering Laboratory (INEL) and Pennsylvania State University, while TRAC-G is a GE-proprietary version. TRAC-BF1 has one-dimensional neutron kinetics capability, while TRAC-G has full 3D neutron kinetics capability. Typically, TRAC runs are very expensive in computational time and most runs are limited to the minimum number of thermalhydraulic regions (about 20). RETRAN is a time-domain transient code developed by the Electric Power Research Institute (EPRI). It has one-dimensional and point-kinetics capability and is a relatively fast running code because it models a single, radial, thermalhydraulic region and uses the so-called three-equation approximation. SABRE is a time-domain code developed and used by Pennsylvania Power and Light for transient analyses that include BWR instabilities. SABRE uses point kinetics for the neutronics and a single thermalhydraulic region. TOSDYN-2 has been developed and used by Toshiba Corporation. It includes a 3D neutron kinetics model coupled to a five-equation, thermalhydraulic model and models of multiple parallel channels as well as the balance of plant. STANDY is a time-domain code used by Hitachi Ltd. It includes 3D neutron kinetics and parallel channel flow across at most 20 thermalhydraulic regions. DITENU is a time and frequency-domain code used by the National Commission on Nuclear Safety and Safeguards of Mexico, which was developed by Autonomous Metropolitan University, and that includes the reactor vessel, reactor recirculation system, and nuclear thermalhydraulics with a five-equation model and point neutron kinetics model.

Stability in boiling water reactors: Models and digital signal processing

11

The latest generations of BWRs are designed with natural circulation as the operation mode under both normal and abnormal conditions. The economic simplified boiling water reactor (ESBWR), designed by GE, and the advanced heavy water reactor (AHWR), being developed in India, are natural circulation BWRs (NCBWRs). In NCBWRs, the heat removal from the core takes place by natural circulation during the rated full power operating condition as well as the start-up and accidental conditions. Prasad et al. (2007) present the state of the art on flow instabilities in natural circulation boiling loops, including numerical codes. The trend in numerical nuclear thermalhydraulic codes is in the development of multiphasic and multiscale computational platforms. These platforms are currently being developed in Mexico, the United States, and Europe, to name a few. These trends include high-performance computational platforms with fine mesh simulations and multichannel approaches that include nuclear thermalhydraulic behaviors in time and in three dimensions. The coupling of multiphysics phenomena is a very complex subject with several possible combinations. An extensive description of this is given in (D’Auria, 2004), and some details and key points are found in different works (e.g., Ivanov & Avramova, 2007; Ud-Din, Peng, & Zubair, 2011). In the past, thermalhydraulic analysis used simplified neutronics models, such as point kinetics, including, when needed, the balance of plant. The results of such simulations provided the necessary boundary conditions for the core, such as mass flow and temperature distribution of the coolant at the core inlet, together with the time functions for pressure, which could be analyzed with detailed 3D neutronics models in order to get more information. However, in reality, these boundary conditions are functions of the power generation in the reactor core. The application of these models is, therefore, limited by the consideration of proper core thermalhydraulic interface conditions and it may lead to very unrealistic accident conditions if all uncertainties are taken into account by demanding conservative boundary conditions. The coupledcode calculation approach constitutes the normal evolution of these methods. This is especially true for cases in which strong feedback between the core neutronics behavior and the plant thermalhydraulics is present, as well as in situations in which neutron flux distortions and excursions are important and the spatial distribution changes during the transient. The coupled-code calculation approach constitutes the normal evolution of these methods. In the case of system codes coupled with 3D neutron kinetics models, six basic components of the coupling methodologies have been identified in order to be able to couple two codes. Thus, the way of coupling (internal or external), the coupling approaches (serial integration or parallel processing coupling), the spatial mesh overlays (fixed or flexible), the coupled time steps algorithms (synchronization of the time steps), the numeric coupling (explicit, semiimplicit and implicit), and the coupled convergence schemes must be considered and implemented. The development of coupled codes in the past years has been fostered by safety analysis requirements. Although the most important part of a nuclear reactor is the core, several accidents originate in the primary or secondary loops or even in some other component, such as the turbine. The coupling can be achieved in three different ways: internal, external, and combined.

12

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

1.4

Digital signal processing and nuclear reactors

In 1998, the Signal Processing Society of the IEEE published an article entitled, 50 years of signal processing: The IEEE Signal Processing Society and its technologies 1948-1998 (Nebeker, 1998). It is important to mention that IEEE stands for the Institute of Electrical and Electronics Engineers. However, as the world’s largest technical professional association, IEEE’s membership has long been composed of engineers, scientists, and allied professionals. These include computer scientists, software developers, information technology professionals, physicists, medical doctors, and many others in addition to IEEE’s electrical and electronics engineering core. For this reason, the organization no longer goes by the full name, except on legal business documents, and is referred to simply as IEEE. Why this article and why this clarification? The answer is simple: for the importance that both represent in the past and present of digital signal processing, as we will explain in the following paragraphs, but also for the future, as it is conceived now, based on the interaction of many disciplines to solve the big problems of humankind. First, what exactly does digital signal processing (DSP) mean? As is mentioned in the IEEE-SPS article above, DSP is an immense and diverse field that remains mysterious or quite unknown to most people. Indeed, it is not a simple task to give a formal, concise, and precise definition of this domain. However we will try to explain it in a simple way using the scheme shown in Fig. 1.3. Here, DSP is depicted as the result of the interaction of three big blocks: theoretical bases, methods, and instrumentation. Digital filtering

Stochastic processes Quantum mechanics

FFT

Number theory Linear algebra

Probability

Physics of waves Functional analysis

Hardware

EMD

Compressed sensing

Time-frequency analysis

Instrumentation

Astronomy Music Geophysics Seismology Vibrations

Adaptive filtering

Methods

Theoretical bases Systems theory

Wavelets

Neural networks

Optics

Applications

Software

Nuclear energy Telecommunications Biology Biotechnology Transport Chemical industry

Radio astronomy Archeology

Fig. 1.3 Scheme for Signal Processing: theory and methods.

Electrical industry Medicine

Multirate filtering

Stability in boiling water reactors: Models and digital signal processing

13

All three of these work together with only one goal: the application. Possibilities based on this paradigm are huge. To illustrate this, we consider a recent application: the speech communication in a smartphone. The information (digital speech) sent to the receptor (destination) is a signal traveling through free space as an electromagnetic wave. Nevertheless, strictly speaking, the transmission of this signal is not considered to be signal processing, but rather the change made to this signal to fit (or improve) the transmission. What kind of improvement is needed for this specific signal? The problem in a digital transmission is the bits per seconds (bits/s) necessary to represent speech correctly (64,000 bits/s is a standard quantity). For real-time communications this issue becomes a challenge. The solution was a signal compression method, which is basically a sophisticated technique involving an autoregressive model to synthesize the speech, and to code, in a smart way, the coefficients of this model and the model error, which are sent to the receptor instead of the complete real speech signal. This signal compression method needs fast algorithms to be implemented in real time, achieving all necessary operations in a time window of 20 ms. The final transmission rate, achieved in a recent smartphone, is around 5000 bits/s, and this is twelve times less than the original transmission rate. We can mention some improvements that can be applied to the signals, which are considered to be signal processing: filtering, coding, estimating, detecting, analyzing, recognizing, synthesizing, recording, and reproducing. All of these types of processing can be applied in both the analog and digital domains, however, recent technology is dominated by the digital world. It is for this reason that we focus on digital signal processing in this book. Now, if we think about a nuclear reactor application and digital signal processing, we can begin with a simple method to estimate the power spectral density of the neutron flux inside the core using a classic fast Fourier transform (FFT) in order to observe the appearance of the density wave, through a significant frequency around 0.5 Hz, in order to detect an instability event in the core of a BWR. Here, the application could be the detection of the density wave, the method is the power spectral density, and the instrumentation is the fast Fourier transform, a fast algorithm that permits us to observe the appearance (or not) of the frequency around 0.5 Hz. However, we can use more complex methods to accomplish the same goal but getting a better and more reliable estimation. That is the principal objective of this book: to look for new methods to detect an instability event in a BWR. Before explaining in detail the methods and why they were chosen, we want to present the most relevant historical moments of signal processing over more than 60 years. In Fig. 1.4 we represent a timeline with the events that were a watershed in the domain of signal processing (analog and digital). In this timeline, we include three divisions corresponding to people or institutions, the events themselves, and the most important applications developed in that context. The year 1948 is preceded by some important contributions, considered as influential in the development of DSP: the Fourier theory (Fourier, 1822), the stochastic processes theory (Wiener, 1928), the pulse code modulation (PCM) theory (Rainey, 1926; Reeves, 1938), and others. The PCM theory corresponds to the theoretical base of analog to digital conversion, impossible to implement at that time. This year 1948 marked a first big step toward DSP and it is considered an annus mirabilis.

Shannon’s article, 1948

1940s

1950s

1960s

Radar Fast Speech processing, filters: sonar convolution and Transversal, adaptive, digital filtering multirate

JPEG 2000

2000s

1970s

MPEG

1990s

Empirical mode decomposition, 1998

Green signal processing

Fig. 1.4 Principal events in the history of signal processing, bases, and tendencies.

GSM

JPEG-DCT

First toy with DSP, 1978

Mandelbrot set

1980s

Neural networks and fractals, 1980 Wavelets theory: Multiresolution analysis, 1988 and compact support orthogonal basis, 1989

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

Multiprocessor architectures with embedded systems: Cloud computing Parallel processing System virtualization

FFT, 1965

Error Invention Audio engineering correcting of transistor, Stereo recording codes 1948 Perceptron

Tendencies DSP with artificial intelligence: Big data, Internet of things, data Science, etc. complex systems

Linear predictive coding, 1970; Image processing, 2D-filtering; Markov chains

14

Influences and bases Fourier theory: 1829, Fourier Stochastic processes: 1928, Wiener Wavelets: 1910, Haar PCM: 1926, Rainey 1939, Reeves

Stability in boiling water reactors: Models and digital signal processing

15

Indeed, Claude E. Shannon published his seminal article (Shannon, 1948) in which he introduces the concept of information and its (average) measure, the entropy, and the basic ideas underlying the coding of a discrete source without noise (decoding error). He also published together with Oliver, Pierce, and Shannon (1948) the “Philosophy of PCM” indicating how this technique would be important in the field of communications. Indeed, now all DSP, in all application fields, is practically implemented in a digital platform. In this same year, Maurice Bartlett and John Tukey developed digital methods of spectrum estimation that have remained in use ever since. These methods, applied especially in radar signals, attempted at that time to estimate a continuous spectrum based on sampled data: thus, the modern spectral estimation domain was born. Also, in 1948, the Bell Telephone Labs announced the invention of the transistor, the solid-state device that would change the field of electronics in the next decades. This device was the catalyst of DSP development in the decade of the 1960s, in which some important events happened. The most important event of this decade was, without a doubt, the discovery of a fast algorithm that allowed the number of operations in the calculation of the discrete Fourier transform to be reduced in an incredible way. This algorithm, discovered by James Cooley and John Tukey in 1965, was called the fast Fourier transform, or simply FFT. This discovery, together with the fast convolution invention (Helms, 1967; Stockham, 1966), permitted the implementation of a digital filtering (Gold & Rader, 1969; Kaiser, 1966), which in turn allowed the implementation of autoregressive models under a scheme of a linear prediction coding (Atal & Hanauer, 1971; Itakura & Saito, 1970). Digital filtering had an immense impact in the 1970s, especially in speech processing: transversal filters, adaptive filters, multirate filters, etc. The image processing boom began with 2D filters and many other applications integrated DSP into their analysis. The Markov chains theory began also to be applied. In addition, hardware development led to the appearance, at the end of this decade, of the first toy, developed by Texas Instruments (1978) and called Speak and Spell, that integrated an advanced algorithm for speech synthesis that was implemented with inexpensive integrated circuits. The 1980s, with the development of integrated circuits, became the era of image processing as a recognized technical specialty. The need for automated image recognition for many military and civil applications permitted the boom, or the rediscovering, of certain techniques, such as neural networks, which began in the 1940s with the idea to emulate the way the human brain works. A neural network consists of a large number of processing units, called neurons, working simultaneously and that are interconnected by multiple links. The perceptron of Rosenblatt (1958) was revived in the 1980s together with an adaptive linear element introduced by Widrow and Winter (1988). Another important event was the creation of the JPEG (Joint Picture Expert Group), which introduced the standard in image compression, still used in the present day. The JPEG standard is based on a real 2D-FFT, the discrete cosine transform (DCT). JPEG motivated the creation of another important group, the MPEG, integrated by specialists in moving pictures, and which has also created an important standard (of the same group name) in our day. Indeed, the layer 3 of the standard MPEG 2, corresponding to the audio layer for moving pictures, became

16

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

the most important standard in audio compression: the MP3. Besides, in this same decade, a mathematical decomposition technique appeared, which was called wavelets from the French word ondelette, used for the first time by Morlet in 1948 (Goupillaud, 1997) and considered as an extension of the Fourier analysis. The wavelets theory is marked by the appearance of some elements along the timeline of the 20th century: 1910, Haar families; 1981, Morlet’s wavelet concept; 1984, Morlet and Grossman’s wavelet; 1985, Meyer’s orthogonal wavelet. However, in 1988 and 1989, with the works of Stephan Mallat and Ingrid Daubechies, two fundamental elements were introduced into this theory: the multiresolution analysis and its implementation with fast algorithms, which allowed a compact support orthogonal (or bi-orthogonal) basis to be obtained. The use of wavelets theory, especially of the multiresolution analysis, has never stopped since then. Based on this, a new standard for image compression appeared, the JPEG 2000, and applications based on wavelets became uncountable. In this same decade appeared fractals, objects in geometry that have nonintegral dimension. DSP integrates the concept of fractals, together with chaotic models and wavelets to improve certain applications, especially in image compression. In the 1990s, some techniques, considered as a part of artificial intelligence, were employed in many areas of signal processing. Two examples of this, among others, are fuzzy logic and the genetic algorithms used as an optimization tool. It was the beginning of many disciplines converging on the solution of a problem. At the end of this decade, in 1998, Norden E. Huang et al. published an important article introducing a new technique called empirical mode decomposition, which seeks to decompose a signal in its intrinsic oscillating modes, forming an inherent base of the signal, that is, a data-driven signal decomposition method. This nonlinear technique and its improved and multidimensional variants have an important place in the DSP domain, which has lasted into the present and spread its influence over many application fields. In 1965, Gordon Moore, an Intel founder, established his famous law: the Moore’s law, which holds that approximately every 18 months the number of transistors in a dense integrated circuit (IC) would be duplicated. In 2015, Intel admitted that this advancement had started to slow down, which means the future of Moore’s law is currently undetermined. In 1988, Gene Frantz, a Texas Instruments Principal Fellow, interested in how to make products more portable, considered that power dissipation must be an important part in the design of ICs, in order to transform, at that time, the cellular world from analog to digital (Frantz, 2012). Based on his experience in the DSP domain, he took the basic function in the DSP hardware domain, the multiplyaccumulate (MAC) function, and plotted along time the power consumption, measured in milliwatts per millions of MACs, necessary to perform this MAC instruction. The result was amazing, with this power dissipation being reduced at the rate of half every 18 months. He called this result the Gene’s law. This law, in conjunction with the high integration, allowed the creation of specialized digital signal processors for cellular phones, with the promise of low-cost, ubiquitous, portable communications. Nowadays, due to the uncertainty of high integration in a few nm (Moore’s law), many hardware designs are oriented to a green design with a low power consumption in every device in order to produce a minor carbon footprint, that is, there is more interest in

Stability in boiling water reactors: Models and digital signal processing

17

reaching the Gene’s law. Based on these paradigms, high integration and low power consumption, the principal theories of digital signal processing, are now widely implemented in an embedded processor world. How about DSP, BWRs, and the stability issue? The decade of the 1980s saw a boom in studies about the stability problem in BWRs, when practically all of the proposed methods, as will see in Chapters 5 and 6, were focused to determine (estimate) the stability parameter linked to a classic linear second order model, the damping ratio (or the decay ratio, DR). This linear model (described by a linear differential equation with constant coefficients) and its stability parameter can be linked first to the Laplace domain, in order to get poles and zeroes, and hence to the Fourier domain to get the frequency response of the BWR system. Both linear inverse transformations, from the Laplace or Fourier domain, allow the impulse response of the system to be obtained. The damping ratio (or the decay ratio) is simply the parameter, in this linear modeling, representing how much this impulse response is damped. Then the goal is to determine when this parameter is going to an undamped (instability) case. For this, the usual measure is simply the ratio of taking the two first crests of this impulse response as it is shown in Fig. 1.5. This ratio corresponds to the decay ratio of the model and the system is undamped (unstable) when the DR > 1. Regarding the BWR stability issue, as we mentioned before, the instability phenomenon is due to the appearance of the density wave and it is marked by a growing oscillation around 0.5 Hz. To detect this oscillation, a system integrated by sensors measuring the neutron flux is deployed axially and longitudinally inside the core. These local detectors, called LPRMs, provide the only signals available inside the core to track a possible instability event. If we are using a second order model to describe the dynamics of the neutron flux in a BWR, we need to match this model with the available signals that stem from the LPRMs. In practically all of the linear methods developed until now, this is accomplished in the following way. First, the LPRM (or an averaged ensemble of these, called an APRM) is fitted by an autoregressive model (AR). This kind of linear model is completely described by the autocorrelation function (ACF) of the signal to be modeled (as we will see in Chapter 6), in this case the Impulse response or ACF A1 A2 DR =

A2 A1

t

Fig. 1.5 Decay ratio measured from the impulse response of a linear second order model.

18

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

LPRM or the APRM signal. The Fourier transform of the ACF corresponds to the power spectrum density (PSD) of the modeled signal, that is, the signal power content along frequency. Now, the squared magnitude of the Fourier transform of the impulse response of the linear second order model (presented in Chapter 5) also corresponds to the PSD, linking in this way this model to a real signal stemming from a BWR, that is, the ACF and the impulse response, for practical purposes, are the same. In general, in all practical linear methods developed until now, the idea is simply to fit the ACF (obtained or not by an AR model) and to get the DR measuring the ratio of its two first peaks (see Fig. 1.5). The DR obtained in this way is considered as the global stability parameter of the BWR system. All possible techniques of digital signal processing mentioned before can be applied to the LPRM or APRM signals to detect the oscillation around 0.5 Hz and to infer its growth. Indeed, the PSD will indicate the appearance of the density wave showing a resonant frequency around 0.5 Hz. A simple FFT is necessary to implement this. However, it is important to mention that to get the real meaning of this transformation the concept of signal stationarity must be fulfilled. This condition also must be accomplished by the ACF. In Chapter 6, we begin introducing the most important aspects of stochastic processes theory to understand some important concepts involving these linear models: stationarity, ergodicity, and the Wiener-Kinchin theorem. Then we present in detail the theory about the AR models and the aspects linked to estimation of the ACF. We examine the classic linear methods based on the AR model, used to estimate the ACF or the PSD in order to determine the DR and the resonant frequency, and also new methodologies involving linear methods also, but based on different digital signal methods, such as the wavelets theory (continuous and discrete), time-frequency representations, etc., corresponding to around 30 years of DSP in the BWR stability issue.

1.5

Toward a new paradigm in BWR stability analysis: Models and digital signal processing, a nonlinear approach

The use of the DR as a feasible BWR stability measure has been widely accepted, nonetheless, it has been observed that a BWR working at an operating point with a small DR can be close to instability. Also, the DR often jumps discontinuously from the well-stable to the far-unstable region. So, according to these issues, the default DR might not be a reliable stability parameter after all. Besides, in regular operating conditions, the need for linear and stationary signals might be a handicap for the DR estimates, as we will see in Chapter 7. Thus, it is important to explore new alternative methodologies and stability parameters adapted to accommodate real BWR behavior, which happens to be nonlinear (possibly chaotic) and nonstationary. In Chapter 5 we focus our efforts on introducing and explaining in detail nonlinear models fitting the dynamic behavior of a BWR, specifically reduced order models (ROM). We also present the Lyapunov exponent to study the BWR stability, but

Stability in boiling water reactors: Models and digital signal processing

19

now from a nonlinear point of view, determining whether it might become a suitable stability parameter that could replace or accompany the DR in the task of assessing unstable scenarios. We will be interested by the chaos indicator most used in practice, the most positive Lyapunov exponent, called the largest Lyapunov exponent (LEE). Based on this, we present two techniques to estimate the LLE for the studied BWR stability scenarios. One is the orbit separation method (OS), which requires a priori knowledge of a differential equation of a dynamical system (such as the ROM). So, the OS only works in the case that an ordinary differential equation (ODE)-based model of the phenomena of interest is known. The second method is a practical LLE, known as the Rosenstein estimation method. This method computes the LLE of a time series after reconstructing the associated trajectory of the system (the attractor reconstruction). No a priori knowledge of an ODE model of the dynamical system of interest is needed for the Rosenstein technique to work. However, there are a series of parameters that are not known a priori that are needed to use this method properly, so such parameters must be guessed by rule of thumb. As in the case of the linear models, we need to ask how to link these models with real signals to determine the associated stability parameter. That is not an easy task for some important reasons that we explain next. For example, in the case of the Rosenstein estimation method, this requires large data sets to compute reliable LLE estimates. Thus, the computation time for this technique is quite demanding. Besides the unknown operation parameters mentioned before, this technique is not suitable for any real-time stability monitoring of BWR nuclear power plants, as we will show in Chapter 7. Based on these facts, we can only estimate a reliable LLE with the BWR data extracted from real nuclear power plants, processed with a powerful informatics laboratory and completely offline. However, this stability parameter can be used as a reference to validate any other possible proposal (stability indicator), permitting a real-time monitoring system to be implemented. At this point, and based on all these facts, we decided to give ourselves the task of studying sophisticated methods of digital signal processing and applying them to real signals from a BWR. The methods chosen require a tradeoff: to have a reliable nonlinear stability indicator with possible real-time implementation. We clearly make a distinction between a stability parameter (global or local; e.g., DR, LLE) of a stability indicator, which is not associated directly to a linear or nonlinear model but that it is capable of establishing when a system is stable or not. Of course, we are aware that any proposal we make must necessarily be compared with the results obtained with the global stability parameter in linear systems (DR) and later with that used in nonlinear systems (LLE). The challenge is detecting the density wave oscillation, tracking in real time this oscillation, and computing (estimate) a nonlinear stability indicator. All these operations must be implemented in a certain time, the minimum possible, in order to offer a reliable alarm to the operator in a nuclear power plant. The task of detecting the oscillation around 0.5 Hz could be considered as an easy one, but in real signals of neutron flux that is not the case. Indeed, in a stable state, there is not an oscillation, and in a completely unstable state this oscillation will be totally developed. Using a global FFT we can detect a resonant frequency around 0.5 Hz, but we cannot decide when, in time, this effect appears due to the density wave.

20

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

The short time Fourier transform (as we will see in Chapter 6) is a technique used to study the time dependence of the frequency content when a signal is nonstationary, which is the case of a real BWR signal. This technique (and other time-frequency representations) segments the nonstationary signal is small time windows in which the signal can be considered as stationary and performs a classic FFT, giving a representation of the power content of this signal over time. In each time window (signal segment) we can apply the theory of stochastic processes to get an AR model, the ACF, to finally determine a local stability parameter, the DR, in the linear case. Actually, this technique is widely applied in other domains, such as telecommunications for instance, in which the speech is coded in cellular phones and it is considered stationary every 20 ms and processed with an AR model in this time window. In a real stability monitor, this can be a plausible solution with a local DR estimated at each certain time and applying an adequate decision rule after confirmation of an incipient density wave oscillation. Besides in this same chapter, we also applied the wavelet theory to explore new alternatives for instability analysis. Nonetheless, in general, BWR signals are nonstationary and nonlinear, thus Fourier-based or wavelet-based methods might lead to biased stability studies. In considering how to overcome issues like linearity and stationarity, we decided to explore nonlinear methods with the goal of detecting the density wave oscillation and an incipient instability event, and thus permitting us to raise an opportune alarm. First, we focused on a nonlinear method to detect the oscillation around 0.5 Hz over time. The chosen method was the technique introduced in 1998 by Huang et al., the empirical mode decomposition (EMD). This simple but powerful technique is widely used in many disciplines with excellent results and consists of extracting, in time, oscillation modes contained in a signal with no a priori knowledge and completely independently from signal stationarity. The extracted modes are called intrinsic mode functions (IMFs). In Chapter 7, we introduce this technique and, as we shall see, the results are excellent; we can perfectly track the oscillation around 0.5 Hz, even one that is not completely developed. For this, we used the Hilbert transform coupled with EMD to get an instantaneous frequency of the extracted IMF, corresponding to the density wave oscillation. This concept of instantaneous frequency (explained first in Chapter 6 and implemented with EMD in Chapter 7) is better adapted than the global idea of a resonant only frequency because, in practice, the oscillation around 0.5 Hz can appear weakly, rendering the FFT or the short time Fourier transform incapable of detecting it adequately. We coupled in a first approach, the EMD with the ACF of modes extracted from a LPRM (or APRM) signal (synthetic or real) and we estimated a ratio based on the two first crests of this ACF, similar to DR estimation, in time windows of 15 s. In order to test the performance of this methodology in detecting the oscillation around 0.5 Hz and the beginning of an instability event, we used real BWR signals from controlled instabilities in the Swedish reactors conforming to the Forsmark and Ringhals stability benchmarks (see Chapter 8) and a real instability event that occurred in 1995 in the Laguna Verde nuclear power plant (LVNPP). Results, shown in Chapter 7, clearly are fruitful. However, an important issue was recently found in EMD: the mixing problem. In the case of BWR signals, this means that the mode around 0.5 Hz sometimes is not

Stability in boiling water reactors: Models and digital signal processing

21

completely extracted in one IMF but is partially contained in some different IMFs, making the detection biased or completely wrong. Many efforts, including ensemble EMD, complete EMD with assisted noise (CEEMDAN), improved CEEMDAN, and more recently, noise assisted-multivariate EMD (NA-MEMD), have been developed to avoid this issue. The NA-MEMD, as we will see in Chapter 7, is no doubt the better technique. The other big issue is how to link this nonlinear mode decomposition with a nonlinear model, in a first step, with a nonlinear stability indicator at the end. EMD does not have a general mathematical framework for its construction (even though some recent efforts are moving toward this), making this practically impossible. Based on these facts, we focused on the research of a nonlinear stability indicator, in a separate way, of the resonant frequency of 0.5 Hz. The LLE would be the natural nonlinear stability parameter, however, as we mentioned before, this requires large data sets to get a reliable estimate. Based on this, we present two nonlinear stability indicators: the Shannon entropy and the fractal dimension. In Chapters 7 and 8, we develop some ideas, including two nonlinear measures capable of indicating when an oscillation due to the density wave appears and to infer its growth, permitting us to perfectly differentiate the stable state from the unstable one. The first functional introduced is the entropy, defined by Claude E. Shannon in 1948 in his seminal work: A mathematical theory of communication. This work was the beginning of the new domain of information theory. Shannon proposed this functional, that he called entropy, to measure the average information in a discretetime source having a finite number of amplitude values, each with a specific probability of appearance. The idea is simple but sophisticated at the same time: to give a unique average measure of randomness of a discrete-time source. The maximum value of the Shannon entropy occurs when all amplitude values of the source have a uniform distribution, that is, they have exactly the same probability of appearance. Any other amplitude distribution has minor Shannon entropy. Observing how real signals from a BWR behave, we consider the possibility to use this measure as a stability indicator. Indeed, in a stable state, the BWR operates at nominal power rate (a constant value). Applying the concept of Shannon entropy to a perfectly known value, this measure is simply zero. When the density wave oscillation appears and grows until it become a perfect sinusoid of 0.5 Hz (a completely developed oscillation), the amplitude distribution of this signal tends toward a uniform one (but not completely) and the Shannon entropy tends to be maximal. For the case of real signals, acquired from detectors placed inside the core, these have continuous amplitude. To apply the concept of Shannon entropy, it is necessary to get an estimate of the amplitude distribution of these signals. This is necessary to implement a method that permits us to get the best histogram, in which the bins number will be the number of discrete amplitudes of the BWR signal, to compute the Shannon entropy. To have a measure ranging between 0 and 1, we normalize this measure dividing by the bins number. For the stable state, a real BWR signal behaves as a Gaussian noise with a specific variance around the nominal power rate (Gaussian distribution mean). The normalized value of the Shannon entropy of this Gaussian distribution can be established as the lowest threshold for detecting an instability event. When the oscillation around 0.5 Hz

22

Linear and Non-Linear Stability Analysis in Boiling Water Reactors

is detected, the Shannon entropy will increase, approaching the maximum value of 1 (an instability event completely developed). In Chapter 8, we present the stability monitoring system of a BWR, emphasizing the process computer that integrates the algorithms capable of detecting an instability event in real time. We begin with a brief presentation of two available BWR stability benchmarks that consist of controlled instability events implemented in two Swedish reactors Ringhals and Forsmark, the final reports of which were delivered in 1996 and 2001, respectively. The goals of these benchmarks were different. While Forsmark was designed to test algorithms based on linear models and obtain the corresponding DR and resonance frequency, Ringhals was designed to validate codes and the stability issue. In this book we use both benchmarks to validate our proposals. We continue with the monitoring system, beginning with an introduction of the BWR core instrumentation, describing the local detectors, the LPRMs. The most important part of this chapter is the presentation of the stability monitors. For this, we explore some possibilities, including the EMD, for detecting the oscillation around 0.5 Hz and implementing both linear (DR) and nonlinear (Shannon entropy and fractal dimension) stability indicators. These stability monitors are designed to be implemented in real time, as will be explained in this chapter. We focus on implementing the monitoring system, considering all local detectors (LPRMs) and not only the averaged ones (APRMs). This monitoring system permits also the detection, in real time, of a possible appearance of an out-of-phase oscillation. For this, we propose to use the bivariate EMD (BEMD) and the cross-correlation function (see Chapter 7) as a first method, and the instantaneous phase (see Chapter 8) as a second method. The decision rules are adapted in correspondence with every stability indicator and the EMD used: EMD, BEMD with DR, and NA-MEMD, instantaneous phase with Shannon entropy and fractal dimension. To close this chapter, we want to show how the convergence of three big disciplines physics, mathematics, and computing - has become necessary. In Fig. 1.6A we show a classical path beginning with a physical phenomenon modeled by the heat conduction equation, leading to the mathematical theory developed by Fourier in 1829 and implemented in its discrete version via a fast algorithm (the FFT) in 1965, as we Heat equation Density wave

Complex system: Nonlinear models Physics

Density wave Physics

Physics

Computing FFT FFT-based methods

(A)

Mathematics Fourier transform AR modeling

Computing Parallel processing

(B)

Mathematics M-EMD with Shannon entropy Fractal dimension

Computing Multiprocessing architectures Embedded hardware

Mathematics Nonlinear DSP methods with Artificial Intelligence

(C)

Fig. 1.6 Some examples in how disciplines of physics, mathematics, and computing converge.

Stability in boiling water reactors: Models and digital signal processing

23

mentioned before. In the same way, we can mention an earthquake (physical phenomenon) in which recordings of the acceleration of the ground (accelerograms) during the phenomenon were analyzed by J. Morlet in 1948 through little waves dilated and translated in time. This analysis laid the foundations of the mathematical framework now known as wavelet theory, which was formalized in the 1980s and implemented using fast algorithms (à trous algorithm) at the end of this decade. In the BWR stability analysis, we can establish a similar path: the density wave (physical phenomenon) producing an oscillation around 0.5 Hz, fitted by a linear AR model (mathematics) and implemented by FFT-based methods to determine a stability parameter, the decay ratio. However, now more possibilities are open if we consider the convergence of these disciplines to solve a big issue (Fig. 1.6B) and not just in the same way as a classical path. For instance, for the BWR stability issue, we can begin considering what real signals are available in practice, namely, many LPRMs placed inside the core. Why not think about what we can do with all of them simultaneously? The answer is straightforward: practically anything. Indeed, now we have, with multiprocessing architectures embedded in hardware systems and capable of solving the most complex algorithms in incredible times, the possibility to execute, in real time and in well-designed systems, even algorithms, such as multivariate EMD. For instance, we can use parallel processing implemented in cheap systems, such as graphics processor units (GPUs), grouping a certain number of LPRMs conforming channels to get just an APRM signal, as is now implemented in a real BWR. But now the methods (linear or not, sophisticated or not) can be applied directly, separately in each LPRM, or considering the complete ensemble of LPRMs, to get a more reliable stability parameter. In this way of thinking, the mathematics can be open to any method capable of offering a reliable stability indicator, but also capable of being implemented in a short time. Now systems are viewed as complex systems, in which a simple method of DSP is not enough to solve problems. Thus, many solutions using DSP are linked to artificial intelligence and its implementation is multiprocessors based (Fig. 1.6C). This book is intended for both researchers and engineers interested in exploring new horizons in the BWR stability issue, focusing on classical modeling and digital signal processing, but also on tendencies in these disciplines. This background is not limited to BWRs and can be applied to other kinds of reactors, such as PWRs, or in other disciplines. Possibilities are now open.