Spatial Harmonic Decomposition as a tool for unsteady flow phenomena analysis

Hydropower is already the largest single renewable electricity source today but its further development will face new deployment constraints such as large-scale projects in emerging economies and the growth of intermittent renewable energy technologies. The potential role of hydropower as a grid stabilizer leads to operating hydro power plants in "off-design" zones. As a result, new methods of analyzing associated unsteady phenomena are needed to improve the design of hydraulic turbines. The key idea of the development is to compute a spatial description of a phenomenon by using a combination from several sensor signals. The spatial harmonic decomposition (SHD) extends the concept of so-called synchronous and asynchronous pulsations by projecting sensor signals on a linearly independent set of a modal scheme. This mathematical approach is very generic as it can be applied on any linear distribution of a scalar quantity defined on a closed curve. After a mathematical description of SHD, this paper will discuss the impact of instrumentation and provide tools to understand SHD signals. Then, as an example of a practical application, SHD is applied on a model test measurement in order to capture and describe dynamic pressure fields. Particularly, the spatial description of the phenomena provides new tools to separate the part of pressure fluctuations that contribute to output power instability or mechanical stresses. The study of the machine stability in partial load operating range in turbine mode or the comparison between the gap pressure field and radial thrust behavior during turbine brake operation are both relevant illustrations of SHD contribution.


Introduction
This article presents a mathematical method that combines information coming from a set of sensors in order to give a spatial description of the measurement. This method called spatial harmonic decomposition is first exposed from a theorical point of view. Then, two examples will be developed. The first one explores load influence on draft tube pulsations in a Francis turbine. The second one links spatial pattern of vaneless area pressure pulsations to the increase in radial thrust during turbine brake operation.

Spatial Patterns
Standard frequency analysis gives an idea of the time coherence of the phenomena but no information about the spatial coherence. Thus, a distinction between asynchronous and synchronous pulsations has long been introduced by Nishi et al. in [4] to describe pressure field around a draft tube cross section. Synchronous pulsations are seen in phase from one sensor to another whereas asynchronous pulsations are not. According to this definition, a rotating pressure field gives rise to asynchronous pulsations. A simple mathematical model of a rotating pattern p rot (α, t) = cos (ωt − kα + ϕ rot ) depending on the spatial azimuth α combining with a synchronous pattern p sync (t) = cos (ωt + ϕ sync ) at the same pulsation ω is: p tot (α, t) = p sync (t) + p rot (α, t) = real A sync e iϕsync + A rot e i(−kα+ϕrot) e iωt = real A tot (α)e iϕtot(α) (1) This helps to explain why such a combination with the same frequency leads to a deviation in the amplitudes between sensors (A tot (α 1 ) = A tot (α 2 )) and phase lags that are no more directly related to the geometrical angular distance between sensors (ϕ tot (α 1 ) − ϕ tot (α 2 ) = α 1 − α 2 ).

Short review of existing decomposition methods
A least three methods have been proposed to separate synchronous from rotating patterns: (i) With two sensors seperated by 90 • , Nishi et al. still in [4] proposed to extract the rotating part assuming the equation p rot (π/2, t) = p rot (0, t − T /4), where T is the rotation period. (ii) Angelico and Muciaccia in [2] proposed a method called Vectorial Decomposition that takes advantage of the complex vectorial representation as illutrated in Eq. 1. Indeed the signal at one given frequency may be seen as vectors in the complex plane whose absolute values are the amplitude of the signal and whose arguments are their phase lag to a reference. In this complex plane, a pure rotating pattern is represented on different sensors by vectors describing a circle whereas a plane wave is a constant vector independent of the azimuthal position. Thus, Vectorial Decomposition consists basically in finding the circle of the rotating pattern in the complex plane. (iii) In the peculiar case where pairs of sensors are diametrically opposed sensors (for example two pairs at α = {0, π} and {π/2, 3π/2}), Doerfler and Ruchonnet implicitly assume in [1] that the rotating part verifies p rot (0, t) + p rot (π/2, t) + p rot (π, t) + p rot (3π/2, t) = 0.
Spatial Harmonics Decomposition is an alternative to these methods. The leading idea is to consider a closed curve that links each sensor position.

Mathematical considerations
The main purpose of this section is to introduce the mathematical framework necessary to compute and give the correct interpretation of the decomposition into spatial harmonics of a phenomenon. The first paragraph deals with abstract considerations. This is to establish such a decomposition in the most general case: a continuous problem in both space and time dimensions. Then, these equations are discretized to match with the experimental real life. Then, a quick overview of the consequences of spatial sampling leads to some basic recommendations on instrumentation and some keys to understand aliasing.

Mathematical framework
The key point to introduce the decomposition is to consider that the scalar field we want to study is defined (or known) along a closed curve. Some mathematical reminders are also necessary.

Prerequisites and reminders:
The space L 2 τ of the τ -periodic and square-integrable functions on [0 τ ] is a Hilbert space. It is fitted with an inner product given for any two elements f and g by Eq. 2: The set of functions {e n (x) = e 2iπnx τ , n ∈ Z} is an orthonormal basis of L 2 τ . Each element f of L 2 τ can be written as linear combination of e n (x), called complex trigonometric polynomial: 2.1.2. Defining a closed curve: Let's introduce the closed curve along which the scalar field is known. C is a C 0 closed curve in an N-dimensional Euclidian space E . C is defined by the given of C 0 periodic function r from R to E . Let's call τ ∈ R the period of r. C is then defined by parametrically, with x being the parameter. Note that r is absolutely not unique. Fig. 1 represents such a curve and an associated parametrisation. To simplify the expression of e n (x), it could be useful to subsitute x by a polar equivalent α = 2πx/τ .

2.1.3.
Defining a scalar field along C : On the other hand, let U be a differentiable real scalar field defined on E . It means that there is a function U , C 1 from E to R that associates a real number u to each point p from E : u = U (p). Taking an interest of the restriction U |C of U along the curve C means considering the function U • r. Fig. 1 gives a representation of what could be such a distribution, considering the altitude from C to the red curve as an image of U |C . Considering the τ -periodicity of r, U • r is also τ -periodic.
Writing u as a trigonometric polynomial allows a geometrical interpretation of the spatial repartition of the field U |C . As U is a real scalar field, a k and a −k are complex conjugates. Then, the sum a k · e k (x) + a −k · e −k (x) belongs to R and can be written as follow: ,where α = 2πx/τ (6) This sum describes an harmonic signal whose order is k. The complex coefficient a k contains both phase and magnitude information. This harmonic is also called a k nodal diameter mode.

Interpretation:
It is time to focus on geometrical interpretation of this decomposition. U |C is now decomposed as a sum of pure harmonic modes whose magnitude and phase informations are contained in the corresponding a k coefficient. Fig. 2 illustrate this interpretation describing the two first modes contained in the field represented on Fig. 1.
(i) a 0 : The coefficient a 0 belongs to R. It matches the average value of U • r(x) and then corresponds to constant amplitude field along C . (ii) a 1 : The coefficient a 1 and its conjugate complex a −1 describes the shape (magnitude and position) of the first spatial harmonic mode.  Figure 2. Interpretation of the two first orders 2.1.6. Time projection: Now considering U as a time-dependent scalar field, the spatial harmonic decomposition leads to the determination of a set of complex signals a k (t) also timedepedendant. The field U |C can still be written as a trigonometric polynomial in which space and time variable are separated.
Each of the time-dependent complex coefficients a k (t) contains the phase and magnitude information regarding the behavior of the harmonic spatial mode that presents k nodal diameters. Assuming that U is also continuous with respect to the time variable, a k (t) signals are now compatible with a frequency analysis thanks to continuous Fourier transform. It yields: This next step into the decomposition process leads to a new reduction of U |C into simple elements. Indeed, the frequency analysis of a k (t) coefficients reduces the behavior of each k nodal diameter mode into pure rotating or pulsating components. Let's rewrite the equation of the synchronous mode a 0 (t) using the fact thatâ 0,ω =â 0,−ω (see eq. 9).
Thus, it clearly appears that each elementâ 0,ω of the spectrum of a 0 (t) represents the magnitude of a synchronous mode pulsating with the frequency |ν = ω/2π| as described in [1] . In addition, the peak amplitude of this synchronous pulsating mode is given by 2|â 0,ω |. A spectral representation of a 0 (t) with a strong emergence at ω 1 as shown in Fig. 3 indicates that the behaviour of the synchronous puslation is dominated by a phenomenon that pulses at ω 1 .
For example, a spectral representation of a 1 (t) with a strong emergence at ω 1 as shown in Fig. 4 indicates that the behaviour of the first harmonic mode is composed by a scheme rotating at ω 1 rad/s. Through this writing, each elementâ k,ω of the spectrum of a k (t) contains the magnitude of a k nodal diameters mode moving along C with the frequency |ν = ω/2kπ| and the phase shift ϕ k,ω . The peak amplitude of this rotating mode is given by 2|â k,ω |. As a k (t) ∈ C, its complex spectrum does not present any symmetric property. Therefore, considering both positive and negative pulsations is necessary. Negative pulsations deals with modes that move in the positive direction of C (x growing). Positive pulsations deals with modes that move in the negative direction.
2.1.7. Energetic considerations: From an energy point of view there is an equivalence between what could be calculated integrated U |C during a time T, its decomposition a k (t) and the spectrum of this decompositionâ k,ω :

Application to experimental processes
The application of the spatial harmonic decomposition on measured datas is relatively simple. It supposes the implementation of a set of sensors measuring the same quantity and the adaptation of the mathematical theory developed in section 2.1 to sampled signals. Let's consider a set of N sensors placed at different positions c n and call x n the unique antecedent of c n by r on [0 τ ]. We call a regularly spaced set of sensor an implantation that satisfies a constant step dx = x n − x n−1 between each sensor position. U |C is no longer continuously defined but known on a finite number N of positions.
2.2.1. Defining C : Most of the time, the experimenter chooses the most logical curve with respect to the topology of the problem provided that it passes through all of the sensors. Regardless of this, the shape of C does not have an impact on the calculation of coefficients, but modifies the way they are geometrically interpreted. On the other hand, the function r that parameterizes the curve C has a direct impact on a k (t) = U (r(x), t), e k (x) . This is due to the fact that r could deform the image in E of a pure harmonic mode. Generally speaking, r is chosen to conserve linearity between [o τ ] and a curvi-linear abscissae of C .
Given the number N of sensors, only a limited number of independent coefficients could be calculated. This is a classical result regarding sampled data and it will be further discussed in the section concerning aliasing. This been said, it is important to notice that the spatial decomposition is a reversible process. It means that, knowing N independent coefficients a k (t), k ∈ [0 N − 1] allows an exact reconstruction of the signal of each sensor using Eq. 4.
2.2.3. Calculatingâ k,ω :â k,ω coefficients are nothing but a Fourier transform of a k (t). Then the best way to compute them is to use a FFT algorithm as long as it is compatible with complex signals. Note that asâ k,−ω =â −k,ω , keeping negative pulsations in FFT result is not necessary.

Energetic considerations:
The transposition of eq.11 into a discrete form leads to the following formulation (14) that links RMS (Root Mean Square) estimation of sensors with RMS of |a k (t)| or |â k,ω | coefficients.
2.2.5. Spatial aliasing: Using a finite number N of time-sampled signals raises the problem of space and time resolution. Regarding impacts of time resolution, experimenters are well informed today and they care to use a suitable sample frequency and anti-aliasing filters. Regarding spatial sampling of U |C using N sensors, aliasing problem can only be removed by increasing the number of sensors. But this is often not realistic. Then, the alternative consists in accepting the problem and studying aliasing consequences in order to separate phenomena during the analysis. In the case of a regularly spaced set of sensors, spatial aliasing comes from the periodicity of e k (x n ) = e 2iπkxn τ when x n = nτ /N, n ∈ [0 N − 1]. Then it comes easily an equivalence between aliasing and congruence relation of orders modulo N : Under these conditions a k (t) = a k (t). The spatial harmonic decomposition cannot differentiate a k nodal diameter mode from a k' nodal diameters mode. It means that the signal a k (t) contains information of all modes whose order k' is congruent to k modulo N . The best way to differentiate them is to compute the frequency componentsâ k,ω of a k signal. Thus relying on his perspicacity, the experimenter must be able to discriminate most of the time if the strong emergence at ω is due to a k nodal diameter mode rotating at ω/k rad/s or a k' nodal diameter mode rotating at ω/k rad/s. 2.2.6. Nyquist spatial frequency: This phenomenon occurs only in the case of an even number of sensors. In both sensors configuration cases, regularly spaced or non-regularly spaced, the calculation of a k (t) when k ≡ N/2[N ] leads to an indetermination: it is not possible to compute correctly both phase and magnitude. This constitutes a real leak of information whereas aliasing is only a mix of information. Fig. 6 shows a situation in which a regularly spaced set of 4 sensors cannot be used to dertermine phase and magnitude of a second order harmonic scheme. Figure 6. Second order harmonic indetermination using 4 regularly spaced sensors 2.2.7. Summary and recommendations: In section 2 we established the possibility to project measurements to a set of pure harmonic modes. It leads to a decomposition of these measurements into rotating and pulsating elements that allows the identifications of main phenomenon driving the measured signals. On top of that, further mathematical developments help to provide some recommendations for the installation of sensors: • Preferably use an odd number of sensor in order to avoid Nyquist spatial frequency.
• Preferably use a regularly spaced set of sensors to simplify aliasing estimation.

Draft Tube Pressure Pulsations
The case considered here is a medium head Francis turbine with 15 blades. The SHD has been applied to experimental records of dynamic pressure from 4 sensors located at the standard cross section beneath the runner outlet. Data has been processed for different operating points of a load variation to the head of best efficiency. Furthermore, the turbine model has been operated at a low Thoma Number and without anti-resonance devices in order to emphasize resonances and instabilities. The discharge will be expressed with respect to the discharge of best efficiency Q opt while frequencies are reported relative to the runner rotation frequency f 0 .
3.1.1. Observed phenomenon : On Fig. 7, RMS values for each sensor are reported alongside with the waterfall of one sensor. From this waterfall, at least four phenomena can be identified: (i) The Part Load Vortex Rope (PLVR) with three harmonics. The frequency of the first harmonic is f P LV R . It is independent of the discharge and equal to approximately 1 4 f 0 . (ii) The Upper Part Load Resonance (UPLR) with one or two (pseudo-)harmonics. Its frequencies f U P LR,i strongly depends on the discharge. Additional terms with frequencies of the form f U P LR,i ± f P LV R related to the vortex rope frequency are also present. f U P LR,i + f P LV R term has a bigger amplitude than f U P LR,i − f P LV R . (iii) A High load Instability (HLI) with one harmonic occuring on a narrow range of discharges. (iv) A Very High Load Instability (VHLI) developed at a main frequency f V HLI . PLVR, UPLR, HLI and VHLI are already known and IEC Code [3] gives examples of such frequency patterns. Besides the PLVR that is a common feature for all Francis turbines, other phenomena strongly depend on the runner design and the operating condition.
3.1.2. Application of SHD: SHD coefficients were computed along this load variation. For the reasons explained in section 2.2, four sensors give only access to P 0 , P −1 and P 1 , P −2 and P 2 . From an energetic point of view, this decomposition leads to the equivalence: Then Fig. 8 shows the evolution of SHD components along this load variation. A comparison with sensor's RMS 2 average is also made. From a mathematical point of view this figure illustrates two results explained in section 2: Figure 8. RMS values of the signal before any frequency decomposition, of P 0 , P −1 , P 1 and P 2 = P −2 • The conservation of energy, as described by Eq. 16, imposes that the average energy measured by the sensors is equal to the sum of the energy of each SHD component. This equivalence is met from 0.6 Q/Qopt until the maximum discharge. • Nyquist spatial frequency leads to an undervaluation of the amplitude of the second order harmonics scheme coefficients P 2 and P −2 using four regularly spaced sensors. That explains the difference between the average energy measured by sensors and the energy contained in SHD coefficients for discharges lower than 0.6 Q/Qopt. Regarding pulsating phenomenon described in section 3.1.1, SHD tells us additional information. Moreover, a frequency analysis on P k allows a complete description of each of them in terms of pure rotating schemes: . Colored waterfall of the magnitude of P 0 , P 1 , P −1 and P 2 as a function of the frequency and of the discharge. Colors refer to the same scale. Waterfalls of P 0 , P −1 , P 1 and P −2 = P 2 are pictured with the same scale on Fig. 9: (i) PLVR is mainly caught by antisymmetric scheme P −1 at the frequency f P LV R . It means that this phenomenon is driven by a one nodal diameter mode rotating the same direction as the runner at the angular velocity 2πf P LV R k = ω P LV R . The residual emergence at the same pulsation on P 0 component represents the induced synchronous pulsation as described in [4]. A much smaller second harmonic is seen on P −2 = P 2 , which gives it the same angular velocity as the first harmonic.
(ii) UPLR first and second harmonics are captured by P 0 . Terms with frequency f U P LR,i + f P LV R and f U P LR,i − f P LV R (the latter is hardly visible) are registered on respectively P −1 and P 1 that links them to the rotating pattern of the vortex rope. (iii) HLI is exclusively captured by P 0 . On this model, HLI is actually characterized by an axial pulsating rope (iv) VHLI is mainly captured by P 0 for its main frequency f V HLI . It has a smaller part on P 1 with a band of frequency below 2f 0 . The fact that VHLI is composed of both a synchronous and a rotating pattern is an interesting difference with respect to HLI.
The low P 0 contribution in the PLVR should be noted. Indeed, this feature should be correlated with the fact that, usually, PLVR does not generate high torque or power fluctuations contrary to UPLR or HLI/VHLI which have stronger P 0 components.

Gap Pressure Pulsations
In continuation of what has been done with draft tube wall pressure sensors, SHD can be applied for detecting occurrences of coherent structures in the vaneless area. Applied to pump-turbines, this study particularly relates these structures with the detection of high radial thrust in the area of turbine brake operation called "S shape". The tested model is a 9 bladed high head pump-turbine instrumented with a non-regularly spaced set of 4 sensors located at the same radius between the stay-vane and guide-vane. The variation follows a constant opening curve from high to low discharge represented in red on Figure 10. Figure 11 shows a strong emergence at 0.6f 0 on radial thrust Fx spectrum.  speed. Calculating SHD highlights the same emergence on P 1 (t) only (see fig. 12 and 13) and describes a pure rotating pressure pattern at the same pulsation. Thus it links the spatial coherence of the gap pressure field with a mechanical effect. Note that a non coherent pressure field does not induce radial thrust that could not be seen in analysing sensors separately. Figure 12. magnitude ofP −1,ω as a function of the frequency and the discharge Figure 13. magnitude ofP 1,ω as a function of the frequency and the discharge

Conclusion and Outlook
The Spatial Harmonic Decomposition is a mathemathical method that allows spatial interpretation of a measured scalar quantity. In the field of hydraulic machinery, it is a more general alternative for existing time-space decompositions applied to a network of pressure sensors since it gives more information about rotating structures. Furthermore, introducing spatial aliasing and Nyquist frequency gives a new understanding about the number of sensors and their relative postition required to properly capture rotating patterns. It has already proved to be very efficient to extend both draft tube and gap pressure pulsations into meaningful patterns. As it ensures energetic conservation, it could be a very suitable way to go beyond a simple pressure fluctuation analysis using the peak to peak estimation of each sensor. Indeed, manufacturers know today how to design machines with low peak to peak pressure fluctuations but they should more focus on the phenomena that directly affect the functioning in terms of power stability or mechanical constraints is more relevant. This leads to a better understanding of operating limits and helps to work on their extension.