© A.P. Daga et al., Published by EDP Sciences 2022

1 Introduction

Hydropower generation units (HGUs) are electromechanical systems meant to transform the potential energy of the flowing water (i.e., a renewable energy source) into electrical energy. In particular, a turbine is used to turn the hydraulic energy into mechanical energy, while a generator finally performs the transformation into electrical energy. The power that can be obtained by the flowing water depends on its hydraulic head and its volumetric flow rate [1] as: (1)

where is the rated power, η is the efficiency of the system, ρ is the water density, g is the acceleration of gravity, h is the difference in height between the water inlet and outlet (i.e., the head) and is the volumetric flow rate.

In general, two main groups of hydropower plants can be distinguished. High-flow, low-head power plants are common for run-of-the-river stations usually involving Kaplan or Francis turbines. Such plants need a continuous water flow and therefore have less ability to provide power on demand. On the contrary, low-flow, high-head plants usually exploit high altitude dams/reservoirs combined to (Kaplan or) Pelton turbines, which allow quick start-up and high manoeuvrability [2].

This peculiarity makes the HGUs mission-critical assets for grid operators (e.g., TERNA for the Italian high voltage electricity grid), which often have the need of compensating sudden increases in the electricity demand [3]. Furthermore, the global energy policy is pushing for an ever more ecological and healthier energy production, so that it is reasonable to believe that in the future HGUs will gain further and renewed interest [4,5]. Hydropower, in fact, despite few concerns about hydro-morphological alterations and aquatic ecosystems issues, which could be easily solved adopting new technologies and practices, remains a key renewable energy source [5].

Condition monitoring of such HGUs becomes then a fundamental task for optimizing their maintenance regime, allowing both safety and economic advantage at the same time.

In this regard, this work focuses on a real hydropower plant of the low-flow/high-head kind located in Valle d’Aosta (Italy), whose HGU features a SCADA (i.e., a Supervisory Control And Data Acquisition) system for vibration monitoring compliant to the ISO20816-5:2018 latest standard. The aim of the work is that of proposing an improvement of the available vibration monitoring system based on orbital analysis, which will be introduced in Section 3. A description of the particular power plant on the contrary, is given in the next section.

2 The Signayes hydroelectric power plant from C.V.A. S.p.A.

The Signayes plant from C.V.A. S.p.A. (Compagnia Valdostana delle Acque − Compagnie Valdôtaine des Eaux) is a hydroelectric power plant located upstream of Aosta at 604 MSL and uses water from Gran San Bernardo and Valpelline valleys in Valle d’Aosta, Italy.

The power plant is built in a cave at the bottom of a 120m well and features a gross head of 351 m, a nominal flow rate of 16 m³_/s and a nominal power of 42 MW. It is capable of producing an average of 220 GWh/year. The Signayes hydropower plant is shown in Figure 1. Object of the present work is the first HGU.

The Signayes HGU is made by a horizontal shaft with twin Pelton turbines featuring z_b = 20 buckets, supported by two pedestal bearings on rigid foundations. The left bearing is able to give both radial and axial thrust (i.e., thrust bearing), while the right one can only carry radial load (i.e., guide bearing). The Signayes HGU features a SCADA system for vibration monitoring compliant to the ISO20816-1:2016 and the ISO20816-5:2018 latest standard [6,7]. In particular, as shown in Figure 2, it is equipped with absolute bearing housing vibration sensors on both the supports for monitoring the root-mean-square (RMS) vibration velocity in mm/s (i.e., integrated accelerometers signal) and relative shaft vibration sensors (i.e., two proximity 90° apart) on the right support for acquiring the peak-to-peak shaft displacement in μm. The location of the sensors complies with the ISO standard and is shown schematically in Figure 3. Further information is given in the next section.

Auxiliary data such as the produced power (i.e., the load), the rotational speed and the oil temperatures are also recorded by C.V.A. S.p.A.

Fig. 1 A picture of the Signayes Hydropower Plant from C.V.A. S.p.A.

Fig. 2 A picture of the HGU bearing housing with sensors.

Fig. 3 Sensors’ location for absolute bearing housing vibration (a) and for relative shaft vibration (b) according to ISO20816-5:2018 latest standard [7].

2.1 Simple HGU bearing hydro-dynamic model

In order to validate the here proposed analysis, a simple dynamic model of the shaft motion within the journal bearing was created on the basis of the real HGU system at nominal work conditions.

Under the assumption of infinitely long journal bearing, the shaft motion can be proved to generate a load-supporting pressure profile in the (incompressible) lubricant film described by the Reynolds equation [8]: (2)

where p is the pressure profile, μ is the dynamic viscosity, r is the nominal radial dimension of the shaft, ω is the (stationary) angular shaft speed and h is the lubricant film thickness, which can be described in terms of the relative angular coordinate ϑ, highlighted in Figure 4: (3)

where c is the radial shaft-support clearance and ε (t) is the absolute value of the displacement of the centre of the shaft with respect to the centre of the support, described by the vector .

Once equations (2) and (3) are established, it is easy to produce a numerical scheme to double-integrate the pressure profile and compute the supporting forces from the lubricant film. In particular, using the relaxation method [9], derivatives are replaced by finite difference approximations leading to a linear system of equations which can be solved considering Full-Sommerfeld periodic boundary conditions (i.e., null initial and final pressures). Nevertheless, the Full-Sommerfeld solution does not consider the occurrence of cavitation, so that the Half-Sommerfeld approach will be used in this simplified model (historically accepted journal bearing design approximation [10]).

Integrating the so-computed Half-Sommerfeld hydro-dynamic pressure profile on the shaft surface (i.e., summing all the in Figure 4 over the whole angular coordinate ϑ), it is straightforward to find the overall pressure force at present time step . Such a force, summed to the reaction of the support on the shaft (i.e., ), is supposed to balance the inertia forces, so that the eccentricity at the next time step can be found by finite difference approximation of the second order differential equation of motion obtained from the free body diagram in Figure 5 under the assumption of half the shaft mass M concentrated in its centre of gravity: (4)

Despite the strong assumptions and simplifications, the here described model can be used to derive theoretical orbits for nominal working conditions. For example, considering an additional harmonic contribution synchronous to the rotation speed ω, the shaft centre point will move according to the trajectory reported in Figure 6. Assuming a defective operation featured by an increase in the static unbalance , it is also possible to derive an unhealthy orbit to be compared to the previously found reference one.

From Figure 6, several considerations arise:

• The model, despite simple, is able to produce physically consistent, stable orbits.

• The order-one orbits induced by synchronous perturbations of the equilibrium position can be reliably modelled as ellipses.

• Increasing the perturbation amplitude, the orbit centre is shifted, the orbit shape may slightly deviate from the perfect ellipse, while the orbit size and orientation undergo variations.

Fig. 4 Graphical problem statement and definition of the reference frames for the bearing hydro-dynamic model.

Fig. 5 Free body diagram of the half-shaft.

Fig. 6 Orbits due to a synchronous harmonic forcing − the black dots represent truly elliptical orbits.

3 ISO20816-based signal processing and the proposed improvement

The fundamental scope of vibration monitoring is that of inferring the state of health of a machine given the available acquisitions. In the simplest case, the problem simplifies to a binary classification: a “healthy” or a “damaged” label is to be assigned to a chunk of data. Indeed, it is common for industrial continuous monitoring equipment to store data in temporary buffers, aggregate the data by extracting the selected feature (e.g., the RMS) and compare the so computed feature to a given threshold which separates the “healthy” from the “anomalous” levels. Usually, in normal conditions, just the aggregated information is stored (i.e., the feature), while the vibration signal over time is saved only when some level is exceeded, to allow higher level analyses (e.g., spectra, envelopes, etc.) [11]. In principle, this corresponds to the setting of a threshold through a hypothesis testing for normality [12], which can be considered a Novelty Detection (or Anomaly Detection) task, as “Novelty” (i.e., a deviation from the normal behaviour), can be used as a damage indicator in absence of confounding influences [13–17].

In this regard, the ISO standard [7] actually identifies four different zones for the assessment of a machine health state:

• Zone A and B: the machine is new or anyway in acceptable condition.

• Zone C: the machine condition is unsatisfactory.

• Zone D: the vibration level is anomalous and may lead to damage.

From zone A&B to zone C a first action limit is defined. An alarm should be then triggered upon exceeding this threshold, and further analysis is recommended (i.e., spectral analysis to identify unbalance, misalignment, oil instability etc., shaft orbit analysis, shaft centreline shift analysis, bearing temperatures check, visual inspection, etc.).

A second action limit, on the contrary, separates zone C from zone D. When a threshold of 1 ÷ 1.25 times this second limit is reached, a trip system is supposed to either reduce the load or shut down the machine, while a complete visual inspection and further checks are required (i.e., tightness of all the mounting bolts, stress evaluation, modal analysis, Non-Destructive Testing, etc.).

It is important to point out that the ISO standard, in its annex A, derived precise action limits for all the main HGU types, based on a statistical analysis of an international vibration database containing more than 7000 datasets of hydraulic machines all over the world. A Burr distribution was fitted on the data, and the 75-th and 90-th percentile critical values were computed to be used as action limits. For a horizontal Pelton HGU with pedestal bearings, for example (i.e., the Signayes case), the generator drive-end bearing action limits are 1, 8mm/s and 2,9mm/s for the RMS velocity of the absolute vibration and 145μm and 225μm for the peak-to-peak relative shaft displacement.

These values hold for an HGU at 100% of the rated power, with bearing and oil temperatures stabilized at normal values (i.e., roughly around 50° for a conventional guide Babbitted bearing and 70° for a conventional thrust Babbitted bearing). Load, speed and temperature are in fact the main confounding influences for a correct diagnostic, and their variation should be taken into account with a more complex signal processing or simply avoided by comparing only pre-filtered data at nominal, stable values for the confounders.

It is important to point out that, according to the ISO standard, an alarm should be triggered even in case of a sudden change of at least 25% of the normal value (i.e., a stable value for at least three months prior to the detected change), even if the first action limit is not reached. Hence it can be useful to keep track of the trends of the average vibration levels.

Finally, the standard also prescribes the main acquisition parameters:

• The acquisitions should represent a minimum of 10 revolutions of the machine shaft

• The frequency range of interest is 0,1 ÷ 3 ⋅ z_b times the rotational frequency, where z_b is the number of blades of the runner (i.e., the number of buckets for a Pelton turbine)

• The sampling frequency should be at least 2,56 times the maximum frequency. Higher sampling rates are advisable, but if it is reasonable to believe that the high frequency components cannot induce significant stress levels, it is suggested to filter them out (e.g., with a low-pass filter).

Hence, it is clear that the maximum rotational frequency is the main design parameter for the SCADA acquisition system. According to [5], the HGUs rotational speed is synchronous with the grid frequency, as the speed governor works to control the water flow so that the produced power matches the need of the grid.

The Signayes plant works at a nominal speed of 500 rpm (i.e., 8.33 Hz), hence, every second of acquisition contains 8 full revolutions. The acquisitions duration was then set to three seconds. Also, the frequency range of interest corresponds to 0.8 ÷ 500Hz, so that the selected sampling frequency was 2560Hz (i.e., twice the 2.56 ⋅ 500Hz requirement).

3.1 Weaknesses and possible improvements

As it is clear from the previous section, the main features acknowledged by the ISO standard are the RMS velocity of the absolute vibration (in mm/s) and the peak-to-peak relative shaft displacement (in μm). No other feature is accepted by the standard, as in section B.2 of [7] it is said that, despite the various claims made in support of features different from the RMS (i.e., the peak value, the crest factor, etc.), no one can be trusted in all the conditions. Even if one feature outperforms the others in some cases, in fact, no general rule can be derived, so that the RMS-velocity of relative vibration and the peak-to-peak absolute vibration remain the reference features. In this regard, a multivariate Novelty Detection able to merge the information from multiple features was proposed by the authors in [11–17] to improve the basic univariate approach and proved successful in many different applications.

However, on exceeding the first action limit, spectral analysis is suggested to identify the cause of an abnormal dynamic behaviour in the absolute vibration, while orbit analysis is proposed for studying anomalies of the relative vibration. The standard, unfortunately, does not elaborate on these diagnostic tools, but focusing on [18], it appears clear that such analyses are founded on additional features:

• Frequencies of the main relative vibration components and vibration amplitudes relative to the first orders

• Orbital paths and orbits of absolute vibration:

• Amplitude and Phase angle of the orbits.

• Shaft centreline position.

Similar considerations arise in [19–24], where orbit shape (i.e., circular vs elliptical), directivity (i.e., forward or backward) and inclination (i.e., angle of the major axis of the ellipse with respect to one of the sensors) are recognized as relevant indicators for relative vibration analysis, as defects or faults can cause changes to the orbits [19].

Hence, the goal of the present paper is to widen the ISO-suggested normal monitoring process (i.e., based only on RMS velocity and peak-to-peak displacement), including additional features from spectral analysis of relative vibration and from orbital analysis of absolute vibration. Unfortunately, if it is easy to find the main spectral contributions (e.g., with a simple Fast Fourier Transform (FFT)), this is not the same for orbit-related features. In this regard, the present work is meant to propose a methodology, which is described in the next section, based on a fitting of the filtered orbits via Principal Component Analysis.

3.2 Orbit-related features and the proposed extraction methodology

As specified in [19,22] and proven by the here proposed simulations in Section 2.1, defects in rotating machines cause a change in the filtered orbits, which in general can be considered as Lissajous figures (i.e., parametric equations describing complex harmonic motion). When the x (t) and y (t) displacements from the two orthogonal probes are considered, if the signals are band-filtered around a particular harmonic (i.e., 1x, 2x of the rotational frequency are of interest), the resulting orbit can be generally assumed elliptic.

The main features describing an elliptic shape are:

• The first semi-axis r_x;

• The second semi-axis r_y;

• The inclination of the first semi-axis from the x axis (i.e., the angle ϑ).

Indeed, the standard parametric equation of a discrete rotated ellipse can be written as a particular Lissajous curve rotated by an angle ϑ (see annex A for derivation): (5)

where n = 0 ÷ N − 1 is the sample index, while N is the total number of samples.

Anyway, it is common to summarize the two semi-axes r_x and r_y using the major semi-axis R and a measure of eccentricity e (i.e., a measure of dissimilarity of the elliptic shape to a circular shape, answering the question “how far is the ellipse from a perfect circle?”) depending on R = max(r_x, r_y) and r = min(r_x, y_y) as: (6)

The aim is then to find a suitable estimate of R, e and ϑ given a noise-affected measure of the trajectory .

In order to perform the function fitting, a Principal Component Analysis (PCA) is proposed in this work. A mathematical investigation of PCA applied to a theoretical (i.e., noise free) elliptic trajectory is derived in annex A to prove the effectiveness of the methodology. In the demonstration, the covariance matrix S of the discrete trajectory of one entire revolution is defined given the triplet of variables r_x, r_y and ϑ. According to PCA [25], the eigenproblem Sv = λv (where S is a 2 × 2 covariance matrix, v is a 2 × 1 eigenvector and λ is a scalar eigenvalue) is then solved to prove that the eigenvalues correspond to half the squared semi-axes, while the inclination angle of the first eigenvector is exactly the inclination angle ϑ of the first semi-axis from the x axis: (7)

This perfectly works for the theoretical discrete trajectory of an entire orbit (i.e., Eq. (2)). Anyway, PCA can easily deal with noisy acquisitions and multiple revolutions to find the best fit, proving to be the ideal tool for estimating the triplet of orbital features R, e and ϑ from measured data. This shown also in Figures 7 and 8, on a synthetic orbit, and proved via Monte Carlo repetition in the next section.

Fig. 7 PCA on ideal elliptic orbit (synthetic data).

Fig. 8 PCA on noisy elliptic orbit (synthetic, noise std = 0.15).

3.3 Validation via Monte Carlo repetition

In order to test the robustness of the here proposed PCA ellipse fitting, a Monte Carlo simulation was set up. Normal random noise was added to a theoretical orbit having r_x = 5.1μm, r_y = 2.8μm and ϑ = 15°, similar to the ones derived in Section 2.1. The standard deviation of the random normal noise was set to vary in a range 0 ÷ 2. A thousand repetitions for increasing values of the noise standard deviation were run. The average results are shown in Figure 9, from which the algorithm robustness can be appreciated, as it is able to find satisfactory estimates (i.e., with errors lower than 5%) even if the orbit is corrupted by a noise with SNR up to 4.5 (i.e., roughly 2.8/0.6 for r_y and 5.1/1.1 for r_x).

Fig. 9 Monte Carlo repetition for noisy elliptic orbit: PCA estimated rx, ry and ϑ Vs random normal noise standard deviation. In red, a 5% threshold for the estimation error.

4 Results and discussion

In order to increase the confidence in the here proposed algorithm, this was tested both on modelled orbits (Sect. 2.1) and on real acquisitions from the Signayes HGU.

The PCA fitting of orbits modelled according to Section 2.1 led to the results reported in Figure 10. In particular, it can be noticed that:

• the PCA ellipse-fitting algorithm is effective in summarizing the entire orbit with just three features (i.e., R, e and ϑ),

• in stationary working conditions, (i.e., when the only possible source of variability is damage), the so produced features are effective in recognizing damage.

Real acquisitions from the Signayes HGU were then treated with the here proposed algorithm. At first, a FIR filter was optimized to extract a single order from the absolute vibration channels, here named x and y. An order of 500 samples was selected for the bandpass filter, whose spectral properties are shown in Figure 11. In Figure 12, the resulting filtered time signals are shown together with their composition in a 2D plane where the orbit can be visualized.

The proposed algorithm was then run on such filtered data to automatically identify the ellipse-characteristic parameters R, e and ϑ. Repeating the PCA ellipse-fitting algorithm over time, the stability of the 1× and 2× orbits and of their estimated characteristic parameters encourage their use as diagnostic features. This is highlighted in Figure 13, where the results of several months of monitoring are shown. The here considered measurements refer to nominal, regimated conditions (i.e., nominal speed, nominal bearing temperatures, almost nominal power load), so as to limit the confounding influences. Figure 13 shows that the estimates of R, e and ϑ are stationary (i.e., wide sense stationarity), and apart from some natural variability, they feature a constant average value and variance. This holds also when monthly data are aggregated (see Tab. 1) highlighting the absence of long-term trends.

Fig. 10 Modelled orbits and their elliptical fits (dotted).

Fig. 11 Amplitude spectra of the x and y relative displacements with highlighted the selected filter for the 2× orbit.

Fig. 12 PCA on real signal from Signayes HGU (nominal, regimated condition) filtered around the second order. The 2× Orbit is visualized in red and features an eccentricity of almost 1. The major semi-axis (i.e., the first principal component) is represented in yellow and is inclined of 138.9°.

Fig. 13 Time evolution of the PCA estimates of R, e and ϑ on real signal from Signayes HGU (nominal, regimated condition). In blue the 1× orbit, in red the 2× orbit.

5 Conclusions

The present work proposes the integration of higher-level features in ISO20816 HGUs vibration monitoring systems with minimum computational effort. Standard vibration monitoring SCADA, indeed, are commonly equipped with accelerometers for recording absolute bearing housing vibration and with proximity sensors for measuring the relative shaft vibration. Hence, in order to take full advantage of the sensors, it can be wise not to focus only on standard features such as the RMS-velocity of relative vibration and the peak-to-peak absolute vibration, but to include also additional features from spectral analysis and orbit analysis, which usually are reserved for offline additional monitoring if particular action limits are exceeded. Anyway, if it is true that the main spectral contributions can be simply highlighted by tried-and-tested algorithms (i.e., FFT), this is not the same for orbit-related features.

It is then the scope of this paper to propose a simple yet effective algorithm for extracting orbit-characteristic features from filtered, relative vibration signals (i.e., bandpass filtered around the first spectral orders). In particular, a PCA ellipse fitting algorithm is proposed in this work to identify parameters such as the major semi-axis (i.e., the first principal component R), its inclination (i.e., the angle ϑ) and the eccentricity e of the filtered orbit.

The parameters identification was tested both on synthetic data (without and with added noise) that on real acquisitions from the HGU of the Signayes hydropower plant located in Valle d’Aosta (Italy). The analysis proved the reliability and stability of the indicators, which are known to be very sensitive to damage or faults that can cause changes to the orbits (as reported in [19,22] and proven by the here proposed hydro-dynamic simulations in Sect. 2.1). Given the low computational burden, such an analysis could be easily implemented online, fostering the monitoring and diagnostics of HGUs.

Future works may involve different filtering strategies, different fitting strategies, as well as the integration of different sensors (e.g., temperature, absolute and relative vibration sensors fusion or machine vision [26–28]).

Acknowledgments

Thanks to C.V.A. S.p.A. for their availability and support.

Annex A PCA for ellipse fitting

• 1)

  Definition of the trajectory of an entire orbit having the shape of a rotated ellipse, starting from the equation of a discrete parametric equation of a Lissajous figure: (A1)

  where n = 0 ÷ N − 1 is the sample index, while N is the total number of samples.

  In matrix notation, using the symbol “’_” for transpose: (A2)

• 2)

  Counter-clockwise rotation of an angle ϑ in the range 0 ÷ π/2 while setting : (A3) (A4)

• 3)

  Arranging the trajectory in a single data-matrix and solving for the corresponding covariance matrix S: (A5) (A6)

  Remembering the Euler identities: (A7)

  it is possible to derive: (A8) (A9) (A10)

  Considering that the exponential sum has a well-known solution: (A11)

  if x = 2π/N, i.e., an entire orbit is taken, the exponential sums are all null, so that the covariance S formulation simplifies to: (A12)

• 4)

  Eigendecomposing the covariance matrix S, it is possible to write: (A13)

  with: (A14)

• 5)

  Eigenvalues can be found given: (A15)from which a characteristic equation can be easily derived: (A16)leading to: (A17)

• 6)

  Eigenvectors can be found as: (A18)

The system of equations from an eigenvalue/eigenvector pair is underdetermined, so that, if the first pair is considered, the following identity can be found: (A19)

So that: (A20)

References

All Tables

All Figures

	Fig. 1 A picture of the Signayes Hydropower Plant from C.V.A. S.p.A.
In the text

	Fig. 2 A picture of the HGU bearing housing with sensors.
In the text

	Fig. 3 Sensors’ location for absolute bearing housing vibration (a) and for relative shaft vibration (b) according to ISO20816-5:2018 latest standard [7].
In the text

	Fig. 4 Graphical problem statement and definition of the reference frames for the bearing hydro-dynamic model.
In the text

	Fig. 5 Free body diagram of the half-shaft.
In the text

	Fig. 6 Orbits due to a synchronous harmonic forcing − the black dots represent truly elliptical orbits.
In the text

	Fig. 7 PCA on ideal elliptic orbit (synthetic data).
In the text

	Fig. 8 PCA on noisy elliptic orbit (synthetic, noise std = 0.15).
In the text

	Fig. 9 Monte Carlo repetition for noisy elliptic orbit: PCA estimated rx, ry and ϑ Vs random normal noise standard deviation. In red, a 5% threshold for the estimation error.
In the text

	Fig. 10 Modelled orbits and their elliptical fits (dotted).
In the text

	Fig. 11 Amplitude spectra of the x and y relative displacements with highlighted the selected filter for the 2× orbit.
In the text

	Fig. 12 PCA on real signal from Signayes HGU (nominal, regimated condition) filtered around the second order. The 2× Orbit is visualized in red and features an eccentricity of almost 1. The major semi-axis (i.e., the first principal component) is represented in yellow and is inclined of 138.9°.
In the text

	Fig. 13 Time evolution of the PCA estimates of R, e and ϑ on real signal from Signayes HGU (nominal, regimated condition). In blue the 1× orbit, in red the 2× orbit.
In the text