- Email: [email protected]
Conditioning of regression matrices for simultaneous estimation of the residual unbalance and bearing dynamic parameters
ARTICLE IN PRESS
Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 19 (2005) 1082–1095 www.elsevier.com/locate/jnlabr/ymssp
Conditioning of regression matrices for simultaneous estimation of the residual unbalance and bearing dynamic parameters Rajiv Tiwari Department of Mechanical Engineering, Indian Institute of Technology, Guwahati North Guwahati, Guwahati 781 039, India Received 18 November 2003; received in revised form 16 August 2004; accepted 27 September 2004 Available online 21 November 2004
Abstract A general identification algorithm has been presented for multi-degree-of-freedom rotor-bearing systems. A specific identification algorithm for simultaneous estimation of the residual unbalance and bearing dynamic parameters in rotor-bearing systems has been developed. The identification algorithm is found to be highly ill-conditioned for the conventional measurement techniques. The measurement of responses, while rotating the rotor alternatively in clockwise or counter clockwise direction leads to a well-conditioned identification algorithm. Numerical examples have been used to illustrate the developed algorithms. The algorithm is tested against the measurement noise and found to be robust. r 2004 Elsevier Ltd. All rights reserved.
1. Introduction The residual unbalance estimation in rotor-bearing system is an age-old problem. From the state of the art of the unbalance estimation, the unbalance can be obtained with fairly good accuracy [1–5]. Now the trend in the unbalance estimation is to reduce the number of test runs required, especially for the application of large turbogenerators where the downtime is very expensive [6]. Tel.: +91 0361 2582 667; fax: +91 0361 2690 762.
E-mail address: [email protected] (R. Tiwari). 0888-3270/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2004.09.005
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1083
Bearing dynamic parameters play an important role in the prediction of the dynamic behaviour of rotor-bearing systems. Due to the difficulty in obtaining actual test conditions the well-developed theoretical analysis to obtain bearing dynamic parameters by using the Reynolds equation fails to give satisfactory results [7]. This forces designers to look for the experimentally obtained dynamic bearing parameters. Exhaustive work has been done on the experimental identification of the eight linearised bearing dynamic parameters [8–11]. Most of the work is based on the dynamic force, which is given to the system, and corresponding response is measured. At least two-independent force–response relationships are used to obtain bearing dynamic parameters. Researchers have reported that often the regression matrix of the estimation equation becomes ill-conditioned and this leads to scattering of the estimated bearing dynamic parameters [12]. Advances in the sensor technology and the increase in computing power in terms of the amount of data can be collected/handled and the speed at which it can be processed leads to the development of methods that could be able to estimate residual unbalance along with bearing dynamic parameters simultaneously [13–16]. These methods generally estimate the residual unbalance accurately but the estimation of bearing dynamic parameters suffers from scattering due to the ill-conditioning of the regression matrix of the estimation equation [1,16]. In the present paper, a general identification algorithm to obtain unknown parameters of the multi-degree of freedom rotor-bearing systems with gyroscopic effects has been presented. Lots of researches were carried out on bearing dynamic parameters estimation [9] by treating the rotor as rigid. Because of high level of uncertainty in the bearing dynamic parameters estimates, still researchers have been concentrating on rigid rotor models, while making dedicated experimental set-up for the bearing dynamic parameters estimation [11]. Residual unbalances gives rise to some amount of error in estimates of bearing dynamic parameters. The present paper addresses the issue of improving bearing dynamic parameters estimation by simultaneously estimating the unknown residual unbalances. Gyroscopic couples are significant when the rotor speed is very high and the diametral mass moment of inertia is large. In laboratory experimental set-up researchers do go for high speed, however, the diametral mass moment of inertia is not significant. Moreover, gyroscopic effects give rise to rotational degree of freedom in equations of motion and accurate measurement of the same is still a challenging problem [17] (unless they could be eliminated by condensation methods). The effect of gyroscopic effects is ignored for the present case illustrations. In the present work, an analysis of the ill-conditioning of the regression matrix, while estimating both the residual unbalance and bearing dynamic parameters of rotor-bearing systems, has been presented. The conditions in which ill-conditioning can occur have been also examined. Theoretical and experimental methods are suggested, as how to improve the condition the regression matrix, in the form of different formulations and types of measurements, respectively. The improvement in the conditioning of the regression matrix has been illustrated by using simulated examples. The effect of noise in the measurement responses has also been tested. 2. Development of a general identification algorithm in rotor-bearing systems The linearised model for any rotating machine (assuming no internal rotor damping) has the form M€q þ ½C þ o1 G1 þ o2 G2 þ _q þ Kq ¼ f;
(1)
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1084
where M is the mass matrix, K is the stiffness matrix, C is the damping matrix and G1, G2, etc. are gyroscopic matrices (skew-symmetric) associated with the different parts of the ‘‘rotational system’’ spinning at the different speeds o1, o2, etc. Transforming Eq. (1) into the frequency domain produces this form DðOÞPðOÞ ¼ ½O2 M þ jO½C þ o1 G1 þ o2 G2 þ þ KpðOÞ ¼ gðOÞ;
(2)
where g(O) is non-zero only for the discrete frequencies O ¼ o1 ; O ¼ o2 ; O ¼ o3 ; etc. (assuming that the only forcing on the system is due to rotating unbalances). Here, D(O) is the dynamic stiffness matrix of the complete rotating machine system at forcing frequency o: The general form of the g(O) may be written gðOi Þ ¼ O2i S i ui ;
(3)
where matrices S1, S2, S3, etc. are ‘‘selection’’ matrices which reflect not only that the unbalance at certain different speeds can exist only at different ‘‘shaft-stations’’ but they also contain the information that for positive Oi ; the forcing in the y-direction lags forcing in the x-direction by 901. Now, partition p into a part which is measured (pm) and a part which is not measured (pu). pm p¼ (4) pu and partition D(O) accordingly thus Dmm ðOÞ Dms ðOÞ DðOÞ ¼ : Dsm ðOÞ Dss ðOÞ
(5)
Note that D itself contains some parts which are known with good confidence and some parts which are quite unknown (the bearing parameters would be in this category). In order to avoid confusion with subscripts, these two separate contributions to D are denoted E and F, respectively, where E is the known part and F is the unknown part. DðOÞ ¼ EðOÞ þ FðOÞ: Matrices E and F are partitioned in the same way as D. Now combine Eqs. (2)–(6). Fmm ðOi Þ Fmu ðOi Þ Pm ðOi Þ Emm ðOi Þ Emu ðOi Þ 2 2 Sim ui : þ ¼ Oi Si u i ¼ Oi Eum ðOi Þ Euu ðOi Þ Fum ðOi Þ Fuu ðOi Þ Pu ðOi Þ Siu
(6)
(7)
Observe that the Si are known. Now, in those cases where Fmu and Fuu are both zero, Eq. (7) can clearly be set out as a linear relationship between the known quantities and the unknowns. The detail of how the unknown quantities are arranged in a vector is irrelevant. The ordering neither helps nor hinders the subsequent processing of this information. After choosing an appropriate ordering, the problem manifests itself as Ax ¼ b;
(8)
where x is the vector of all unknowns (to be determined) and b is a vector of known quantities. In subsequent sections, the issue of what to do to make sure that A is well-conditioned will be discussed in detail.
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1085
3. A rigid rotor mounted on flexible bearings 3.1. Equations of motion and response For a rigid rotor supported on two identical flexible bearings with an independent unbalance excitation unit [9], a schematic diagram is shown in Figs. 1 and 2, equations of motion can be written as jot jOt mx€ þ cxx x_ þ cxy y_ þ kxx x þ kxy y ¼ o2 U Res þ O2 U Tri x e x e ;
(9)
jot jOt þ O2 U Tri my€ þ cyx x_ þ cyy y_ þ kyx x þ kyy y ¼ o2 U Res y e y e ;
(10)
and
where m is the rotor mass per bearing, cij and kij (with i; j ¼ x; y) are the linearised damping and stiffness coefficients of the bearing, U is the unbalance (i.e. multiplication of the unbalance mass and its eccentricity; in general it is a complex quantity and it contains the magnitude and phase information), x and y are the rotor linear displacements in the vertical and horizontal directions, respectively, o is the rotational speed of the rotor, O is the frequency of excitation from the pffiffiffiffiffiffiffi auxiliary unbalance excitation unit, t is the time instant and j ¼ 1: Superscripts: Res and Tri represent the residual and trial unbalance masses, respectively. Taking the solution of Eqs. (9) and (10) in the following form: x ¼ X Res ejot þ X Tri ejOt
(11)
y ¼ Y Res ejot þ Y Tri ejOt ;
(12)
and where the X and Y contain the magnitude and phase information of the respective displacements and in general they are complex quantities (i.e. X ¼ X r þ jX i ; etc.; where subscripts r and i Disk
m
c
Independent unbalance excitation unit
m Bearing k
c
Rigid Foundation Fig. 1. A schematic diagram of a test rig rotor-bearing model.
k
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1086
Fig. 2. An anti-synchronous excitation by an independent unbalance excitation unit.
represent real and imaginary, respectively). Since these equations of motion are linear in nature, on substituting Eqs. (11) and (12) into Eqs. (9) and (10), we get ) ( ) " #( Res U ðkxy þ jocxy Þ ðkxx mo2 þ jocxx Þ X Res x ¼ o2 (13) U Res ðkyx þ jocyx Þ ðkyy mo2 þ jocyy Þ Y Res y and
"
#(
ðkxx mO2 þ jOcxx Þ
ðkxy þ jOcxy Þ
ðkyx þ jOcyx Þ
ðkyy mO2 þ jOcyy Þ
X Tri Y Tri
)
( ¼O
2
U Tri x
)
U Tri y
(14)
For known unbalance information (both the residual and trial unbalances) and rotor-bearing parameters, Eqs. (13) and (14) can be used to obtain the displacement amplitude and phase components. The Eqs. (13) and (14) can be combined as 3 2 ðDxx ðoÞ mo2 Þ Dxy ðoÞ 0 0 7 6 7 6 ðDyy ðoÞ mo2 Þ 0 0 Dyx ðoÞ 7 6 7 6 2 7 6 ðOÞ mO Þ D ðOÞ 0 0 ðD xx xy 5 4 0 0 8 Res 9 8 2 Res 9 o Ux > X > > > > > > > > > > > > > > > > > > Res 2 Res < Y = < o Uy > =
¼ > X Tri > > > > > > O2 U Tri > > > > x > > > > > > > > : Tri ; : 2 Tri > ; Y O Uy
Dyx ðOÞ
Dyy ðOÞ mO2 Þ
ð15Þ
with Dij ðoÞ ¼ kij þ jocij ;
i; j ¼ x; y;
(16)
Dij ðOÞ ¼ kij þ jOcij ;
i; j ¼ x; y:
(17)
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1087
For developing the identification algorithm for the simultaneous estimation of the residual unbalance and bearing dynamic parameters, the frequency domain equations of motion in form of Eq. (15) is used. For the case when the auxiliary excitation unit is not used, however, the trial mass are attached to the rotor itself, then the solution of Eqs. (9) and (10) could be assumed in the following form: x ¼ X ejot
(18)
y ¼ Y ejot ;
(19)
and
where the X and Y contain the magnitude and phase information of displacements and in general they are complex quantities. On substituting Eqs. (18) and (19) into Eqs. (9) and (10), we get ( ) " # Res Tri U þ U X Dxy ðoÞ ðDxx ðoÞ mo2 Þ x x (20) ¼ o2 Tri : U Res Y Dyx ðoÞ ðDyy ðoÞ mo2 Þ y þ Uy For the known unbalance information (i.e. both the residual and trial) and rotor-bearing parameters, Eq. (20) can be used to obtain the displacement amplitude and phase at any spin speed of the rotor. For developing the identification algorithm for simultaneous estimation of the residual unbalance and bearing dynamic parameters, the frequency domain equations of motion in the form of the Eq. (20) will also be used.
3.2. Simultaneous estimation of the residual unbalance and bearing dynamic parameters The Eq. (20) can be rearranged such that all the unknown quantities (i.e. the residual unbalance and bearing dynamic parameters) are in the left-hand side of the equation and the right-hand side contains all the known information, in the following form: 2 ¼ o2 U Tri Dxx ðoÞX þ Dxy ðoÞY o2 U Res x x þ mo X
(21)
2 ¼ o2 U Tri Dyx ðoÞX þ Dyy ðoÞY o2 U Res y y þ mo Y :
(22)
and
Noting that the unbalance force components (both for the residual and trial unbalances) in two orthogonal directions are related as U y ¼ jU x ;
(23)
where the first sign (i.e. negative) is for the case when the direction of rotation of the rotor is in accordance with the positive axis direction conventions and the second sign (positive) is for the case when the direction of rotation of the rotor is in the negative direction of the right-hand axis
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1088
conventions. Noting Eqs. (21)–(23) can be combined as 8 9 Dxx ðoÞ > > > > > > > > > " #> ( ) D ðoÞ > > xy = 2 2 < o2 U Tri X Y 0 0 o x þ mo X Dyx ðoÞ ¼ : 2 Tri 2 > o U þ mo Y 0 0 X Y jo2 > > > y > Dyy ðoÞ > > > > > > : Res > ; Ux
(24)
For the case when the auxiliary unbalance excitation unit is used, it is convenient to express the estimation equation from the Eq. (15) in the following form: A1 ðoÞ q1 ðoÞ fbg ¼ (25) A2 ðOÞ q2 ðOÞ with 2
X Res r 6 X Res 6 ½A1 ðoÞ ¼ 6 i 4 0 0
Y Res r Y Res i
0 0
0 0
oX Res i Res oX r
oY Res i Res oY r
0 0
0 0
o2 0
0
X Res r
Y Res r
0
0
oX Res i
oY Res i
0
0
X Res i
Y Res i
0
oX Res r
oY Res r
o2
0
3 0 o2 7 7 7; o2 5 0
(26) 2
X Tri r
6 X Tri 6 ½A2 ðOÞ ¼ 6 i 4 0 0 fbg ¼ f kxx
kxy
3
Y Tri r
0
0
OX Tri i
OY O i
0
0
0
0
Y Tri i
0
0
OX Tri r
OY Tri r
0
0
0
0
X Tri r
Y Tri r
0
0
OX Tri i
OY Tri i
0
07 7 7; 05
0
X Tri i
Y Tri i
0
0
OX Tri r
OY Tri r
0
0
U Res xr
T U Res xi g ;
kyx
fq1 ðoÞg ¼ o2 f ðmX Res r Þ
kyy
cxx
cxy
cyx
Res ðmX Res i Þ ðmY r Þ
Tri fq2 ðOÞg ¼ O2 f ðmX Tri r þ U xr Þ
cyy
(28)
T ðmY Res i Þg ;
Tri Tri Tri ðmX Tri i þ U xi Þ ðmY r þ U yr Þ
(27)
(29) Tri T ðmY Tri i þ U yi Þ g :
(30)
In this section, two form of the estimation equation have been developed i.e. Eqs. (24) and (25). These are the required form of the estimation equation in which all the unknowns i.e. the residual unbalance and bearing dynamic parameters are stacked in a column vector. These equations can be used to obtain these unknown parameters with the knowledge of the response, the rotational frequency of the rotor/excitation, the rotor mass and the trial unbalance information. From both forms of estimation equations, it can be seen that the number of unknowns are more than the number of equations, which is a case of underdetermined system of linear simultaneous equations. The entire unknown can be obtained with the help of sets of independent force–response measurements such that the number of equations increases at least equal to or more than unknowns. However, the consideration of the condition of the regression matrix to be inverted has important role in obtaining the better estimate of parameters.
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1089
3.3. Condition of the regression matrix The present section deals with various possibilities of measurements of independent responses while keeping the speed of the rotor constant. The following three types of measurements are discussed (a) With sets of trial unbalances. (b) For the above case with rotating the rotor in clockwise and counter clockwise direction, alternatively. (c) Rotor with an independent unbalance excitation unit. For each of the above cases the condition of the regression matrix is also discussed. 3.3.1. Method 1 By taking measurements with three different sets of trial masses (one without trial mass and other with two different trial masses; it is assumed that the rotor has always some amount of residual unbalance) at the same rotational speed of the rotor, six equations can be obtained to solve for five unknowns. Noting Eq. (24), the form of regression equation would be (31)
½AðoÞfbðoÞg ¼ fqðoÞg with 2
X1
6 0 6 6 6 X2 ½AðoÞ ¼ 6 6 0 6 6 4 X3
Y1
0
0
0
X1
Y1
Y2 0
0 X2
0 Y2
Y3
0
0
0
X3
Y3
0
o2
3
jo2 7 7 7 27 o 7 ; jo2 7 7 7 o2 5
(32)
jo2
fbðoÞg ¼ f Dxx ðoÞ
Dxy ðoÞ Dyx ðoÞ
fqðoÞg ¼ o2 f mX 1
mY 1
T Dyy ðoÞ U Res x g ;
1 ðU Tri þ mX 2 Þ x
1 2 ðU Tri þ mY 2 Þ ðU Tri þ mX 3 Þ y x
(33) 2 ðU Tri þ mY 3 Þ gT : y
(34) It should be observed that the vector, {b(o)}, is dependent on the rotor speed. The Eq. (31) is a standard form of the regression equation and it can be used to obtain the residual unbalance and bearing dynamic parameters with the help of the least squares fit from the following equation: fbðoÞg ¼ ð½AðoÞT ½AðoÞÞ1 ½AðoÞT fqðoÞg:
(35)
It should be noted that the solution of the above equation requires inversion of the square matrix ð½AðoÞT ½AðoÞÞ: The matrix condition can be judged by obtaining the condition number (see Appendix B) of the matrix. It has been found that the matrix ð½AðoÞT ½AðoÞÞ is highly ill-conditioned (however, masked by measurement noise while actual responses are processed)
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1090
since the determinant of the matrix is nearly zero i.e. j½AðoÞT ½AðoÞj 0:
(36)
This poses difficulty in identifying the required parameters uniquely. It has been found that even while increasing the number of trial runs to four or even more, does not help in improving the condition of the above matrix. Moreover, as it could be seen from Eq. (32) by scaling of the last column (i.e. by 1=o2 ) the condition of the regression matrix may be improved. However, the condition of Eq. (36) holds after the scaling also. 3.3.2. Method 2 In this method measurements are taken by rotating the rotor first in the clockwise direction and then in the counter clockwise direction. It is assumed that bearing parameters do not change with the change in the sense of rotation (i.e. clockwise or counter clockwise). For symmetric and most commonly used bearings e.g. cylindrical fluid film bearings, rolling element bearings, etc., it may be the case. This approach is applied to the method described above. The measurements with only one trial unbalance is sufficient and the form of equations would be (37)
½AðoÞfbðoÞg ¼ fqðoÞg with 2
o2
X1
Y1
0
0
6 0 6 6 6 X2 ½AðoÞ ¼ 6 6 0 6 6 4 X3
0 Y2
X1 0
Y1 0
0
X2
Y2
Y3
0
0
0
X3
Y3
0 fbðoÞg ¼
Dxx ðoÞ
fqðoÞg ¼ o2 f mX 1
Dxy ðoÞ mY 1
3
7 7 7 7 7; jo2 7 7 7 2 5 o jo2 o2
(38)
jo2
Dyx ðoÞ Dyy ðoÞ mX 2
mY 2
U Res x
T
2 U Tri þ mX 3 x
; 2 U Tri þ mY 3 gT ; y
(39) (40)
where (X1, Y1) are measurements when the rotor having only the residual unbalance is rotated in the positive axis direction, (X2, Y2) are measurements when the rotor with the residual unbalance is rotated in the negative axis direction and (X3, Y3) are measurements when a trial mass is also added on the rotor and it is rotated in the positive axis direction. On comparing the matrix [A(o)] of the Eq. (38) with that of the Eq. (32), it can be seen that the matrix element (4, 5) gets its sign changed because of the negative axis direction of rotation of the rotor (refer the Eq. (23)). The condition number of the matrix ð½AðoÞT ½AðoÞÞ improves drastically. The fourth measurement can also be incorporated with the trial unbalance when the rotor is rotated in the opposite direction. More measurements can also be incorporated with different trial unbalances. 3.3.3. Method 3 For the case when the rotor cannot be rotated in either direction (i.e. clock wise and anti-clock wise) or the bearing parameters change with the direction of rotation of the rotor, an independent
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1091
unbalance excitation unit is generally used. For such case it is assumed that rotor is always rotated in the positive axis direction. The condition of the regression matrix for identification of the bearing dynamic parameters improves when such arrangement is used especially for the case when two excitation frequencies are anti-synchronous [12]. The three-independent measurements could be obtained, first corresponding to the rotor rotational frequency and two measurements corresponding to two excitation frequencies (O1 and O2 i.e. the negative sign indicates the sense of rotation is opposite to the positive axis direction). It is assumed that for all the three measurements the rotor speed remains constant. From the Eq. (25) the regression equation takes the following form: (41)
½AðoÞfbg ¼ fqðoÞg with ½AðoÞ ¼ ½ A1 ðoÞ
A2 ðO1 Þ A2 ðO2 Þ T ;
fqðoÞg ¼ f q1 ðoÞ
q2 ðO1 Þ
(42)
q2 ðO2 Þ gT ;
(43)
A2 ðO2 Þ ¼ A2 ðO2 Þ;
(44)
where matrices [A1(o)] and [A2(O)] are defined in Eqs. (26) and (27), respectively, and vectors {q1(o)} and {q2(O)} are defined by Eqs. (29) and (30), respectively. The vector {b}, as defined in Eq. (28), contains residual unbalances in the last two rows. Noting Eq. (23), accordingly matrices [A1(o)] and [A2(O)] are affected only in the last two columns when the direction of residual unbalance changes. Since for the present case the direction of the residual unbalance does not change, Eq. (44) holds good. Moreover, it can be observed that Eq. (27) has last two columns as zero entries. With the above arrangement it can be seen that the regression matrix remains illconditioned since on rotating the trial mass in two different direction does not affect the regression matrix as the trial unbalance are contained in the vector {q} (i.e. Eq. (30)).
4. Results and discussions A rotor-bearing system as shown in Fig. 1 is considered for the numerical simulation. The rotor mass is 0.447 Kg. Both the bearings are assumed to be identical and its dynamic properties are given in Table 1 [12]. Residual and trial unbalance magnitudes and angular positions are tabulated in Table 2. The rotor response is generated at different unbalance configurations at a particular rotor speed by using Eq. (20). The simulated response is fed into the identification algorithm of Method 1 (i.e. Eq. (31)) to estimate the residual unbalance and bearing dynamic parameters. The estimated Table 1 Bearing dynamic parameters Bearing parameters
kxx (N/m)
kxy (N/m)
kyx (N/m)
kyy (N/m)
cxx (Ns/m)
cxy (Ns/m)
cyx (Ns/m)
cyy (Ns/m)
Values
16,788.0
1000.0
1000.0
18,592.0
6.0
0.0
0.0
6.0
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1092
Table 2 Residual and trial unbalances S.N.
Type of unbalance
Magnitude (kg)
Radius (m)
Angular location (deg.)
1 2 3
Residual unbalance Trial unbalance 1 Trial unbalance 2
0.001 0.002 0.003
0.03 0.03 0.03
0.0 90 60
Table 3 Comparison of the identified unbalance and bearing dynamic parameters Parameters chosen for simulation of responses
kxx (N/m) kxy (N/m) kyx (N/m) kyy (N/m) cxx (Ns/m) cxy (Ns/m) cyx (Ns/m) cyy (Ns/m) Residual unbalance (kgm) Residual unbalance phase (deg.) Condition number of regression matrix without scaling Condition number of the regression matrix after column scaling
Parameters estimated by the identification for the unbalance and bearings dynamic parameters Method 1 (Section 4.1) Three runs Three runs without with 1% noise noise
Three runs with 5% noise
Method 2 (Section 4.2) Three runs Three runs without with 1% noise noise
Three runs with 5% noise
16788.0 1000.0 1000.0 18,592.0 6.0 0.0 0.0 6.0 3.0E–5
16,679.0 0.0 699.8 0.0 8.1 0.0 262.2 0.0 3.0E–5
16,693.0 0.0 700.6 0.0 8.1 0.0 262.5 0.0 3.0E–5
16,753.0 0.0 703.9 0.0 8.1 0.0 263.5 0.0 3.0E–5
16,788.0 1000.0 1000.0 18,592.0 6.0 1.9E–11 6.8E–12 6.0 3.0E–5
16,803.0 999.4 1001.0 18,582.0 6.01 3.8E–11 3.9E–12 5.996 3.0E–5
16,863.0 997.0 1005.0 18,542.0 6.03 2.42E–11 1.64E–11 5.982 3.0E–5
0.0
1.57E–13
1.04E–13
1.38E–13
2.9E–14
3.85E–14
1.47E–13
—
1.62E9
1.63E9
1.66E9
0.667E9
0.67E9
0.68E9
—
41.38
41.58
42.42
17.03
17.11
17.46
parameters are tabulated in Table 3 along with the condition number of the regression matrix and results suggest that it suffers from ill-conditioning. However, it should be noted that the residual unbalance could be estimated accurately and this has been reported in [16]. To improve (i.e to reduce) the condition number of the regression matrix (i.e. Eq. (32)) the last column is scaled by a factor of o2 1E4, so that the last column entries are of the same order as that of the rest. The corresponding condition number of the matrix is also tabulated in Table 3. Drastic improvement in the condition number occurred due to the scaling, however, for the present case no improvement in the estimated has been found. The similar identification exercise is performed by contaminating simulated responses by the random noise of 1% or 5%. The corresponding
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1093
estimated parameters and the condition number of the regression matrix are also tabulated in Table 3 and most of the bearing dynamic parameters have ill-conditioning effects. By rotating the rotor alternatively in the clockwise or counter clockwise direction the rotor response is generated at different unbalance configurations at a particular rotor speed by using Eq. (20). The simulated response is fed into the identification algorithm of Method 2 (i.e. Eq. (37)) to estimate the residual unbalance and bearing dynamic parameters. The estimated parameters are tabulated in Table 3 and results suggest that the regression equation is now well conditioned. The corresponding condition number of regression matrix (i.e. Eq. (38)) is also tabulated in Table 3 and as compared to Method 1 the condition number of the regression matrix is much lower. For the present case, however, it is observed that there is no improvement in the estimates of the parameters. The same exercise is performed by contaminating simulated responses by the random noise of 1% or 5%. The corresponding estimated parameters and condition numbers are also tabulated in Table 3 and they suggest the robustness of the present algorithm against measurement noise. With the auxiliary unbalance unit (i.e. Method 3) the identification algorithm suffers from illconditioning similar to the Method 1 and result trends are identical to the Method 1.
5. Conclusions The present method is aimed at improving the existing method of estimating the bearing dynamic parameters by simultaneously estimating the unknown residual imbalance, which might cause error in the estimates of bearing parameters, inadvertently. For a linear multi-degreeof-freedom rotor-bearing system a general identification algorithm of systems’ unknowns has been presented. Then, an identification algorithm for simultaneous estimation of the residual unbalance and bearing dynamic parameters has been developed. The problem is found to be ill-posed. By taking the response measurements, while rotating the rotor alternatively in clockwise or counter clockwise gives identification algorithm as well-conditioned. The algorithm is tested against the measurement noise and found to be robust. However, with an auxiliary unbalance unit (in which excitation frequency can be chosen independent of rotor speed both in magnitude and sense) does not help in improving the estimates of bearing dynamic parameters while estimating residual unbalance simultaneously as against while estimating the bearing dynamic parameters alone [12]. The present method proposes essentially three methods by which to improve the condition of matrix A. The first involves multiple trial masses, which was noted in the introduction as being uneconomical and in Section 3.3.1 as ineffective. The second involves running the rotor in opposite directions; however, this is rarely possible in real rotors. Moreover, it is possible that the bearing coefficients also change. It will be more appropriate to call bearing stiffness and damping coefficients in this paper as bearing support parameters. The third measure involves an auxiliary unbalance unit. This is useful educational and research experimental set-ups. The estimation of bearing dynamic parameters is still thought to be relatively reliable only in laboratory test rigs [11]. However, it is unlikely to fit one to a working rotor. The development given in Section 2 potentially takes this idea of the present paper to the stage where it might be applied to a multiple-rotor machine (like a gas-turbine engine).
ARTICLE IN PRESS 1094
R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
Acknowledgement The author acknowledges anonymous referees for their useful contribution in improving the presentation of the paper especially Section 2. Appendix A. Nomenclature c the bearing damping coefficient C the damping matrix D k+jc D the dynamic stiffness matrix G theffiffiffiffiffiffigyroscopic matrix p ffi j 1 k the bearing stiffness coefficient K the stiffness matrix m the rotor mass per bearing M the mass matrix t the time instance U the unbalance (mass eccentricity) x, y the orthogonal axis system X, Y the complex displacements o the rotor rotational frequency O the excitation frequency from an independent unbalance excitation unit Subscripts r the real part i the imaginary part Superscripts Res the residual unbalance Tri the trial unbalance
Appendix B. condition number The condition number of a matrix (with respect to inversion) measures the sensitivity of the solution of a system of linear equations to errors in the data. It gives an indication of the accuracy of the results from matrix inversion and the linear equation solution. Values of the condition number near 1 indicate a well-conditioned matrix. The 2-norm condition number, which is used in the present case, is defined the ratio of the largest singular value of regression matrix to the smallest [18].
References [1] W. Kellenburger, Should a flexible rotor be balanced in N or N+2 planes?, Transactions of the American Society of Mechanical Engineers: Journal of Engineering for Industry 94 (1972) 560–584.
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1095
[2] J. Drechsler, Processing surplus information in computer aided balancing of large flexible rotors, Institution of Mechanical Engineers Conference on Vibrations in Rotating Machinery, Cambridge, UK, 1980, pp. 65–70. [3] P. Gnilka, Modal balancing of flexible rotors without test runs: an experimental investigation, Journal of Vibration 90 (1983) 157–172. [4] J.M. Krodkiewski, J. Ding, N. Zhang, Identification of unbalance change using a non-linear mathematical model for rotor bearing systems, Journal of Vibration 169 (1994) 685–698. [5] M.S. Darlow, Balancing of High-Speed Machinery, Springer, Berlin, 1989. [6] S. Edwards, A.W. Lees, M.I. Friswell, Experimental identification of excitation and support parameters of a flexible rotor-bearings-foundation system from a single run-down, Journal of Sound and Vibration 232 (5) (2000) 963–992. [7] J.W. Lund, Review of analytical methods in rotor-bearing dynamics, Tribology International 13 (5) (1980) 233–236. [8] P. Ramsden, Review of published data and their application to the design of large bearings for steam turbines, Proceedings of Conference in Lubrication and Wear: Fundamentals and Application to Design, IMechE 182 (3A) (1967–68) 75–81. [9] M.J. Goodwin, Experimental techniques for bearing impedance measurement, Transactions of the American Society of Mechanical Engineers: Journal of Engineering for Industry 113 (1991) 335–342. [10] E.E. Swanson, R.G. Kirk, Survey of Experimental Data for Fixed Geometry Hydrodynamic Journal Bearings, ASME Journal of Tribology 119 (4) (1997) 704–710. [11] R. Tiwari, A.W. Lees, M.I. Friswell, Identification of dynamic bearing parameters: a review, Shock & Vibration Digest 36 (2) (2004) 99–124. [12] R. Tiwari, A.W. Lees, M.I. Friswell, Identification of speed-dependent bearing parameters, Journal of Sound and Vibration 254 (5) (2002) 967–986. [13] J.H. Chen, A.C. Lee, Estimation of linearised dynamic characteristics of bearings using synchronous response, International Journal of Mechanical Science 37 (2) (1995) 197–219. [14] A.C. Lee, Y.P. Shih, Identification of the unbalance distribution and dynamic characteristics of bearings in flexible rotors, Proceedings of IMechE, Part C: Journal of Mechanical Engineering Science 210 (5) (1996) 409–432. [15] R. Tiwari, N.S. Vyas, Parameter estimation in imbalanced non-linear rotor-bearing systems from random response, Journal of Sound and Vibration 208 (1) (1997) 1–14. [16] J.K. Sinha, M.I. Friswell, A.W. Lees, The identification of the unbalance and the foundation model of a rotating machine from a single run-down, Mechanical System and Signal Processing 16 (2–3) (2002) 255–271. [17] N. Dharmaraju, R. Tiwari, S. Talukdar, Development of a novel hybrid reduction scheme for identification of an open crack model in a beam, Mechanical System and Signal Processing 19 (3) (2005) 633–657. [18] G.H. Golub, C.F. Van Loan, Matrix Computations, Third ed, The Johns Hopkins University Press, London, 1996.
Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 19 (2005) 1082–1095 www.elsevier.com/locate/jnlabr/ymssp
Conditioning of regression matrices for simultaneous estimation of the residual unbalance and bearing dynamic parameters Rajiv Tiwari Department of Mechanical Engineering, Indian Institute of Technology, Guwahati North Guwahati, Guwahati 781 039, India Received 18 November 2003; received in revised form 16 August 2004; accepted 27 September 2004 Available online 21 November 2004
Abstract A general identification algorithm has been presented for multi-degree-of-freedom rotor-bearing systems. A specific identification algorithm for simultaneous estimation of the residual unbalance and bearing dynamic parameters in rotor-bearing systems has been developed. The identification algorithm is found to be highly ill-conditioned for the conventional measurement techniques. The measurement of responses, while rotating the rotor alternatively in clockwise or counter clockwise direction leads to a well-conditioned identification algorithm. Numerical examples have been used to illustrate the developed algorithms. The algorithm is tested against the measurement noise and found to be robust. r 2004 Elsevier Ltd. All rights reserved.
1. Introduction The residual unbalance estimation in rotor-bearing system is an age-old problem. From the state of the art of the unbalance estimation, the unbalance can be obtained with fairly good accuracy [1–5]. Now the trend in the unbalance estimation is to reduce the number of test runs required, especially for the application of large turbogenerators where the downtime is very expensive [6]. Tel.: +91 0361 2582 667; fax: +91 0361 2690 762.
E-mail address: [email protected] (R. Tiwari). 0888-3270/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2004.09.005
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1083
Bearing dynamic parameters play an important role in the prediction of the dynamic behaviour of rotor-bearing systems. Due to the difficulty in obtaining actual test conditions the well-developed theoretical analysis to obtain bearing dynamic parameters by using the Reynolds equation fails to give satisfactory results [7]. This forces designers to look for the experimentally obtained dynamic bearing parameters. Exhaustive work has been done on the experimental identification of the eight linearised bearing dynamic parameters [8–11]. Most of the work is based on the dynamic force, which is given to the system, and corresponding response is measured. At least two-independent force–response relationships are used to obtain bearing dynamic parameters. Researchers have reported that often the regression matrix of the estimation equation becomes ill-conditioned and this leads to scattering of the estimated bearing dynamic parameters [12]. Advances in the sensor technology and the increase in computing power in terms of the amount of data can be collected/handled and the speed at which it can be processed leads to the development of methods that could be able to estimate residual unbalance along with bearing dynamic parameters simultaneously [13–16]. These methods generally estimate the residual unbalance accurately but the estimation of bearing dynamic parameters suffers from scattering due to the ill-conditioning of the regression matrix of the estimation equation [1,16]. In the present paper, a general identification algorithm to obtain unknown parameters of the multi-degree of freedom rotor-bearing systems with gyroscopic effects has been presented. Lots of researches were carried out on bearing dynamic parameters estimation [9] by treating the rotor as rigid. Because of high level of uncertainty in the bearing dynamic parameters estimates, still researchers have been concentrating on rigid rotor models, while making dedicated experimental set-up for the bearing dynamic parameters estimation [11]. Residual unbalances gives rise to some amount of error in estimates of bearing dynamic parameters. The present paper addresses the issue of improving bearing dynamic parameters estimation by simultaneously estimating the unknown residual unbalances. Gyroscopic couples are significant when the rotor speed is very high and the diametral mass moment of inertia is large. In laboratory experimental set-up researchers do go for high speed, however, the diametral mass moment of inertia is not significant. Moreover, gyroscopic effects give rise to rotational degree of freedom in equations of motion and accurate measurement of the same is still a challenging problem [17] (unless they could be eliminated by condensation methods). The effect of gyroscopic effects is ignored for the present case illustrations. In the present work, an analysis of the ill-conditioning of the regression matrix, while estimating both the residual unbalance and bearing dynamic parameters of rotor-bearing systems, has been presented. The conditions in which ill-conditioning can occur have been also examined. Theoretical and experimental methods are suggested, as how to improve the condition the regression matrix, in the form of different formulations and types of measurements, respectively. The improvement in the conditioning of the regression matrix has been illustrated by using simulated examples. The effect of noise in the measurement responses has also been tested. 2. Development of a general identification algorithm in rotor-bearing systems The linearised model for any rotating machine (assuming no internal rotor damping) has the form M€q þ ½C þ o1 G1 þ o2 G2 þ _q þ Kq ¼ f;
(1)
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1084
where M is the mass matrix, K is the stiffness matrix, C is the damping matrix and G1, G2, etc. are gyroscopic matrices (skew-symmetric) associated with the different parts of the ‘‘rotational system’’ spinning at the different speeds o1, o2, etc. Transforming Eq. (1) into the frequency domain produces this form DðOÞPðOÞ ¼ ½O2 M þ jO½C þ o1 G1 þ o2 G2 þ þ KpðOÞ ¼ gðOÞ;
(2)
where g(O) is non-zero only for the discrete frequencies O ¼ o1 ; O ¼ o2 ; O ¼ o3 ; etc. (assuming that the only forcing on the system is due to rotating unbalances). Here, D(O) is the dynamic stiffness matrix of the complete rotating machine system at forcing frequency o: The general form of the g(O) may be written gðOi Þ ¼ O2i S i ui ;
(3)
where matrices S1, S2, S3, etc. are ‘‘selection’’ matrices which reflect not only that the unbalance at certain different speeds can exist only at different ‘‘shaft-stations’’ but they also contain the information that for positive Oi ; the forcing in the y-direction lags forcing in the x-direction by 901. Now, partition p into a part which is measured (pm) and a part which is not measured (pu). pm p¼ (4) pu and partition D(O) accordingly thus Dmm ðOÞ Dms ðOÞ DðOÞ ¼ : Dsm ðOÞ Dss ðOÞ
(5)
Note that D itself contains some parts which are known with good confidence and some parts which are quite unknown (the bearing parameters would be in this category). In order to avoid confusion with subscripts, these two separate contributions to D are denoted E and F, respectively, where E is the known part and F is the unknown part. DðOÞ ¼ EðOÞ þ FðOÞ: Matrices E and F are partitioned in the same way as D. Now combine Eqs. (2)–(6). Fmm ðOi Þ Fmu ðOi Þ Pm ðOi Þ Emm ðOi Þ Emu ðOi Þ 2 2 Sim ui : þ ¼ Oi Si u i ¼ Oi Eum ðOi Þ Euu ðOi Þ Fum ðOi Þ Fuu ðOi Þ Pu ðOi Þ Siu
(6)
(7)
Observe that the Si are known. Now, in those cases where Fmu and Fuu are both zero, Eq. (7) can clearly be set out as a linear relationship between the known quantities and the unknowns. The detail of how the unknown quantities are arranged in a vector is irrelevant. The ordering neither helps nor hinders the subsequent processing of this information. After choosing an appropriate ordering, the problem manifests itself as Ax ¼ b;
(8)
where x is the vector of all unknowns (to be determined) and b is a vector of known quantities. In subsequent sections, the issue of what to do to make sure that A is well-conditioned will be discussed in detail.
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1085
3. A rigid rotor mounted on flexible bearings 3.1. Equations of motion and response For a rigid rotor supported on two identical flexible bearings with an independent unbalance excitation unit [9], a schematic diagram is shown in Figs. 1 and 2, equations of motion can be written as jot jOt mx€ þ cxx x_ þ cxy y_ þ kxx x þ kxy y ¼ o2 U Res þ O2 U Tri x e x e ;
(9)
jot jOt þ O2 U Tri my€ þ cyx x_ þ cyy y_ þ kyx x þ kyy y ¼ o2 U Res y e y e ;
(10)
and
where m is the rotor mass per bearing, cij and kij (with i; j ¼ x; y) are the linearised damping and stiffness coefficients of the bearing, U is the unbalance (i.e. multiplication of the unbalance mass and its eccentricity; in general it is a complex quantity and it contains the magnitude and phase information), x and y are the rotor linear displacements in the vertical and horizontal directions, respectively, o is the rotational speed of the rotor, O is the frequency of excitation from the pffiffiffiffiffiffiffi auxiliary unbalance excitation unit, t is the time instant and j ¼ 1: Superscripts: Res and Tri represent the residual and trial unbalance masses, respectively. Taking the solution of Eqs. (9) and (10) in the following form: x ¼ X Res ejot þ X Tri ejOt
(11)
y ¼ Y Res ejot þ Y Tri ejOt ;
(12)
and where the X and Y contain the magnitude and phase information of the respective displacements and in general they are complex quantities (i.e. X ¼ X r þ jX i ; etc.; where subscripts r and i Disk
m
c
Independent unbalance excitation unit
m Bearing k
c
Rigid Foundation Fig. 1. A schematic diagram of a test rig rotor-bearing model.
k
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1086
Fig. 2. An anti-synchronous excitation by an independent unbalance excitation unit.
represent real and imaginary, respectively). Since these equations of motion are linear in nature, on substituting Eqs. (11) and (12) into Eqs. (9) and (10), we get ) ( ) " #( Res U ðkxy þ jocxy Þ ðkxx mo2 þ jocxx Þ X Res x ¼ o2 (13) U Res ðkyx þ jocyx Þ ðkyy mo2 þ jocyy Þ Y Res y and
"
#(
ðkxx mO2 þ jOcxx Þ
ðkxy þ jOcxy Þ
ðkyx þ jOcyx Þ
ðkyy mO2 þ jOcyy Þ
X Tri Y Tri
)
( ¼O
2
U Tri x
)
U Tri y
(14)
For known unbalance information (both the residual and trial unbalances) and rotor-bearing parameters, Eqs. (13) and (14) can be used to obtain the displacement amplitude and phase components. The Eqs. (13) and (14) can be combined as 3 2 ðDxx ðoÞ mo2 Þ Dxy ðoÞ 0 0 7 6 7 6 ðDyy ðoÞ mo2 Þ 0 0 Dyx ðoÞ 7 6 7 6 2 7 6 ðOÞ mO Þ D ðOÞ 0 0 ðD xx xy 5 4 0 0 8 Res 9 8 2 Res 9 o Ux > X > > > > > > > > > > > > > > > > > > Res 2 Res < Y = < o Uy > =
¼ > X Tri > > > > > > O2 U Tri > > > > x > > > > > > > > : Tri ; : 2 Tri > ; Y O Uy
Dyx ðOÞ
Dyy ðOÞ mO2 Þ
ð15Þ
with Dij ðoÞ ¼ kij þ jocij ;
i; j ¼ x; y;
(16)
Dij ðOÞ ¼ kij þ jOcij ;
i; j ¼ x; y:
(17)
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1087
For developing the identification algorithm for the simultaneous estimation of the residual unbalance and bearing dynamic parameters, the frequency domain equations of motion in form of Eq. (15) is used. For the case when the auxiliary excitation unit is not used, however, the trial mass are attached to the rotor itself, then the solution of Eqs. (9) and (10) could be assumed in the following form: x ¼ X ejot
(18)
y ¼ Y ejot ;
(19)
and
where the X and Y contain the magnitude and phase information of displacements and in general they are complex quantities. On substituting Eqs. (18) and (19) into Eqs. (9) and (10), we get ( ) " # Res Tri U þ U X Dxy ðoÞ ðDxx ðoÞ mo2 Þ x x (20) ¼ o2 Tri : U Res Y Dyx ðoÞ ðDyy ðoÞ mo2 Þ y þ Uy For the known unbalance information (i.e. both the residual and trial) and rotor-bearing parameters, Eq. (20) can be used to obtain the displacement amplitude and phase at any spin speed of the rotor. For developing the identification algorithm for simultaneous estimation of the residual unbalance and bearing dynamic parameters, the frequency domain equations of motion in the form of the Eq. (20) will also be used.
3.2. Simultaneous estimation of the residual unbalance and bearing dynamic parameters The Eq. (20) can be rearranged such that all the unknown quantities (i.e. the residual unbalance and bearing dynamic parameters) are in the left-hand side of the equation and the right-hand side contains all the known information, in the following form: 2 ¼ o2 U Tri Dxx ðoÞX þ Dxy ðoÞY o2 U Res x x þ mo X
(21)
2 ¼ o2 U Tri Dyx ðoÞX þ Dyy ðoÞY o2 U Res y y þ mo Y :
(22)
and
Noting that the unbalance force components (both for the residual and trial unbalances) in two orthogonal directions are related as U y ¼ jU x ;
(23)
where the first sign (i.e. negative) is for the case when the direction of rotation of the rotor is in accordance with the positive axis direction conventions and the second sign (positive) is for the case when the direction of rotation of the rotor is in the negative direction of the right-hand axis
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1088
conventions. Noting Eqs. (21)–(23) can be combined as 8 9 Dxx ðoÞ > > > > > > > > > " #> ( ) D ðoÞ > > xy = 2 2 < o2 U Tri X Y 0 0 o x þ mo X Dyx ðoÞ ¼ : 2 Tri 2 > o U þ mo Y 0 0 X Y jo2 > > > y > Dyy ðoÞ > > > > > > : Res > ; Ux
(24)
For the case when the auxiliary unbalance excitation unit is used, it is convenient to express the estimation equation from the Eq. (15) in the following form: A1 ðoÞ q1 ðoÞ fbg ¼ (25) A2 ðOÞ q2 ðOÞ with 2
X Res r 6 X Res 6 ½A1 ðoÞ ¼ 6 i 4 0 0
Y Res r Y Res i
0 0
0 0
oX Res i Res oX r
oY Res i Res oY r
0 0
0 0
o2 0
0
X Res r
Y Res r
0
0
oX Res i
oY Res i
0
0
X Res i
Y Res i
0
oX Res r
oY Res r
o2
0
3 0 o2 7 7 7; o2 5 0
(26) 2
X Tri r
6 X Tri 6 ½A2 ðOÞ ¼ 6 i 4 0 0 fbg ¼ f kxx
kxy
3
Y Tri r
0
0
OX Tri i
OY O i
0
0
0
0
Y Tri i
0
0
OX Tri r
OY Tri r
0
0
0
0
X Tri r
Y Tri r
0
0
OX Tri i
OY Tri i
0
07 7 7; 05
0
X Tri i
Y Tri i
0
0
OX Tri r
OY Tri r
0
0
U Res xr
T U Res xi g ;
kyx
fq1 ðoÞg ¼ o2 f ðmX Res r Þ
kyy
cxx
cxy
cyx
Res ðmX Res i Þ ðmY r Þ
Tri fq2 ðOÞg ¼ O2 f ðmX Tri r þ U xr Þ
cyy
(28)
T ðmY Res i Þg ;
Tri Tri Tri ðmX Tri i þ U xi Þ ðmY r þ U yr Þ
(27)
(29) Tri T ðmY Tri i þ U yi Þ g :
(30)
In this section, two form of the estimation equation have been developed i.e. Eqs. (24) and (25). These are the required form of the estimation equation in which all the unknowns i.e. the residual unbalance and bearing dynamic parameters are stacked in a column vector. These equations can be used to obtain these unknown parameters with the knowledge of the response, the rotational frequency of the rotor/excitation, the rotor mass and the trial unbalance information. From both forms of estimation equations, it can be seen that the number of unknowns are more than the number of equations, which is a case of underdetermined system of linear simultaneous equations. The entire unknown can be obtained with the help of sets of independent force–response measurements such that the number of equations increases at least equal to or more than unknowns. However, the consideration of the condition of the regression matrix to be inverted has important role in obtaining the better estimate of parameters.
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1089
3.3. Condition of the regression matrix The present section deals with various possibilities of measurements of independent responses while keeping the speed of the rotor constant. The following three types of measurements are discussed (a) With sets of trial unbalances. (b) For the above case with rotating the rotor in clockwise and counter clockwise direction, alternatively. (c) Rotor with an independent unbalance excitation unit. For each of the above cases the condition of the regression matrix is also discussed. 3.3.1. Method 1 By taking measurements with three different sets of trial masses (one without trial mass and other with two different trial masses; it is assumed that the rotor has always some amount of residual unbalance) at the same rotational speed of the rotor, six equations can be obtained to solve for five unknowns. Noting Eq. (24), the form of regression equation would be (31)
½AðoÞfbðoÞg ¼ fqðoÞg with 2
X1
6 0 6 6 6 X2 ½AðoÞ ¼ 6 6 0 6 6 4 X3
Y1
0
0
0
X1
Y1
Y2 0
0 X2
0 Y2
Y3
0
0
0
X3
Y3
0
o2
3
jo2 7 7 7 27 o 7 ; jo2 7 7 7 o2 5
(32)
jo2
fbðoÞg ¼ f Dxx ðoÞ
Dxy ðoÞ Dyx ðoÞ
fqðoÞg ¼ o2 f mX 1
mY 1
T Dyy ðoÞ U Res x g ;
1 ðU Tri þ mX 2 Þ x
1 2 ðU Tri þ mY 2 Þ ðU Tri þ mX 3 Þ y x
(33) 2 ðU Tri þ mY 3 Þ gT : y
(34) It should be observed that the vector, {b(o)}, is dependent on the rotor speed. The Eq. (31) is a standard form of the regression equation and it can be used to obtain the residual unbalance and bearing dynamic parameters with the help of the least squares fit from the following equation: fbðoÞg ¼ ð½AðoÞT ½AðoÞÞ1 ½AðoÞT fqðoÞg:
(35)
It should be noted that the solution of the above equation requires inversion of the square matrix ð½AðoÞT ½AðoÞÞ: The matrix condition can be judged by obtaining the condition number (see Appendix B) of the matrix. It has been found that the matrix ð½AðoÞT ½AðoÞÞ is highly ill-conditioned (however, masked by measurement noise while actual responses are processed)
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1090
since the determinant of the matrix is nearly zero i.e. j½AðoÞT ½AðoÞj 0:
(36)
This poses difficulty in identifying the required parameters uniquely. It has been found that even while increasing the number of trial runs to four or even more, does not help in improving the condition of the above matrix. Moreover, as it could be seen from Eq. (32) by scaling of the last column (i.e. by 1=o2 ) the condition of the regression matrix may be improved. However, the condition of Eq. (36) holds after the scaling also. 3.3.2. Method 2 In this method measurements are taken by rotating the rotor first in the clockwise direction and then in the counter clockwise direction. It is assumed that bearing parameters do not change with the change in the sense of rotation (i.e. clockwise or counter clockwise). For symmetric and most commonly used bearings e.g. cylindrical fluid film bearings, rolling element bearings, etc., it may be the case. This approach is applied to the method described above. The measurements with only one trial unbalance is sufficient and the form of equations would be (37)
½AðoÞfbðoÞg ¼ fqðoÞg with 2
o2
X1
Y1
0
0
6 0 6 6 6 X2 ½AðoÞ ¼ 6 6 0 6 6 4 X3
0 Y2
X1 0
Y1 0
0
X2
Y2
Y3
0
0
0
X3
Y3
0 fbðoÞg ¼
Dxx ðoÞ
fqðoÞg ¼ o2 f mX 1
Dxy ðoÞ mY 1
3
7 7 7 7 7; jo2 7 7 7 2 5 o jo2 o2
(38)
jo2
Dyx ðoÞ Dyy ðoÞ mX 2
mY 2
U Res x
T
2 U Tri þ mX 3 x
; 2 U Tri þ mY 3 gT ; y
(39) (40)
where (X1, Y1) are measurements when the rotor having only the residual unbalance is rotated in the positive axis direction, (X2, Y2) are measurements when the rotor with the residual unbalance is rotated in the negative axis direction and (X3, Y3) are measurements when a trial mass is also added on the rotor and it is rotated in the positive axis direction. On comparing the matrix [A(o)] of the Eq. (38) with that of the Eq. (32), it can be seen that the matrix element (4, 5) gets its sign changed because of the negative axis direction of rotation of the rotor (refer the Eq. (23)). The condition number of the matrix ð½AðoÞT ½AðoÞÞ improves drastically. The fourth measurement can also be incorporated with the trial unbalance when the rotor is rotated in the opposite direction. More measurements can also be incorporated with different trial unbalances. 3.3.3. Method 3 For the case when the rotor cannot be rotated in either direction (i.e. clock wise and anti-clock wise) or the bearing parameters change with the direction of rotation of the rotor, an independent
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1091
unbalance excitation unit is generally used. For such case it is assumed that rotor is always rotated in the positive axis direction. The condition of the regression matrix for identification of the bearing dynamic parameters improves when such arrangement is used especially for the case when two excitation frequencies are anti-synchronous [12]. The three-independent measurements could be obtained, first corresponding to the rotor rotational frequency and two measurements corresponding to two excitation frequencies (O1 and O2 i.e. the negative sign indicates the sense of rotation is opposite to the positive axis direction). It is assumed that for all the three measurements the rotor speed remains constant. From the Eq. (25) the regression equation takes the following form: (41)
½AðoÞfbg ¼ fqðoÞg with ½AðoÞ ¼ ½ A1 ðoÞ
A2 ðO1 Þ A2 ðO2 Þ T ;
fqðoÞg ¼ f q1 ðoÞ
q2 ðO1 Þ
(42)
q2 ðO2 Þ gT ;
(43)
A2 ðO2 Þ ¼ A2 ðO2 Þ;
(44)
where matrices [A1(o)] and [A2(O)] are defined in Eqs. (26) and (27), respectively, and vectors {q1(o)} and {q2(O)} are defined by Eqs. (29) and (30), respectively. The vector {b}, as defined in Eq. (28), contains residual unbalances in the last two rows. Noting Eq. (23), accordingly matrices [A1(o)] and [A2(O)] are affected only in the last two columns when the direction of residual unbalance changes. Since for the present case the direction of the residual unbalance does not change, Eq. (44) holds good. Moreover, it can be observed that Eq. (27) has last two columns as zero entries. With the above arrangement it can be seen that the regression matrix remains illconditioned since on rotating the trial mass in two different direction does not affect the regression matrix as the trial unbalance are contained in the vector {q} (i.e. Eq. (30)).
4. Results and discussions A rotor-bearing system as shown in Fig. 1 is considered for the numerical simulation. The rotor mass is 0.447 Kg. Both the bearings are assumed to be identical and its dynamic properties are given in Table 1 [12]. Residual and trial unbalance magnitudes and angular positions are tabulated in Table 2. The rotor response is generated at different unbalance configurations at a particular rotor speed by using Eq. (20). The simulated response is fed into the identification algorithm of Method 1 (i.e. Eq. (31)) to estimate the residual unbalance and bearing dynamic parameters. The estimated Table 1 Bearing dynamic parameters Bearing parameters
kxx (N/m)
kxy (N/m)
kyx (N/m)
kyy (N/m)
cxx (Ns/m)
cxy (Ns/m)
cyx (Ns/m)
cyy (Ns/m)
Values
16,788.0
1000.0
1000.0
18,592.0
6.0
0.0
0.0
6.0
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1092
Table 2 Residual and trial unbalances S.N.
Type of unbalance
Magnitude (kg)
Radius (m)
Angular location (deg.)
1 2 3
Residual unbalance Trial unbalance 1 Trial unbalance 2
0.001 0.002 0.003
0.03 0.03 0.03
0.0 90 60
Table 3 Comparison of the identified unbalance and bearing dynamic parameters Parameters chosen for simulation of responses
kxx (N/m) kxy (N/m) kyx (N/m) kyy (N/m) cxx (Ns/m) cxy (Ns/m) cyx (Ns/m) cyy (Ns/m) Residual unbalance (kgm) Residual unbalance phase (deg.) Condition number of regression matrix without scaling Condition number of the regression matrix after column scaling
Parameters estimated by the identification for the unbalance and bearings dynamic parameters Method 1 (Section 4.1) Three runs Three runs without with 1% noise noise
Three runs with 5% noise
Method 2 (Section 4.2) Three runs Three runs without with 1% noise noise
Three runs with 5% noise
16788.0 1000.0 1000.0 18,592.0 6.0 0.0 0.0 6.0 3.0E–5
16,679.0 0.0 699.8 0.0 8.1 0.0 262.2 0.0 3.0E–5
16,693.0 0.0 700.6 0.0 8.1 0.0 262.5 0.0 3.0E–5
16,753.0 0.0 703.9 0.0 8.1 0.0 263.5 0.0 3.0E–5
16,788.0 1000.0 1000.0 18,592.0 6.0 1.9E–11 6.8E–12 6.0 3.0E–5
16,803.0 999.4 1001.0 18,582.0 6.01 3.8E–11 3.9E–12 5.996 3.0E–5
16,863.0 997.0 1005.0 18,542.0 6.03 2.42E–11 1.64E–11 5.982 3.0E–5
0.0
1.57E–13
1.04E–13
1.38E–13
2.9E–14
3.85E–14
1.47E–13
—
1.62E9
1.63E9
1.66E9
0.667E9
0.67E9
0.68E9
—
41.38
41.58
42.42
17.03
17.11
17.46
parameters are tabulated in Table 3 along with the condition number of the regression matrix and results suggest that it suffers from ill-conditioning. However, it should be noted that the residual unbalance could be estimated accurately and this has been reported in [16]. To improve (i.e to reduce) the condition number of the regression matrix (i.e. Eq. (32)) the last column is scaled by a factor of o2 1E4, so that the last column entries are of the same order as that of the rest. The corresponding condition number of the matrix is also tabulated in Table 3. Drastic improvement in the condition number occurred due to the scaling, however, for the present case no improvement in the estimated has been found. The similar identification exercise is performed by contaminating simulated responses by the random noise of 1% or 5%. The corresponding
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1093
estimated parameters and the condition number of the regression matrix are also tabulated in Table 3 and most of the bearing dynamic parameters have ill-conditioning effects. By rotating the rotor alternatively in the clockwise or counter clockwise direction the rotor response is generated at different unbalance configurations at a particular rotor speed by using Eq. (20). The simulated response is fed into the identification algorithm of Method 2 (i.e. Eq. (37)) to estimate the residual unbalance and bearing dynamic parameters. The estimated parameters are tabulated in Table 3 and results suggest that the regression equation is now well conditioned. The corresponding condition number of regression matrix (i.e. Eq. (38)) is also tabulated in Table 3 and as compared to Method 1 the condition number of the regression matrix is much lower. For the present case, however, it is observed that there is no improvement in the estimates of the parameters. The same exercise is performed by contaminating simulated responses by the random noise of 1% or 5%. The corresponding estimated parameters and condition numbers are also tabulated in Table 3 and they suggest the robustness of the present algorithm against measurement noise. With the auxiliary unbalance unit (i.e. Method 3) the identification algorithm suffers from illconditioning similar to the Method 1 and result trends are identical to the Method 1.
5. Conclusions The present method is aimed at improving the existing method of estimating the bearing dynamic parameters by simultaneously estimating the unknown residual imbalance, which might cause error in the estimates of bearing parameters, inadvertently. For a linear multi-degreeof-freedom rotor-bearing system a general identification algorithm of systems’ unknowns has been presented. Then, an identification algorithm for simultaneous estimation of the residual unbalance and bearing dynamic parameters has been developed. The problem is found to be ill-posed. By taking the response measurements, while rotating the rotor alternatively in clockwise or counter clockwise gives identification algorithm as well-conditioned. The algorithm is tested against the measurement noise and found to be robust. However, with an auxiliary unbalance unit (in which excitation frequency can be chosen independent of rotor speed both in magnitude and sense) does not help in improving the estimates of bearing dynamic parameters while estimating residual unbalance simultaneously as against while estimating the bearing dynamic parameters alone [12]. The present method proposes essentially three methods by which to improve the condition of matrix A. The first involves multiple trial masses, which was noted in the introduction as being uneconomical and in Section 3.3.1 as ineffective. The second involves running the rotor in opposite directions; however, this is rarely possible in real rotors. Moreover, it is possible that the bearing coefficients also change. It will be more appropriate to call bearing stiffness and damping coefficients in this paper as bearing support parameters. The third measure involves an auxiliary unbalance unit. This is useful educational and research experimental set-ups. The estimation of bearing dynamic parameters is still thought to be relatively reliable only in laboratory test rigs [11]. However, it is unlikely to fit one to a working rotor. The development given in Section 2 potentially takes this idea of the present paper to the stage where it might be applied to a multiple-rotor machine (like a gas-turbine engine).
ARTICLE IN PRESS 1094
R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
Acknowledgement The author acknowledges anonymous referees for their useful contribution in improving the presentation of the paper especially Section 2. Appendix A. Nomenclature c the bearing damping coefficient C the damping matrix D k+jc D the dynamic stiffness matrix G theffiffiffiffiffiffigyroscopic matrix p ffi j 1 k the bearing stiffness coefficient K the stiffness matrix m the rotor mass per bearing M the mass matrix t the time instance U the unbalance (mass eccentricity) x, y the orthogonal axis system X, Y the complex displacements o the rotor rotational frequency O the excitation frequency from an independent unbalance excitation unit Subscripts r the real part i the imaginary part Superscripts Res the residual unbalance Tri the trial unbalance
Appendix B. condition number The condition number of a matrix (with respect to inversion) measures the sensitivity of the solution of a system of linear equations to errors in the data. It gives an indication of the accuracy of the results from matrix inversion and the linear equation solution. Values of the condition number near 1 indicate a well-conditioned matrix. The 2-norm condition number, which is used in the present case, is defined the ratio of the largest singular value of regression matrix to the smallest [18].
References [1] W. Kellenburger, Should a flexible rotor be balanced in N or N+2 planes?, Transactions of the American Society of Mechanical Engineers: Journal of Engineering for Industry 94 (1972) 560–584.
ARTICLE IN PRESS R. Tiwari / Mechanical Systems and Signal Processing 19 (2005) 1082–1095
1095
[2] J. Drechsler, Processing surplus information in computer aided balancing of large flexible rotors, Institution of Mechanical Engineers Conference on Vibrations in Rotating Machinery, Cambridge, UK, 1980, pp. 65–70. [3] P. Gnilka, Modal balancing of flexible rotors without test runs: an experimental investigation, Journal of Vibration 90 (1983) 157–172. [4] J.M. Krodkiewski, J. Ding, N. Zhang, Identification of unbalance change using a non-linear mathematical model for rotor bearing systems, Journal of Vibration 169 (1994) 685–698. [5] M.S. Darlow, Balancing of High-Speed Machinery, Springer, Berlin, 1989. [6] S. Edwards, A.W. Lees, M.I. Friswell, Experimental identification of excitation and support parameters of a flexible rotor-bearings-foundation system from a single run-down, Journal of Sound and Vibration 232 (5) (2000) 963–992. [7] J.W. Lund, Review of analytical methods in rotor-bearing dynamics, Tribology International 13 (5) (1980) 233–236. [8] P. Ramsden, Review of published data and their application to the design of large bearings for steam turbines, Proceedings of Conference in Lubrication and Wear: Fundamentals and Application to Design, IMechE 182 (3A) (1967–68) 75–81. [9] M.J. Goodwin, Experimental techniques for bearing impedance measurement, Transactions of the American Society of Mechanical Engineers: Journal of Engineering for Industry 113 (1991) 335–342. [10] E.E. Swanson, R.G. Kirk, Survey of Experimental Data for Fixed Geometry Hydrodynamic Journal Bearings, ASME Journal of Tribology 119 (4) (1997) 704–710. [11] R. Tiwari, A.W. Lees, M.I. Friswell, Identification of dynamic bearing parameters: a review, Shock & Vibration Digest 36 (2) (2004) 99–124. [12] R. Tiwari, A.W. Lees, M.I. Friswell, Identification of speed-dependent bearing parameters, Journal of Sound and Vibration 254 (5) (2002) 967–986. [13] J.H. Chen, A.C. Lee, Estimation of linearised dynamic characteristics of bearings using synchronous response, International Journal of Mechanical Science 37 (2) (1995) 197–219. [14] A.C. Lee, Y.P. Shih, Identification of the unbalance distribution and dynamic characteristics of bearings in flexible rotors, Proceedings of IMechE, Part C: Journal of Mechanical Engineering Science 210 (5) (1996) 409–432. [15] R. Tiwari, N.S. Vyas, Parameter estimation in imbalanced non-linear rotor-bearing systems from random response, Journal of Sound and Vibration 208 (1) (1997) 1–14. [16] J.K. Sinha, M.I. Friswell, A.W. Lees, The identification of the unbalance and the foundation model of a rotating machine from a single run-down, Mechanical System and Signal Processing 16 (2–3) (2002) 255–271. [17] N. Dharmaraju, R. Tiwari, S. Talukdar, Development of a novel hybrid reduction scheme for identification of an open crack model in a beam, Mechanical System and Signal Processing 19 (3) (2005) 633–657. [18] G.H. Golub, C.F. Van Loan, Matrix Computations, Third ed, The Johns Hopkins University Press, London, 1996.