Decay ratio estimation based on time–frequency representations

Annals of Nuclear Energy 37 (2010) 93–102

Contents lists available at ScienceDirect

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

Decay ratio estimation based on time–frequency representations José E. Torres-Fernández, Alfonso Prieto-Guerrero, Gilberto Espinosa-Paredes * División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, México D.F. 09340, Mexico

a r t i c l e

i n f o

Article history: Received 25 April 2009 Accepted 23 November 2009 Available online 23 December 2009

a b s t r a c t A novel method based on bilinear time–frequency representations (TFRs) is proposed to determine the time evolution of the linear stability parameters of a boiling water reactor (BWR) using neutronic noise signals. TFRs allow us to track the instantaneous frequencies contained in a signal to estimate an instantaneous decay ratio (IDR) that closely follows the signal envelope changes in time, making the IDR a measure of local stability. In order to account for long term changes in BWR stability, the ACDR measure is introduced as the accumulated product of the local IDRs. As it is shown in this paper, the ACDR measure clearly reflects major long term changes in BWR stability. Last to validate our method, synthetic and real neutronic signals were used. The methodology was tested on the Laguna Verde Unit 1, two events were reported in the Forsmark stability benchmark. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The estimation of the stability parameters of a boiling water reactor (BWR) has been a hot topic of research. Indeed different techniques have been developed to estimate the fundamental frequency and the decay ratio (e.g. Verdu et al., 2001; Navarro-Esbri et al., 2003; Sunde and Pázsit, 2007). Some of them are only applicable when the signal is stationary (Verdu et al., 2001), and others when the signal is non-stationary using short-time windows that provide reliable results. Navarro-Esbri et al. (2003) studied the time dependence of the natural frequency using two different tools: the short-time Fourier transform (STFT) and the time dependent power spectral density (TDPSD). Both techniques split the signal in short segments with a high degree of overlapping. Their results show that the STFT provides better accuracy than the TDPSD despite the lower sensitivity to noisy neutronic signals of the TDPSD method. It is worth mentioning that an autoregressive method (AR) was used by those authors to estimate the decay ratio (DR). However one of the most important disadvantages in the methods above discussed is the large number of floating operations (multiplications and additions) needed to be implemented. Blazquez and Ruiz (2003) analyzed the Laguna Verde event in the time and frequency domains. In the frequency domain, the transition from stability to instability was detected when the real pole of the power spectrum density (PSD) almost vanished. In the time domain, the AR method does not show a transition, but indicates how the reactor loses its stability progressively. Recently, the wavelets theory (Daubechies, 1992) has been used to explore new alternatives for transient instability analysis

* Corresponding author. Fax: +52 55 5804 4900. E-mail address: [email protected] (G. Espinosa-Paredes). 0306-4549/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2009.11.015

(Espinosa-Paredes et al., 2005). It has been shown that stability depends on several variables such as control rod patterns, void fraction, burnup, inlet mass flow, among others. A key point is that in general, BWR signals are non-stationary, therefore traditional methods such as the Fourier transform (Gabor, 1946) might lead to biased stability parameters. Sunde and Pázsit (2007) proposed an original work using the wavelet transform in combination with the autocorrelation function (ACF). They applied the wavelet ridges technique and the ACF to estimate the DR of the BWR signals. In the last few decades, research has been devoted to the study of power oscillations and the mechanisms that generate them (e.g. Saha and Zuber, 1978; Peng et al., 1984; Lahey and Podowski, 1989; March-Leuba, 1990; March-Leuba and Blakeman, 1991; March-Leuba and Rey, 1993). Several approaches have been taken to address the stability of BWRs, March-Leuba (1986) pioneered the study of reduced-order models for coupled thermal-neutron dynamics, followed by Turso et al. (1997), Van der Hagen (1998) and Muñoz-Cobo et al. (2004), among others. Some of those works were developed to gain insights about the BWRs dynamics, while others were focused on a more complete description of the heat transfer process (e.g. Uehiro et al., 1996; Guido et al., 1997; Podowski and Pinheiro, 1997). The models presented in the mentioned works were able to describe, in a qualitative way, low-frequency oscillations and even instabilities, but neither sustained oscillations of relatively high frequency nor highly non-stationary behavior could be described accurately by such models. In this paper, a method based on time–frequency representations (TFRs) is proposed to estimate linear stability parameters. We apply this methodology using the neutronic power signals of two BWRs: Laguna Verde and Forsmark. The proposed methodology is based on the Wigner TFR, also called Wigner distribution (WD). In the first part of this paper a brief theoretical background of time–frequency analysis (TFA) is presented, in the second part

94

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

we present the methodology based on the bilinear TFRs and at the end the proposed method is applied to a synthetic signal and three real BWR signals for validation purposes. 2. Preliminaries The main objective of the time–frequency analysis (TFA) is to describe how the energy of a signal is distributed in the time–frequency domain (TFD). Therefore with a given one-dimensional signal, the TFA generates a bidimensional time–frequency representation (TFR) that describes how the signal energy is jointly distributed in the time and frequency. As it can be easily inferred, for a given signal there are many solutions (TFRs), not all of them will be non-trivial and useful. Therefore, the TFA problem must be constrained to eliminate the trivial and useless solutions. The constraints are derived from the properties that a TFR ideally should meet. For any square integrable signal, s(t), the properties that the TFR should met are (Cohen, 1995): Non-negativity: Cðt; xÞ P 0 8ðt; xÞ Reality: Cðt; xÞ # R Time shift: sðtÞ ¼ f ðt  t 0 Þ; ) C s ðt; xÞ ¼ C f ðt  t 0 ; xÞ Frequency shift: sðtÞ ¼ f ðtÞ expðjx0 tÞ; ) C s ðt; xÞ ¼ C f ðt; x x0 Þ R  Time marginal: 21p C s ðt; xÞdx ¼ jsðtÞj2  PðtÞ R 2  Frequency marginal: C s ðt; R xÞdt ¼ jSðxÞj

   

xC s ðt;xÞdx

 InstantaneousRfrequency: R ¼ u0 ðtÞ  hxit C s ðt;xÞdx tC s ðt;xÞdt  Group delay: R ¼ t g ðxÞ

 Time support: sðtÞ ¼ 0 for jtj > t c ; ) C s ðt; xÞ ¼ 0 for jtj > t c  Frequency support: SðxÞ ¼ 0; for jxj > xc ; ) C s ðt; xÞ ¼ 0 for jxj > xc  Reduced interference (refer to cross terms, see Jeong and Williams, (1992)). It is worth mentioning that not all TFRs met all the properties listed above and it is not necessary at all. For instance, the most well known and used TFR is the spectrogram, also known as short-time Fourier transform (STFT), which does not satisfy the marginals. Also most TFRs take negative values despite the fact that TFRs are supposed to be densities. Last, desirable TFRs properties are application dependant in general. Given the above list of desired TFRs properties, it might seem that the TFRs design is kind of heuristic or by proof and error. Fortunately, any TFR can be uniquely expressed in terms of the general class equation (GCE) given by Cohen, (1995):

Cðt; xÞ ¼

ð2Þ

where xn is the natural resonant angular frequency, A0 is the residue magnitude, and f is the damping parameter. Here the oscillation term is represented by a complex exponential, and the damping effect is modeled using an exponentially decaying factor. We do use a complex signal for mathematical convenience to avoid cross term interference between the positive and negative frequencies components of the signal spectrum (Cohen, 1995). Next, we will show how to estimate the parameters, f and xn. The Wigner TFR or Wigner distribution (WD) of the signal, given in Eq. (2), is calculated by setting /(h, s) = 1 in the GCE:

  Z Z Z qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x u þ j xn s 1  f2  ht  sx exp 2f n 2 4p  þhu dudsdh ð3Þ

WDðt; xÞ ¼

Rearranging terms and integrating with respect to h, and then with respect to u, we obtain:

WDðt; xÞ ¼

1 expð2fxn tÞ 2p   Z qffiffiffiffiffiffiffiffiffiffiffiffiffi   exp j sx  xn 1  f2 s ds

ð4Þ

Last we integrate with respect to s to get the final result, the WD of Eq. (2):

 qffiffiffiffiffiffiffiffiffiffiffiffiffi WDðt; xÞ ¼ expð2fxn tÞd x  xn 1  f2

C s ðt;xÞdt

Z Z Z

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sðtÞ ¼ A0 expðfxn tÞ exp jxn 1  f2 t

ð5Þ

In this work, we have used the WD for mathematical simplicity but any TFR that satisfies the marginal and the instantaneous frequency could have been used. Despite any TFR satisfying those two conditions could be used, we recommend for mathematical simplicity using a signal-independent kernel leading to the well known bilinear TFRs which are easy to design and implement. If the reader is interested in signal-dependent kernel design, we refer him to Cohen (1995), Baraniuk and Jones (1993). Next it is shown how to estimate, xn, from the WD, taking in consideration that the WD satisfies the instantaneous frequency and the marginal properties (Cohen, 1995). In this case:

R qffiffiffiffiffiffiffiffiffiffiffiffiffi xC ðt; xÞdx R s ¼ u0 ðtÞ  hxit ¼ xn 1  f2 C s ðt; xÞdx

ð6Þ

Observe that u(t) is the phase of Eq. (2). On the other hand the time marginal of Eq. (5) is given by:

PðtÞ ¼

Z

WDðt; xÞdx ¼ expð2fxn tÞ

ð7Þ

s  ðu  s=2Þsðu þ s=2Þ/ðh; sÞ

 expðjht  jsx þ jhuÞdudsdh

ð1Þ

where s(t) is any square integrable signal, and /(h, s) is the TFR kernel,  represents a complex conjugate. The desired TFR properties listed above are translated into kernel constrictions; making the TFR design equivalent to the kernel design. Note that the kernel can be a function of the signal, however in this work we consider signal-independent kernels which lead to bilinear TFRs. The term bilinear is used to denote that the signal enters only twice in the TFR calculation (see Eq. (1)); examples of bilinear TFRs can be found in Cohen (1995). 3. Decay ratio estimation using the TFA In this work, we use the well known impulse response approximation to the reactor signal. Using that approximation the reactor signal can be expressed as the real part of:

Derivating this last equation with respect to t and solving for f, the damping parameter, we have:

1 d ½logðPðtÞÞ 2xn dt

fðtÞ ¼ 

ð8Þ

Substituting Eq. (8) into Eq. (6) and solving for xn, we have:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 1 d xn ðtÞ ¼ ðu0 ðtÞÞ2 þ log PðtÞ 4 dt 0

ð9Þ

where u (t) is the instantaneous frequency (see Section 2). Here it is worth mentioning that in order to use Eq. (9), the WD(t, x) must satisfy the time marginal and the instantaneous frequency property (Cohen, 1995). Therefore, the instantaneous frequency can be calculated as the conditional frequency average for a given time, hxit (see previous section). In general, once that we have estimated the damping parameter, f(t), we can estimate the decay ratio (DR) as:

95

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

0

1

2pfðtÞ C B DRðtÞ ¼ exp @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA ð1  fðtÞ2

ð10Þ

Last, observe that given the signal of interest is s(t), we can calculate its WD. Then by using the WD, we can obtain the time marginal, 0 P(t), and the instantaneous frequency, u (t). This last parameter is calculated as the conditional mean, hxit, over a frequency range that contains the resonant frequency. From these two last parameters, we can calculate the natural resonant angular frequency, xn, using Eq. (9). Then, we can use Eq. (8) to estimate the damping parameter, f, and finally, we obtain the DR from Eq. (10). It is worth mentioning that the TFA will get an instantaneous decay ratio, IDR, that varies in time and closely reflects the following definition of an envelope based IDR (EIDR):

ACDRðtn ; t 0 Þ ¼

n Y

IDRðti ; t i1 Þ

ð13Þ

i¼1

We can rewrite Eq. (13) in a recursive way, i.e., if we have the ACDR at time tn1, we do not need all the n  1 multiplications in Eq. (13) to calculate the ACDR at time tn, we only need to do one multiplication:

ACDRðtn ; t 0 Þ ¼ IDRðt n ; tn1 ÞACDRðtn1 ; t0 Þ

Aðt 2 Þ EIDR ¼ Aðt 1 Þ

ð11Þ

where A(ti) is the instantaneous envelope of the nuclear signal at time ti. Observe that t2 > t1. Here it is worth mentioning that the EIDR and the IDR are two different ways to measure the DR. For the last term, we understand an IDR calculated based on TFA as it is described in this work. On the other hand, EIDR is defined by Eq. (11) as an envelope ratio. As we will see, the TFA based IDR closely follows the EIDR. For practical reasons we need to perform local averaging to declare the system stable or unstable and be able to take appropriate action to control the nuclear reactor. Let us consider what kind of average would be the most appropriate for our cause. If we take a second look at Eq. (11), we realize that due to the multiplicative action of the IDR, it is not the arithmetic mean that we should consider but the geometric mean. Based on this observation, we can estimate the signal envelope at a discrete time tn, using its value at the initial time t1 and the product of the EIDR factors that take us from t1 to tn1: n Y

EIDRðt i ; t i1 Þ

ð12Þ

i¼1

ACDRðt1 ; t 0 Þ ¼ IDRðt 1 ; t0 Þ

ð15Þ

where t0 is our time reference (zero time) and t1 is the next time point of interest, i.e., t i R ftn 8n ¼ 0; 1; 2; 3; . . .g. When using the ACDR we obtain results that confirm our intuition and knowledge about the physical phenomena going on the reactors. 4. Experiments and results First we use a synthetic signal to validate our method before taking the challenge of real signals. Among those, we include the Laguna Verde (LV) instability event and two BWR signals from the Forsmark stability benchmark: Averaged Power Range Monitors (APRMs) from cases 4 and 5. 4.1. Simulation using a synthetic signal The implementation of the model was made in MATLAB to estimate f for the signal given by Eq. (2) with the following parameters: A0 = 1, xn = 2.5p and f = 0.005 that represents a DR based on an autoregressive model of 0.9691. Based on these parameter

IDR Synthetic signal

1.04 1.03 1.02 1.01 1 0.99 0.98 0.97 0.96

0

20

40

ð14Þ

This is a recursive equation with the following initial value:

1.05

Amplitude

Aðt n Þ ¼ Aðt1 Þ

Therefore, due to the multiplicative nature of this last equation, we need to consider a geometric mean instead of the arithmetic mean. However, if we consider the geometric mean we would get an average DR but we are most interested in the maximum accumulated value over a certain period of time. Therefore, consider the accumulated value of the IDR (ACDR), from the first sample to the current one (from time zero to the current time):

60

80

Time [s] Fig. 1. IDR and synthetic signal (real part).

100

120

96

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

values, the IDR calculated using Eq. (10) with an instantaneous damping parameter, f(t), ranges from 0.97 to 1 and it is shown in Fig. 1. The accumulated decay ratio (ACDR) defined in Eqs. (13)–(15), for this case, is shown in Fig. 2. Observe that the ACDR behaves as expected decreasing with time, indicating that the envelope of the signal is decreasing (oscillatory phenomena decreases in amplitude as a function of time). Also observe that in Figs. 1 and 2, only the real part of Eq. (2) has been plotted. Notice that the sig-

nal has been scaled and an offset value has been added to it with purposes of comparison and visualization only.

4.2. Experiments and results with real BWR signals 4.2.1. Case of the instability event of Laguna Verde The proposed methodology was applied in a real instability event where the Averaged Power Range Monitor (APRM) signal

ACDR

1

Synthetic signal

0.8 0.6

Amplitude

0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0

20

40

60

80

100

120

Time [s] Fig. 2. ACDR and synthetic signal (real part).

40

38

Power (%)

36

34

32

30

28 0

100

200

300

400

500

Time [s] Fig. 3. Laguna Verde (LV) signal.

600

700

800

97

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

was obtained from the Laguna Verde Nuclear Power Plant (LVNPP) Unit 1. On January 24 of 1995 a power instability event occurred in LVNPP Unit 1, which is a BWR-5 and is operated since 1990 at a rated power of 1931 MWt. The instability event happened during a Cycle 4 power ascension without fuel damage. When the thermal power reached 37% of the rated power, the recirculation pumps were running at low speed driving 37.8% of the total core flow. The flow control valves were closed to their minimum position in order to shift the recirculation pumps at a high speed. The drop in drive flow resulted in a core flow reduction of 32% and a power reduction of 32%. Two control rods were also partially withdrawn during valve closure. The new low flow operating conditions re-

sulted in growing power oscillations (González et al., 1995). Fig. 3 shows the register of the APRMs, which was obtained via the Integral Information Process System (IIPS) with a sampling rate of 0.2 s for the whole time series. The figure shows that the power signal is non-stationary and the oscillations reached up to 10.5% of the maximum peak to peak amplitude. In order to get a clear interpretation of results, each signal is depicted as a function of time instead of indexing it using samples. For the sampling rate of 0.2 s, the time axis for this signal varies from 0 to 700 s. According to González et al. (1995) and Farawila et al. (1996) the channel ‘‘A” of the APRM trace shows no unstable behavior at 03:28:00 h. The valve closure was initiated at 03:28:20 h. A small core flow reduction was noticeable 40 s later, and the APRM-A

1.2 IDR LV signal

1.15

Amplitude

1.1

1.05

1

0.95

0.9

0.85 60

80

100

120

140

160

180

200

220

Time [s] Fig. 4. LV signal and IDR.

1.025 1.02 1.015

Amplitude

1.01 1.005 1 0.995 0.99 0.985

IDR LV signal

0.98 0.975

480

500

520

540

Time [s] Fig. 5. LV signal and IDR.

560

580

98

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

trace depicts signs of instability although the variations in the magnitude of the signal remained within the noise level. As the valve continued to close, the APRM-A trace shows clear unstable behavior starting at 03:30:30 h. The valve reached the minimum position at 03:31:30 h, and the oscillations continued without a significant increase in their growth rate. The valve reached the minimum position at 03:31:30 h, and the oscillations continued without any significant increase in their growth rate. The operator attempted to stabilize the power level by increasing the core flow opening the valves at 03:33:20 h. As a result of increasing the core flow, the oscillation started to decay at 03:34:40 h. At 03:35:20 h the oscillation reached 3% of amplitude, as a consequence, the reactor was manually scrammed to avoid core damage. Hence this is a very interesting case to validate our algorithm. Indeed, the LVNPP signal exhibits very clear trends (see the three marked segments in Fig. 3). At the beginning the LVNPP signal is

increasing very slowly, in the middle the amplitude increases much faster and in the end the amplitude decreases rapidly, since manual operator intervention was required to avoid a potential dangerous emergency at the LVNPP. We preprocess the LVNPP signal in the following way:  We get rid of any DC constant value and any linear constant trend by applying a high pass filter with cutoff frequency at 0.16 Hz.  We use a low pass filter with the appropriate cutoff frequency to only get the interest signal component related to the stability of the nuclear reactor.  Also in Figs. 4–6 for purposes of comparison only, we multiply the LVNPP signal by an appropriate scalar and add a DC value of one. This does not affect in any way the LVNPP signal but helps to visualize and to interpret our results.

IDR LV signal 1.1

Amplitude

1.05

1

0.95

0.9

0.85 620

630

640

650

660

670

680

690

700

710

720

Time [s] Fig. 6. LV signal and IDR.

1200 1000 800

Amplitude

600 400 200 0 -200 -400 LV signal ACDR

-600 -800

0

100

200

300

400

Time [s] Fig. 7. ACDR for the LV signal.

500

600

700

99

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

is decreasing locally even if later it continues its climb towards higher values. This climb continued until cautious operators decided to apply manual controls to avoid a potential dangerous situation at LVNPP. This action resulted in an envelope reduction, fact that is clearly picked up by the IDR measure (see Fig. 6). In conclusion we have an excellent measure of the IDR as opposed to classical methods that consider a constant DR over the analysis window length (usually 2000 samples at 12.5 Hz). In Fig. 7 it is depicted the ACDR for the LVNPP signal. Note that LVNPP signal has been scaled for purposes of visual comparison only.

Let us consider Fig. 4 where we have plotted a very slowly amplitude increasing section of the LVNPP signal. Observe that the IDR closely follows the signal amplitude variations reflecting closely the EIDR defined by Eq. (11). Note that the IDR is >1 when the signal envelope increases locally and it is <1 when the amplitude locally decreases (see Fig. 5). In that figure the amplitude of the LVNPP signal increases much faster than the signal section plotted in Fig. 4. As it was expected the IDR value is most of the time >1 showing that the signal envelope is increasing. Also observe that sometimes the IDR is <1, due that the IDR is a local measure and the signal envelope

10 4 ACDR Threshold for 50% power increase

Logarithmic Scale

10

10

10

10

10

3

2

1

0

-1

0

100

200

300

400

500

600

700

Time [s] Fig. 8. ACDR for LV signal in a logarithmic scale.

73

72

71

Power (%)

70

69

68

67

66

65

64

0

50

100

150

200

250

300

Time [s] Fig. 9. The APRM 1 signal from the case 4 of the Forsmark stability benchmark.

100

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

As expected the ACDR increases sharply and then it decreases when manual controls went off. In Fig. 8, the ACDR for the LVNPP signal is depicted in a logarithmic scale. Also observe that we have placed a threshold that indicates when the LVNPP signal increases its power in 50% with respect to time zero (the time reference).

tween a global oscillation mode and a regional (half core) oscillation. Results of the proposed method are shown in Figs. 10 and 11. Here it is worth mentioning that we preprocess the signal by applying a band pass filter to isolate the signal components of interest. Also as in the LVNPP case, we scaled the APRM signal and add a DC offset for purposes of facilitating the interpretation of results. From Fig. 10, we observe that the IDR obtained using TFA closely reflects the envelope changes in the filtered nuclear signal. Thus the TFA probes capable of reflecting the instantaneous changes in the DR as opposed to the well known AR modeling which considers a constant DR over the window analysis length, i.e., the DR obtained using an AR model is piecewise constant over a long period of time. However we can see from Fig. 10 that the

4.2.2. ‘‘Case 4” of Forsmark stability benchmark The TFA method was applied to neutronic signals of Forsmark stability benchmark, specifically to the case ‘‘c4.aprm” (Verdu et al., 2001). The benchmark is based on data from several measurements performed in the Swedish BWR reactor Forsmark 1 and 2, in the period 1989–1997. The ‘‘c4.aprm” signal corresponds to the average power range monitor (APRM) register during the stability event showed in Fig. 9. This case contains a mixture be-

IDR APRM 1

2.5 2 1.5

Amplitude

1 0.5 0 -0.5 -1 -1.5 -2 -2.5 0

10

20

30

40

50

60

70

80

90

Time [s] Fig. 10. IDR and filtered nuclear signal. IDR closely reflects the signal amplitude changes.

8

ACDR APRM1 C4 7

Amplitude

6

5

4

3

2

1

0

50

100

150

200

Time [s] Fig. 11. ACDR for the APRM C4 signal.

250

300

100

101

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

signal envelope is better modeled by a Gaussian-like function since first goes up and then it decays. Indeed, from Fig. 10 the system produces a continuous rise in the amplitude signal followed by a decay that will eventually go to zero if the system was left in repose (zero inputs). Here it is worth mentioning that the DR can be measured on the output of the system for any given excitation, there is no need to know the system transfer function or its excitation. From Fig. 11 we observe that there are no major trends in the ACDR measure that indicate the onset of an unstable state. Therefore, the IDR as well as the ACDR measures are oscillating around one, reflecting the natural rise and fall of the signal envelope. Thus we could argue that the reactor is in a stationary state and, that the

ACDR can be best used to detect major changes in the nuclear reactor stability. A long term steady increase in the ACDR will indicate the onset of an unstable state in the nuclear reactor. 4.2.3. ‘‘Case 5” of Forsmark stability benchmark Last, the TFA was applied to the neutronic signal of Forsmark stability benchmark, specifically to ‘‘case 5” (Verdu et al., 2001). The signal corresponds to the APRM registered during the instability event, which is shown in Fig. 12. The time axis, considering the sampling time of 0.08 s, goes from 0 to 320 s. This event corresponds to a situation where the neutronic power of the reactor displays abnormal and apparently unstable oscillations. That is, it

85

Power (%)

80

75

70

65

60

55 0

50

100

150

200

250

300

Time [s] Fig. 12. The APRM 1 signal from the case 5 of the Forsmark stability benchmark.

25

ACDR APRM1 C5 20

Amplitude

15

10

5

0

-5

-10

-15

0

50

100

150

200

Time [s] Fig. 13. ACDR for the APRM C5 signal.

250

300

350

102

J.E. Torres-Fernández et al. / Annals of Nuclear Energy 37 (2010) 93–102

seems that a medium-frequency unstable oscillation dominates over the stable high-frequency oscillations usually found under normal operation conditions (see Fig. 12). Thus the corresponding time series represent a non-stationary case where the autoregressive methods have a limited validity according with Verdu et al. (2001). The reported plant conditions corresponding to the event were: 64.7% of the rated power and 4014 kg/s. We also preprocess this signal by applying a band pass filter to isolate the interest components from the signal and, for purposes of comparison only, we scale adequately the APRM signal. The ACDR and the signal for this case are plotted in Fig. 13. In this case, at the beginning the ACDR value is oscillating around one without showing a clear increase or decrease trend. Thus we could argue the reactor is at a stationary state. However at 80 s the ACDR starts showing an increase in the magnitude of the oscillations going clearly above a DR = 1. Thus we can argue the reactor is becoming unstable. This tendency continues for around 70 s, and then the ACDR steadily decreases until reactor stability is achieved. Therefore the ACDR measure clearly reflects our understanding of this case. 5. Conclusions In this work, we propose a novel approach to analyze the stability of BWRs based on TFRs. For mathematical convenience the WD is used but any bilinear TFRs that satisfy the marginals and the instantaneous frequency property could have been used instead of the WD. Also notice that non-bilinear TFRs can be employed but they involve signal-dependent kernel design as opposed to the signal-independent kernels used by all bilinear TFRs. Thus if a bilinear TFR is used, there is no need to do a dedicated signal design kernel for every signal to be analyzed. Besides that significant results were obtained using bilinear TFRs in this work. Using TFRs to estimate the decay ratio (DR) allows to track its variations in time as opposed to commonly signal processing techniques which consider a constant DR over a period of time (analysis window duration), i.e., a piecewise constant DR is obtained over a long period of time. Hence we can only talk about stability on the time scale of the analysis window function. However a DR based on TFRs is a function of time, we would rather refer to it as an IDR to differentiate it from a constant window DR. Indeed the IDR closely follows the signal envelope variations making it an excellent local DR measure. Also it makes possible to talk about local signal stability due to its local nature. Despite the excellent local properties of the IDR, we need a measure that allows us to track the long term stability trend in the BWRs. This is achieved by the ACDR measure defined as the accumulative multiplication of the local IDRs (see Eqs. (13)–(15)). As it is shown in this work the ACDR closely follows the long term stability trends in the BWRs, especially in the LVNPP case where the signal stability trends are very clear (see Fig. 3). Regarding the APRMs cases 4 and 5 of the Forsmark stability benchmark, the ACDR clearly reflects major changes in the signals amplitude indicating instability, specially for the case 5. Also in these cases,

the oscillation of the ACDR around one indicates stationary states in the BWRs. By stationary states we mean that there is not a long term increasing or decreasing trend in the signal. Therefore the ACDR closely reflects the long term trends in the BWRs stability, while the IDR is a local measure of the DR. References Baraniuk, R.G., Jones, D.L., 1993. A signal-dependent time–frequency representation: optimal kernel design. IEEE Trans. Signal Process. 41, 1589– 1602. Blazquez, J., Ruiz, J., 2003. The Laguna Verde BWR/5 instability event. Prog. Nucl. Energy 43, 195–200. Cohen, L., 1995. Time Frequency Analysis. Prentice Hall, Englewood Cliffs, NJ. Daubechies, I., 1992. Ten Lectures on Wavelets. SIAM, Philadelphia. Espinosa-Paredes, G., Prieto-Guerrero, A., Nuñez-Carrera, A., Amador-García, R., 2005. Wavelet-based method for instability analysis in boiling water reactors. Nucl. Technol. 151, 250–260. Farawila, Y.M., Pruitt, D.W., Smith, P.E., Sánchez, L., Fuentes, L.P., 1996. Analysis of the Laguna Verde instability event. In: Proc. National Heat Transfer Conf., Houston, Texas, August 3–6, vol. 9. American Nuclear Society, pp. 198– 202. Gabor, D., 1946. Theory of communication. J. IEE 93, 429–457. González, V.M., Amador, R., Castillo, T., 1995. Análisis del evento de oscilaciones de potencia en la CNLV: Informe Preliminar. CNSNS-TR-13, REVISION 0. Comisión Nacional de Seguridad Nuclear y Salvaguardias, MEXICO. Guido, G., Converti, J., Clause, A., 1997. Density wave oscillations in parallel channels, an analytical approach. Nucl. Eng. Des. 125, 121–139. Jeong, J., Williams, W.J., 1992. Kernel design for reduced interference distributions. IEEE Trans. Signal Process. 40, 402–412. Lahey, R.T., Podoswski, M.Z., 1989. On the analysis of various instabilities in twophase flow. Multiphase Sci. Technol. 4, 183–371. March-Leuba, J., 1986. A reduced-order model of boiling water reactor linear dynamics. Nucl. Technol. 75, 15–22. March-Leuba, J., 1990. LAPUR Benchmark Against In-Phase and Out-of-Phase Stability Test. NUREG/CR-5605, ORNL/TM-11621. March-Leuba, J., Blakeman, E.D., 1991. A mechanism for out-of-phase power instabilities in boiling water reactors. Nucl. Sci. Eng. 107, 173–179. March-Leuba, J., Rey, J.M., 1993. Coupled thermal-hydraulic neutronic stabilities, in boiling water reactor a review of the state of the art. Nucl. Eng. Des. 145, 97– 111. Muñoz-Cobo, J.L., Chiva, S., Sekhri, A., 2004. A reduced order model of BWR dynamics with subcooled boiling and modal kinetics: application to out of phase oscillations. Ann. Nucl. Energy 31, 1135. Navarro-Esbri, J., Ginestar, D., Verdu, G., 2003. Time dependence of linear stability parameters of a BWR. Prog. Nucl. Energy 43, 187. Peng, S.J., Podowski, M.Z., Lahey, R.T., Becker, M., 1984. NUFREQ-NP: a computer code for the instability analysis of boiling water nuclear reactors. Nucl. Sci. Eng. 88, 404–411. Podowski, M.Z., Pinheiro, M., 1997. Modeling and numerical simulation of oscillatory two phase flow with application to boiling water nuclear reactor. Nucl. Eng. Des. 177, 179–188. Saha, P., Zuber, N., 1978. An analytical study of the thermally induced two-phase flow instabilities including the effects of thermal non-equilibrium. Int. J. Heat Mass Transfer 21, 415–426. Sunde, C., Pázsit, I., 2007. Wavelet techniques for the determination of the decay ratio in boiling water reactors. Kerntechnik 72, 7–19. Turso, J.A., March-Leuba, J., Edwards, M., 1997. A modal-based reduced-order model of BWR out-of-phase instabilities. Ann. Nucl. Energy 12, 921–934. Uehiro, M., Rao, Y.F., Fukuda, K., 1996. Linear stability analysis on instabilities of inphase and out-of-phase modes in boiling water reactors. J. Nucl. Sci. Technol. 33, 628–635. Van Der Hagen, T.H.J.J., 1998. Fuel heat transfer modelling in reduced-order boiling water reactor dynamics models. Ann. Nucl. Energy 25, 1287–1300. Verdu, G., Ginestar, D., Muñoz-Cobo, J.L., Navarro-Esbri, J., Palomo, M.J., Lansaker, P., Conde, J.M., Recio, M., Sartori, E., 2001. Forsmark 1 and 2 Stability Benchmark. Time Series Analysis Methods for Oscillations During BWR Operation. Final Report, NEA/NSC/DOC. 2.