#### Abstract

Rotordynamic stability is crucial for high performance centrifugal compressors. In this paper, the weighted instrumental variable (WIV) based system identification method for rotating machinery stability is investigated based on a sine sweep forward excitation with an electromagnetic actuator. The traditional multiple input multiple output (MIMO) frequency response function (FRF) is transformed into a directional frequency response function (dFRF). The rational polynomial method (RPM) combined with WIV is developed to identify the rotor’s first forward mode parameters. This new approach is called the COMDYN method. Experimental work using the COMDYN method is carried out under different rotating speeds, oil inlet temperatures, and pressure conditions. Two sets of bearings with preloads 0.1 and 0.3 are investigated. A numerical rotor-bearing model is also built. The numerical results correlate reasonably well with the experimental results. The investigation results indicate that the new method satisfies the desired features of rotating machine stability identification. Furthermore, the system log decrement was improved somewhat with the increase of oil inlet temperature. The increase of oil supply pressure affects the rotor-bearing system stability very slightly. The results of this paper provide new and useful insights for potentially avoiding instability faults in centrifugal compressors.

#### 1. Introduction

Assuring the rotordynamic stability is of crucial importance for high performance centrifugal compressors. Improved rotordynamic stability assurance will be welcomed by both end users and centrifugal compressor manufacturers. The stability estimation method developed in this work combines complex dFRF, RPM, and WIV and is called the COMDYN method. It contributes to meeting the objective of better stability identification. Over the years, rotordynamicists have successfully followed the dual approach of excitation source reduction and vibration absorption through increased effective damping [1]. Now the increasing introduction of magnetic actuators for rotor stability testing can contribute new methods of implementing stability testing, using a theoretical background that has existed for a long time, but not many test results have been published.

Large numbers of papers have focused on the rotordynamic stability, by both numerical and experimental work, on the component and system levels. On the component level, Childs [2] presented a method used to analyze the dynamic parameters of seals based on a bulk-flow method. Later, a three-dimensional CFD method has become widely used to analyze the stiffness, damping, and mass properties of gas seals [3] and compressor impellers [4]. Also, experimental works were carried out to investigate the gas seal factors affecting these dynamic parameters [5]. Much attention was also focused on the dynamic performance of fluid film bearings. Dimond et al. [6] gave a review of tilting pad bearing theory, including synchronous and nonsynchronous models. Future trends in bearing analysis were discussed and a path for experimental verification is proposed. Childs et al. [7] investigated the frequency dependent properties of the bearing stiffness and damping coefficients. Delgado et al. [8] investigated the dynamic response of direct lubrication for a 5-pad, rocker-back pivot tilting pad bearing in a controlled motion test rig.

On a system level, Baumann [9], Moore et al. [10], Sorokes et al. [11], Bidaut et al. [12], and Takahashi et al. [13] investigated the stability performance of centrifugal compressors under full-scale, full load shop test conditions using the frequency domain method through exciting the rotor with an active magnetic bearing. Cloud [14] presented a time domain method named multiple output backward autoregression (MOBAR) to investigate the stability of a flexible rotor supported by tilting pad bearing. Later, Pettinato et al. [15] investigated shop acceptance testing of compressor rotordynamic stability and the theoretical correlation. The results indicate that the estimated results between the MOBAR method and prediction error methods (PEM) are in very close agreement. Ikeno et al. [16] investigated the stability of a compressor rotor supported on four-pad tilting pad bearings and five-pad tilting pad bearings, respectively. For the five-pad tilting pad bearing case, they also researched the difference between load on pad (LOP) and load between pads (LBP). The outcome of that paper indicated that the LOP style provided better stability performance for the compressor rotor. Squeeze film dampers have been widely used to improve the stability of rotor. Gerbet et al. [17] carried out the rotordynamic evaluation of a full scale rotor on two tilting pad bearings with integral squeeze film dampers. Besides the experimental damping ratio identification to assure the stability of rotor in operating conditions, many researchers including Gupta et al. [18] and Dimond et al. [19] investigated how to predict the stability of the rotor by an analytical method. The results of these research efforts show that the analytical method can predict the rotordynamic stability result fairly well. Although these previous research projects are excellent, there are still many challenges such as forward and backward overlap, uncertainties in the numerical simulation, and high noise in the shop test environment.

The instrumental variable (IV) method, with the reputation of unbiasedness as an estimator, was first introduced into the economic parameters estimation. Later, Kendall and Stuart [20] introduced it into mechanical systems to estimate the mechanical parameters. Fritzen [21] applied the IV method to estimate the mechanical parameters in the mass-damping-stiffness dynamic equations (or the MCK-equations) and obtained a more accurate estimation than the ordinary least square (OLS) technique. In the system parametric identification area, Gao and Xu [22] developed a WIV method to estimate the autoregressive moving average (ARMA) statistical model with noise input and output in time domain, obtaining much better estimation compared to the IV method, especially when the signal to noise ratio (SNR) is low.

In this paper, the traditional MIMO FRF matrix is transformed into the dFRF matrix from the real number field to the complex field using a transformation matrix. This eliminates or depresses the influence of forward and backward modal superposition. The theoretical basis for this complex dFRF method is provided in Lee [23, 24] with an extensive, detailed presentation of the complex form of the equations of motion. However, in the complex dFRF method developed in this paper, pure sinusoidal excitations are used as compared to the pseudorandom excitations recommended by Lee [23, 24]. Finally, the modal parameters are obtained for stability estimation by applying the RPM, combined with WIV method, to fit and identify the dFRF in the complex field. This method is called the COMDYN method of stability identification. The test rig was the same as the one employed by Cloud [14] but the testing method was different.

In the experimental tests, two sets of bearings with preloads 0.1 and 0.3 are tested for their damping effect, as measured by the system log decrement, on the rotor-bearing system. This work also compares rotordynamic stability under different bearing loading configurations (LBP/LOP) both experimentally and theoretically. A numerical model was built to simulate the log decrement of the rotor-bearing system. The numerical results correlate with the experimental results reasonably well. Also, the effect of rotating speed as well as bearing lubrication oil pressure and temperature on the stability is investigated by both experimental and numerical methods. The result indicates that as the rotating speed increases, the log decrement decreases, approximately linearly, as the rotating speed becomes larger than the first natural frequency. It also shows that the system log decrement increases slightly as the supply temperature of the bearing lubricating oil increases. Although the effect of oil inlet pressure also has this tendency, the degree is very slight again. The LOP configuration provides more damping to the rotor system and the oil inlet temperature affects the log decrement in this configuration somewhat more than the LBP configuration. The results of this paper provide some new insights for potentially avoiding instability faults in centrifugal compressors.

#### 2. Weighted Instrumental Variable/Complex dFRF-COMDYN-Identification Method

Figure 1 shows a schematic of a rotor-bearing system with multiple degrees of freedom (DOF). This example includes four discs and a single bearing for illustration purposes. It can be described by the following differential equation: In (1), the external forces represent unbalance forces, external fluid excitations such as net hydraulic load from each stage, and other effects. The force may be either a constant vector or a function of displacement or velocity of rotor vibration. The rotor nodal displacement vector can be described as , and it includes both translational and rotational degrees of freedom. Although the nodal displacement vector has four DOF, it is difficult to measure the rotational and in the process of modal identification for multidegree rotor-bearing systems. Therefore, the rotational DOF are neglected in the testing procedure.

The number of measurement locations for most rotor systems is limited, which results in a limited number of rotor responses that can be used to characterize the system. Equation (1) can be written more compactly as Here, the displacement response corresponding to the single measurement plane in Figure 1 can be described as Similarly, the excitation force applied at a single excitation plane can be described as Following Lee [23, 24], using a Fourier transform of (2), the FRF matrix for the measurement plane and the excitation plane is then found as The FRF matrix contains information concerning the system natural frequencies and damping ratios, but this information is not guaranteed to be complete. For example, if the node of a mode occurs at the measurement plane, then that mode will not be observed. The individual FRFs can be written in terms of a polynomial representation. When properly decomposed, this representation contains the system modes and the natural frequencies. This modal representation is given as In (6), the overbar notation indicates the conjugate of complex variable.

Usually, the PEM for estimating the FRF is first used to give the error function, including the parameters to be identified. The PEM is then used to obtain the minimum conditions of the 2-norm of error vector, which is the familiar OLS techniques. However, the OLS is a biased estimator. This paper applies an unbiased estimator, the WIV method, to estimate the parameters, which yields a more accurate estimate, even in high noise circumstances.

One way to obtain the error function involves using the RPM to establish the FRF model. For the FRF in (6), the RPM can be described as Rewriting the above formula with rational polynomials in the numerator and denominator gives Equation (8) provides the theoretical model for the FRF without any reduction, so it is expected to have high accuracy.

After establishing the model, the error function can be written as where is the FRF value measured for the frequency point and is the total number of measurements.

Substituting (8) into (9) and premultiplying the equation by the denominator function, , then the error function becomes the linearized form Combining all of the single error functions to form the matrix description yields where

Here the footnote represents a complex parameter.

Because the parameter vector is a real number and the other vectors and matrix are complex numbers, (11) needs to be transformed to real functions of the form Then, rewrite (13) as

Traditionally, the parameters in (14) can be estimated using OLS, and the estimation result is .

To obtain an unbiased estimation, the IV matrix is brought into the estimation process. The defined IV matrix satisfies the following conditions [21]: where , is nonsingular, and is the observation time.

Next premultiply (14) by with the result Using (16) and noting that according to the definition of the IV, the unbiased estimation equation can be written as

Since the term is nonsingular, it is invertible. Solving (17) to obtain the estimated parameters vector yields

The main problem is how to choose the IV matrix. Fritzen [21] presented an easy way to obtain the matrix . In addition, the measured FRF values have different SNR at different measured frequency points (). Thus, (14) can be weighted as where . The weighting values are the same for the real and imaginary equations at one frequency point.

Then, the objective function of the 2-norm error vector becomes

Obviously, the OLS estimation for (20) is

Next, the WIV for parameter estimation becomes

Adapting the same process described in Fritzen [21], the parameter vector will be well estimated.

Usually, the frequency point with the greatest amplitude in the amplitude spectrum has the highest SNR. Thus, weighting function values can be defined as where is the variation factor of the weighting function. When , the 2-norm of error vector is just linearly weighted.

#### 3. Stability of Rotor Supported by Tilting Pad Bearing

##### 3.1. Description of the Test Rig

To investigate the stability of a flexible rotor supported by a set of two tilting pad bearings, a test rig was built some years ago in the ROMAC lab [14]. Figure 2 presents the picture of the test rig. The test rig consisted of a 1,549.4 mm long rotor supported by two tilting pad journal bearings and three active magnetic bearings (AMBs). The test apparatus was mounted on somewhat flexible steel rails in a large concrete block. The rotor was separated from the motor drive by a coupling. Two of the magnet actuators, designated as Actuators 1 and 2, were located between the journal bearings for the application of either cross-coupled stiffness or control forces. These locations mirror those on common industrial between-bearing machines, where various destabilizing components, such as oil and labyrinth seals, exist [14]. The third magnetic actuator (called the shaker) was located at the outboard end, where previous industrial stability measurement investigations, as referenced above, applied an external excitation. Actuator 1 was used to provide the sine sweep exciting force in this paper.

Figure 3 shows a picture of the tilting pad bearing assemblies. The spherically pivoted bearings are axially split for easier installation. The bearing’s outer housings were provided with individual nozzles to feed oil to each pad. Oil then leaves the bearing housing through drain holes behind the pads and axially through two floating ring end seals. The effects of bearing preload (0.1 and 0.3) and load configuration (LOP and LBP) on the rotor stability were investigated. The assembled clearance of two bearings was maintained at approximately 172.7 microns.

As Figure 4 indicates, the log decrement measurements were carried out by applying a forward harmonic force with a sweep over a frequency range from 20 Hz to 200 Hz. To minimize any distortion of the measured resonance curves, a slow frequency increase rate of 1 Hz/sec was used. The level of the excitation force applied by the exciter was adjusted in order to ensure high quality responses. Also, some backward harmonic excitations were employed.

##### 3.2. Results of the Experimental Test

Figure 5(a) indicates the full spectrum waterfall plot during the test rig run-up process. The preload of the five tilting pads bearing is 0.1 and the configuration was LOP in Figure 5(a). The sampling rate is 3 kS/sec. The recorded averaged time series of the measured shaft vibrations were transformed into the electromagnet actuator (shaker) pole face axes and then into fast Fourier transform (FFT) signals. Thereby, resonance peaks were determined. To make the acquired rotor sates stable, the rotating speed was increased at a rate of only 10 rpm/second. From this plot we can see that in most of the speed range, the direction of the whirl obit is forward. The amplitude reaches its peak at 5,160 rpm. The emergence of backward whirl near the first natural frequency is mainly because of the bearing stiffness asymmetry in the and directions. This is a clear demonstration of a backward mode response excited by a forward excitation. Few backward mode responses under forward mode excitations have been shown in the literature. Figure 5(b) shows the full spectrum waterfall under forward sine sweep excitation. The parameters of the bearing are same as those in Figure 5(a). Figure 5(c) shows the full spectrum waterfall under a forward sine sweep excitation. The preload of the 5-pad tilting pads bearing is 0.3 and the LOP configuration for this case. Comparing these two diagrams, it can be seen that the amplitude of both backward and forward vibration at 10000 rpm is lower with 0.1 preload bearing than that with 0.3 preload bearing. Thus it can be concluded that the bearing with preload 0.1 has more isotropic properties and damping than that with preload 0.3, as one would expect.

**(a) Waterfall for speed increase process (5-pad, LOP, preload 0.1)**

**(b) Forward sine sweep at 10,000 rpm (5-pad, LOP, preload 0.1)**

**(c) Forward sine sweep at 10,000 rpm (5-pad, LOP, preload 0.3)**

**(d) Backward sine sweep at 9,000 rpm (5-pad, LOP, preload 0.3)**

**(e) Forward sine sweep at 9,000 rpm (5-pad, LOP, preload 0.3)**

Figure 5(d) indicates the full spectrum waterfall of the rotor, supported by 5-pad tilting pad bearings with preload 0.3 and LOP, under backward sine sweep conditions at 9,000 rpm. Figure 5(e) indicates the full spectrum waterfall of the rotor supported by the 5-pad tilting pad bearings with preload 0.3 and load on pad under forward sine sweep conditions at 9,000 rpm. These two plots show that forward excitation leads to high amplitude backward whirl and backward excitation can also lead to prominent forward whirl because of the highly anisotropic stiffness properties of the bearing dynamic performance which is evidenced by the split responses near the first natural frequency. Figure 6 shows the frequency response function (FRF) under sine sweep excitation with 0.3 preload and LOP configuration at 10,000 rpm. A fourth order fitting method is used to identify the natural frequency and log decrement of first forward mode.

#### 4. Rotordynamic Model of the Test Rig

To predict the damped natural frequency and log decrement of the test rig, two analysis packages developed by the ROMAC Lab were used. The bearing software code MAXBRG [25] was used to analyze the dynamic characteristics of bearing based on a thermoelastic-hydrodynamic (TEHD) analysis. Different from some commercial codes, it also included the effects of pad pivot flexibility. FORSTAB [26], based on a finite element beam (FEM), was used to analyze the rotor dynamic performance of the test rig. Figure 7 indicates the model of the test rig rotor which is composed of 63 beam elements, similar to Dimond [19]. The disks and rotor lamination stacks of the AMBs are modeled as lumped masses. Figure 8 indicates the whirl model of the first forward natural frequency when the rotor rotates at 10,000 rpm. The bearing preload is 0.3 and the LOP configuration is considered. The plot in Figure 8 shows that the natural frequency is 5,206 rpm and with the log decrement 0.102. The somewhat flexible steel rails supporting the test rig above the concrete block were not modeled, as they are extremely difficult to characterize.

#### 5. Factors Affecting the Stability of the Rotor System

##### 5.1. Log Decrement at Different Rotating Speeds

Theoretically, the log decrement of the first forward procession mode will decrease with the increase of the rotating speed. This tendency is mainly due to the dynamic performance of the bearing. With the increase of speed, the damping coming from the oil film decreases while the stiffness usually increases. The log decrement of the rotor-bearing system at different rotating speeds from 7,000 to 10,000 rpm was tested, with the results that are shown in Figure 9. Figure 9(a) also provides the simulation result based on the software MAXBRG (for the tilting pad bearings) and FORSTAB (for the rotor response) packages at the rotating speed range from 4,500 to 12,000 rpm. It is found that the log decrement has a small increase with rotating speed below the first natural frequency. This tendency is reversed for the speed above the natural frequency. It is also seen that the numerical results correlate with the experimental ones very well. Figure 9(b) shows the difference between the numerical and experimental result. Some of the differences between the simulated and measured log decrement values are due to the unmodeled flexible steel mounting rail.

**(a) Log decrement of the test rig at different speed**

**(b) Differences between simulated and expermental result**

##### 5.2. System Log Decrement at Different Oil Inlet Temperatures

Under practical situations, maintaining the oil inlet temperature constant in the whole performance test process is difficult, as Pettinato et al. [22] indicated. On the other hand, the change of oil inlet temperature provides one convenient way to somewhat improve rotor stability. This paper investigates the level of the rotor-bearing system stability due to the inlet oil temperature effect. The results are shown in Figure 10. The LOP load orientations were conducted with a bearing preload of 0.1 at the operating speed of 10,000 rpm. The supply temperature points are 43°C, 45°C, 48°C, and 51°C, respectively. It is shown that increasing the bearing oil supply temperature slightly raises the log decrement of the system. This is mainly because the increasing of oil inlet temperature decreases the stiffness of bearing. Also, the increase of oil supply temperature will decrease the viscosity of oil, which will decrease the damping of bearing. It indicates that the rotor-bearing system effective damping will increase somewhat more in the testing temperature range. Furthermore, the numerical and experimental results correlate well when the oil inlet temperature is lower than 45°C.

**(a) Log decrement of the test rig at different oil inlet temperature**

**(b) Derivation between simulated and expermental results**

##### 5.3. Log Decrement at Different Oil Inlet Pressures

The oil inlet pressure is potentially another factor affecting the system stability. The paper presents tests of the system log decrement at four different points. The inlet pressures were 1.0 bar, 1.25 bar, 1.5 bar, and 1.75 bar. Both LBP and LOP load orientation tests were conducted, respectively, with a bearing preload of 0.1. From the results shown in Figure 11, it can be found that the damping ratio increases with the increase of inlet oil pressure but the degree is very slight. This is mainly because the increasing of oil inlet pressure will reduce the cavity region so that the bearing damping coefficient will increase.

Figures 10 and 11 indicate that the rotor stability can be affected by either the oil inlet temperature or pressure. To this point, the stability tendency corresponding to the two factors that can be described by a mesh plot was given in Figure 12. 16 test points are selected which are formed by 4 pressure points and 4 temperature points as indicated above. Figure 12(a) shows the log decrement with oil inlet temperature and inlet pressure of LOP and Figure 12(b) gives results for the LBP.

**(a) LOP**

**(b) LBP**

#### 6. Summary

This paper investigates a new method of determining the damping properties of rotors on fluid film bearings using a frequency domain forward sine sweep excitation, the complex dFRF combined with WIV—the COMDYN method. The complex frequency domain approach provides modal parameters with complex zeros and poles. The approach also determines a clear distinction between forward and backward modal responses including a new method of decomposing the response into forward and backward modal responses. It should be noted that the MOBAR method can also handle both forward and backward modal responses. However, the new COMDYN method differs from the MOBAR method in that it does not use the blocking technique to measure system log decrement during transient vibration decay as well as utilizing both input and output signals. The approach eliminates the influence of forward and backward excitation modal overlap. The new method can potentially be used to determine bearing stiffness and damping coefficients. It was shown in this work that pure sinusoidal excitation gives very good estimations of the system damping which is different from the theoretical literature development using pseudorandom excitation. No test results similar to this new complex dFRF method using the magnetic exciter on a flexible rotor and sinusoidal excitation have been published.

Further, the WIV approach was developed and implemented experimentally as an essential component of this research. Also, the work studied two key relationships between the system log decrement with the increase of rotating speed and oil inlet temperatures by experimental and numerical methods, based on two sets of tilting pad bearing with preloads 0.1 and 0.3. In the experimental research work, the COMDYN method is used to identify the log decrement of the rotor-bearing system at different rotating speeds, oil inlet temperatures, and inlet pressures with both LOP and LBP load orientations, respectively. Under the experimental conditions carried out in this paper, it can be concluded that the system log decrement decreases with the increase of rotating speed almost linearly after the operating speed passes through the first critical speed, similar to the results in [16]. The bearing contributes more system damping when it is in the LOP configuration than when in the LBP configuration when the other conditions are the same. In general, the bearing with preload 0.1 provides more damping to rotor-bearing system than that with preload 0.3. Furthermore, the stability of the rotor-bearing system can be improved by increasing the oil inlet temperature slightly. The outcomes of this paper with the COMDYN approach can contribute to understanding the process of shop tests for centrifugal compressor stability evaluation. It has also shown some limited potential for improving that stability by some relatively simple changes in tilting pad bearing inlet temperatures and inlet pressures.

#### Nomenclature

AMBs: | Active magnetic bearings |

ARMA: | Autoregressive moving average |

CFD: | Computational fluid dynamics |

dFRF: | Directional frequency response function |

DOF: | Degrees of freedom |

FFT: | Fast Fourier transform |

FRF: | Frequency response function |

LBP: | Load between pads |

LOP: | Load on pads |

MIMO: | Multiple input multiple output |

MOBAR: | Multiple output backward autoregression |

OEM: | Original equipment manufacturer |

OLS: | Ordinary least square |

RPM: | Rational polynomial method |

SNR: | Signal to noise ratio |

TEHD: | Thermoelastic-hydrodynamic |

WIV: | Weighted instrumental variable. |

*Arabic*

a: | Numerator in rational polynomials |

B: | Backward precession |

b: | Denominator in rational polynomials |

C: | Damping matrix |

c: | Complex number |

e: | Error vector |

f: | External forces applying on the rotor |

F: | Forward precession |

G: | Gyroscopic matrix |

H: | Frequency response function matrix |

: | Measured frequency response function matrix |

H_{mn}: | Frequency response in the m direction due to an excitation in the n direction |

j: | |

K: | Stiffness matrix |

: | Weighting value matrix |

k: | Weighting function value |

M: | Mass matrix |

S_{N}: | Number of support points |

s: | Laplace frequency point |

: | Parameter vector |

: | Vector |

: | Instrumental variable matrix |

, , : | Rotor displacement, velocity, and acceleration |

x: | Rotor horizontal displacement degree of freedom; modal displacement |

y: | Rotor vertical displacement degree of freedom; adjoint matrix. |

*Greek Symbols*

α: | Rotor rotation about x-axis degree of freedom |

β: | Rotor rotation about y-axis degree of freedom |

λ: | Eigenvalue |

ω: | Angular frequency. |

*Subscripts*

b: | Bearing terms, finite element model |

c: | Cross-coupled stiffness terms, finite element model; complex number symbol |

k: | kth measured point |

m: | Lumped mass terms, finite element model |

s: | Shaft terms, finite element model. |

#### Conflict of Interests

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

#### Acknowledgments

The work described in this paper was supported financially by the Natural Science Foundation of China (51135001) and “973” Program (2012CB026000). These supports are gratefully acknowledged.