Abstract

With the increase of the centrifugal compressor capability, such as large scale LNG and CO2 reinjection, the stability margin evaluation is crucial to assure the compressor work in the designed operating conditions in field. Improving the precision of parameter identification of stability is essential and necessary as well. Based on the time-varying characteristics of response vibration during the sine-swept process, a short-time Fourier transform (STFT) filter was introduced to increase the signal-noise ratio and improve the accuracy of the estimated stability parameters. A finite element model was established to simulate the sine-swept process, and the simulated vibration signals were used to study the filtering effect and demonstrate the feasibility to identify the stability parameters by using Multiple-Input and Multiple-Output system identification method that combines the prediction error method and instrumental variable method. Simulation results show that the identification method with STFT filter improves the estimated accuracy much well and makes the curves of frequency response function clearer. Experiment was carried out on a test rig as well, which indicates the identification method is feasible in stability identification, and the results of experiment indicate that STFT filter works very well.

1. Introduction

With the ever increasing capacity in today’s ethylene, refinery, LNG, and high pressure CO2 reinjection, it is critically important to optimize the design of rotors, impellers, and journal bearings to improve rotor dynamic stability in field. Definitely, it is crucial to ensure that centrifugal compressors are rotor dynamically stable under operating conditions before they are delivered.

Rotor stability identification methods for rotor-bearing systems have been investigated by many researchers, and they estimated rotor stability during the shop test based on both frequency and time domain signal analysis [1, 2]. Takahashi et al. [1], Kocur and Cloud [2], and Pettinato et al. [3] did a plenty of work on shop rotor dynamic testing by using the electromagnetic bearing as a shaker to excite the rotor system. Cloud [4] developed the Multiple Output Backward Autoregressive (MOBAR) method for time domain signals and carried out significant work on the rotor stability investigation. A stability identification method combining the PEM and IV method to estimate the frequency response function (FRF) in frequency domain was developed by Wang et al. [5], and numerical simulation results show that the method can identify the damping ratio of the first forward and backward modes with high accuracy, even in a severe noise environment; furthermore, two example centrifugal compressors (nine-stage straight-through and six-stage back-to-back) were employed to verify the feasibility of identification method in industrial configurations as well. With this method, Wang et al. [6, 7] carried out experiments to investigate two sets of bearings with preloads 0.1 and 0.3 under different rotating speeds, oil inlet temperatures, and pressure conditions; Wang et al. [8] studied the effect of specific load on rotor stability as well. Recently, a new method, Operational Modal Analysis (OMA) method, which was used in structure modal analysis successfully, is very popular for its reputation of being convenient to use and having no need of any artificial excitation to the rotor system. Guglielmo et al. [9] and Carden and Morosi [10] used the method to evaluate the stability of rotating machinery. The only drawback is that OMA requires the rotor running at steady state condition and a sufficient aerodynamic excitation force acting on the compressor as pointed out by Guglielmo et al. [9].

The key of identification method is of high antinoise ability and identification accuracy, and improving quality of the raw vibration data is crucial to the reliability of the identified results. Also, Kocur and Cloud [2] pointed out that the estimation process faced the challenges that data contain additional response from the significant presence of the immeasurable, internal excitations within the machine such as unbalance, misalignment, aerodynamics when dealing with the measurement data. The traditional filtering method, like band-pass filter, can retain the signal just in interest frequency range and decay other frequency component signals, but the noise signal mixing in the interest frequency band cannot be removed. Particularly, for the sine-swept response signal, its sinusoidal frequency changes over time, so a filter with a fixed bandwidth cannot eliminate the effect of noise in the interest frequency band. Focusing on the time-varying characteristic of the response vibration, a time-frequency filter, such as short-time Fourier transform (STFT) [11], was studied to eliminate other uncorrelated interference and retain the useful signal component. To some extent, STFT is essentially the technique of using Compact Harmonics Wavelets (CHW) for time-frequency analysis and system identification; many researchers vested effort to the development of CHW theory and its applications. Coifman and Guido [12] did some pioneer works about CHW. Daubechies [13] provided lectures about the wavelets and its applications in signal processing and other fields, in which the comparison between STFT and wavelets was mentioned. Farge et al. [14] applied CHW ideas into the computation of two-dimensional turbulent flows, where the wavelet packet best-basis seems to distinguish the low-dimensional dynamically active part of the flow from the high-dimensional passive components, providing a hope of drastically reducing the number of degrees of freedom necessary to the computation. Le [15] did a practicable framework to represent the wavelets transform of linear dynamics delivering a formula for projecting linear time-invariant dynamics and stochastic processes on the wavelet spaces of a multiresolution analysis; applications of this formula to system and parameter identification for control and modeling of mechanical systems were demonstrated. Latest, Yao et al. [16] used the STFT to estimate the linear frequency modulation signal parameters and the experimental results verified its validity and feasibility. Liu et al. [17] developed a variable window based on STFT to study the rotating machinery vibrating signal and its relationship to rotating speed in the process of startup.

In this paper, the STFT filter was introduced to increase the signal-noise ratio and improve the accuracy of the estimated stability parameters. Numerical and experimental study verified its feasibility and effect in identification process, and using the STFT filter made the FRFs clearer and the estimated results more reliable. An improved Multiple-Input and Multiple-Output (MIMO) system identification was adapted to estimate the modal parameters. Studying a stable and accurate identification method provides possibility and quantitative criteria to study the influence of bearing, seal, load, or other operation parameters on the rotor system.

2. STFT Theory

The STFT is used to observe and analyze the signal whose sinusoidal frequency changes over time. For a given signal , it can be intercepted to series by a window function that is slid along the time axis and is nonzero for only a short period of time, and then the intercepted signal series can be transferred into complex form that contains the phase and magnitude information of the signal over time and frequency by the Fourier transform, resulting in a two-dimensional representation of the signal. The transform for the continuous-time case is written as

Making , then is essentially the Fourier transform of . In practice, signal that we obtained is discrete, so it can be broken up into chunks. Fast Fourier transform (FFT) is used to transform the chunks to form a complex matrix, which records magnitude and phase for each point in time and frequency. Mathematically, this can be expressed as where and is discrete because of using FFT in computation.

Like Fourier transform, the STFT is also invertible and the most widely accepted way of inverting the STFT is by using the overlap-add (OLA) method. Selesnick [11] provided a signal reconstruction method to obtain the inverse STFT and chose a half-cycle sine window for analysis.

Window function is expressed aswhere is the length of window and “” can be taken as the scaling factor that determines frequency resolution of each intercepted signal.

And the th intercepted signal can be written as

Then the STFT is , and its inverse form is .

For the original signal , it can be reconstructed to be as

Submitting (4) to (5), then

Evidently, once , then ; that is, the signal is a perfect reconstruction. Hence, the condition can be simplified as [11]

The half-cycle sine window in (3) satisfies this condition [11]. Furthermore, many other windows, like rectangular window, can also be designed to do the STFT and inverse STFT. In this paper, a half-cycle sine window was chosen for its less spectral leakage when doing FFT.

The length of the window () determines the scaling factor “” that affects the resolution of resulting FFT. In this paper, the length of time of each intercepted signal was set to be 1 s because the step of increasing frequency was 1 Hz/s, and the resolution was 1 Hz which is enough for analysis and filtering the noise. Because of the linear characteristic of the frequency, the scaling factor was chosen to be constant when decomposing and reconstructing the whole signal.

It is worth pointing out that the STFT idea is essentially the forward CHW using windowed sinusoidal bases for analysis, and that the inverse FFT on such transformed signals is just the typical inverse wavelet transform; thus, many other methods based on CHW can complete the work of STFT as well. But STFT is very simple to use without changing the scaling factor “”.

3. Sine-Swept Response Signal

Logarithm decrement of the rotor system is usually estimated using frequency method by applying the sine sweep force on the rotating rotor system. Because the rotor system is assumed to be a linear system, the frequency of response displacement is theoretically corresponding or equal to that of the exciting force during the test process. Other frequencies appearing in the signal are not needed for the modal identification, and they affect the results that are identified with FRF method.

In practice, the exciting force is not a strict chirp signal, and it is accompanied with little white noise, which can excite the vibration near the natural frequency of rotor system, just like the right rectangular portion in the waterfall graph (see Figure 1(a)). Although the vibration is little in a certain time slice, it can accumulate together to be big enough to affect the identification when doing FFT for the whole time signal. Besides, some interferential and irrelevant frequency vibrations that are not needed are excited as well.

The traditional filtering method, like band-pass filter, can retain the signal just in interest frequency range and decay other frequency component signals, but the noise signal mixing in the interest frequency band cannot be removed. Particularly, for the sine-swept response signal, its sinusoidal frequency changes over time, so a filter with a fixed bandwidth cannot eliminate the effect of noise in the interest frequency band. Thus, a time-frequency filter is needed to eliminate noises in the special time and frequency range or points.

4. STFT Filter

The STFT is used to express the characteristics of time-varying signal. Figure 1(b) is the STFT spectrogram that shows the projection of the vibration at different time and frequency points, and the maximum scale of amplitude is set to 4 μm to see the characteristic line clearer. Evidently, two evident lines exist on the time-frequency coordinate plane. The horizontal one stands for unbalance vibration related to the rotating speed, and the slope line, which is the most useful signal component for identifying the rotor system, stands for the vibrations that are excited by the sine-swept force. The other spots or short lines are regarded as noise. Therefore, in order to obtain high quality response signal, we just need to retain the two lines in the coordinate, or set other portion to zero. Moreover, Figure 2 shows the filtering process.(1)Transfer the noise signal to a complex matrix in time-frequency domain by using STFT.(2)Set filtering regulation according to the STFT spectrogram (e.g., retain the elements near the two lines and set other elements to zero in the complex matrix).(3)Reconstruct the processed complex matrix to the time signal by the inverse STFT.

5. Numerical Study

5.1. Finite Element Rotor Model for Sine-Swept Simulation

In order to simulate the sine-swept process, we established a finite element (FE) model of a rotor that is 1558.6 mm in length and the span of supports is 1079.8 mm as shown in Figure 3. Two active magnetic bearings (AMBs) are mounted on the rotor. The FE model contains 38 elements using Timoshenko beam with gyroscopic effect [8], and the support nodes are Node 7 and Node 33. Another AMB used as a shaker to excite the force on the rotor is located at Node 38. The impeller and balance disk are located at Node 17 and Node 20, respectively. Table 1 shows the support parameters of tilting pad bearings for FE calculation.

The FE equation for rotor [18] is usually expressed aswhere is mass matrix, is the damping matrix associated with rotor speed, is the stiffness matrix, and and are the displacement vector and excited force vector, respectively.

By calculating the eigenvalues of (8), the modal parameters of the rotor system can be predicted, and these modal parameters are taken as theoretical values in this paper. Figure 4 shows the first backward and forward damped modal shape at rotating speed of 7,000 rpm. In the shop test of rotor stability, only the first forward whirl mode is concerned.

In the sine-swept process, a transient integration method, Newmark solution technique [18], was used to calculate the transient vibration response as below:

5.2. Simulation of Identification of Rotor Stability by Sine-Swept Method

The FE model of test rig descripted in Figure 3 and the supported parameters are used to simulate the sine-swept process exciting a forward sine-swept force from 0 Hz to 200 Hz by step 1 Hz/s on the non-drive end at Node 38. In this simulation, an unbalance weight of 1000 g·mm is added on the impeller at Node 17 to simulate the practical unbalance condition. Finally, the - and -direction vibration of TPJP Bearing at Node 7 are obtained by using the transient integral equation (9) with calculation step of 1/2048 s. Figure 5(a) shows the time domain waveforms of and , and resonance region appears in the range from 50 to 100 s apparently. Figure 5(b) shows the full scale waterfall of the vibration during the sine-swept process. The positive frequency represents the forward whirl and the negative frequency represents the backward whirl. Apparently, the forward force stimulates the forward mode because the rotor is assumed as a linear system. The slash spectra with resonance are stimulated by the exciting force and the fix vibration frequency at 116.7 Hz is stimulated by unbalance force.

5.3. Study of SFTF Filter in the Identification of Rotor Stability

In order to show the filter effect, we used the simulated signals and of the left bearing as standard or reference signals. The simulated Gauss noise signals were obtained by random function from MATLAB: where is the unit white noise with standard deviation of 1 and is the noise level.

In our simulation, the noise level was set to be 2 μm, 5 μm, and 10 μm separately to study the filtering effect. Figure 6 shows the waveforms of with and without STFT filtering at Node 7, and the filtered signal and original signal are basically identical. Figure 6(b) zooms in for details, and the agreement indicates that the STFT filter has good filtering performance.

Moreover, the FRFs curves were calculated to show the filtering effect as well. Figure 7 shows the theoretical FRFs curves of the rotor system, and Figure 8 shows FRFs curves calculated from the time waveforms with noise. Comparing Figures 7 and 8, we can find the SNR is very low in Figure 8 and characteristic of FRFs is very terrible. After using the STFT filter on the vibration signals, the calculated FRFs in Figure 8 are much better than those in Figure 8, and this improvement indicates STFT filter can restrain the noise frequency component and keep the interest frequency component to make the frequency response characteristic much clearer. The FRFs are nearly identical in Figures 7 and 9, so there is no doubt that using the filtered signals to identify the rotor system will be much reliable.

In order to show the STFT’s effect of improving identification accuracy, we added three levels of noise—μm, 5 μm, and 10 μm on to the theoretical signals. Figure 3 has shown that the theoretical first forward natural frequency and damped log decrement are 76.88 Hz and 0.2643, respectively. The identification method with IV [6] was used to estimate the modal parameters, and we did some improvements described in the appendix.

It should be pointed out that the selected frequency range of FRFs affects the identified results because the signal simples are different. In order to eliminate the effect of selection of frequency range, we made an average algorithm that follows steps as below:(1)Estimate raw modal parameters with a simple selection of frequency range, and then obtain the approximate damped frequency (called center frequency).(2)Define four frequency points based on the as shown in Figure 10.(3)Estimate the modal parameters at multiple intervals.(a)Select f_start from the interval , which has an increment of 1 Hz.(b)Select f_end from the interval , which has an increment of 1 Hz.(c)Do all estimations at all the combined intervals of [f_start, f_end].(d)Sort all the estimated results by log decrement from small to large.(4)Evaluate the statistic of the middle 50% results, obtaining the mean and variance of the modal parameters.

We simulated ten identification processes and Table 2 shows the estimated results of forward modal parameters at different noise levels with STFT method, and the numbers with bold font show the minimum or maximum of the estimated log decrements. Apparently, the estimated results are stable and basically identical to the theoretical values when the noise levels are 2 μm and 5 μm. When the level is 10 μm, as shown in Figure 11, the vibration signal is covered by noise, and the resonance peak disappears. Obviously, the SNR of the signal is very low under this condition. After using the STFT method, the resonance peak appears again, showing the characteristic of the original signal. Although when μm, the results become undesirable at some simulations (such as number 7 and number 8 in Table 2), the estimated results are acceptable in most simulations. Additionally, the vibration signals cannot be so terrible to be used for estimation in shop test or study experiment, and we can do multiple experiments to increase the SNR to reduce or eliminate these undesirable cases. In a word, increasing the SNR is essential for improving the accuracy of estimated results, and STFT filtering makes the FRFs clearer and estimated results more reliable and accurate. At this point, STFT filtering can handle the noise very well in practice.

6. Experimental Study

We did experiment at rotor dynamic test rig described in Figure 12. A slender steel shaft is supported by two steel oil journal bearings lubricated by ISO VG32 turbine oil. The middle of the shaft holds a steel sleeve coupled with labyrinth seals and with a diameter of 180 mm. Near the bearings, two steel balance disks are used to regulate the original vibration of the rotor and provide unbalanced exciting forces. The magnitudes and angles of unbalanced excitation forces can be changed by the location of unbalance blocks, which are fixed on the balance disks. Six rings of copper labyrinth seals are fixed on the interior wall of the cylinder. The cylinder is hung by springs, which can be exchanged to adjust stiffness, from vertical and horizontal directions. The cylinder is made sure to be easily excited only by seal force which has approximately the same frequency as the rotor rotation. A 15 KW DC motor drives the shaft through a 4.5 : 1 ratio speed-increasing gearbox with a disc-type coupling. The rotor speed can be adjusted from 500 rpm to 6000 rpm. An AMB as a shaker is used for applying dynamic load to the cylinder with loads up to 500 N. The dynamic load applied to the cylinder is measured with the load cell located between the stinger and shaker frame. In order to investigate positive preswirl test to explore the potential impact of swirl brakes, four air inlet pipes which feed in the middle of the cylindrical are equally spaced on the circumference. The angle between each inlet pipe and the tangential direction of rotor surface is 45 degrees, and this layout can provide positive preswirl that the inlet gas has the same rotating direction as the rotor.

In our experiment, we applied the sine-swept force from 30 Hz to 120 Hz by step of 1 Hz onto the non-drive end of the test rig and collected the vibration and exciting force data. In order to reduce phase lag, we did not use antialiasing filter to cut off high frequency signals in the acquisition system, and we set a high sample rate of 30 KHz to reduce influence of aliasing of high frequency noise on the interest frequency range (40–100 Hz). The time waveforms were full of high frequency noise because of high sample rate, but the spectra at low frequency port were clear; moreover, what we focus on are the spectra of low frequency portion of the signals, so we filtered the high frequency noise to make the time waveform clearer and resampled the signals at a rate of 3 KHz to reduce the amount of computation when calculating FRF. After obtaining the raw vibration and exciting force data, we did STFT filtering for them to get higher SNR signals for identification of log decrement. A moving half-cycle sine window with 3000 points was used when doing STFT filtering.

Figures 13 and 14 show the time waveforms and STFT spectra of . The STFT spectrum in Figure 14(b) is much clear than that in Figure 14(a), and the characteristic lines become evident. Furthermore, the SNR of the FRF is much better with SFTF filtering than without SFTF filtering as shown in Figure 15. Undoubtedly, FRF with higher SNR provides higher accuracy for parameter estimation and increases the reliability of the estimated results. Thus, the estimated results in Figure 15(b) are more reliable than that in Figure 15(a), even if they have almost the same estimation values. From another perspective, the identification method has high antinoise ability when comparing the estimated results in Figures 15(a) and 15(b).

7. Conclusions

(1)Adapting the moving frequency range for modal parameter estimation and averaging the middle 50% results sorted by log decrement can eliminate the error caused by random selection of range.(2)STFT filter can remove the noises in the special time and frequency range or points, obtaining high SNR signals used for system estimation.(3)Using STFT makes the FRF characteristic clearer and the estimated results more stable and reliable. Moreover, the identification method with STFT can handle the low SNR case very well, and the COMDYN method itself has high antinoise ability according to the comparison between the estimated results with and without STFT filter in experiment.(4)Experimental results verified that the identification method with STFT filter could estimate the log decrement very well.

Appendix

The individual FRFs can be written in terms of a polynomial representation. When properly decomposed, this representation contains the modes and the natural frequencies. The modal representation is given as

The FRF can be can be written as [6]

Use the rational polynomial method to fit this damping system transfer function; then the th FRF can be expressed as

The identification process involves simultaneously curve fitting all the FRFs, producing one identified model of the entire system. Mathematically, all the FRFs have the same denominator .

The error function between original and fitted data can be defined as where is the fitted data, is the test data, and is the number of test points of th FRF.

Multiply at both side of (A.4); then the error function of th FRF can be

Changing (A.5) to matrix form yieldswhere , , and . Also,Assume that there are FRFs, and combine all the error functions of them asSimply, the entire error vector is expressed as whereMaking , , , yields

Here the footnote represents a complex parameter.

Because the parameter vector is a real number and the other vectors and matrix are complex numbers, (A.11) needs to be transformed to real functions of the formRewrite (A.12) as

We can introduce an Instrumental Variable [6] to estimate the parameter vector asOnce we obtain the of this system, we can calculate the modal parameters, frequency and damping ratio (or log decrement), from roots of the polynomial.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

Acknowledgment

The work described in this paper was supported financially by the Natural Science Foundation of China (51135001, 11302133). These supports are gratefully acknowledged.