Parametric Identification of a Dynamic Foundation Model of a Rotary Machine Using Response Data Due to Unknown

The estimation of a model of the foundation of a rotary machine has been recently attempted by using the difference between two sets of response data at some of the bearing locations from two consecutive rundowns of the machine, with and without known unbalance weights at certain positions on the two balance discs of each rotor respectively. However, it would be a great advantage to be able to perform the estimation with a single rundown. Due to practical restrictions in performing such tests (accessibility, costs etc.), there are cases in which data for only one rundown are available. In this case, the unbalance configuration is unknown and has, therefore, to be estimated, in addition to the unknown foundation model. Due to the special form of the unbalance force, this overall inverse problem can be solved by eliminating the unbalance configuration from the model estimation process. The remaining equation to estimate the foundation model consists of the projection of the response data, where the associated projector depends on the foundation model parameter. First results using the method, applied to a laboratory test rig and to a commercial turbo-generator, are presented.


INTRODUCTION
The successful condition monitoring of generators in modern power stations can be significantly enhanced by reliable mathematical models of the complete machines.Although sufficiently reliable models of the rotor and the bearing are well established, the influence of the foundation on the machine dynamics is not yet fully understood.In recent years several attempts have been made to model the foundation by finite elements, but due to the complexity of the foundation those attempts revealed unsatisfactory results (Lees and Simpson,  1983).454 U. PRELLS AND A.W. LEES Mathematical modelling is always purpose- orientated (Natke, 1995).In case of foundation modelling the purpose is not focussed on estimat- ing foundation mass and stiffness but to establish a foundation model which reproduces the contribu- tion of the foundation to the dynamic of the entire system with sufficient accuracy.The criterion for the quality of the model estimate is the fit between calculated model response and measurement.In order to determine the contribution of the foundation to the rotor's dynamic behaviour, the real- valued symmetric matrices Ar of a frequency filter model* (for example Mottershead and Stanway,  1986)   N F(co) (jco) Ar E C nxn   (1) r=0 have to be estimated using response data at the bearing locations during a machine rundown cover- ing a frequency range c E 2. Note that the matri- ces Ar have to be symmetric because Maxwell's Theorem of Reciprocity must hold true.Although F results analytically from dynamic condensation (see Appendix A), Ar can be assumed to be real- valued because the damping contribution of the foundation is negligible.Assembling the np'= n(n + 1)(N + 1)/2 independent entries of the matri- ces Ar in one vector x Rn, the dynamic stiffness matrix F(c, x) can be understood as a function of the model parameter vector (see Appendix B).For given model parameters x, F(c,x) maps the response u(co) C to the forces fB(co) C at the bearing locations.The latter can be calculated as fB(co) Q(co)u(co) + C(co)p, (2) since the matrices Q(co) and C(co) depend only on the models of the rotor and the bearings, as shown in Appendix A. The real-valued vector p R d repre- sents the unbalance configuration and consists of the masses, the eccentricities and the angles between the positions of the masses on the balance discs and the shaft marker.Using for u(co) the difference of the responses of two consecutive rundowns, with and without balance weights, the parameter vector x can be estimated from F(co, x)u(co) fB(co). (3) This method of estimating a foundation model has been discussed in several papers, for instance Lees (1988), Zanetta (1992), Feng and Hahn (1995), Vania (1996), Smart et al. (1996), Lees and Friswell  (1997)   (1997,  1998).However, the applicability of this method is based on two major requirements: (1) the unbalance configuration p has to be known, (2) the vertical and horizontal responses at all bear- ing locations have to be measured.
These criteria are often not fulfilled, due to equip- ment failure, costs or accessibility.The objective of this paper is to explore the possibility of estimating the model parameters x in the case where only an incomplete set of data from one rundown is avail- able, i.e. the unbalance configuration p is unknown and only n' components of the n-dimensional re- sponse vector have been measured, i.e.
where H0 E 1R ' is the measurement matrix which selects the measured part of the complete response vector (see Appendix A).The main problem here is that the unbalance configuration p is unknown.One possibility consists in estimating p in addition to the model parameter vector x. Background to that problem can be found in Lees and Friswell   (1997).However, the method presented in that paper requires response data at all bearing loca- tions, which are often not available.The method presented in this paper aims to separate the esti- mation of the model parameters x R n" from the estimation of the unknown unbalance configurationp R d * Note that in the case N 2 the matrices A0, A1, A2 correspond to the contribution of stiffness, damping and inertia respectively.
In the first section, the basic estimation equation is introduced and extended to the case of unknown unbalance configurations.In the second section, the numerical handling of the resulting inverse problem is discussed and countermeasures to regu- larise the problem are suggested.x) Q(co))-I x)p, the equation error can be calculated as g-V(x)V+(x)g (hmn' P(x))g-N(x)g.( 8) The symmetric and idempotent matrix P(x) is the projector into the subspace spanned by V(x) and N(x) is the projector into the orthogonal comple- ment.Since N(x) varies with x this method is some- times called variable projection method (see for instance Golub and Pereyra, 1973; Krogh, 1974;  Kaufman, 1975).A cost function can be formulated using the squares of the relative Euclidean norm of Eq. ( 8) J(x) "-IlN(x)gllZ/llgll 2 g-CP(x)g/llgll 2 (9) Note that 0 <_J(x)_< due to the orthonormal decomposition g/llg[I P(x)g/llg[I + N(x)g/l[gll. (lO)   where all terms are complex.For a set of discrete frequencies ft := {cl,...,COm} Eq. ( 5) can be ex- tended and written in real-valued form by doubling the order (see Appendix C) g V(x)p. (6) Here g E I 2mn' contains the measured responses and V(x)E]2mn'xd is the real equivalent of the frequency response matrix Z(co, x).Note that Eq. ( 6) is linear in the rotor unbalance configuration p.
Equation (6) states that there exists a linear com- bination of the columns of V(x) that generates g, which means that the vector g, the generalised response vector, is contained in the subspace spanned by V(x).Although p and x are unknown this statement is independent of the value of p.To realise this one may solve Eq. ( 6) in the least-squares sense for p yielding where V + is the Moore-Penrose inverse of V.
Inserting the normal solution/5 of p into Eq. ( 6 If x has been estimated by minimising J(x) then Eq. ( 7) delivers an estimate for p.Although the formulation (9) of the estimation problem has the advantage of being independent of the unknown unbalance configuration p, it has the disadvantages (existence, uniqueness and stability) typical of out- put residual-related methods for solving non-linear inverse problems (Natke et al., 1995).In the next section, the existence, uniqueness and stability of the solution are explored, and the numerical and computational aspects of the estimation procedure are explained.

THE ESTIMATION PROBLEM
The existence of a solution depends essentially on two conditions" (1) the projectors must be non-trivial, i.e.P(x) I2m,,', or equivalently N(x) O, (2) the degree ofthe chosen filter model (see Eq. (1)) must reflect the most important peaks of the available data.
A sufficient condition for the non-triviality of the projectors is 2mn' > d, which is the dimension of the unbalance vector p.If one assumes that there are two balance discs on each shaft, then d= 4ns (in the case of an ns-shaft rotor).Since each shaft is supported by two journal bearings and the mea- surements are taken at each bearing in horizontal and vertical directions, the maximum number of responses is 4nsn _> n'.From this one finds rn > n/2n', which means, in principle, two frequencies are suf- ficient in the case that only half of all possible horizontal and vertical displacements .havebeen measured.However, this limit case is far from real- istic.Even if the structure only has a single reso- nance, rn 2 frequencies are insufficient to sample that peak.An optimum discretisation in the fre- quency domain for the purpose of model updating has been discussed by Cottin (1991).It was stated that at least three (better five) frequencies per reso- nance are necessary, where the step size depends on the damping ratio of each individual peak.Moreover the usual identifiability conditions (see for instance Natke, 1992) have to be satisfied, which requires a sufficiently large number of resonances within the frequency range to ensure a small con- dition number of V(x).
The uniqueness of the estimates depends on the chosen degree of the frequency filter model.It is well known from polynomial approximations that the error between model and data can be made arbitrarily small by increasing the degree.The determination of an adequate degree for the filter model is a non-trivial task, because the peaks of the data may be (1) due to measurement errors or (2) due to resonances of the subsystem.

rotor/bearing
Those resonances of Z, which are sensitive with respect to resonances of the foundation, are diffi- cult to distinguish.Investigations using simple test models reveal that often peaks of the frequency response matrix Z cannot be controlled by adjust- ing the resonances of F. The following two limit cases are helpful to find those peaks of Z which are sensitive with respect to the peaks of F: lim Z(o, x) -H-Q-1 ()C(), (12)   lim Z(o,x) HF-I(o,x)C(), (13)   Ilxll q with the positive scalar q >> Q(w)[I, Vco E f.In the first limit case, the frequency response function reflects only the contribution of the rotor/bearing model, whilst in the second case the influence of the rotor/bearing model contained in Q is negligible, and the contribution of the rotor/bearing model to the response is only due to the matrix C. Plotting the peaks of the response function Z(co, x) for some cases Ilxll [0, q], using a simple initial foundation model is useful in deciding the appropriate model degree N. Of course, this choice is closely related to the stability of the estimates, i.e. the sensitivity of the estimates with respect to data errors.
The stability of the estimates depends on three points: (1) the moderate choice of the degree of the filter model, (2) the condition of the matrix V(x), (3) the dimension of span(V(x)) relative to the dimension of x.
These are related in their effect on the foundation parameter estimates.As already mentioned, one has to analyse the data in order to choose an appro- priate degree of the filter model.However, data are always corrupted by noise, thus deciding which peaks are real and which are due to noise is rarely simple.In the case that foundation and noise- related peaks are of comparable magnitude, fre- quency filter models of different degrees will lead to the same order of equation error.To avoid non- unique estimates, those data which are most noise corrupted have to be excluded from the estimation process.This can be done by introducing a moder- ate weighting, or simply by using only those data with relatively large magnitudes.Assuming that the same type of sensors are used, data with small relative norms ri max lui()l ( 14) will be more noise polluted.For a given threshold r0 > 0 the selection of data .(CO){,/(&): ri ro} (15) for the estimation process leads to a reduction of the dimension of the parameter space.A reduction in the number of parameters to be estimated not only stabilises the estimates but also reduces the degree of non-uniqueness.A large number of parameters is related to a high degree frequency filter model.The larger the chosen frequency filter model degree, the more likely are non-unique estimates, as different parameter vectors xx may generate the same space, span(V(x)) span( V(x ')).The physical explanation is that some parameters are mainly related to resonances of Z outside the frequency range f.A change in those parameters will have little effect on the subspace spanned by V. A method to eliminate parameters related to those modes a posteriori has been reported elsewhere (Friswell et al., 1997).However, this method assumes that a high-dimensional parameter vector has been estimated, which is often not feasible due to numerical and computational limitations.The problem of local minima of output-related methods is well known (Oeljeklaus, 1999).Depend- ing on the size of the problem, the cost function J(x) defined in Eq. ( 9) has several local minima.Those local minima are due to several levels of fit between peaks of the data and the model response.If, for instance, the model response for x matches the data at one peak whilst the model response for x does not match any peak of the data, then the cost func- tion value J(x) is lower than J(x').Since there are in general model parameters, which lead to model responses that match the data at several peaks, the associated cost function values are local minima.
An iterative solution procedure, as for instance Newton-Raphson, will, in general, lead to a local minimum, depending on the quality of the chosen initial parameter vector.To avoid the use of locally convergent non-linear solvers, a controlled random search has to be applied which evaluates the cost function many times (several thousand evalua- tions).The time necessary for one evaluation of J(x) depends essentially on the inversion of the dynamic stiffness matrix (see Eq. ( 5)).A modal re- presentation of the dynamic stiffness matrix is, in general, not possible, due to the speed-dependence of the model of the journal bearings.Thus, the inversion has to be performed at each frequency step, which is time-consuming.Already for rela tively small problem sizes, the evaluation of the cost functions at a given parameter vector x can take several seconds.Even on a fast computer this can lead to several hours of execution time.The over- all computation time of the estimation method depends (1) on the number of parameters to be estimated, (2) on the number of frequency points and (3) on the quality of the initial parameter vector.
The choice of the initial model parameter is a difficult task, since, in general, the foundation model estimation is a 'black-box problem', i.e. no a priori model is available.In order to find an appropriate initial parameter vector Xini, it is sug- gested to start with a 2nd order diagonal frequency filter model F(,Xini).This simple model is then optimised in a pre-process by maximising the num- ber ofpeaks of V(xini) within the frequency range f, using the limit cases defined by Eqs. ( 12) and ( 13).
This strategy tries to minimise the risk of finding a local minimum.Using this initial parameter vector, the cost function J(x) is minimised by a modified Nelder-Mead simplex algorithm (for details see Nelder and Mead, 1965) using only those response vector components which have magnitudes above the threshold r0.If the minimum of the cost function is insufficient, the degree of the frequency filter model is increased.
In the next section, the results of two applications are presented to demonstrate the capability of the method described above.

APPLICATIONS
The first application is an experimental rotor rig located in the laboratory of the Department of Mechanical Engineering at the University of Wales Swansea.The second application is a commercial turbo-generator.The estimation routines have been coded in MATLAB and executed on a PI1450 MHz PC and with 256 MB RAM under WINDOWS'98.

Experimental Rotor Rig
A schematic of the the rotor rig is depicted in Fig. 1.A detailed description of the test rig can be found in Edwards et al. (1999).It consists of a solid steel shaft of length 740 mm and diameter 10 mm.At positions A and B the shaft is connected via two brass-bush bearings to a foundation, hence the foundation model has dimension n 4. At position M, the rotor is connected to an electric motor, via a flexible coupling.The bearings are assumed stiff and free of damping.At positions A and B, piezo-electric accelerometers are attached in the horizontal and vertical directions.The responses at the bearing locations due to an unknown residual unbalance distribution (pl ps) along the shaft have been measured, during a slow rundown of the rotor covering a frequency range between 15 and 45 Hz.Thus, in this case, n' n 4 and d 8. Since the number of frequency points was m 1343 the size of V(x) is 10744-by-8.The dimension of the parameter vector x dim(x) 10. (N-t-1) (16) depends on the chosen degree N of the filter model.
The time necessary to evaluate J(x) was about 9 s and was independent of the chosen degree of the filter model.Diagonal frequency filter models of various degrees have been calculated in a pre- process in order to maximise the number of peaks of V(x) within the running range.Using these initial models, the cost function J(x) has been minimised.
The minimum values of the cost function for different degrees of foundation frequency filter  FIGURE 2 Measured (dotted) and calculated responses (solid) using a foundation frequency filter model of degree 3. 459 models are shown in Table I.For frequency filter models of degrees 3 and 4, the cost function minimum is about 0.07.In contrast to the model of degree 3 the additional matrix A4, that has to be estimated for the model of degree 4, is almost zero.Thus, both models lead to the same estimates.
Using the foundation model estimate of degree 3, the residual unbalance vector p can be estimated by U. PRELLS AND A.W. LEES Eq. ( 7) yielding 10 .4. (-0.547, -0.099, 1.752, 0.26, 1.848, -0.322, 0.511,0.128)T, ( 17) which appears to be physically realistic.Using the estimates of the foundation model and of the unbalance, the response can be calculated.In Fig. 2, the model response (solid) is compared to the mea- sured data (dotted).The calculated displacements match the data at the two groups of peaks in the upper and lower frequency range.Although the match about 32 and 28 Hz is not quite right, the overall fit is satisfactory.This lack of fit is due to the fact that the data signal is relatively low at those frequencies.
Application to a Power Turbo-Generator   The method outlined above has been applied to a commercial turbo-generator.The rotor, in this case, consists of 6 shafts which are connected via 12 journal bearings to a steel foundation; hence n 24.
A finite element beam model of the rotor was statically condensed from over 2000 to n-196 degrees-of-freedom (DoF).The non-linear model of the journal bearings was linearised at each frequency point to deliver a frequency dependent damping and stiffness matrix.Twenty-two hori- zontal and vertical responses have been measured at all bearing locations during one rundown, covering a frequency range between 5 and 50 Hz, with rn 275 frequency points.In Fig. 3, the relative maximum magnitudes of the measured response vector components are plotted.The first component corresponds to the horizontal DoF at the first journal bearing of the HP rotor, and the last component correspond to the vertical DoF of the 2nd journal bearing of the exciter.In order to reduce computational effort, a frequency filter model of degree N= 2 has been chosen to model the dynamic contribution of the foundation.Moreover, the data used for the estimation have been restricted to those components having relative maximum magnitudes larger than r0 0.5.The n= 4 selected data components are circled in Fig. 3.A residual unbalance distribu- tion was assumed along the entire rotor.Taking into account all translational DoF, except those at the bearing locations, the dimension of the unbal- ance vector p was d 72.In the N--2 case consid- ered here, the parameter vector x is of dimension np 24(25 + 1)/2 900, and the matrix V(x) (see Eq. ( 6)) has the size 2200-by-72.Evaluation of the cost function J(x) takes about 26s.Using an undamped diagonal initial model, estimated in a pre-process to maximise the number of peaks F() within the frequency range, the controlled random j search procedure has been initiated to minimise the N cost function J(x).The value of the cost function at Ar the initial parameter vector was 0.86.After 467 improvements and a total computation time of more than 48 h, the value decreases to 0.12. Figure 4 shows the data (dotted) and the calculated model C responses (solid).
Although the model response does not fit pre- cisely, the overall shapes of the spectrum are reproduced.Due to the relatively low degree of the frequency filter model, no better fit was expected, u() The necessity of using a frequency filter model of higher degree was confirmed by the corresponding unbalance estimate: the values were too large to be physically realistic.It is also possible that the algo- rithm has found a local minimum.The use of a fre- p quency filter model of higher degree would lead to an increase in the overall computation time.The first attempt to estimate np= 2100 parameters of a Q() frequency filter model of degree N 6 was cancelled C(a) after 4 days of execution time because no significant improvement of the cost function value could be ei found.

CONCLUSIONS
A method has been applied to a high-dimensional parameter estimation problem for the estimation of a rotor foundation model in the case of incomplete measurements and unknown unbalance regimes.The usability of the method has been tested by two practical applications: an experimental test rig and an industrial power turbo-generator.Although the results are so far encouraging, the computational problems encountered require further investigation.

NOMENCLATURE z(,x)
frequency in rad/s frequency filter complex unit (-1) 1/2 degree of filter model symmetric and real-valued matrix for each r 0,..., N set of all complex-valued n-by-n matrices complex vector space of dimension n n-dimensional complex-valued force vector at the bearing locations n-dimensional complex-valued response vector n'-dimensional complex-valued vector of measured responses real vector space of dimension d real-valued d-dimensional vector of unbalance configuration complex-valued n-by-n matrix complex-valued n-by-d matrix k-by-k identity matrix ith column vector of the identity matrix real-valued n-by-n selecting matrix real-valued np-dimensional In general let Ax():= -co2Mx-b-jcoC2c-b Kx denote the dynamic stiffness matrix at frequency , where the subscript X E {R,B,F} refers to the rotor, the bearing and the foundation respectively.Mx, Cx and Kx denotes the matrices of inertia, damping and stiffness.The dynamic stiffness matrices are partitioned according to the internal DoF, indi- cated by an additional subscript/, and by those DoF connected to the bearings, indicated by the additional subscript B: A FUF fF >> AFro A FII U FI In Eqs. ( 18)-( 20) the two rules of model synthesis are involved: (1) Kinematic compatibility requires the responses at the connection DoF to be the same.
(2) Due to Newton's principle of actio et reactio the forces at the connection DoF differ only in sign.
Moreover it is assumed that besides the unbal- ance force ft0 of the rotor no other (relevant) forces act on the system.Of course, the non-zero entries of the force vector fRi are the components of the vector of unbalance forces ft0, i.e. there exists a selecting matrix S := [ei,..., ei] such that fm= Svfv. (21) The unbalance force has the form fU(OO) C02Ip, where the real-valued d-dimensional vector p of the unbalance configuration contains the eccentricity, the mass and the angle to the shaft marker of each balance disc.The matrix represents the harmonic excitation in the frequency domain and has a block- diagonal structure j 0 (23) In addition, the special structure of the bearing model is already taken into account.It can be shown that the matrix B B(o) is block diagonal, i.e.where in general the matrix BiE C 22 for the ith bearing is given by Bi Ki(&) +j&Di(), { 1,..., rib}. (25) The stiffness matrices and the damping matrices, Ki, Di 22, result from linearisation and are in general non-symmetric and non-singular, and depend on the excitation frequency.The dynamic stiffness matrix F used in Eqs.
=:f Note that AFII is non-singular because the entire system is grounded.Coupling the models of the rotor from Eq. ( 18) and the model of the bearings from Eq. ( 19) leads to the input/output equation of for symmetric real-valued matrices.Each matrix At, r 0,..., N, of the frequency filter of degree N can be written uniquely as where Z(co, x) H(F(co, x) + Q(co))- For m excitation frequencies Eq. ( 44) can be extended to give gc-Vc(x)p, where .() / I z(, x) g <(x) :: .

FIGURE 3 FIGURE 4
FIGURE 3 Relative maximum magnitudes of the response vector components.
vector of model parameters complex-valued n-by-d frequency response matrix set of m discrete frequencies Because of commercial confidence the depicted values are scaled to unity.
and real-valued adjustment parameters xr.Collecting all parameters in one vector x of dimension np := (N+ 1)n(n + 1)/2 the frequency filter can be written as np F(co, x) Z Ra(co)Xa"

F
in Eq. (40) one finds u M (co) Z(co, x)p,

E
EN NE ER RG GY Y M MA AT TE ER RI IA AL LS S Materials Science & Engineering for Energy SystemsEconomic and environmental factors are creating ever greater pressures for the efficient generation, transmission and use of energy.Materials developments are crucial to progress in all these areas: to innovation in design; to extending lifetime and maintenance intervals; and to successful operation in more demanding environments.Drawing together the broad community with interests in these areas, Energy Materials addresses materials needs in future energy generation, transmission, utilisation, conservation and storage.The journal covers thermal generation and gas turbines; renewable power (wind, wave, tidal, hydro, solar and geothermal); fuel cells (low and high temperature); materials issues relevant to biomass and biotechnology; nuclear power generation (fission and fusion); hydrogen generation and storage in the context of the 'hydrogen economy'; and the transmission and storage of the energy produced.As well as publishing high-quality peer-reviewed research, Energy Materials promotes discussion of issues common to all sectors, through commissioned reviews and commentaries.The journal includes coverage of energy economics and policy, and broader social issues, since the political and legislative context influence research and investment decisions.S SU UB BS SC CR RI IP PT TI IO ON N I IN NF FO OR RM MA AT TI IO ON N Volume 1 (2006), 4 issues per year Print ISSN: 1748-9237 Online ISSN: 1748-9245 Individual rate: £76.00/US$141.00Institutional rate: £235.00/US$435.00Online-only institutional rate: £199.00/US$367.00For special IOM 3 member rates please email s su ub bs sc cr ri ip pt ti io on ns s@ @m ma an ne ey y. .cco o. .uuk k E ED DI IT TO OR RS S D Dr r F Fu uj ji io o A Ab be e NIMS, Japan D Dr r J Jo oh hn n H Ha al ld d, IPL-MPT, Technical University of Denmark, Denmark D Dr r R R V Vi is sw wa an na at th ha an n, EPRI, USA F Fo or r f fu ur rt th he er r i in nf fo or rm ma at ti io on n p pl le ea as se e c co on nt ta ac ct t: : Maney Publishing UK Tel: +44 (0)113 249 7481 Fax: +44 (0)113 248 6983 Email: subscriptions@maney.co.uk or Maney Publishing North America Tel (toll free): 866 297 5154 Fax: 617 354 6875 Email: maney@maneyusa.comFor further information or to subscribe online please visit w ww ww w. .mma an ne ey y. .cco o. .uuk k C CA AL LL L F FO OR R P PA AP PE ER RS S Contributions to the journal should be submitted online at http://ema.edmgr.comTo view the Notes for Contributors please visit: www.maney.co.uk/journals/notes/emaUpon publication in 2006, this journal will be available via the Ingenta Connect journals service.To view free sample content online visit: w ww ww w. .i in ng ge en nt ta ac co on nn ne ec ct t. .cco om m/ /c co on nt te en nt t/ /m ma an ne ey y Friswell et al. (1997) and Prells et al.