MULTIDIMENSIONAL DECOMPOSITION OF TIME-VARYING VIBRATION RESPONSE SIGNALS IN ROTATING MACHINERY

Mechanical Systems and Signal Processing (1997) 11(4), 577–601

MULTIDIMENSIONAL DECOMPOSITION OF TIME-VARYING VIBRATION RESPONSE SIGNALS IN ROTATING MACHINERY I. B Faculty of Mechanical Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel

 D. J. E Department of Mechanical Engineering, Head Dynamics Section, Imperial College, London SW7 2BX, U.K. (Received October 1996, accepted February 1997) A new method for decomposition of speed-dependent vibration signals measured on rotating machinery is presented. The method uses a number of measured signals from an array of sensors to decompose the vibration response into: (1) the direction of the progression of motion; (2) its frequency content; (3) the various wavelengths; and (4) speed of rotation (or time). For the simple case of shaft vibration, two sensors spaced at 90° are required to separate the vibration pattern into forward- and backward-travelling components. For the case of a rotating disc, more sensors are required to decompose the response further into different spatial wavelengths and their directions of progression. This decomposition of the response allows monitoring and identification of potentially dangerous vibration patterns which exhibit large backward- or forward-progressing components. The new method is shown to provide better resolution in the time–frequency (speed–frequency) representation where usually overlapping patterns are separated. Several analytical and experimental results that demonstrate the usefulness of the proposed method are shown. 7 1997 Academic Press Limited

1. INTRODUCTION

This paper describes a new method for the decomposition and visualisation of vibration signals measured on rotating structures. This decomposition allows one to analyse the vibration patterns and to identify their generating mechanism with greater insight than before. The separation provided by this new method has four elements: (i) the spatial domain, (ii) the temporal domain, (iii) the direction of progression, and (iv) the frequency domain are all treated. The effects of dynamic forces acting upon a rotating structure are numerous and thus one seeks a tool that provides a qualitative as well as a quantitative appreciation of various parameters affecting the response. The procedure described here provides both qualitative information in terms of vibration patterns and quantitative information in the lines showing the values of the natural frequencies and by the computed response levels. The time or speed dependency of this information can be evaluated directly by the display produced. The main parameters that affect the nature of the response which can be examined using the suggested method are: (a) speed of rotation; (b) external loading; and (c) asymmetries in both rotating and non-rotating parts of the structures. Although the 0888–3270/97/040577 + 25 $25.00/0/pg970092

7 1997 Academic Press Limited

578

.   . . 

interpretation of the decomposed results still requires much effort and physical understanding, the separation provided by the new method greatly reduces this effort. When an axi-symmetric rotating structure (such as a disc) is being studied, generally a mode shape has to be assumed [1] and so the lack of spatial knowledge obtained from just a few measured locations can be overcome by using the proposed method. Separation into the different directions of rotation and the different spatial wavelengths can be used to check those assumptions and to give a better graphical estimate of the actual mode shape. Customarily, the vibration behaviour during run-up or run-down operations of rotating machinery is recorded and plotted on a time–frequency distribution diagram [2] or a Campbell diagram [3] or, more recently, using wavelet maps [4]. Both time–frequency and Campbell diagrams provide important information about the evolution of the frequency spectrum of the response with time or with speed of rotation. This information assists considerably in detecting changes in structural resonance frequencies, useful for detecting structural changes (or possibly damage). The method introduced here is an enhancement to any of the existing time–frequency distributions, as long as the phase information is retained. Combining several such distributions derived from several sensors (here the usual Campbell diagram is used which is particularly applicable to rotating machinery) one can better understand the measured response of rotating machinery parts having complex vibration patterns (mode shapes). The combination of rotation effects and a complex deformation pattern yields many spectral lines in the measured response which, in turn also result in a very dense Campbell diagram. Such a diagram is observed [3] whenever there are discs, bladed discs or rotational asymmetries on the foundation (bearings) or the rotating parts commonly found in high-speed machinery like jet engines and steam-turbines [3, 5]. The complex spectra measured on rotating machines are the result of propagating waves having different wavelengths (caused by different numbers of nodal diameters for discs) travelling in a co-rotating or counter-rotating direction. Even in the case of a simple rotating shaft mounted on bearings, forward and backward orbits are observed [6]. The existence of a crack [6], or other forms of asymmetry in the rotating shaft, manifests itself in the form of periodically varying stiffness or inertia which also results in additional harmonics in the response [3]. Following the Introduction, Section 2 provides a theoretical overview which briefly covers time–frequency distributions and Campbell diagrams. Also given in this section are some examples for rotating systems whose response patterns can be described in terms of waves travelling in different directions. These examples provide the motivation for the development of the main algorithms. In Section 3, the decomposition procedure for a rotating shaft and a rotating disc is developed and discussed. A method to transform the results to the disc (rotating) co-ordinates is shown as well. Some examples using analytical and experimental data are given in Section 4. The examples are discussed and certain features of the presented information are outlined. The final part, Section 5, offers the conclusions and recommendations for further research. 2. THEORETICAL BACKGROUND

This section outlines some assumptions, results and tools which are embedded in the proposed method. A comprehensive comparison and review of all available time–frequency computation methods is not attempted, but a brief introduction to Campbell diagrams and to the computation of time-varying spectra is given. The method presented in this paper can be applied to any type of time–frequency or speed–frequency distribution, wavelet maps or short-time Fourier transforms, as long as the phase

   -  

579

information is retained. The phase information is also required for the spatial decomposition performed by the proposed method. In terms of wavelet analysis this means, for example, that the harmonic wavelets [4] need to be used rather than simply taking the customary amplitude map. The operating conditions of rotating machinery vary with time, especially during start-up and rundown. During start-up, the speed of rotation is increased until the machine reaches its operating speed. While in the transient regime, all excitation forces (e.g. due to unbalance) change both in amplitude and frequency. Thus, in order to capture the time dependency of the spectrum, one needs a tool to compute and present such behaviour. This analysis tool is briefly presented below. 2.1. -     It is a well-established fact that the frequency content of the response of rotating machinery is a function of (amongst other factors) the speed of rotation. The monitored vibration signal changes its characteristics in a systematic way when the speed of rotation is being changed. For example, the frequencies at which the various maxima occur vary when the speed of rotation is changed and the amplitude of vibration is high when the frequency of rotation (or one of its integral multiples) coincides with one of the structural natural frequencies. The estimation of a continuously varying spectrum from sampled response data has recently received much attention [2], and from these distributions diagrams describing the evolution of the spectrum with time can be produced. Usually, the distribution of the response levels for all frequencies in the measured range (spectrum) is plotted against time (spectrogram, time–frequency diagram) or against speed of rotation (Campbell diagram, speed–frequency diagram). In the present work, the so-called short-time Fourier transform (STFT) is used to construct the time–frequency diagram. The Campbell diagram is computed from the time–frequency plot by ensuring a linear relation, t = aV + b (where a and b are constants) or by repositioning the spectra at equal speed increments. An estimate for the time-frequency distribution (TFD) is commonly computed by using the STFT by application of the following integral: S(v, t) =

g

a

s(t)W(t − t) e−ivt dt

(1)

−a

where s(t) is the measured signal at time t, S(v, t) is the time–frequency distribution, W(t − t) is the localisation window centred around t. The role of the window is described in Figure 1.

W(t – τ) window

Signal

Time t

Figure 1. Localisation around a short-time/speed record for the short-time Fourier transform using a sliding window.

580

.   . . 

The underlying assumption behind the STFT method is that the measured spectrum does not change too rapidly [2], e.g. the time record within the selected window is quasistationary. An important definition in time–frequency analysis is the one of the instantaneous frequency of a time function. The instantaneous frequency, fi (t), is related to the instantaneous phase, 8(t) by [7]: fi (t) =

1 d8(t) 2p dt

(2)

Equation (2) is used throughout this paper in the various examples to define the instantaneous phase of the signals. 2.2.          In this section, the two basic systems treated in this paper are briefly described. Some of the structural properties are attributed to the expected response patterns thus assisting in the interpretation of the examples and the analysis shown later. 2.2.1. Rotors—rotating shafts with rigid discs A very important feature of the response of rotating machinery is the speed dependency of the response spectrum. Changes of the response with speed of rotation are commonly observed for rotating shafts with discs mounted on them. The unbalance response resulting from the combined effect of very small imperfections and structural asymmetries manifests itself in both the amplitude and frequency of the response. Frequencies associated with the natural frequencies of the structure are typically measured on rotors. These frequencies often vary with speed of rotation due to gyroscopic and centrifugal effects altering the natural frequencies of the structure. The gyroscopic forces have a stiffening effect on forward (co-rotating) whirl of shafts and a softening effect on backward (counter-rotating) whirling rotors. The whirl motion of shafts can be described conveniently in terms of travelling waves. 2.2.1.1. Travelling waves in rotating shafts. In this paragraph two physical systems are described for which the response can be described in terms of forward and backward whirl. Whirling motion of shafts causes the centre of gravity of a shaft cross-section to move in an elliptical orbit (or a combination of elliptical orbits at different frequencies) around the centre of rotation. Whirl can also be described in terms of travelling waves progressing in a co-rotating direction (forward) or counter-rotating (backward). Backward whirl is usually encountered when the bearings and the support structure are anisotropic. This was demonstrated analytically by several researchers [5, 6] for a simple Jeffcott rotor model (Fig. 2). This analysis is given in Appendix A for reference and to assist in the derivation of the main results. It can be shown that backward whirl is excited for such a system by ordinary, forward-rotating unbalance (see Appendix A) [6]. Hence, separation of the whirling response into forward and backward components provides valuable information about the asymmetry of the bearing stiffness in the x- and y-directions. Under steady-state conditions when the speed of rotation, V, is constant, and when the only excitation is due to unbalance, the rotor will whirl at a single frequency and the response measured by s1 (t) and s2 (t) will have the shape of an ellipse. As shown in Fig. 3, this ellipse can be decomposed into two opposite directional whirl circles. When there are several sources of excitation or when transient conditions are in effect, each response signal contains many frequencies. For every frequency, there exist forward and backward complex amplitudes.

   -  

581

Figure 2. A dual sensor system measuring rigid vibrations.

A rigid disc mounted on a rotating shaft experiences similar motion to that of a shaft, this motion being characterised by one nodal diameter mode (see Figure 4). In spatial terms, the two wave functions, C1 (r0 ) eiu and C−1 (r0 ) e−iu (considering here a polar r0 , u system) suffice to describe the axial or lateral motion of every point in the system. As we consider the dynamic response we treat a time-varying combination of these two wave functions at every moment of time. The response of the shaft at a single frequency at any point along the shaft, can be described by means of a displayed elliptical orbit. This elliptical orbit can be described by means of two phasors (and a constant displacement), one phasor is rotating in a co-rotating (forward) direction and the other in a counter-rotating (backward) direction

Figure 3. Elliptical orbit showing the major and minor axes as a function of the forward (Af ) and backward (Ab ) whirl amplitudes.

582

.   . . 

Figure 4. Elliptical orbit decomposed into forward and backward circular whirling orbits (phasors).

(see Fig. 5), i.e. for a constant whirl frequency, the angle equals u = vt and we have: s˜ = C1 eivt + C−1 eivt

(3)

Where C1 (C−1 ) is the constant complex forward (backward) whirl amplitude. 2.2.2. Rotating shafts with flexible discs Flexible discs exhibit a spatial vibratory response having more than a single nodal xdiametral component. At a constant radial position r = r0 , one can express the spatial displacement at any angular location by superposition of wave functions [8]: Cn (r0 ) einu

and C−n (r0 ) e−inu

(3a)

where r0 and u measure the radial and angular locations on the disc respectively and n represents the number of nodal diameters on the disc. In the disc co-ordinate system, a pictorial representation of the first nine possible wave functions (i.e., n = 1, 2, . . . , 9) is shown in Fig. 6, for a specific radial location r = r0 . Since we are dealing with a rotating structure, the combinations of all the mentioned patterns (wave functions) can be time-varying, so that a stationary sensor located at (r0 , u)

Figure 5. Different mode-shapes (having different wavelengths) of a cyclically symmetric structure.

   -  

583

Figure 6. One nodal diameter motion of a shaft with a rigid disc.

in a polar co-ordinate system will measure a time-varying superposition of the spatial forms: a

su (t) = s pn (t) einu + p−n (t) e−inu

(4)

n=0

Here pn (t) (p−n (t)) are the time-varying forward (backward) n-nodal diameter related amplitudes. su (t) is the signal measured by a sensor located at an angle u. A method allowing for the separation of the forward and backward related amplitudes for every number of nodal diameter is described in Section 3. 3. SPATIAL DECOMPOSITION OF TIME OR SPEED VARYING SPECTRA

In this section, the proposed decomposition method is developed for two basic systems: a flexible shaft and a flexible disc mounted on a shaft. 3.1. /        The lateral motion of the shaft and the axial motion of the disc can be described by means of equation (4), where only a range of nodal diameter (and hence a range of wavelengths) is considered to contribute. This range is n = nmin · · · nmax . nmax

su (t) = s pn (t) einu + p−n (t) e−inu

(5)

n = n min

When only a shaft is considered, the effective range (as shown in Appendix B) is nmax = nmin = 1 where for a rigid disc, the axial motion described by n = 0 needs to be considered as well. We can write explicitly the expression for a signal measured on a shaft, as shown in Fig. 2 or on a disc in the z-direction, as shown in Fig. 4. su (t) = p1 (t) eiu + p−1 (t) e−iu + C

(5a)

A rigid disc attached to a flexible shaft is depicted in Fig. 4. In this figure, the intersection of the disc with a plane parallel to the x–y plane is seen to create a single nodal line. In what follows, two methods allowing us to extract the forward and backward whirl with amplitudes are outlined. 3.1.1. Method I: using a combined signal to separate forward and backward whirl for shafts Combining the response measured from s1 (t) and s2 (t) (see Fig. 2 and 7), to create a complex signal,

584

.   . .  s˜ (t) = s1 (t) + is2 (t)

(6)

we obtain, for a single frequency v: s˜ (t) = Af eivt + Ab e−ivt + C

(7)

where Af , Ab are the forward and backward vibration amplitudes respectively and C equals (x0 + iy0 ) and is constant. It is shown in Appendix B that the complex signal created in equation (7) can be substituted into equation (1) to obtain a double-sided spectrogram for which positive frequencies represent the forward whirl and negative frequencies represent backward whirl. Here, we prefer a more systematic approach that allows us to solve for the wave amplitudes by using the spatial discrete Fourier transform (DFT). This approach is later generalised for the multi-wavelengths case. 3.1.2. Method II: a spatial approach for separating forward and backward whirls for shafts We can write three equations for the three unknowns (f1 (t), f−1 (t), C) (see Appendix D). Here we prefer to use the Fourier (spatial) transform in order to obtain the desired coefficients and for that purpose we use a hypothetical arrangement consisting of four sensors, as shown in Fig. 7. Separating the static part (x0 , y0 ) for the dynamic part (f1 (t), f2 (t)) it can be observed that the following relationships hold for the various sensors: s1 (t) = x0 + f1 (t)

s2 (t) = y0 (t)y0 + f2 (t)

s3 (t) = −x0 − f1 (t)

s4 (t) = −y0 − f2 (t)

(8)

From equation (8) it can deduced that only two sensors are required, but for convenience we use all four sensors in the DFT formula, as follows (see also appendix D): 3

pn (t) = s sk (t) e−i(pkn/4)

(9)

k=0

where k, n = 1, 2, 3 and pn (t) is the term associated with the nth nodal diameter motion. Due to the cyclic repetition property of the DFT and since p/2 is the Nyquist wavelength, p3 (t) = p−1 (t) is the 1-nodal diameter backward-propagating wave term. This fact will be used later in Section 3.2 and is also discussed in Appendix D.

Figure 7. (a) Three-dimensional view of the four sensors measuring the shaft whirl orbit. (b) Planar view of the four sensors measuring the shaft whirl orbit indicating the undeflected position.

   -  

585

Figure 8. Rotating shaft with a flexible disc and an array of sensors to measure disc vibrations.

Applying the linear integral transformation in equation (1), we obtain: 3

Pn (v, t) = s Sk (v, t) e−i(pkn/4).

(10)

k=0

Here Pn (v, t) and Sk (v, t) are the time–frequency distributions of pn (t) and sk (t) respectively. Pn (v, t) is the time–frequency distribution associated with the n nodal diameter related motion. It can be shown (Appendix D) that the forward and backward time–frequency distributions can be expressed in terms of the individual time–frequency distributions (each sensor processed independently): forward whirl P+1 (v, t) = S1 (v, t) − iS2 (v, t)

(11a)

backward whirl P−1 (v, t) = S1 (v, t) + iS2 (v, t)

(11b)

a where Sk (v, t) = f−a sk (t)W(t − t) e−ivt dt is the time–frequency distribution measured from the kth sensor.

3.2. /        A measurement system is depicted in Fig. 8 capable of separating flexible discs’ vibration into different wave patterns (in fact, the method is applicable for any type of rotating structure and not only discs). A stationary sensor, sk , will measure the response which the disc experiences modulated by the rotating mode shapes, so that by using an array of sensors we are able to decompose the different wavelengths present. In this paper, most derivations assume that the sensors measuring the vibration signals on the surface are equally spaced (circumferentially, Du), but this is not a requirement. For non-equally spaced sensors, a general set of linear equations need to be solved [9] instead of using a standard DFT.

586

.   . . 

In Fig. 8, a vector of measurements at time t, [s1 (t) s2 (t) · · · sm (t)], is constructed from all m sensors. The decomposition of the different wavelengths necessitates the extraction of the coefficients in equation (5). Since the sensors are assumed to be spaced equally at an angle Du, the shortest wavelength that can be extracted, according to the Nyquist sampling criterion [10, 11], (expressed in terms of the equivalent number of nodal diameters) is: nnyq =

p Du

(12)

For example, in the system depicted in Fig. 2, the two sensors spaced at 90° can detect up to nnyq = 180°/90° = 2. Given the m discretely distributed sensors, a system of linear equations can be set up: nmax

su (t) = s pn (t) einu + p−n (t) e−inu

(13)

n=0

for every m+1 F m odd G 2 u = kDu, k = 1 . . . m and nmax = g m m even G f2 The solution for equation (13) consists of m coefficients pn (t) which are generally complex. When the number of wavelengths present in the response nmax is larger than nnyq (practically nmax is restricted by the number of sensors), spatial aliasing will occur and the estimate of the coefficients will be erroneous. This apparent difficulty can be resolved by using band-limiting constraints (wavelength-limited, when applicable) or by using a sufficient number of sensors. In any case the terms from the spatial aliasing can be identified since they appear at the ‘wrong’ frequencies so that physical understanding of the expected behaviour can be used to detect them. Wavelength-limited interpolation can be demonstrated for the shaft depicted in Fig. 2. As stated before, the two sensors spaced at 90° can theoretically detect up to nnyq = 180°/90° = 2, but in this case we have only n = +1 and −1. Thus, we can force C2 = C−2 = 0 which leads to the result described in equation (8) and in full in Appendix D. Therefore, we need use only two sensors (instead of four) while no information is lost. A similar procedure can be taken when the contribution of short wavelengths is very small. Then short wavelengths are neglected without greatly affecting the results. 3.2.1. Combined decomposition of wavelengths and direction of progression in a multidimensional time–frequency or Campbell diagram plots In the previous sections, the operations and transformations leading to the time- or speed-dependent spectra, as well as spatial decomposition of the response, were described. Both operations use a one-dimensional Fourier transform (or alternatively a curve-fitting procedure) to obtain the spatial decomposition. For the kth sensor the time (speed) frequency distribution is computed according to equation (1), a Sk (v, t) = f−a sk (t)W(t − t) eivt dt. The spatial pattern of the displacement for each time and frequency instant, is a combination of forward and backward travelling waves. Thus, the time–frequency distribution measured by the kth sensor, (located at the angular

   -  

587

location, u = kDu), can be expressed as: nmax

Sk (v, t) =

s P+n (v, t) e+inkDu + P−n (v, t) e−inkDu

(14)

n = −n max

where P+n (v, t)(P−n (v, t)) is the time–frequency distribution of the forward (backward) progressing n-nodal diameter wave. The solutions for equation (14) are P+n (v, t) and P−n (v, t) which are functions of three parameters: (1) n, the number of nodal lines—determines the spatial wavelength; (2) t, the time instant—determines the time or the speed computed from V = V(t); (3) v, the frequency—determines the frequency for the specific wavelength and time or speed. This function P+n (v, t) or P−n (v, t) can be plotted separately for each wavelength (n) to reveal a multitude of information, as will be shown in Section 4. 3.2.2. Tranformations of the multidimensional spectrograms to rotating co-ordinates In many applications, it is desired to determine the travelling waves with respect to the rotating system. In this case, the above results can be transformed using the instantaneous phase angle of the disc:

g

c(t) = V(t) dt

(15)

The phase angle in inertial co-ordinates u can be related to the phase angle in rotating co-ordinates, 8, by: u = 8 − c(t)

(16)

where 8 is the angle in the disc local co-ordinates. Practically, the phase can be measured directly from a device called a resolver or by using an absolute encoder. Since different wavelengths travel at different speeds, it is most convenient to transform each wavelength, namely P+n (v, t), separately. A simple way of understanding the shift in the frequency scale each term undergoes is to apply the Doppler effect (a more rigorous derivation can be found in [8, 9]. Thus, for the nth nodal diameter related travelling wave on the disc, the shift in frequency would be: forward wave frequency shift v − nc (t) = v − nV(t)

(17a)

backward wave frequency shift v + nc (t) = v + nV(t)

(17b)

The time–frequency distribution experienced by a sensor attached to the rotating body is simply: R P2n (v, t) = P2n (v 3 nV(t), t)

(18)

R where P2n (v, t), P2n (v 3 nV(t), t) are in the rotating and stationary co-ordinates respectively. The shift in frequency is more easily implemented in the transformed domain, using the shift property of the Fourier transform. The result must be transformed back to the

.   . . 

588

frequency domain as shown below: R P2n (v 3 nV(t), t) =

1 2p

g$g

%

P2n (v, t) e2ivt dv e+infV(h) dh e−ivt dt

(19)

or, using equation (15), R (v 3 nV(t), t) = P2n

g$g

1 2p

%

P 2 n (v, t) e−ivt dv e2inc(t) e−ivt dt

(20)

Equation (20) allows us to transform the time–frequency distribution measured in stationary co-ordinates to rotating, body-fixed ones. The only extra information needed is the instantaneous phase angle c(t) of the rotating disc. The dynamic behaviour in the body-fixed co-ordinates is important to estimate the stresses in the rotating part. It is important to note that this transformation is not possible unless the separation into a different number of nodal diameters has been performed. This is due to the fact that different wavelengths travel at different speeds and have different phase-shifts due to the rotation of the disc. 4. EXAMPLES

In this section, some simulated and experimental examples are shown. Various aspects of the proposed method are demonstrated. 4.1.   4.1.1. Separation of forward and backward components for a rotating shaft: Assume that two signals s1 (t) and s2 (t) are measured from the system in Fig. 2, and that the two measurements show the following time behaviour: s1 (t) = Af cos (8f (t)) + Ab cos (8b (t)) s2 (t) = Af sin (8f (t)) − Ab sin (8b (t))

(21)

Substituting equation (21) in (6), one obtains the forward–backward representation: s˜ (t) = Af e+i8f (t) + Ab e−i8b (t)

(22)

The phase arguments, 8f (t) and 8b (t) are defined via their time derivatives [instantaneous frequencies computed according to equation (2)], as follows: 2pff (t) = vf (t) =

2pfb (t) = vb (t) =

0

1

d8f (t) t = 2p 2 + 110 dt tmax

0

0 11

d8b (t) t = 2p 190 − 174 dt tmax

2

(23)

Processing the signals provides the results shown in Fig. 9. The example shown here consists of two waves, one travelling in a forward direction while the other is in the reverse direction. While Fig. 9(a) provides the information about the frequency content at each time instant, Fig. 9(b) provides the direction in addition and clearly separates the forward (positive frequency) and backward (negative frequency)

   -  

589

components. It is clear that the proposed method provides much better separation and does not display the smearing effects due to the overlap of the two frequencies of (a). The separation into forward and backward components, apart from providing more information, also yields better resolution. 4.1.2. Example: rotating disc This example simulates an array of nine sensors, fixed in space, (as in Fig. 8) all positioned at the same radial location r = r0 , but having different angular locations. The sensors are measuring the displacement of the disc in the axial direction. The disc is assumed to rotate at a constant speed V. The qth sensor is located at an angle uq = q · 12°, q = 1 . . . 9; each sensor measures the following time function: 4

sq (t) = s Afn cos (n(uq − 8f (t))) + Abn cos (n(uq + 8f (t)))

(24)

n=1

where the forward amplitudes are: {Ab 1 , Ab 2 , Ab 3 , Ab 4 } = {0.0, 2.0, 0.0, 1.0} the backwards amplitudes are: {Af 1 , Af 2 , Af 3 , Af 4 } = {0.7, 0.0, 2.0, 0.0} the instantaneous phase is: 8(t) = p(13.33 + 26.6667t2 ) and n is the number of nodal diameters. The resulting STFT (time–frequency distributions) are presented in Figs 10 and 11. The threefold decomposition demonstrated in this example shows a dramatic increase in the understanding and separation of the vibration patterns in a rotating disc. The ordinary Campbell diagram plotted in Fig. 10 does not change from one sensor to another, while the directional and spatial decomposition in Fig. 11 clearly shows much more information and provides better resolution by separating the measurements both into forward–backward and different wavelengths components. 200 (b) 150

200 180 160 140 120 100 80 60 40 20 0

(a)

Frequency (Hz)

f (Hz)

100 50 0 –50

–100 –150 0

10

20

30 40 Time (s)

50

60

–200

10

20

30 40 Time (s)

50

60

Figure 9. (a) One-sided STFT using s1 (t) alone. (b) Double-sided STFT using s1 (t) and s2 (t).

.   . . 

590

1000 900 800

Frequency (Hz)

700 600 500 400 300 200 100 0

0

0.5

1

1.5

2 Time (s)

2.5

3

3.5

Figure 10. Campbell diagram using any one of the sensors.

4.2.  4.2.1. Rotating shaft mounted in bearings In this section, a brief background concerning the response of rotating shafts mounted in linear bearings is provided. Prior knowledge of the behaviour of such a system assists in the interpretation of an experimentally obtained Campbell diagram shown in the second part. A rotating and vibrating shaft exhibits several phenomena unique to a rotating structure. Among these are the forward and backward whirl and vibration response due to unbalance inevitably existing is such a system. The unbalance response manifests itself in the form of a speed-dependent harmonic excitation (see Appendix A). A rotating system could also have speed-dependent properties resulting from gyroscopic for centrifugal forces acting on the structure. In fact, the coexistence of multiple response phenomena in most rotating structure makes the distinction and separation quite difficult. The decomposition method presented in this paper helps in some cases to gain a better visualisation and separation of the various phenomena. One of many phenomena observed in the experimental example shown later is the intersection of the so-called first engine order lines (thick lines in Fig. 12) with some of the structural natural frequencies (thin lines). At these points, the unbalance excitation possesses a frequency identical to one of the structural natural frequencies which results in a resonance condition, e.g. high vibrations. When some degree of anisotropy in the foundations exists, the backward whirling modes exhibit similar resonance behaviour even though the unbalance excitation force is co-rotating. One expects to observe high amplitudes of response at these intersection points (circles in Fig. 12), because at those speeds, the excitation due to unbalance coincides with a natural frequency. Intersection in the positive frequency region means the excitation of a co-rotating mode. In the negative frequency region, the intersection describes a backward whirl resonance.

Frequency (Hz)

   -  

1000

1000

n=1

500

500

0

0

–500

–500

–1000

–1000

Frequency (Hz)

20

1000

40

1000

500

500

0

0

–500

–500

–1000

–1000 20

40 Time (s)

n=2

60

n=3

60

591

20

40

60

20

40 Time (s)

60

n=4

Figure 11. Directional and spatial decomposition of the Campbell diagrams for n = 1, 2, 3, 4 in the frequency range f = −1000 . . . 1000 (Hz).

4.2.1.1. Measured response and the computed Campbell diagram. In this section, experimentally obtained results measured on a rotating shaft are presented. The experimental rig consists of a DC motor coupled to a shaft supported in bearings. Also, a rigid disc is fixed to free end of the shaft (see Fig. 15). Two proximity probes (90°

Figure 12. Natural frequencies of clamped shaft vs speed and first engine orders.

.   . . 

592 200 180

Frequency (Hz)

160 140 120 100 80 60 40 20 0

500

1000

1500 2000 Speed (rpm)

2500

3000

Figure 13. Conventional diagram measured on a rotating shaft.

apart) are arranged as shown in Fig. 2 and the only excitation forces are due to the inherent unbalance in the system. The conventional and the two-sided directional Campbell diagrams are shown below in Figs 13 and 14, respectively. In this example, an experimental rig was used to demonstrate the value of the proposed method. Two plots are shown, the first one is an conventional Campbell diagram where

200

150

Frequency (Hz)

100

50

0

–50

–100

–150

–200 500

1000

1500 2000 Speed (rpm)

2500

3000

Figure 14. Forward–backward (two-sided) Campbell diagram measured using two sensors.

   -  

593

Figure 15. Shaft–disc system and the electromagnetic exciter.

all the information is overlaid on the positive frequency range. The second plot shows separated positive-frequency (forward) and negative-frequency (backward) vibration components and the superior resolution of the directional Campbell diagram is clearly demonstrated. It can be seen when the unbalance excitation is exciting the forward or the backward modes by inspection of the magnitudes (darker colours). The increase in the amplitude of the backward vibration is noticed when the excitation due to the first engine order exceeds the first (backward) resonance. 4.2.2. Rotating disc response measured with an array of sensors In the following section, the proposed method is applied to data measured on a rotating disc attached to a shaft. This is an initial implementation of the method to a more

500

Frequency (Hz)

400

300

200

100

0

5

10 Time (s)

15

20

Figure 16. Conventional spectrogram using sensor on the disc; log amplitude.

594

.   . . 

complicated system. The experimental and display routines are constantly being improved, and thus the presented example is an intermediate result and is likely to be refined in the future. In many real-world situations, a fixed-in-space excitation force is exerted on a rotating disc. This force can result from the interaction between the stator and the rotating disc, or in the laboratory it can be supplied by an electromagnetic exciter depicted in Fig. 15. A fixed-in-space, force, F0 eivt, was applied by an electromagnetic exciter. More information about the experimental set-up can be found in [11]. The experimental rig described here was used to test the proposed method for disc vibration. In this case, the disc modes for which the number of nodal diameters is greater than 1, cannot be excited by unbalance in the system, and so an external excitation has to be provided to generate a response in the required mode. The force experienced by the disc is modulated by the disc modes, thus the excitation of a specific modal pattern is achieved by the constant (DC) part of the excitation inevitably existing in the magnetic exciter. The time–frequency distribution for a single sensor and for the individual numbers of nodal diameters are depicted in Figs 16 and 17. The STFT can be displayed in a three-dimensional plot which reveals some more details. In Fig. 18, some distinctive features are outlined on the plots. Using an electromagnetic exciter, several modes were excited simultaneously. The pure axial motion and the umbrella mode are represented by the n = 0 case in Fig. 17. This type of motion seems to be very significant compared with n q 0, but other modes (e.g. n q 0) are also clearly visible from Figs 17 and 18. It is worth mentioning that for n = 0 and for n = nnyq = 4 no phase information can be obtained (the DC and the Nyquist points have arbitrary phase or are real) thus the forward/backward separation cannot be performed (see Figs 17 and 18 for n = 0 and n = 4). The spatial separation of the vibration patterns reveals that much of the spectral information in Fig. 16 is overlapped. On the other hand, the spatial decomposition shown in Figs 17 and 18 allows one to assess the contribution of each wavelength (number of nodal diameters). The ability to detect the direction of progression (forward–backward) for each wavelength provides valuable information assisting in the interpretation of dynamic response of the rotating disc. A large part of the response is due to the pre-bending of the disc (wobble). This pre-bending appears as engine orders (straight lines at f = kV, k = 1, 2, . . .) and this makes it difficult to observe some of the dynamics phenomena. This apparent difficulty is likely to be improved in future set-ups. The procedure which was described in Section 3.2.2 and equation (20), was applied to the measured disc response. The result is shown only for n = 2, due to space limitation, but the same transformation can be applied to every n. By inspection of Fig. 19 where the time–frequency distribution for n = 2 is shown in the disc co-ordinates, it can be observed that the forward and backward frequencies seem to increase (in absolute value). The physical interpretation for that effect is that the centrifugal stiffening of the disc is more dominant than the gyroscopic effect (that tends to reduce the backward frequency). All fixed patterns on the disc (wobble) which appear as engine orders in Fig. 17 are mapped into one line for which f = 0. The frequency shift due to the rotation (Gllilean tranformation), which is clearly noted in Fig. 17, has disappeared. Making use of the instantaneous phase angle, the measured response signals were resampled to obtain equal angular spacing. This type of processing results in new time spacing, yielding an equal number of response samples per one shaft revolution. For this type of processed signal the frequency axis is replaced replaces by harmonics of the rotation (engine orders), as can be seen from Figs 20 and 21.

   -  

Frequency (Hz)

–500

–500

n=0

–400

–400

–300

–300

–200

–200

–100

–100

0

0

100

100

200

200

300

300

400

400

500 –500

0

5

15

10

20

500

–400

–400

–300

–300

–200

–200

–100

–100

0

0

100

100

200

200

300

300

400

400

500

0

n=1

0

–500

n=2

5

10 15 Time (s) –500 n=4 –400

20

500

595

10

5

15

20

10 15 Time (s)

20

n=3

0

5

Frequency (Hz)

–300 –200 –100 0 100 200 300 400 500

0

5

15 10 Time (s)

20

Figure 17. Time–frequency distributions (STFT) decomposed into forward and backward and for n = 0, 1, 2, 3, 4 nodal diameters; log amplitude.

.   . . 

596 (a)

(b)

Amplitude

+f (Hz) forward

0

Amplitude

+f (Hz) forward

–f (Hz) backward

0

–f (Hz) backward

Time Time First engine order

Aliasing (spatial)

n=2 bwd n=2 fwd

n=1 fwd

Second engine order

(c)

(d)

Amplitude

+f (Hz) forward Time

Third engine order (bwd)

0

Time

Interference (fwd)

–f (Hz)

Symmetric

Interference (fwd/bwd)

n=3 bwd

Third engine order

Amplitude

+f (Hz)

0 –f (Hz) backward

Interference (bwd)

Fourth engine order

Figure 18. Three-dimensional plot of the STFTs in Fig. 17 for (a) n = 1, (b) n = 2, (c) n = 3, (d) n = 4; log amplitude.

500 400 300 200 100 0 –100 –200 –300 –400 –500

0

50

100 Time (s)

150

Figure 19. Spectrogram of the n = 2 component transform to the disc (rotating) co-ordinates.

   -  

597

12

Engine order

10 8 6 4 2 20

40

60

80 100 Time (s)

120

140

Figure 20. Engine order—ordinary time distribution using a single sensor.

Applying the wavelength separation method introduced in Section 3.2 separates the different engine orders as each of them belongs to a different number of nodal diameters. In the example shown in Fig. 21, which goes one step further and separates the n = 2 order into forward and backward components, the second engine order is clearly separated and all other engine orders visible in Fig. 20 have disappeared. Another feature of the

10 8 6

Engine order

4 2 0 –2 –4 –6 –8 –10

20

40

60

100 80 Time (s)

120

140

Figure 21. Forward/backward engine order—time distribution for the second engine order (n = 2).

598

.   . . 

separation method is the distinction between forward and backward progression which can clearly be seen in Fig. 21. 5. CONCLUSIONS

In this paper a new tool for directional and spatial decomposition of vibrations is presented. The suggested method proves very useful for analysing the vibration response of rotating structures where the forward–backward as well as the different wavelengths can be decomposed. The presented results show some initial experiments with the new method, and it is expected that other aspects will be learnt while testing and refining the method. The method has the disadvantage of using several sensors in contrast to one in the usual Campbell diagram; in return, a much better decomposition of the vibration patterns is achieved and thus our understanding of the analysed response increases. Further refinement of the method is expected in the experimental procedure where static deflections (pre-bending) will be removed by using signal processing methods. Better tracking procedures of engine orders in both rotating and stationary co-ordinates are currently under development. REFERENCES 1. D. J. E 1969 Journal of Sound and Vibration 9, 65–79. The effects of detuning upon the forced vibrations of bladed disks. 2. C L. 1989 Proceedings of the IEEE 77, 941–981. Time–frequency distributions—a review. 3. G. G 1992 Journal of Sound and Vibration 155, 385–402. A fast modal technique for the computation of the Campbell diagram of multi-degree-of-freedom rotors. 4. D. E. N 1993 An Introduction to Random Vibration and Spectral and Wavelet Analysis, 3rd edn. Longman Group. 5. A. M 1996 The 6th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Honolulu, Hawaii. Dynamics of anistropically supported rotors. 6. E. K 1993 Dynamics of Rotors and Foundations. 7. J. C. M and J. K. H 1994 Mechanical Systems and Signal Processing 8, 243–258. A comparison between the modified spectrogram and the pseudo-Wigner–Ville distribution with and without modification. 8. Y. H, H. M and S. S 1985 Journal of Sound and Vibration 102, 457–472. Modal response of a disk to a moving concentrated harmonic force. 9. I. B, D. J. E and P. S 1995 ASME, Winter Annual Meeting, October, Boston. Multi-dimensional directional spectrograms and Campbell (Zmod) diagrams for rotating machinery. 10. V. C. C 1988 Signal Processing—The Modern Approach. McGraw-Hill. 11. P. A. J 1984 Deconvolution. San Diego: Academic Press. 12. I. B, P. S, D. A. R and D. J. E 1994 International Conference on Vibration Measurement by Laser Techniques, October, Ancona, Italy. A laser-based measurement system for measuring the vibration on rotating discs.

APPENDIX A: UNBALANCE RESPONSE OF A LONG RIGID ROTOR WITH ANISOTROPIC BEARINGS

The rigid rotor in Fig. 2 behaves according to the equation of motion: IT b + IpVa˙ + Kb b = Mb (t) IT a¨ − Ip Vb + Ka a = Ma (t) where: Kb = 12 Kx L 2, Ka = 12 Ky L 2. A conical unbalance, uu manifests itself as two moments:

(A1)

   -   Ma (t) = uu V (It − Ip ) sin Vt

Mb (t) = uu V (It − Ip ) cos Vt

2

2

599 (A2)

Therefore, neglecting the transient solution, we obtain: a(t) = A(V) sin Vt =

uu V 2((Ip − It )Kb + V 2(It2 − Ip2 )) sin Vt D

b(t) = B(V) cos Vt =

uu V 2((Ip − It )Ka + V 2(It2 − Ip2 )) cos Vt D

(A3)

where: D = (It2 − Ip2 )V 4 + (Ka + Kb )It V 2 + Ka Kb . Backward whirl occurs when A(V)B(V) Q 0. For this system assuming Ka Q Kb , the speed range in which backward whirl occurs is: Ka Kb QVQ It + Ip I t + Ip

(A4)

APPENDIX B: FORWARD/BACKWARD SEPARATION USING A COMPLEX COMBINATION OF TWO MEASUREMENTS

Two sensors, s1 (t) and s2 (t), spaced 90° apart, can be combined to create, s˜ (t) = s1 (t) + is2 (t), in order to compute the double-sided spectrogram directly.  :   For this case we have: s1 (t) = Af cos Vt,

s2 (t) = Af sin Vt

(B1)

The Fourier transform of F[s˜ (t)], has one line at v = V F[s˜ (t)] = 2pAf d(v − V)

(B2)

 :   For this case we have: s1 (t) = Ab cos Vt,

s2 (t) = −Ab sin Vt

(B3)

The Fourier transform of F[s˜ (t)], has one line at v = −V F[s˜ (t)] = 2pAf d(v + V)

(B4)

 :  /  For this case we have: s1 (t) = A1 cos Vt + A2 sin Vt,

s2 (t) = A3 cos Vt + A4 sin Vt

F[s˜ (t)]/p = [i(A3 − A2 ) + (A1 + A4 )]d(v − V) + [i(A2 + A3 ) + (A1 − A4 )]d(v + V)

(B5) (B6)

   The standing wave, usually used in antenna transmission applications, measures the ratio between the travelling and the standing waves. This ratio is defined as: SWR =

=Af = + =Ab = =Af = − =Ab =

where Ab = (A2 + A3 ) + (A4 − A1 ) and Af = i(A2 − A3 ) + (A1 + A4 ).

(B7)

.   . . 

600

The SWR is the ratio between the major and minor principle axes of the ellipse shown in Fig. 3. For a pure standing wave, this ratio tends to infinity (=Af = = =Ab =); a pure travelling wave (forward or backward) yields SWR = 1.       (  ) A different separation than forward/backward can be performed. As mentioned before the travelling waves and standing waves have a different effect on the structure, hence we can separate the response into a standing wave part: As = >Af = − =Ab >

(B8)

At = =Af = + =Ab = − >Af = − =Ab >

(B9)

and a travelling wavepart:

which simplify into: At =

6

2=Ab =, =Af = q =Ab = 2=Af =, =Af = Q =Ab =

(B10)

APPENDIX C: FORWARD/BACKWARD SEPARATION OF SHAFT WHIRL

It can be shown that the forward and backward whirl amplitudes can be computed by means of solving three equations for the three unknowns, p1 (t), p−1 (t), C: Therefore, we write three equations by substituting three values for t in equation (5a). for u = 0

s0 (t) = p1 (t) + p−1 (t) + C

p 2

sp/2 (t) = ip1 (t) − ip−1 (t) + C

for u = p

sp (t) = −p1 (t) − p−1 (t) + C

for u =

(C1)

Similarly to equation (D2) (apart from the scaling factor), we obtain: p1 (t) = 12 (s0 (t) − isp/2 (t)) p1 (t) = 12 (s0 (t) + isp/2 (t))

(C2)

APPENDIX D: FORWARD/BACKWARD SEPARATION OF SHAFT WHIRL

A convenient way to derive the spatial Fourier series is via the DFT. This option can be used for the case where the sensors are equally spaced and cover the whole space. The DFT in its matrix form for the four-sensor case, is obtained by using equation (9). We substitute equation (8) for the measured signals to obtain: p0 (t) 1 F G p1 (t) J G K G1 G G=G G p2 (t) G G 1 f p−1 (t) j k 1

1 −i −1 i

1 1 −1 i 1 −1 −1 −i

x0 + f1 (t) L GF G y0 + f2 (t) J G GG G GG −x0 − f1 (t) G lf −y0 + f2 (t) j

(D1)

   -  

601

Evaluating equation (D1), we obtain: p0 (t) = 0, 1 2 1 2

p2 (t) = 0

p1 (t) = f1 (t) − if2 (t)

p2 (t) = 12 p−1 (t) = f1 (t) + if2 (t)

(D2)

Here: forward component p1 (t) = Af backward components p−1 (t) = Ab Application of equation P21 (v, t) = S1 (v, t) 2 iS2 (v, t).

(1)

to

(D2),

yields

(D3) equation

(11a, 11b):