MULTIPLE FAULT IDENTIFICATION METHOD IN THE FREQUENCY DOMAIN FOR ROTOR SYSTEMS

Fault identification in rotor systems has been studied by many authors, but the considered malfunction is one single fault only, generally an unbalance. Real machines can be affected by several different types of faults; moreover sometimes also two different faults may develop simultaneously. A model based method for identifying multiple faults acting simultaneously on a rotor system in the frequency domain is briefly described and its robustness with regards to measuring and modelling errors is evaluated, by means of numerical simulations performed on the models of two typical power plant machines: a steam turbogenerator and a gas turbogenerator.


Introduction
In recent years the topic of the fault identification, or diagnosis, in general technical processes has been the object of a huge number of papers in different fields of engineering science.Some of them deal with general procedures that can be applied on different physical systems and a quite complete review is reported in [1].If we focus our interest on fault identification in rotating machines under a mechanical point of view, also in this case several recent studies have appeared and indicated some preferred research path.
Generally speaking, an identification procedure can be performed by means of causality correlations of measurable symptoms to the faults.Two main approaches are commonly used that can roughly be divided in qualitative and quantitative methods.In qualitative approach, the symptoms can be defined using qualitative information, based on human operators' experience, which creates a knowledge base.Then faultsymptom matrices, fault-symptom trees, if-then rules or fuzzy logic classifications are used to indicate in a probabilistic approach the type, and sometimes also the size and the location of the most probable fault.Also artificial neural networks (ANN) can be used for creating the symptom-fault correlation.A recent contribution is given in [2]: an expert system can be built up in which different diagnostic reasoning strategies can be applied.The advantage of the qualitative approach is to furnish a probabilistic likelihood classification of the impeding fault type, whilst the quantitative approach considers initially a fault type at once.Some kind of rough qualitative classification is generally offered in standard or advanced machinery or process supervisory systems, but the amount of measured data and the availability of suitable models of the process or of the machines could be conveniently exploited, allowing more precise and reliable fault identification.
As regards the quantitative approach, it is normally a model based fault detection method and a reliable model of the system or of the process is used for creating the symptom-fault correlation, or the input-output relation.However this method has many different ways of applications, someone of them merge also qualitative and quantitative approaches.An example of this is given [3].In [4] a modal expansion of the frequency response function of the system, on both numerical model and experimental results is used to identify the unbalance distribution on a rotor.In [5,6] a model is presented in which equivalent loads due to the faults (rub-bing and unbalances) 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.The identification is then performed by least square fitting in the time domain.In [7] a model based identification in the frequency domain is employed to identify unbalance on a test-rig, while in [8] a model based procedure exploiting analytical redundancy is used to detect the faults in a single-shaft gas turbine.A rather common aspect of these effective analysis is their "ad hoc" application.A systematic approach has been introduced by the authors in [9] to identify several different types of faults and to discriminate among faults with similar harmonic components.This method has been experimentally validated on different test-rigs and real machines (see also [10][11][12][13]) for many types of faults, such as unbalances, rotor permanent bows, rotor rubs, coupling misalignments, cracks, journal ovalization and rotor stiffness asymmetries.
Here a generalization of the method, in order to take into account several simultaneous faults, is presented along with a numerical analysis about the robustness as regards measuring and modelling errors.
The method requires the definition of the models of the elements that compose the system, i.e. the rotor, the bearings and the foundation, along with the models of the faults, which can be represented by harmonic components of equivalent force or moment systems.The identification of the multiple faults is made by a least square fitting in the frequency domain, by means of the minimization of a multidimensional residue between the measured vibration in some measuring planes on the machine (usually, but not necessarily, the bearings) and the calculated vibrations due to the acting faults.

Fault modelling and identification
A complete description of the system modelling is far from the scope of the present paper and also many paper are available in literature about rotor modelling by means of finite elements.Similar considerations can be made about the modelling of the bearings and of the foundation.On the contrary, it is useful to briefly recall some concepts about the fault models and the identification procedure for multiple faults, while it is possible to refer to [14], in which a complete description is reported.
In order to understand the force models, it is necessary to introduce the reference systems used in the 2D f.e.beam model of the rotor.Each node of the model has 4 d.o.f.If we consider the two subsequent nodes, the jth and the j + 1th, they define the element jth, as shown in Fig. 1.
Indicating vectors and matrices with bold letters and scalar quantities with italic letters, the generalized displacements of node jth can conveniently be arranged in a vector x(j): which can be merged with the other ordered displacement vectors of the rotor nodes to form the complete displacement vector of the rotor: x = [. . .x j ϑ xj y j ϑ yj x j+1 (2) ϑ xj+1 y j+1 ϑ yj+1 . ..]T Thus, considering the degrees of freedom along whose a force, or a moment, acts, it has the following representation: in which the vector is the "localization" vector, F (k)  and M (k) the modulus, φ (k) the phase and nΩ the frequency.Among the different model based identification techniques, which can be roughly classified (see [1]) as parameter estimation, state estimation and parity equations, the first has been chosen as more suitable and successfully tested in many practical cases.Anyway, rather than the identification of the changes due to the fault in the system parameters, which influence generally the complete mass, stiffness and damping matrices of the system, the identification of the equivalent external force that produces the same effect is an easier task.This can be shown in mathematical terms starting from the standard matrix equation of the system without fault: where the right hand side (r.h.s.) is composed by the weight W (which is known) and by the original unbalance U and bow M u (which are unknown).If dM, dD and dK are the changes in mass, damping and stiffness matrices due to the arising fault, then Eq. ( 5) becomes: Assuming linearity in the system, the total vibration x t can be split in two superposed parts: the vibration vector x 1 due to the r.h.s.forces, both known and unknown, of Eq. ( 6), which satisfy again Eq. ( 5) and the vibration x due to the fault only.The component x may be obtained by calculating the vector differences of the actual vibrations x t minus the original vibrations x 1 measured, in the same operating conditions (rotation speed, flow rate, power, temperature, etc.): Combining Eq. ( 8) with Eq. ( 5), we obtain finally:  The r.h.s. of Eq. ( 9) can be considered as a system of equivalent external forces, which force the fault-free system to have the change in vibration defined by x that is due to the developing fault only: In Eq. ( 10) the system parameters are known and the fault identification is reduced to a force identification.If Eqs (3) and ( 4) are taken into account, it is evident that few elements of the unknown fault forcing vector are in reality different from zero.This fact makes the equivalent external forces approach more convenient than a plain parameter estimation approach.Due to the considered linearity, the method is convenient again also in the case of multiple simultaneous faults since they act on few d.o.f. of the system.Moreover, if also a steady-state situation is assumed, the harmonic balance criteria from Eq. ( 10) can be applied.This assumption can be justified even if experimental data of real machines are usually available from run-down transient.In big turbogenerators of power plants, due to the high inertia of the system, the transient occurs with slowly changing speed, so that actually it can be considered as a series of different steady state conditions.This allows to use these data in the frequency domain and the following equations, for each n-th harmonic component, are obtained: where the force vector F fn , which can be composed by several vectors Bow or rub Couple, 1x rev. [F has to be identified.This force vector could be function of Ω or not, depending on type of the fault.As regards the number of the harmonics to be considered, in field experience, generally not more than 3 components represent completely the periodical vibration time history.
Introducing the system dynamic stiffness matrix for the speed Ω and for the nth harmonic component, Eq. ( 11) can be rewritten as: Equation ( 13) corresponds actually to a measured vibration at a certain rotating speed.Considering now that, among the set of all the available measured vibrations at the q rotating speed, a subset corresponding to p rotating speeds is used for the fault identification, both speeds and vibrations that can be organized as a vectors: with p q and Eq. ( 13) becomes: If the system dynamic stiffness matrix [E(nΩ)] is inverted, as α n (Ω), the vibration amplitudes can be obtained as: The lines in Eq. ( 16) are rearranged, by partitioning the inverse of the system dynamic stiffness matrix, and omitting from α n and F fn the possible dependence on Ω for conciseness, in order to split the complex amplitude vector Ξ Bn corresponding to the d.o.f of the measured absolute vibrations in the measuring planes from the vector Ξ An of the remaining d.o.f. of the rotor system model: Using the first set of Eq. ( 17), it is possible to define, for each harmonic component, the vector differences δ n , between calculated vibrations Ξ Bn and measured vibrations Ξ Bmn : The problem of identifying the forces F fn that minimize the differences between the calculated and the measured vibrations, can now be solved by means of a least square method.In fact, in Eq. ( 18) the number of equations n m (number of measured d.o.f.) is lower than the number n d (number of d.o.f. of the complete system model) which is also the number of elements of F fn , but the vector of the equivalent fault forces has few non-zero elements (see Eqs (3) and ( 4)) even if the fault is not one only.Therefore a scalar difference, called "relative residue" is defined as the root of the ratio of the squared δ n , divided by the sum of the squared measured vibration amplitudes X Bmn : and the least square approach is used in order to find the solution (identified faults) that minimize the differences which are calculated for all the different rotating speeds which are taken into consideration.By means of the hypothesis of localisation of the fault, the residue is calculated for each possible node of application of each defect.This fact implies that, if we indicate with z k the abscissa along the rotor in correspondence to the kth fault among m faults, the relative residue in Eq. ( 19) is a surface in a R m+1 space, in other terms: Where the residue reaches its minimum, i.e. the minimum of the surface in Eq. (20), there is the most probable position of the fault.Figure 2 shows a sample of the residue surface for m = 2.It is worth to stress that both the location, the module and the phase of the fault are identified even if those fault characteristics could  be redundant for field use: often it is important only to localize the fault position along the shaft, such as in case of a rub, whilst in case of unbalances, coupling misalignments it is also important to know the modulus and the phase.Moreover other considerations should be made about the computational efficiency of the proposed method and the identification of more than two simultaneous faults, even if it is theoretically possible.
First of all, in actual machines it is very unusual the occurrence of more than two simultaneous faults.Secondly, the calculation time needed for the identification can become very large for more than two faults.This can make impossible an on-line identification.In a first approximation calculation time grows linearly with the product between the number p of the rotating speeds Ω and the number n d to the power of the number m of the faults.
In case of two faults, an useful graphical tool to quickly localize the faults along the rotor is the "residue map", a contour map of the residue surface with the rotor model along the x and y axes, as shown in Fig. 3.

Numerical simulations
The robustness of the proposed method as regards modelling and measuring errors has been tested on different types of machines, with different types of simul- taneous faults.Due to the limited space available for the paper, in the following some numerical cases of two faults are presented on two machine types.The first machine model is a 320 MW steam turbogenerator composed by a HP-IP turbine, a LP turbine and a generator.The overall length of the machine is about 28.7 m, the mass is about 135000 kg and seven oil-film bearings, of which those on HP-IP turbine are bi-lobed and the others lemon-shaped, sup-port the group.The model of the rotor is composed by 136 elements (Fig. 4), the 1st critical speed of HP-IP turbine is about 1300 rpm and that of LP turbine about 2700 rpm.The bearing stiffness and damping coefficients are available for rotating speeds in the range 500-3000 rpm.The foundation is modelled by seven 2 d.o.f.pedestals (mass, spring and damper systems) with constant mass, stiffness and damping coefficients.
The second machine is an 125 MW turbogas gener- ator composed by a gas turbine and a generator.The overall length of the machine is about 20.3 m and the mass is about 82000 kg.Two tilting pad bearings support the turbine section, while two oil-film plain circular bearings support the generator.The model of the rotor is composed by 124 elements (Fig. 5).Also in this case, the supporting structure is modelled by four 2 d.o.f.pedestal (mass, spring and damper systems) with constant mass, stiffness and damping coefficients.
In order to limit the number of the test cases presented, only two types of faults are considered: the unbalance and the local bow.The two fault models  L ] is the localization vector, A (k) (Ω) the complex amplitude of the kth equivalent force system, (mr) (k) the product of mass and distance from the rotating axis of the kth unbalance and M (k) the modulus of one of the couple moments.A general discussion on fault causes and modelling can be found in Bachschmid and Pennacchi [9].
Seven different cases of two simultaneous faults have been analysed on the steam turbogenerator model and three on the turbogas generator, with different combinations of types and location of faults.
The two complete models were used to generate several reference data sets of the system responses to the different simultaneous errors.
Then "corrupted" models are defined, introducing errors in the bearing models because the oil film linearized stiffness and damping coefficients are generally believed to be affected by highest modelling errors with respect to the other components of the system.The tests were thereby carried out modifying the characteristics of stiffness and damping of the bearings.Since the errors in stiffness and damping coefficients of the bearings influence also shaft displacements in the bearings, which are the measured quantity in actual machines, the introduced errors also simulate measuring errors.Two kinds of errors were introduced, one which will be referred to as "random" and the other which will be referred to as "systematic".An example is reported in Fig. 6.
In the "random" error case, all coefficients of all bearings are modified adding a random quantity between ±30% of the value of the corresponding coefficient for each rotating speed.The obtained values are then interpolated with a spline in order to avoid discontinuities in the slope.If the "random" error is intended as a measuring error, it can represent for instance noise on the sensors.
The "systematic" error is instead obtained increasing or decreasing bearing coefficients by the same amount of between ±30% all over the speed range.In this case the error is systematic on the single characteristic but it is randomly spread over different bearings.If the "systematic" error is intended as a measuring error, it can represent for instance an offset on the sensor measure.
The changes in bearing characteristics obviously lead to a different response of the models to imposed faults.An example is reported in Fig. 7 for the turbogenerator.Note that the resonance frequency shift cannot be considered as negligible.
The results of the sensitivity analysis are reported in Table 2 to Table 11 for two simultaneous faults on both turbogenerator and turbogas.In the columns labelled with ∆l the percentage error in the identification of the position of the fault is reported, which is calculated comparing the abscissa of the position of the fault with the length of the f.e.beams that compose the model of the single shaft in the shaft line.Columns labelled with ∆m report the percentage error on the module of the fault in the case of the unbalance, while in the case of the local bow the column is labelled ∆φ and the error is calculated on the relative rotation of the nodes that the identified bending moment causes on the element where the fault is applied in the reference case.This allows to consider the different stiffness of the element where the fault is identified.Finally, columns labelled with ∆∠ report the percentage error on the phase in degrees referred to a phase of 180 • .Errors are not reported for the first faults in Tables 6 and 10, since they are located in a wrong position.Residue maps, for conciseness, are reported for the case in Table 2 only.
From the analysis of the results reported in Table 2 to Table 11 it is possible to draw some general conclusions: i) a quite good agreement between the imposed and identified defects is obtained even when errors are introduced in the model definition.The value of the residue is generally rather low.In particular, systematic errors on the bearing coefficients cause the identification to be less effective; ii) the identification procedure for two simultaneous faults can fail in identifying one of the faults in presence of errors if the faults are of the same kind and very close each other.This occurs in the case of unbalances in Table 4, where one fault is identified correctly in position, with a module and phase that approximately represent the vector sum of the reference faults, while the other is located in a wrong position with a little module.In case of local bows, as in Tables 6  and 10, only one of the faults is identified with a good accuracy; iii) the identification procedure for two simultaneous faults proves to be effective in identifying the faults in presence of errors if they are of different kind, even if they are very close each other (see the case in Table 8).

Conclusions
A model based method for the identification of simultaneous faults in rotor systems has been presented in this paper and all the analytic details are explained.Since in model based identification one of the most common cause of errors in the identified faults is due to of lack of accuracy in the fully assembled machine and to noise in the experimental data, in order to test the method effectiveness in presence of modelling or measuring errors, a sensitivity analysis has been performed on different models of real machines, where the errors were introduced by means of "corrupted" models of the bearings.The obtained results have shown that the proposed method is rather effective in identifying the faults also in presence of errors and it is suitable to be applied on full size real machines.

Fig. 2 .
Fig. 2. Residue surface in case of simultaneous identification of two faults.The location of the faults is in the rotor nodes corresponding to the minimum of the surface.

Fig. 7 .
Fig. 7. Comparison of turbogenerator FRF in bearing #3 in reference, random and systematic error case.

Table 1
Models of the considered faults

Table 2
Turbogenerator with an unbalance on HP-IP turbine and an unbalance on LP turbine ∠∆∠Node ∆l Module (kgm) ∆m ∠

Table 3
Turbogenerator with a local bow on HP-IP turbine and a local bow on LP turbine

Table 4
Turbogenerator with two unbalances on HP-IP turbine

Table 6
Turbogenerator with two local bows on HP-IP turbine

Table 7
Turbogenerator with an unbalance on LP turbine and a local bow on HP-IP turbine

Table 8
Turbogenerator with an unbalance on LP turbine and a local bow on LP turbine

Table 9
Turbogas with an unbalance on turbine and an unbalance on generator

Table 10 Turbogas
with a local bow on turbine and a local bow on generator

Table 11
Turbogas with an unbalance on turbine and a local bow on generator

Table 2
case.considered are reported in Table1, where [F