A Super-Harmonic Feature Based Updating Method for Crack Identification in Rotors Using a Kriging Surrogate Model

Abstract

Dynamic model updating based on finite element method (FEM) has been widely investigated for structural damage identification, especially for static structures. Despite the substantial advances in this method, the key issue still needs to be addressed to boost its efficiency in practical applications. This paper introduces the updating idea into crack identification for rotating rotors, which has been rarely addressed in the literature. To address the problem, a novel Kriging surrogate model-based FEM updating method is proposed for the breathing crack identification of rotors by using the super-harmonic nonlinear characteristics. In this method, the breathing crack induced nonlinear characteristics from two locations of the rotors are harnessed instead of the traditional linear damage features for more sensitive and accurate breathing crack identification. Moreover, a FEM of a two-disc rotor-bearing system with a response-dependent breathing crack is established, which is partly validated by experiments. In addition, the associated breathing crack induced nonlinear characteristics are investigated and used to construct the objective function of Kriging surrogate model. Finally, the feasibility and the effectiveness of the proposed method are verified by numerical experiments with Gaussian white noise contamination. Results demonstrate that the proposed method is effective, accurate, and robust for breathing crack identification in rotors and is promising for practical engineering applications.

1. Introduction

Rotors are one of the most important components of rotating machines, which are widely used in many engineering fields. Cracks in rotors are a very dangerous kind of fault that may lead to a sudden and catastrophic failure of a rotating machine. Therefore it is of vital significance to monitor and detect a crack, in order to reduce maintenance cost and avoid failure of a rotating machine. In addition, crack parameter identification can provide a foundation for life prediction of a cracked rotor.

A crack in a rotor introduces an additional local flexibility and reduces the overall stiffness. Usually, a fatigue crack in a rotor is a breathing type, which makes the stiffness of such a rotor changing with time and causes the system to become nonlinear. Because of these, a cracked rotor possesses complex dynamic characteristics. Much work on the dynamic behavior of cracked rotors has been performed to identify cracks in rotors, and a number of reviews on this topic [1,2,3,4,5,6,7] have been published. From the literature, vibration-based crack identification methods can be broadly classified as model-based methods, signal-based methods, modal parameters, and their derivatives-based methods and pattern recognition-based methods.

The approach based on equivalent crack forces is well known as a model-based method for crack identification in rotors, considering the effects of cracks as equivalent forces applied in the intact systems. It has been utilized by Sekhar [8] to identify the crack locations and depths of a rotor with two cracks. Pennacchi et al. [9] further used the approach for cracks identification of a relatively large experiment rig. There are also some other model-based methods. Söffker et al. [10] presented a model-based method based on a proportional-integral observer to identify a crack in an operating rotor. Singh et al. [11] proposed a model-based method to identify a switching crack in a four-degrees-of-freedom Jeffcott rotor utilizing the vibration signal in conjunction with the controller current of the active magnetic bearing. Model-based methods were adopted by Saridakis et al. [12], and Xiang et al. [13] to minimize the difference between real outputs and model outputs to determine the location and depth of a crack in a rotor-bearing system with genetic algorithm (GA). Extended Finite Element methods were developed for cracks and flaws detection and quantification by He et al. [14] and Waisman et al. [15]. Cavalini Jr et al. [16] took nonlinear combinational resonances induced by cracks under external diagnostic forces at certain frequencies as an objective function of a differential evolution optimization method to determine the crack location and depth by minimizing the difference between the measured and modelled rotor systems.

Short Time Fourier Transform (STFT) [17,18], Wavelet Transform (WT) [19,20,21] and Hilbert-Huang Transform (HHT) [22,23] as signal-based methods, were used to detect a crack from acceleration response of a rotor by Chandra and Sekhar [24]. Approximated Entropy based on information theory was used by Sampaio and Nicoletti [25] to detect rotor cracks considering the interference of misalignment. Full-spectrum method was used to detect a breathing crack in a rotor by Shravankumar and Tiwari [26]. Yang et al. [27] utilized the modulation between crack breathing and vibration response, by using squeezing time-frequency transform (STFT) to extract the periodic oscillation in instantaneous frequency for breathing crack detection. Xie et al. [28] proposed a phase-based spectrum analysis method to extract weak harmonics induced by cracks, to help realize crack detection.

For modal parameters and their derivatives based methods, empirical mode decomposition (EMD) and Laplace wavelet finite element combined method was proposed by Dong et al. [29], to obtain the accurate modal parameters to identify a crack in a rotor. Iterative signal extraction method was utilized to reduce the noise interference in modal parameters identification by Liu et al. [30]. Changes in resonant and anti-resonant frequencies were used by Rubio et al. [31] to detect crack locations in a two-crack torsional shaft. Rahman et al. [32] utilized the changes in phase angle of frequency response function to identify the location and depth of an open crack in a rotor. Seo et al. [33] proposed a method for open crack localization of a shaft by comparing the map of the modal constants of the reverse directional frequency response functions with the reference map of the un-cracked model.

The support vector machine (SVM) is a common pattern recognition-based method with many applications [34], and was used to realize crack identification with features extracted by wavelet transform in [10]. A notch in a stationery rotor was identified by the first four natural frequencies using artificial neural network (ANN) by Zapico-Valle et al. [35].

As a summary for the literature review, there are more crack detection methods than crack identification methods for rotors, and among the crack detection methods, signal based methods are in the majority. These signal based methods can realize crack detection conveniently, while they are difficult in crack parameters identification. For modal parameters based methods, the crack will definitely affect the modal parameters, but the modal parameters are not sensitive to shallow cracks. And for the pattern recognition-based methods, a large amount of crack samples are required to train the models, which limits the application of such kind of methods. For model based methods, they could perform well in crack identification, only if the models were accurate enough. Therefore, the accuracy of a modelling method and efficiency of an identification approach should be the focuses.

Based on the above short review, there is a trend to transform the crack identification problem into an optimization problem, and through updating iterations to find the crack parameters minimizing the difference between measured features and calculated features, which can be classified as model updating method [36]. The model updating method has been widely used in structural damage identification [37,38,39,40], and what matters the most are the model construction and identification efficiency. For the static beam or plate, an accurate model is far more easily to obtain, which will be not the case for rotating rotors, especially when a breathing crack is there. Therefore, how to establish an efficient and reasonable updating model is quite important for crack identification in rotors.

Super-harmonic features are good indicators of a crack in a rotor, which emerge since the presence of a crack leads to asymmetric behaviour and a breathing crack introduces nonlinearity. Although super-harmonic features can also be caused by other faults, such as coupling misalignment which also induces super-harmonics in vibration response of a rotor [41], it is much easier to identify and correct [42]. So in this work only the crack is considered. Energy variation of the 3× component (where X represents the frequency corresponding to the rotating speed) extracted based on Wavelet Packet Transform [43,44,45] was used to quantify the crack depth of a crack with known crack location in a rotor by Gómez et al. [46], and the analytical Jeffcott rotor model and the corresponding experimental rotor with a saw-cut crack were studied to validate the method. EMD method [47,48,49] was applied to steady-state responses generated from a Jeffcott rotor to extract the 3× and 2× components in the neighborhood of 1/3 and 1/2 of the critical rotating speed by Guo et al. [50], results showed that the variation of averaged amplitudes of super-harmonic components provided clear and robust signatures of early cracks in rotating rotors. However, from the literature, few researches have been carried out to identify both the location and depth of a crack in a rotor with these super-harmonic features.

For a cracked rotor system, various crack parameters correspond to various amplitudes of super-harmonic components. For crack localization there must be some spatial information in the corresponding features, therefore, if there is just one measurement point, it is not possible to distinguish the cracks at different locations with varied depths. Hence, more measurement points are needed. In view of the difficulties in taking measurement in rotating rotors, only two points are measured for crack identification. It is believed that two measurement points are enough, only if accurate model information is available. Hence the key to identify a crack in rotors with two measurement points is to obtain a crack model as accurate as possible (it should be noted that the model here is not limited to a structural model but also a mapping between any two variables). However, to obtain an accurate model between crack parameters and features, many samples are required, which will be quite costly and the efficiency will be quite low. Therefore, how to obtain the accurate relationship with much fewer samples and to boost its efficiency are challenges.

A Kriging surrogate model is a statistics-based interpolation method [51]. It has been increasingly utilized in optimization problems in which large computing resources are required [52], and it has also been used in parameter identification [53] and crack identification in plates [54]. However, few applications of Kriging surrogate model have been reported in crack identification in rotors.

In this work, a novel Kriging-based FEM updating method for crack identification for rotors with a breathing crack is proposed using nonlinear super-harmonic features from two “measurement” points, integrating with an improved intelligent optimization algorithm of Particle Swarm Optimization (PSO). The proposed method is effective, accurate and robust for breathing crack identification in rotors and can identify all the crack parameters at the same time, moreover, it is promising for practical engineering applications.

The structure of rest of the paper is organized as follows. Section 2 establishes a model of a two-disc rotor-bearing system with a breathing crack and validates this model by experiment and previous published work. Section 3 carries out the sensitivity analysis of crack parameters to super-harmonic characteristics. Section 4 constructs the Kriging surrogate model that relates crack parameters and super-harmonic features. Section 5 describes the detail of the proposed crack identification method. Section 6 presents the numerical simulations and crack identification results. Finally, Section 7 concludes the paper.

2. Modelling and Validation of a Two-Disc Rotor with a Breathing Crack

In order to analyze the characteristics of a two-disc rotor with a breathing crack, the finite element method is adopted to establish the model, in which the breathing crack is simulated by the CCLP (crack closure line position) breathing model proposed by Darpe et al. [55]. Since the method has been elaborated in the previous work by the authors [56], here, for the conciseness of the paper, the modelling procedure will just be briefly introduced. Further, validation of the established model is carried out by modal experiment and comparison with previous published work.

Two-node Timoshenko beam elements with six degrees-of-freedom (three translations and three rotations) per node are used to model the rotor-bearing system. Figure 1 shows the two-disc cracked circular rotor-bearing system, the cracked shaft element, and the finite beam element model. To establish the model of the cracked rotor, the key point is to model the crack appropriately and calculate the stiffness matrix of the crack element that can be calculated by the strain energy release rate approach [56]. After that, through assembling the cracked and un-cracked elements, the finite element model of the rotor can be obtained. Though the method is universal for various crack angles, only the transverse crack is considered in this paper in order to simplify the problem.

For the rotor in Figure 1a, two discs are considered as rigid bodies, which have three translational and three rotational inertias, and they are added to the mass matrix elements at the corresponding degrees-of-freedom. The gyroscopic effect of the two discs is also included. The bearings are simplified as stiffness and damping systems, both constrain 2 lateral degrees-of-freedom. The rotating speed of the rotor is $\Omega$. By assembling the system matrix of the cracked element and un-cracked elements, the finite element beam model shown in Figure 1c can be established.

Denote $q_i$as displacement vector of node i having 6 degrees-of-freedom:

$$q_i={\left\{x_i ,y_i ,z_i ,{\theta }_{xi} ,{\theta }_{yi} ,{\theta }_{zi}\right\}}^{\mathrm{T}}$$

(1)

The equations of motion in the stationary coordinate system can be written as:

$$M\ddot{q}+\left(C+\Omega G_{\mathrm{y}}\right)\dot{q}+K\left(t\right)q=F_{\mathrm{e}}+F_{\mathrm{g}}+F_{\mathrm{ex}}$$

(2)

$$q={\left\{q_1,q_2,\dots q_i,\dots ,q_n\right\}}^{\mathrm{T}}$$

(3)

where $M$ is the system mass matrix;$C=aM+bK$ is system damping matrix considering the Rayleigh damping; $G_{\mathrm{y}}$ is system gyroscopic matrix; $K$ is system stiffness matrix that will be updated as the crack breathes; n is the total node number (in the following numerical calculation, and the rotor is evenly divided into 60 elements, hence n is equal to 61); $F_{\mathrm{e}}$ is the excitation due to eccentricity of discs; $F_{\mathrm{g}}$ is excitation due to the gravitational force; and $F_{\mathrm{ex}}$ is external excitation during operation.

Parameters of the rotor-bearing system are shown in

Table 1. Parameters of the two-disc rotor.

Feature	Value (Units)	Feature	Value (Units)
Shaft length (l0)	0.6 m	Mass of disc 1	0.759 kg
l1	0.46 m	polar moment of inertia of disc 1	5.758 × 10−4 kg·m2
l2	0.07 m	diametrical moment of inertia of disc 1	3.177 × 10−4 kg·m2
l3	0.08 m	Mass of disc 2	0.770 kg
l4	0.1 m	polar moment of inertia of disc 2	5.843 × 10−4 kg·m2
Shaft diameter (D)	0.01 m	diametrical moment of inertia of disc 2	3.232 × 10−4 kg·m2
density (40Cr)	7.87 × 103 kg·m−3	Gravitational acceleration	9.8 m/s2
Yang’s modulus	2.11 × 1011 Pa	Rayleigh damping coefficient (a)	0.68
Poisson’s ratio	0.277	Rayleigh damping coefficient (b)	2.8 × 10−5

, where $a$ and $b$ are calculated by assuming modal damping ratios of the first two modes being 0.005 and 0.01.

To validate the established finite beam element model which called beam model below, on one hand, the intact rotor is established by ABAQUS software, on the other hand, modal experiment by an impact hammer is carried out. The tested rotor is shown in Figure 2, in which two orthogonal eddy current transducers are used to obtain the response excited by the hammer.

The first three natural frequencies are obtained and listed in

Table 2. Comparison of the natural frequencies from numerical calculation and experiment.

Order	Beam Model	ABAQUS	Measured	Error between Beam Model and ABAQUS	Error between Beam Model and Measured Results
1st vertical bending	42.1 Hz	40.9 Hz	42.5 Hz	1.2 Hz (2.85%)	0.4 Hz (0.95%)
1st lateral bending	42.1 Hz	41.0 Hz	42.9 Hz	1.1 Hz (2.61%)	0.8 Hz (1.90%)
2nd vertical bending	107.9 Hz	112.4 Hz	107.4 Hz	4.5 Hz (4.17%)	0.5 Hz (0.47%)
2nd lateral bending	107.9 Hz	112.4 Hz	108.7 Hz	4.5 Hz (4.17%)	0.8 Hz (0.74%)
1st torsion	158.8 Hz	160.2 Hz	-	1.4 Hz (0.88%)	-
3rd vertical bending	349.6 Hz	358.0 Hz	349.6 Hz	8.4 Hz (2.40%)	0.0 Hz (0.00%)
3rd lateral bending	349.6 Hz	357.2 Hz	350.7 Hz	7.6 Hz (2.17%)	1.1 Hz (0.31%)

. It should be noted that the results from the beam model and ABAQUS software were acquired after a few rounds of model updating iterations by adjusting bearing parameters to make the errors of the first three mode frequencies between ABAQUS and measured ones as small as possible. The corresponding bearing stiffness was found to be 960,000 Nm⁻¹, and bearing damping is 100 Nsm⁻¹.

As one can see from Table 2, the maximum error between the beam model and ABAQUS was 4.17%, and the error between the beam model and measured results was even less than 2%. The above results show that the established beam model is validated and can be used for further dynamic research. It will be used to study the dynamic behaviors of rotors with cracks, provided that the crack model is accurate.

Based on the validated model, the Newmark integration method [57] was adopted to compute the rotor unbalance response after assembling a crack into the rotor by the afore-mentioned crack breathing model. To illustrate the dynamic behavior of the cracked two-disc rotor, a typical waterfall diagram of the cracked rotor (with a crack depth of 0.35D, located in the 30th element) under different rotating speeds below the first critical speed of the rotor is shown in Figure 3.

As one can see from Figure 3, 2×, and 3× super-harmonic components appear and do not exist in an un-cracked rotor, because of crack breathing, which is a nonlinear feature. When the rotating speed reaches 1/3 and 1/2 of the critical speed, the amplitudes of 3× and 2× components reach their maximum, respectively. The evolution process of whirl orbits when the rotating speed is at the neighborhood of 1/3 and 1/2 of the critical speed is shown in Figure 4.

As for the case near 1/3 critical speed in Figure 4, the orbits contain two inner loops corresponding to 3× component, and the loops first grow and become dominant, then after passing through 1/3 critical speed, the loops disappear gradually. In contrast, there is only one inner loop in the orbit corresponding to 2× component. It almost shares the same evolution trend as the 3× component. When the rotor passes through 1/3 or 1/2 sub-resonance region, the orbit phase rotates nearly 90°. The behaviors agree with the experiment results in [50] and numerical results in [58], which can verify the correctness of the model and its solution in this paper.

As the model has been established and validated, to study the feasibility and effectiveness of the identification method in theory without interference of other uncertainties in real machines, the crack identification method will be developed and verified based on the validated numerical model of the cracked rotor.

3. Investigation of the Effects of Crack Parameters on Super-Harmonic Characteristics

After validating the model of the cracked rotor, to select the most sensitive features to establish the Kriging surrogate model for crack identification, the effects of crack parameters on the super-harmonic components at 1/3 and 1/2 critical speed of the rotor were investigated.

Responses at node 12 were obtained by the established finite element model. The crack location was made to vary from element 10 to 52, which covers the part of the shaft between the two bearings, and crack depth ranges from 0 to 0.45D with step size of 0.05D. The numerical simulation results are illustrated in Figure 5 below.

From Figure 5a,c above, one can see that the amplitudes of 2× component at 1/3 and 1/2 critical speeds increased monotonously with the increase of crack depths, and they were also distinct between different crack locations. In contrast, the amplitudes of 3× component at 1/3 and 1/2 critical speeds were not monotonous with the increase of crack depths, which can be seen from Figure 5b,d. If the features are monotonous with a crack parameter, then it will be easier to distinct a similar feature value, which will be helpful to reduce the identification errors. Therefore, in view of crack quantification, 2× components at 1/3 and 1/2 critical speed are better features for crack parameters identification than the 3× ones.

As the purpose of this paper is to propose a crack identification method for rotors and investigate the feasibility of the method based on a Kriging surrogate model, to simplify the problem, and reduce the demand of input features, only 2× component at 1/3 critical speed was utilized for crack identification. However, it is difficult to identify a crack from the rotor response at just one measure point, because such a solution may not be unique. To tackle the problem, more measurement points are needed. Considering the difficulties of setting up sensors for rotating rotors, two measurement points were used to determine the crack parameters. Two measurement points should be enough, when the model that maps crack parameters to amplitudes of 2× components from two measurement points is accurate enough. In view of the accessibility of sensors in real applications, the two measurement points A and B were chosen to be node 12 and node 50 respectively in the numerical experiment, which were near the two bearings. Here, it can be said that it is possible to identify the crack in rotors by using the 2× components from two measurement points. However, the cost of experiment or computation is another important issue which will affect the crack identification efficiency. As the nonlinear breathing crack modelled in this work was response-dependent, computation cost of the simplified beam model was still heavy. Therefore, it is not efficient to utilize the beam model directly combined with an optimization method, as the computationally expensive FE models will be called repeatedly during the optimization process. To reduce the cost of experiment or computation in real applications and to avoid the repeated analysis of computationally expensive FE models during the optimization process, a Kriging surrogate model was adopted to establish the accurate relationship with fewer samples, which will be elaborated in the following sections.

4. Construction of the Crack Kriging Surrogate Model

A Kriging surrogate model is a statistics-based interpolation method and has been elaborated by Forrester et al. [59]. In this part, the crack Kriging surrogate model will be established for crack identification using the sensitive super-harmonic features. From the sensitive analysis above, the 2× components at 1/3 critical speed obtained from measurement points A and B will be used for the crack Kriging model construction. The relationship between crack parameters and 2× amplitudes at 1/3 critical speed can be written as: $Y=f\left(l,d\right)=f\left(x\right)$, where, $l$ and $d$ are vectors for crack locations and crack depths respectively and $x$ is the crack parameters set. Specifically, the relationships for measurement points A and B are $Y^1=f^1\left(l,d\right)=f^1\left(x\right),Y^2=f^2\left(l,d\right)=f^2\left(x\right)$. Since two measurement points are required, two surrogate models corresponding to measurement points A and B should be constructed, respectively. Details for the crack Kriging surrogate model construction are as follows.

For measurement point A, given n samples of crack parameters and the corresponding 2× amplitudes:

$$x=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_1\\ x_2\end{array}\\ x_3\\ ⋮\end{array}\\ x_n\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{l}_1\\ l_2\end{array}\\ l_3\\ ⋮\end{array}\\ l_n\end{array} \begin{array}{c}\begin{array}{c}\begin{array}{c}{d}_1\\ d_2\end{array}\\ d_3\\ ⋮\end{array}\\ d_n\end{array}\right],Y^1=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{y}_1\\ y_2\end{array}\\ y_3\\ ⋮\end{array}\\ y_n\end{array}\right]$$

(4)

their relationship can be written as a polynomial regression model $F\left(\mathsf{\beta },x\right)$ and a normally distributed random process $ϵ\left(x_i\right)$ with mean value $\mu$ and standard deviation $\sigma$:

$$y_i\left(x_i\right)=F\left(\mathsf{\beta },x\right)+ϵ\left(x_i\right)$$

(5)

The correlation function between the i-th and j-th samples is $R\left[ϵ\left(x_i\right),ϵ\left(x_j\right)\right]$, which is the Kriging function in this paper. For a k-variable design space, a Kriging function can be expressed as:

$$R\left[ϵ\left(x_i\right),ϵ\left(x_j\right)\right]=\exp \left(-{\displaystyle \sum }_{h=1}^k{\theta }_h{\left|x_i^h-x_j^h\right|}^{p_h}\right)$$

(6)

where ${\theta }_h$ indicates the sensitivity of the variable h in a sample (as only the location and depth of a crack are involved, therefore h = 1, 2.), and $p_h$ determines the smoothness of the correlation function.

The likelihood to eliminate the modelling error can be expressed as a Probability Density Function:

$$L\left(y_1, y_2,\cdots , y_n|\mu ,\sigma \right)=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{n/2}{\left|R\right|}^{1/2}}\exp \left(-\frac{{\left(Y^1-\mathbf{1}\mu \right)}^{\mathrm{T}}{R}^{-1}\left(Y^1-\mathbf{1}\mu \right)}{2{\sigma }^2}\right)$$

(7)

and its logarithmic form is:

$$\ln \left(L\right)=-\frac{n}{2}\ln \left(2\pi \right)-\frac{n}{2}\ln \left({\sigma }^2\right)-\frac{1}{2}\ln \left(\left|R\right|\right)-\frac{{\left(Y^1-\mathbf{1}\mu \right)}^{\mathrm{T}}{R}^{-1}\left(Y^1-\mathbf{1}\mu \right)}{2{\sigma }^2}$$

(8)

where $R$ is the correlation matrix of observed samples and $|·|$ is the operator of determinant:

$$R=\left(\begin{array}{ccc}R\left[ϵ\left(x_1\right),ϵ\left(x_1\right)\right]& \cdots & R\left[ϵ\left(x_1\right),ϵ\left(x_n\right)\right]\\ ⋮& \ddots & ⋮\\ R\left[ϵ\left(x_n\right),ϵ\left(x_1\right)\right] & \cdots & R\left[ϵ\left(x_n\right),ϵ\left(x_n\right)\right]\end{array}\right)$$

(9)

Hence, an essential step to construct a surrogate model is to maximize the likelihood. By setting the derivatives of Equation (8) with respect to $Y^1$ and ${\sigma }^2$ to zero, one can obtain the most likely values of $\widehat{\mu }$ and $\widehat{\sigma }$:

$$\widehat{\mathsf{\mu }}=\frac{1^{\mathrm{T}}{\mathrm{R}}^{-1}{\mathrm{Y}}^1}{1^{\mathrm{T}}{\mathrm{R}}^{-1}1}$$

(10)

$${\widehat{\sigma }}^2=\frac{{\left(Y^1-\mathbf{1}\widehat{\mu }\right)}^{\mathrm{T}}{R}^{-1}\left(Y^1-\mathbf{1}\widehat{\mu }\right)}{n}$$

(11)

Then the concentrated likelihood function can be obtained by substituting $\widehat{\mu }$ and $\widehat{\sigma }$ into Equation (8) and removing the constant term of $-\frac{n}{2}\ln \left(2\pi \right)$:

$$\ln \left(L\right)\approx -\frac{n}{2}\ln ({\widehat{\sigma }}^2)-\frac{1}{2}\ln \left|R\right|$$

(12)

The value of this function depends on the unknown parameters $\mathsf{\theta }$ and $p$corresponding to the parameters of the Kriging function in Equation (6), called sensitivity vector of crack parameters and smoothness coefficient vector of correlation function respectively, therefore the problem of constructing the surrogate model is transformed into finding the $\mathsf{\theta }$ and $p$, which maximize the likelihood. The global search method, such as GA, can be used.

Once $\widehat{\mu }$, $\widehat{\sigma }$, $\mathsf{\theta }$ and $p$ are obtained, the 2× amplitude from measurement point A at an untried point $x^{*}$ can be predicted based on,

$$\widehat{y}\left(x^{*}\right)=\widehat{\mu }+r^{\mathrm{T}}{R}^{-1}\left(Y^1-\mathbf{1}\widehat{\mu }\right)$$

(13)

where $r$ is the correlation vector between the untried point $x^{*}$ and the observed data $x$:

$$r=\left(\begin{array}{c}R\left[ϵ\left(x_1\right),ϵ\left(x^{*}\right)\right]\\ ⋮\\ R\left[ϵ\left(x_n\right),ϵ\left(x^{*}\right)\right]\end{array}\right)$$

(14)

Better performance of a Kriging surrogate model requires the fewer prediction error in the entire design space, and the prediction error can be described by [52],

$$\mathrm{MSE}\left({\mathbf{x}}^{*}\right)={\widehat{\sigma }}^2\left[1-{\mathbf{r}}^{\mathrm{T}}{\mathbf{R}}^{-1}\mathbf{r}+\frac{{\left(1-{\mathbf{1}}^{\mathrm{T}}{\mathbf{R}}^{-1}\mathbf{r}\right)}^2}{{\mathbf{1}}^{\mathrm{T}}{\mathbf{R}}^{-1}\mathbf{r}}\right]$$

(15)

To evaluate the prediction error more reasonable, the relative MSE defined by the ratio between MSE and the corresponding amplitude of 2× component is adopted:

$$\mathrm{RMSE}\left({\mathbf{x}}^{*}\right)=\mathrm{MSE}\left({\mathbf{x}}^{*}\right)/\mathbf{Y}\left({\mathbf{x}}^{*}\right)$$

(16)

Here, the Kriging model for measurement point A is finished, and the same process is for B.

5. Crack Identification Method Based on the Kriging Surrogate Model

Crack identification is a typical inverse problem. The main idea of the proposed method is to transform the inverse problem into an optimization problem, based on the surrogate model established between crack parameters $x$ and super-harmonic components $Y$ from two measurement points. Therefore, crack identification can be stated as:

$$\{\begin{array}{c}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d} x^{*}\\ \min \Vert Y_{*}\left(x^{*}\right)-Y_{\mathrm{Target}}\Vert \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o} \mathbf{l}\mathbf{b}<x^{*}<\mathbf{u}\mathbf{b}\end{array}$$

(17)

where, $Y_{*}$ is the vector of predicted super-harmonic components from two measurement points at the new set of crack parameters $x^{*}$; $Y_{\mathrm{Target}}$ is the vector of amplitudes of the measured super-harmonic components; and $lb$ and $ub$ are vectors of the lower and upper bounds for $x^{*}$.

As the surrogate model established by the initial samples will not guarantee the accuracy in the entire parameter range, the surrogate model should be updated for the more accurate parameters identification. Therefore, to estimate the accuracy of the surrogate model and identified parameters the following two equations are used:

$$C_1=1-\frac{\sqrt{{\sum }_i{\left(Y_{*}^i\left(x^{*}\right)-Y_{\mathrm{FE}}^i\left(x^{*}\right)\right)}^2}}{\sqrt{{\sum }_i{\left(Y_{\mathrm{FE}}^i\left(x^{*}\right)\right)}^2}}$$

(18)

$$C_2=1-\frac{\sqrt{{\sum }_i{\left(Y_{*}^i\left(x^{*}\right)-Y_{\mathrm{Target}}^i\right)}^2}}{\sqrt{{\sum }_i{\left(Y_{\mathrm{Target}}^i\right)}^2}}$$

(19)

where, $C_1$ is the first criterion to evaluate the accuracy of the current surrogate model; and $C_2$ is the second criterion to evaluate the matching degree of the identification results by the current surrogate model. The two criterions are together used to decide whether the current surrogate model needs to be updated or not. $Y_{\mathrm{FE}}^i$ is the ith element of $Y_{\mathrm{FE}}$ with crack parameters $x^{*}$, where $Y_{\mathrm{FE}}$ is the super-harmonic features of crack parameters $x^{*}$calculated by the finite element model (i.e., the beam model). When the surrogate model and identified parameters do not satisfy the pre-set criterions, the surrogate model should be updated by adding the identified optimal solution in the current iteration, until the criterions are satisfied. Figure 6 illustrates the flowchart of the identification process. It should be noted that the surrogate model is composed of the sub-surrogate models from points A and B, and during the identification process, surrogate models from points A and B are operating individually to predict the features of untried parameters, and the two models will be merged when evaluating the fitness and the criterions.

To search for the optimal crack parameters based on the surrogate model, global optimization methods, such as GA and Particle Swarm Optimization (PSO), can be used. Compared with GA, PSO utilizes the real coding in optimization and has many merits, such as good robustness and fast convergence rate [60]. In order to enhance the efficiency of PSO, an improved approach based on multi-point adding is proposed. Its main idea is to predict possible solutions in current iteration that the final optimal solution will be nearby, then to add the possible solutions and the optimal solution in the current iteration into the initial swarms for next iteration. The possibility is measured by the distance between the possible solution and the optimal solution in the current iteration. The nearer the distance is, the larger the possibility is. All the added points include three points, the optimal solution, the nearest and the second nearest solutions in the current iteration.

6. Numerical Experiments

In order to investigate the effectiveness and efficiency of the proposed method, numerical experiments were carried out based on the established cracked rotor finite element model. A sampling plan, as shown in Figure 7, with 20 initial samples generated by Latin Hypercube Sampling (LHS) [52] with crack locations from the 10th element to the 52th element and crack depths in [0, 0.45D] was utilized to construct the surrogate model and identify the crack parameters. Cracks in 10 cases (shown in Figure 7) which cover the typical situations generated by LHS were utilized to test the proposed method.

6.1. Prediction Error in the Entire Design Space

Since the number of initial samples is limited, there will be prediction errors when the current established surrogate model is used to predict new samples, which will affect the accuracy and efficiency of prediction. Therefore, it is important to evaluate the prediction error of the established Kriging model first and to select the proper sample plan. The distribution of prediction error in the entire crack parameters space from the Kriging model of measurement point A is shown in Figure 8, under different sample plans.

Figure 8a corresponds to the prediction error by the Kriging model established by the initial 20 samples, as shown in Figure 7. From Figure 8a, one can see that the RMSE (the ratio between MSE and the corresponding amplitude of 2× component) is quite large at the boundary when the crack depth is zero, which means there is a large error when the Kriging model is used to predict sample cases without a crack. In order to improve the prediction performance for no-crack cases, the samples with crack depth zero and crack location ranging from the 10th element to the 52th with a step of 2 elements were added to the initial samples (it should be noted that the added samples are the same case with no crack, however different in mathematic form). The prediction error by the Kriging model re-established is shown in Figure 8b. As one can see from Figure 8b, the prediction error was reduced after adding the samples. However, the error was still large at the area where the crack location was near the 52th element. Then, nine new samples with crack location in the 52th element with crack depth varying in the range [0.05D: 0.05D: 0.45D] were added to improve the Kriging model, leading to results in Figure 8c. As shown in Figure 8c, the RMSE in the entire crack parameter space had been much reduced, and the error at the boundary had also dropped considerably. Though the error at the boundary remains large relatively to the inside area, the Kriging model is thought to be acceptable as a reasonable model for crack identification. The prediction error is the major source of identification error, which can be reduced by a bigger sample size. Therefore, there will be a trade-off between the cost and accuracy.

6.2. Results of Crack Identification

After trial and errors, PSO was adopted with the swarm size of 50, the maximum swarm of 100, and acceleration constants of 1.5, and criterions of $C_1\ge 0.99$ and $C_2\ge 0.95$ were used to determine if the updating process should be stopped. Figure 9 and

Table 3. Crack identification results of the under-test samples.

Under-Test Sample Number	Actual Parameters		Identified Results		Updating Times	Criterion Values for the Final Updated Model $({\mathit{C}}_1,{\mathit{C}}_2)$
	Crack Location (Element)	Crack Depth (D)	Crack Location (Element)	Crack Depth (D)
1	34	0.1	32	0.1	1	(1, 0.96)
2	46	0.05	46	0.05	4	(1, 1)
3	50	0.35	50	0.35	5	(1, 1)
4	28	0.25	28	0.25	2	(1, 1)
5	20	0.4	20	0.4	2	(0.99, 0.99)
6	12	0.3	12	0.3	5	(1, 1)
7	42	0.2	42	0.2	3	(1, 1)
8	38	0.45	38	0.45	2	(1, 1)
9	16	0.15	16	0.15	3	(1, 1)
10	24	0	28	0	1	(1, 1)

are the update process and identification results of the under-test samples in Figure 7. Especially, in Figure 9b, sample 6 was selected as a representative case to show the $Y_{*}\left(x^{*}\right)$ evaluated at the optimal solution x* in current iteration compared with $Y_{\mathrm{Target}}$ of actual parameters.

As the samples in the initial sampling plan and the samples added to update the surrogate model have the resolution of 2 elements in crack location and 0.05D in crack depth, when the crack location and depth identification errors fall in [−2,2] and [−0.05,0.05] respectively, the identification results are considered to be accurate. From Figure 9, one can see all the samples were identified quite effectively. The samples located in the inside area were identified through updating for no more than 3 times, while for samples at the boundary, numbered 2, 3, and 6, needed updating for 4 to 5 times. As for the identification accuracy, from Table 3, one can see that all the crack parameters of the under-test samples were accurately identified, except the crack location of sample 1, which has been identified as the 32th element. The error is thought to be acceptable, considering the resolution of crack locations. In addition, only after one updating step, the update process was stopped. The identification accuracy can be improved if more severe criterions are set, but it will lead to more updating steps. So, there is a trade-off between identification accuracy and efficiency. In view of all the identification results and the error factors, it can be concluded that the proposed crack identification method based on Kriging surrogate model performs quite well for the above unpolluted samples.

To study the robustness of the identification method, each sample was now added with noise. The noise-polluted response $y_{\mathrm{N}}$ can be expressed as [61]:

$$y_{\mathrm{N}}=y+L_{\mathrm{N}}{\mathsf{\sigma }}_{\mathrm{y}}s$$

(20)

where ${\mathrm{L}}_{\mathrm{N}}$ is a number within (0, 1) that represents the noise level, which was selected as 5% in this work; ${\mathsf{\sigma }}_{\mathrm{y}}$ is the standard deviation of original response y after noise is added; and $s$ is an N-length vector of normally distributed random numbers with zero mean and unit variance. Figure 10 and

Table 4. Crack identification results of the under-test samples with 5% noise.

Under-Test Sample Number	Actual Parameters		Identified Results		Updating Times	Criterion Values for the Final Updated Model $({\mathit{C}}_1,{\mathit{C}}_2)$
	Crack Location (Element)	Crack Depth (D)	Crack Location (Element)	Crack Depth (D)
1	34	0.1	34	0.1	2	(1, 1)
2	46	0.05	50	0.1	8	(1, 0.95)
3	50	0.35	50	0.35	9	(1, 1)
4	28	0.25	28	0.25	1	(0.99, 0.99)
5	20	0.4	20	0.4	2	(1, 1)
6	12	0.3	12	0.3	10	(1, 1)
7	42	0.2	42	0.2	4	(1, 1)
8	38	0.45	36	0.45	1	(1, 0.96)
9	16	0.15	16	0.15	2	(1, 1)
10	24	0	24	0	1	(1, 1)

show the update process and identification results of the under-test samples with 5% noise.

Compared with the update process in Figure 9, from Figure 10 one can see that much more updating times are required for the samples at the boundary numbered 2, 3 and 6 with 5% noise, while the updating times for the non-boundary area are more or less the same. On the identification accuracy of the results in Table 4, there were some errors in samples 2 and 8. For sample 8, the reason for the error was the same as sample 1 without noise shown in Table 3, which was the relaxed threshold of the stop criterions. While for sample 2, the reason for the identification error was the prediction error at the boundary. Therefore, noise mainly affects the identification efficiency and accuracy for samples at the boundary, which can be improved by adding more samples at the boundary to reduce the prediction error. For the non-boundary cases, the proposed method is robust. To conclude, after proper sampling plan and adding samples, the proposed method is efficient, accurate and robust.

7. Conclusions

A new crack identification method is proposed for rotors with a breathing crack based on the Kriging surrogate model using the crack induced super-harmonic responses from two measurement points. Numerical experiments were carried out to investigate the performance of the proposed method based on a validated finite element rotor model with a breathing crack. The work in this paper validates the feasibility of the proposed method for crack identification in rotors. The results clearly show its effectiveness and accuracy. In addition, the presented method is robust to noise. Only a small number of samples are needed to construct the Kriging model, which reduces the cost of experiments or simulations. Furthermore, only two measurement points are required and this will be very useful for real applications. The two measurement points can be arranged at the bearing pedestals, which is quite convenient. However, only 2× components are considered to establish the surrogate model. To improve the reliability of the method for real rotors, more crack features should be applied. In addition, only one transverse crack with two crack parameters is considered. Multiple cracks with complex crack parameters may appear in real applications. Though the identification process is almost the same, much more effort is required for multi-crack situations, which will be studied in the near future combined with other efficient optimization methods.