Model Based Fault Identification in Rotor Systems by Least Squares Fitting

In the present paper a model based method for the on-line identification of malfunctions in rotor systems is proposed. The fault-induced change of the rotor system is taken into account by equivalent loads which are virtual forces and moments acting on the linear undamaged system model to generate a dynamic behaviour identical to the measured one of the damaged system.


INTRODUCTION
Rotating machinery installed in power plants are normally equipped with vibration sensors which continuously provide large amounts of vibration data during operation, (Bruel and Kjaer, 1986).For diagnosing the state of a machine, usually signal based monitoring systems are used as a good tool, although they do not fully utilise the information contained within the vibration data, (Bruel and Kjaer, 1986); (Glendenning et al.,  1997).These approaches to machinery diagnostics are generic rather than machine specific and the interpretation of the data is based on qualitative rather than quantitative information.
Contrary to signal based monitoring systems, model based diagnostic systems developed in recent years utilise all information contained in Corresponding author.Tel." + + 49-(0)6151-16 2785, Fax: + + 49-(0)6151-16 4125, e-mail: markert@mechanik.tu-darmstadt.de the continuously recorded vibration signals by taking into account models of the specific machine and possible faults as well as the recorded data of the intact machine.The new method may be used together with or alternatively to conventional sig- nal based monitoring systems.
Model based monitoring systems give more accurate and faster information than conventional signal based systems, since a-priori information about the rotor system is systematically included in the identification process.Therefore, the type, position and severity of a fault can be estimated with more reliability and in most cases during operation without stopping the machine.
In the BRITE--EURAM project MODIAROT differ- ent model based identification methods for finding malfunctions in rotor systems of power plants have been developed, (Penny et al., 1995); (Mayes  and Penny, 1998); (Bachschmid and Dellupi,  1996).These methods work either in the time or in the frequency domain depending on the mal- function type and the operating state for which the vibration data are available.In this paper a new time domain identification method is presented based on least squares fitting in the time domain.
For a model based identification of faults, reliable mathematical models of the rotor system and of possible faults are required.The mathematical models have to be so simple that fast on- line calculations are possible, but must be exact enough to reproduce the measured vibrations of a real rotor-bearing-foundation system.
Exact models of faults are in many cases non- linear.Therefore, time integration algorithms are needed to solve the non-linear equations of the damaged rotor system directly, although the undamaged system is linear.Such a direct identifi- cation process is extremely time consuming, espe- cially for models having many degrees of freedom.On-line identification of faults would not be possible with this approach, because the time con- suming integration has to be carried out for every possible position and extent of each fault type.
The new model based method for identification of malfunctions avoids the non-linearities, which are normally brought into the system equations by fault models.The method is based on the idea that the faults can be represented by virtual loads acting on the linear model of the undamaged system.Equivalent loads are fictitious forces and moments which generate the same dynamic be- haviour as the real non-linear damaged system does.System models being originally linear remain linear and unchanged, even if non-linear faults occur.Therefore, to handle the system equations, only simple and fast mathematical procedures are necessary.The identification of the nature, position and severity of faults can be carried out on-line during operation, even if the number of degrees of freedom is very high.

IDENTIFICATION METHOD Mathematical Description of the Method
Normally, the model of the undamaged rotor system is linear and described by linear equations of motion: Moi:o (t) + Bofo(t) + Koro(t) Fo(t). ( 1)   The vibrations of the N degrees of freedom due to the load F0(t) acting during normal operation are described by the N-dimensional vector r0(t).M0, B0 and K0 are the mass, damping and stiffness matrices of the undamaged system consisting of rotor, bearings and foundation.Some of these ma- trices may change with the rotor speed f due to gyroscopic effects and journal bearings.
The fault-induced change of the dynamic be- haviour depends on the malfunction type, posi- tion and severity, (Bach et al., 1998).These fault parameters are put together in the fault parameter vector [.For example, the fault vector of a single unbalance includes the position, size and phase angle of the unbalance.
In reality a fault changes the system properties which results in a system behaviour r(t) different from the normal vibrations r0(t).The differences between the vibrations r(t) of the damaged system and the normal vibrations r0(t) of the undamaged system are residual displacements, velocities or accelerations, respectively.These are given by Ar(t) r(t) ro(t), (2) zxe(t) e(t) eo(t), (3) On the other side, one can imagine that the residual vibrations Ar(t) could be caused by additional time-varying forces acting on the undamaged system, Eq. ( 1), instead of being caused by a fault, Mob(t) + Bof(t) + Kor(t) Fo(t) + AF([, t). ( 5) Subtraction of the equations of motion for the damaged and the undamaged system yields MoAi:(t) + BoAr(t) + KoAr(t) AF(I, t), (6)   which relates the measured residual vibrations directly to the equivalent forces of the faults present in the system.For calculating the equivalent load, which re- presents the fault arising in the rotor system, simply the residual vibrations resulting from the actual measurements and the corresponding nor- mal vibrations have to be put into Eq.(6).
vibrations to the time scale of the current measurement.Differences in the sampling fre- quencies can be eliminated by interpolating the recorded normal vibrations to the sampling times of the current measurement.Finally, phase shifts are avoided by recording a trigger signal.
On every measurement position of the rotor system either displacements, velocities or accelera- tions are measured using appropriate types of transducers.But, for the identification process, all three signal types are needed for every measure- ment position.Acceleration, velocity and displace- ment signals are obtained by either differentiation or integration of the direct measured signals.This can be done either in the time domain or via FFT in the frequency domain.

Modal Expansion
In practice, vibration transducers are installed only at a few positions of the rotating machinery.
Therefore, residual vibrations AM(t) are available only for a few degrees of freedom of the model.The number M of measurement locations is much smaller than the number N of the model's degrees of freedom.The data of the non-measurable de- grees of freedom have to be estimated from the measured signals.
The measurable part At(t) of the residual vector is related to the full residual vector A(t) by the measurement matrix C,

Signal Processing
For calculating the residual vibrations, measured vibration data for both the undamaged and the damaged rotor system have to be available for the same operating and measurement conditions.These are, for example, load, rotor speed, phase reference and sampling frequency.Because directly matching data are usually not available, some signal processing has to be done to achieve the same conditions.For example, small differences of rotor speeds can be compensated in some cases by adjusting the time scale of the recorded normal A4(t CA(t). (7) On the other hand, the full residual vector can be approximated by a set of mode shapes f which are put together in the reduced modal matrix Logically, the number K of mode shapes used may not exceed the number M of independently measured vibration signals contained in /(t), K <_ M. The vector of modal co-ordinates Aq(t) is estimated by combining the measurement Eq. ( 7) with the modal representation of the full residual vector and minimising the equation error by the least squares method.
Eventually, the full residual vector at all degrees of freedom is estimated by QA/(t), (10) where the constant matrix Q can be calculated in advance.
For illustrating the relation between fault parameters and the corresponding equivalent load, a single unbalance un with the phase angle 6n acting on the rotor at position number n is considered.The fault parameters of the single unbalance are given by [ [n, u,, 6n] r.The equivalent forces act- ing on node n are AF,h([, t) gt2un sin (f2t + cos (at + (13) in horizontal and vertical direction, respectively.The equivalent forces on all other nodes are zero.

Determination of Equivalent Loads
The equivalent force A'(t) characterising the unknown faults is calculated by putting the residual accelerations A(t), residual velocities A(t) and residual displacements A(t) of the full vibrational state into Eq.( 6), A'(t) MoQA4(t) + BoQAM(t) + KoQAM(t). (11) Only simple matrix multiplications and additions have to be carried out on-line for calculating the equivalent loads A'(t) from the measured vibra- tion signals.

Fault Models
For the proposed model based fault identification method, each fault has to be represented by a mathematical model describing the relation be- tween the fault parameters [ and the equivalent force AF(t).In this context, AF(Ig, t) is a math- ematical expression for the time history of the forces acting on each degree of freedom of the model.These fictitious forces depend on the fault parameters but in many cases also on the rotor speed f as well as on the vibrational state of the rotor system as described for a lot of faults in (Penny et al., 1999); (Bach et al., 1998) and (Bach  et al., 1998).
Least Squares Fitting in the Time Domain In practice, measurement noise and measurement errors, the limitation of the number of sensor sig- nals available as well as modelling inaccuracies of the rotor system generate errors in the measured equivalent forces.Therefore, a least squares algo- rithm is used to achieve the best fit between the measured equivalent forces A(t) and the theore- tical ones AFi([i, t) produced by the fault models by adjusting the fault parameters Ii for all faults taken into account, 2 dt Min.

(14)
Starting from an initial guess, the least squares algorithm varies the fault parameters [i of each fault model until the deviation is minimal.

Probability Measures
The quality of the fit achieved can be used to estimate the probability of the different identified faults.Two probability measures based on corre- lation functions have been developed and success- fully tested.
The first probability measure Pl, called coher- ence, is the normalised correlation of the identified equivalent forces zXFi(i, t) of a particular fault with the measured equivalent forces Al(t)for lag r 0, qS/XF"/X (0) (15) Pl V/AFi,AFi (0) A, A, (0) Due to the normalisation by the auto-correlation functions of AFi(i,t) and A'(t), the coherence takes values between -1 < Pl < 1, where pl means that AFi(li, t) matches A(t) perfectly.
The second probability measure p2, called intensity, measures the contribution of the particular fault to the measured total equivalent forces /Fiden AFi, ()AFi'AFi (0) (16) P2 (/)/kFidentAFident (0)" The intensity measure takes values in the interval 0 <p2 < 1, where P2-1 signifies that the identi- fied fault is the only one present in the rotor system.
Simulation studies showed that both measures Pl and P2 should be used simultaneously to evalu- ate the probability of a fault.Suitable threshold val- ues are Pl > 10% and p2 > 20% indicating that the specific fault is present in the rotor system.

Graphical Overview
In Figure the whole identification procedure is outlined.The measured vibration signals M(t) are the input and the fault parameters Ii for each fault type are the output.The normal vibrations, the model parameters of the undamaged rotor system as well as the fault models are taken from prior measurements and pre-calculations.

NUMERIC EXPERIMENTS A Simple Rotor-Stator-Model
The effects of different modelling, measuring and identification parameters as well as disturbances and errors in the measuring signals on the identification results have been investigated by means of numeric experiments with a simple rotor system.This simple rotor system consists of a rigid foundation body and a flexible rotor carrying four rigid disks (Fig. 2).The foundation as well as the rotor are supported by linear springs and dampers.Gyroscopic effects as well as cross coupling terms between the two radial directions are neglected, so that the dynamic behaviour is characterised by only six degrees of freedom all being in the same direction, two for the foundation and four for the rotor.The first natural frequencies and the corresponding modal damping coefficients are listed in Table I.The rotor speed was 9.5Hz, which is approximately 30% below the third natu- ral frequency of the rotor system.During normal operation of the undamaged system two residual unbalances luo l-50gmm and [uoNI--100grnm   excite the rotor at the middle disks rn2 and m3.Two faults were investigated exemplarily: additional unbalances at disks m2 and m3 and rubbing between disk m2 and foundation ms.

Identification of Unbalances
In different case studies the effects of various parameters on the identification results were examined independently.The additional unbal- ances ]b/2l--30 gmm and lu31 600 gmm represent the fault which has to be identified by the de- scribed algorithm.The first one is smaller and the second one is larger than the respective initial unbalances.The time window of the measurements covers about 5 rotor revolutions.
In the reference case 1, the time-histories of displacements, velocities and accelerations were exactly measured at all six degrees of freedom of the rotor stator system.Therefore, no modal expansion was necessary and applied for the reference case 1.The identification algorithm was constrained to search only for unbalances at the two middle disks.The identified unbalances the corresponding probability measures pj(uk) as well as the relative errors ek are listed in Table II.Figure 3   Natural Frequency fk 1.5 Hz 4.0 Hz 13.7 Hz 31.0Hz Modal Damping Dk 0.8% 2.1% 6.0% 27.0% displacements and Figure 4 the estimated equiva- lent forces at the six degrees of freedom.In this case, the unbalances were identified exactly and the estimated equivalent forces match exactly the ones of the two additional unbalances which exert on the rotor in reality.
In the following test cases, numbered 2 to 10, various influences on the results of unbalance identification are examined, which occur more or less strongly in practical applications: In Test case 2, only accelerations are measured at all degrees of freedom.Velocities and displace- ments are calculated by numeric integration.
In Test case 3, faults were searched for at all degrees of freedom of the rotor stator system.This means, the algorithm was allowed to identify also unbalances at positions, where none can occur.
However, no equivalent forces and unbalances were identified at the stator and at the two balanced rotor disks and 4.
In Test case 4, two different kinds of faults are searched for simultaneously: unbalance and rubbing.
In Test case 5, the time window for measurements was shortened to only two rotor revolutions.
In Test case 6, the forced vibrations excited by unbalances are superimposed by free oscillations decaying with time.
All influences described so far do hardly affect the identification result.The real unbalances at disks 2 and 3 were identified very exactly.
By contrast, the identification results are affected much stronger by measurement noise, calibration errors and reduction of the number of measuring locations.These influences are examined in the following four Test cases 7 to 10.
In Test case 7, all measuring signals were falsified by band-limited noise.The standard deviation of the disturbances was about 2% of the amplitude of the correct signal.
In Test case 8, a calibration error of 20% in the sensor mounted at disk m2 was simulated.But the dominating unbalance was identified exactly.
In Test case 9, the accelerations were measured at only 5 locations.The full vibrational state was estimated by means of modal expansion.The identification result is very accurate because the exact deformation can be approximated very well by the first five eigenvectors which are used.
The approximation of the full vibrational state becomes more inaccurate if measured signals are available only at a few locations.In the case 10, the accelerations were measured only at locations 2, 3 and 4. Therefore, only three eigenvectors could be taken into account for the modal ex- pansion procedure.Figure 5 shows the residual displacements.Three degrees of freedom were directly measured (thick solid lines), the other three were approximated by modal expansion (dashed lines).These are compared to the correct ones (thin solid lines).
Nevertheless, the dominating unbalance at disk 3 was identified quite well, even its phase angle was identified exactly.The modal expansion has the major influence on the identification results.The estimated equivalent forces are smeared over all degrees of freedom, Figure 6.But still, the highest force amplitudes occur at the faulty nodes.

Identification of Rubbing
In the cases 11 and 12 the gap S 2 between the disk m2 and the foundation ms has been reduced, so that the rotor touches the foundation with impacts and rotor to stator rub occurs.The criteria for identifying rub is NEWTON'S law of reaction: At the contact location, the force exerted on the rotor must coincide with the force exerted on the stator at any time.
The total equivalent force consists of the super- position of the equivalent forces caused by the unbalances at m2 and rn3 and those caused by rubbing between disk m2 and foundation ms.In Figure 7 it can be clearly distinguished between these two components.In  In all cases, the identification method deter- mined correctly the location of the fault by finding the largest equivalent force at the position where the fault was.However, substantial equivalent forces were also estimated for other nodes not affected by the faults.In other words, the equivalent forces were smeared over the entire system.The reason is the model error caused by the fact that due to the small number of measuring loca- tions accessible, only a few mode shapes can be taken into account to approximate the full vibra- tional state.This restriction, which can not be overcome in most practical applications, leads to larger deviations of the estimated fault severity than in the numeric experiments described before.
But, shifting all equivalent forces to the position where the largest one occurred, and adding them up yield in all cases acceptable results also for the fault severity.

CONCLUDING REMARKS
As the investigation showed, the identification method is able to localise the fault position, determine the type of fault and specify the fault parameters.
For reliable identification results as much information as possible should be available about the rotor-stator system.This means, that the mathematical model of the undamaged rotor system has to be as accurate as possible.To fulfill this requirement, the number of degrees of free- dom for the mathematical model can be large, because only simple mathematical operations have to be carried out with the linear model and time- consuming numerical integrations are avoided.
On the other hand, the success of the identifica- tion depends essentially on the number of meas- ured locations.

FIGURE 5
FIGURE 5 Test case 10: displacements (exact and reconstructed by modal expansion using 3 transducers).
shows the measured residual

TABLE II
case 1: estimated equivalent forces.
Table III the identifica- tion results are compared to those of test case 4, where no rub was present but the identification procedure searched for rub.
FIGURE 7 Test case 12: estimated equivalent forces for rubbing and unbalance.