Crack Parameters Identification Based on a Kriging Surrogate Model for Operating Rotors

Parameters identification of cracked rotors has been gaining importance in recent years, but it is still a great challenge to determine the crack parameters including crack location, depth, and angle for operating rotors.*is work proposes a newmethod to identify crack parameters in a rotor-bearing system based on a Kriging surrogate model and an improved nondominated sorting genetic algorithm-III (NSGA-III). A rotor-bearing system with a breathing crack is established by the finite element method and the superharmonic components are used as index to detect the cracks, the Kriging surrogate model between crack parameters and the superharmonic component amplitudes of the vibration response for rotors are constructed, and an improved NSGA-III is proposed to obtain the optimal crack parameters. Numerical experiments show that the proposed method can identify the crack location, depth, and angle accurately and efficiently for operating rotors.


Introduction
Rotors are one of the most important parts of rotating machinery, which are widely used in power stations, aircraft engines, motors, and so on.Transverse fatigue crack is easily produced on rotors under bending excitation, and slant fatigue crack is generated under torsional excitation due to the harsh operation conditions.e crack propagation speed might be extremely slow at the initial stage, as the crack has propagated to a certain depth, the crack propagation speed increases sharply and the shaft may malfunction within a few hours of operation, causing a catastrophic accident [1].In order to reduce the maintenance costs of equipment and avoid the occurrence of serious accidents, it is of enormous significance to continuously monitor the health of rotors during operation.
Many researchers have studied crack monitoring for rotors based on vibration signals [2][3][4][5][6], the methods can be split into two groups depending on whether the rotors are rotating or not.e crack parameter identification methods under nonoperating state of rotors are usually based on modal parameters of the cracked rotors.Dong et al. [7] established a rotor system based on wavelet finite element method and determined the position and depth of the crack by the intersection of contours of the first three natural frequencies.Haji and Olutunde Oyadiji [8] put forward the concept of orthogonal natural frequencies (ONFs) and achieved the location of cracks by calculating the ONFs at different positions of a disc.Zapico-Valle et al. [9] and Yu et al. [10] used the modal parameters of rotors corresponding to different crack positions and depths as the neural network input to identify rotor crack parameters.Xiang et al. [11] estimated the crack parameters by minimizing the errors of natural frequencies between simulation and experimental results.But the modal parameters used in crack identification cannot be easily obtained in operating conditions, and it is more meaningful to identify the crack parameters for rotors under operating conditions.
Generally, crack parameters identification methods for operating rotors need to test the dynamic response of the rotor.Prabhakar et al. [12] detected rotor crack during startup and shutdown of a rotor-bearing system qualitatively with continuous wavelet transform.Singh and Tiwari [13] detected the slope discontinuity in the elastic line of the shaft due to cracks and achieved localization of multiple cracks in a stepped shaft.Saravanan and Sekhar [14] identified the crack in a rotorbearing system utilizing the concept of operational deflection shape (ODS) and kurtosis of vibration response.Zhang et al. [15] detected the singularity of the ODS in rotors to identify the crack location and used the approximate waveform capacity dimension (AWDC) to make the singularity more prominent.Ramesh Babu and Sekhar [16] proposed a concept called amplitude deviation curve method based on ODS and realized the two cracks localization in a rotating rotor.Lu et al. [17] proposed a multicrack location method based on proper orthogonal decomposition (POD) using fractal dimension (FD) and gapped smoothing method (GSM).Lees et al. [18] reviewed the rotor crack parameter identification method based on equivalent crack forces, which treated the effects of cracks as equivalent forces exerted on a healthy rotor system.Some researchers used equivalent crack forces method to obtain the crack position and depth of rotors, such as Sekhar [19][20][21], Markert et al. [22], Pennacchi and et al. [23].e second group of studies concentrates on crack detection and location for operating rotors, but it is difficult to quantitatively identify the shape and depth of the crack.
As the crack angle affects the stress, the speed of crack propagation, and fatigue life of rotors, some researchers studied the dynamic behavior of the rotor due to the slant crack [24][25][26][27][28].However, there is no relevant research on crack parameters identification if crack angle is considered for operating rotors at present, transverse crack and slant crack have similar stiffness, and crack detection results could be misleading.Focusing on this challenge, the problem of crack parameters identification is transformed into multiobjective optimization problem.en, NSGA-III [29,30] is improved to obtain the optimal crack location, depth, and angle at the same time, and the Kriging surrogate model [31,32] is applied to increase the optimization speed.
In this paper, a new crack parameters identification method is proposed for operating rotors using Kriging surrogate model and an improved NSGA-III.In Section 2, the model of a rotor-bearing system with breathing cracks is constructed by the finite element method, and the influence of different crack parameters on the amplitude of superharmonic components is discussed.In Section 3, the Kriging surrogate model between crack parameters and the superharmonic component amplitudes is constructed.In Section 4, using an improved NSGA-III to predict the optimal crack parameters based on the Kriging surrogate model by minimizing the multiobjective function values related to the superharmonic components amplitudes.Numerical experiments show that the proposed method can identify the crack location, depth, and angle for operating rotors accurately and efficiently.

Superharmonic Components Analysis of
a Rotor-Bearing System with a Breathing Crack

Equations of Motion of a Cracked Two-Disc Rotor-Bearing
System.As shown in Figure 1, a two-disc rotor-bearing system is established by the finite element method, and the physical parameters are given in Table 1.e shaft is discretized into 60 Timoshenko beam elements with six degrees of freedom at each node, the two discs are assumed as rigid bodies with three translational and three rotational inertias, and the gyroscopic effects of which are considered.e ball bearings are considered as linear springs and dampers.
By assembling the cracked shaft element, the uncracked shaft element, the rigid discs, and the bearings, the equation of motion of the cracked rotor system in a fixed coordinate system can be written as where M s and M r are the mass matrices of the shaft and the two discs, respectively, C s � αM s + βK s is the shaft Rayleigh damping matrix, C b is the damping matrix of the bearings, Ω is the rotational speed of the shaft, G s is the shaft gyroscopic matrix, G r is the gyroscopic matrix of the two discs, q is the displacement matrix of the node, K s is the stiffness matrix of the shaft, K b is the stiffness matrix of the bearing, K c (t) is the stiffness matrix of the crack element, F u is the excitation due to static unbalance of discs, F g is the gravitational force, and F ex is the external excitation during operation.
Taking into account the bending and torsional coupling excitation caused by the eccentricity of the discs, the excitation acted upon the discs is When rotors are rotating at a constant speed Ω, the elements in F u can be expressed as [33] 2 Shock and Vibration where θ x (t) is the torsional angle, β is the unbalance orientation angle of the disc, and e is the disc eccentricity.

Model of the Shaft Element with Breathing Crack.
Figure 2(a) shows a cracked shaft element of length L and radius R, a is the crack depth, 2b is the crack width, P 1 − P 12 are resultant forces acted upon the 12 nodes, x − y − z is the global coordinate system, and x ′ − y ′ − z ′ is the local coordinate system.x L is the distance between the center of the crack and the left end of the element, and θ is the crack angle between the cross section and the shaft centerline.e most probable angle is 45 °for a slant crack and 90 °for a transverse crack in practical engineering [25,26,34], and so the crack angles studied in this paper are 45 °and 90 °. e flexibility matrix G 0 of the uncracked element and the additional flexibility matrix G c of the cracked element can be derived based on SERR theory [4], and the detailed expressions can be obtained from [24].
G ce is the total flexibility matrix of the crack element, given by e stiffness matrix of the cracked element K ce and uncracked element K uce can be represented as where T is the transformation matrix written as So as to simulate the breathing behavior of the crack, the stiffness of the crack element is calculated based on the classical theory of crack closure line (CCL) [24], which is a hypothetical line separating the closed and open parts of the crack (Figure 2).e crack front is subdivided into m parts, the opening mode stress intensity factor K I is calculated by Equation ( 7), and a positive K I corresponds to the open crack state and a negative K I to the closed state.us, the location of the CCL is determined, and the stiffness matrix of the crack element can be obtained as e Newmark method [35] is utilized to solve the equations numerically, setting the Newmark parameters δ � 0.5 and α n � 0.25 to ensure numerical stability.

Superharmonic Components Analysis of the Cracked Rotor
System.
e superharmonic components in subcritical speed region can be used as index to detect the cracks in the rotorbearing systems [36,37].e vertical vibration response of the rotor-bearing system is collected at four measuring points (illustrated as A, B, C, and D in Figure 1), and the rotating speed is 1/3 of the first critical speed (840 r/min) of the system.e 1X, 2X, and 3X superharmonic components amplitudes of the vibration response at the measuring point A are obtained through fast Fourier transform (FFT).
Figure 3 depicts the amplitude curves of superharmonic components of the vertical vibration response with different crack depths as the crack is located at the 21st element of the shaft.It shows that the 2X amplitude increases with increasing crack depths, the 1X and 2X amplitudes resulted from transverse crack are greater than those from slant crack, and the 1X, 2X, and 3X amplitudes by transverse crack are close to those by slant crack as the crack depth is small.Figure 4 depicts the amplitudes curves of superharmonic components of the vertical vibration response  As the 1X and 2X amplitudes of the vertical vibration response for the cracked rotor-bearing system are more sensitive to crack parameters, then the 1X and 2X amplitudes are selected as the index vectors for crack identi cation, and the superharmonic features of the four measuring points are utilized to avoid the multisolution   Shock and Vibration problem in the process of crack parameters identification.

Construction of the Kriging Surrogate Model
e problem of crack parameters identification can be transformed to an optimization problem.Using optimization methods to estimate the optimal parameters via minimizing the objective function relating to output deviations such as dynamical response and modal parameters [38,39].However, the optimization process requires a large number of iterative steps, and each iterative step needs to repeat calculation of complex finite element models, making the problem very complicated.us, it is necessary to establish an uncomplicated relationship between the crack parameters and the superharmonic features for crack identification.
e Kriging surrogate model provides explicit functions to represent the relationships between inputs and outputs of a linear or nonlinear system [39], which is a statistics-based interpolation method and is not affected by random error [31,32].Construction of the Kriging model requires small amount of samples to obtain high estimation accuracy, significantly reducing calculation time.
Given initial crack parameters samples and the corresponding 1X and 2X amplitudes are is the ith component of sampling crack parameters (l i is the crack position, d i is the crack depth, and θ i is the crack angle), represented as a three-dimensional variable vector, l * ∈ [5,55], l * ∈ Z, d * ∈ [0.05, 0.5], and where y i is the ith set of the output Y vector, amp mX n (x i ) is the superharmonic components amplitudes at the nth measurement point, and n � 1, 2, 3, 4 and m � 1, 2.
And their relationship can be written as where f(x i ) is a vector of a linear combination of selected functions, β is the matrix of regression coefficients, ε l (x i ) is the stochastic process with a variance of σ 2 l and a mean of zero, and q is the dimension of the Y vector.e covariance matrix is expressed by R is the Gauss correlation function matrix and θ is the correlation parameter, given by When existing an untested crack parameter x * , the corresponding superharmonic component amplitudes can be estimated as the linear predictor, which is shown as where c is a n × 1 coefficient vector and y l is the lth column of the matrix Y. e correlation vector r(x) can be written as where x i is the initial crack samples and x * is the new sample.e predicted error is in which To keep the predictor unbiased, it is required that Under this condition, the mean squared error of Equation ( 13) is Shock and Vibration e Lagrangian function for the problem of minimizing the mean squared error with the constraint of Equation ( 17) is e gradient of Equation ( 19) with respect to c is According to the first order necessary conditions for optimality Substituting the above equation into Equation ( 10) gives Hence, the relation between crack parameters and the corresponding 1X and 2X amplitudes values has been deduced.
After the Kriging surrogate model is established, the squared multiple correlations (SMC) and the empirical integrated squared error (EISE) criterion [40] can be used to evaluate the preciseness of the Kriging surrogate model.
where  y l is the lth set of the Y vector of the Kriging surrogate model, y l is the lth component of the superharmonic component amplitudes calculated by finite element analysis,  y is the mean of the true values, and q is the length of the vector  y l . is paper adopts the multipoint adding strategy to update the Kriging surrogate model if the predicted superharmonic components amplitudes cannot meet SC and EISE [41].

Crack Parameters Identification Based on an Improved NSGA-III
After the Kriging surrogate model between crack parameters and the superharmonic components amplitudes is established, crack parameters identification problem is transformed into a multiobjective optimization problem, which can be stated as follows: where amp mX n (x * ) and amp mX n (x Target ) are the m× superharmonic component amplitudes at the nth measurement point.amp mX n (x * ) is calculated via the Kriging surrogate model and amp mX n (x Target ) is measured on actual cracked rotor, and n � 1, 2, 3, 4 and m � 1, 2.
It is a great challenge for the classical multiobjective optimization algorithms when the optimization problem involves four or more objectives.To get the global optimal solution, the NSGA-III recently proposed by Deb and Jain [29,30] may be a good choice to optimize 8 objective functions at the same time.NSGA-III is an evolutionary many-objective optimization algorithm based on reference point, which can ensure the diversity of population members via adaptively updating a number of well-spread reference points.
e crack angles considered are 45 °or 90 °, and the crack is located in the l * th element in this paper.If the NSGA-III is directly used for crack parameters identification in the iterative process, the crack parameters of the population members are mostly decimal numbers, causing crack parameters identification results to be inaccurate.In addition, the decimal crack parameters reduce the convergence rate of the objective function and increase the optimization time.
erefore, the NSGA-III needs to be improved for crack identification of the rotor-bearing system.Preprocessing of population members is added before members' selection in every iterative process to obtain crack parameters accurately and efficiently, and the improved part of the NSGA-III is as follows.
(1) Preprocessing of crack position: where l * t is the crack position obtained by preprocessing of l t at the tth generation and Round is the rounding function.

6
Shock and Vibration (2) Preprocessing of crack angle: where θ * t is the crack angle obtained after preprocessing of θ t at the tth generation.
e preprocessing of the crack angle can make the crack angle close to 45 °or 90 °change to 45 °or 90 °.
Replace l t and θ t with l * t and θ * t to form a new population.
ere are many Pareto optimal solutions by the improved NSGA-III.As it is very di cult to have a set of crack parameters that minimize all the objective function values, an index is required to select the optimal crack parameters in the Pareto optimal solutions.By de ning Mobw as the selection index of the crack parameters identi cation results, Mobw can be obtained by weighing each objective function value: where k i is the weight coe cient, and k 1 k 2 k 3 k 4 0.15 and k 5 k 6 k 7 k 8 0.1 in this paper.e corresponding identi cation result of the minimum value of Mobw is taken as the nal crack parameters.

Identi cation Results Based on a Kriging Surrogate Model.
Numerical experiments are carried out to verify the e ectiveness of the method.e detailed parameters of the rotorbearing system are given in Table 1, 300 sets of samples of crack parameters are generated by Latin hypercube sampling [42], and then corresponding vertical vibration response of the rotor-bearing system assuming a 5% random white noise is collected at four measuring points.
en, the Kriging surrogate model between crack parameters and corresponding 1X and 2X amplitudes is established.
e improved NSGA-III is employed to estimate the crack parameters found on the Kriging surrogate model by minimizing the multiobjective function values given in Equation (25), and the initial parameters of the improved NSGA-III are as follows: the population size, maximum number of iterations, crossover percentage, and mutation percentage are 150, 100, 0.5, and 0.5, respectively.e owchart of the method for identifying crack parameters based on the improved NSGA-III and Kriging surrogate model is shown in Figure 5.
is paper considers eight cases (Table 2) to evaluate the performance of the proposed method and selects case 1 to analyze the identi cation process of crack parameters for rotors in detail.
Figure 6 illustrates the mean value of objective function varies with evolutionary generations, and it shows that the mean value of objective function optimized by the improved NSGA-III has already converged in the 30th generation, while the mean value of the objective function with the original NSGA-III not only has higher values but also has a slow convergence speed.
e crack parameters results identi ed with the improved NSGA-III are compared with results with the original NSGA-III (Figure 7), and it can be seen that the crack parameters obtained by the improved NSGA-III are closer to the actual crack parameters in Pareto optimal sets.Figures 6 and 7 indicate that the improved NSGA-III can obtain Pareto optimal solutions accurately and e ciently.
e values of the objective function corresponding to each group of crack parameters in Pareto optimal solution are shown in Figure 8, and it is necessary to select a best set of solutions based on the index Mobw in Equation (29).
As shown in Figure 9, the selection index Mobw achieves the minimum value in the 90th population member, and the corresponding crack parameter is x * i 20, 0.375066, 90 { }.Compared with other crack parameters in Pareto optimal solutions, the 90th set of parameters is closer to the actual values (marked as blue circle in Figure 7), which indicates that it is reasonable to select the optimal crack parameters by the index Mobw.
e identi ed crack parameters are listed in

Shock and Vibration
near the boundary (crack location is close to the 55th element and crack depth is close to 0.1 cm), so the 1X and 2X amplitudes predicted by the Kriging surrogate model are less accurate.
e accuracy of crack parameter identification results is improved by increasing the number of samples in the boundary area.

Comparison with the Results by Artificial Neural
Networks.Artificial neural network (ANN) models are widely used in engineering for regression and classification, which have been proved to be effective in identifying the crack location and depth using modal parameters of static rotors [9,10].In order to explain the accuracy and efficiency of the proposed method in identifying the crack parameters, ANN is used to construct the relationship between superharmonic components amplitudes and crack parameters of a rotating rotor-bearing system.e configuration of ANN used in this paper is backpropagation algorithm (BPA), and the input and target output for training are expressed, respectively, as , where P � Y T , Y is the superharmonic components matrix in Equation ( 9) and Τ� X T , X is the crack parameters matrix in Equation (10).A single hidden layer containing 20 neurons is designed for training process, and the training progress of ANN uses Levenberg-Marquardt backpropagation algorithm.e 300 sets of samples same as that in Section 5.1 are used to construct the neural network, assuming the 5% random white noise.Table 3 shows crack parameters identified by ANN, and it is imprecise compared with identified results obtained by the proposed method (Table 2) in this paper.
In addition, the crack parameters identification results obtained by the training neural network with 600 sets of samples are shown in Table 4, and it can be found that increasing the samples number trained by neural network can improve the identification accuracy of crack parameters.
However, increasing the number of samples requires more time for finite element calculation, and identified results are still imprecise compared with that obtained by the proposed method (Table 2) in this paper.

Figure 1 :
Figure 1:e two-disc rotor-bearing system with a breathing crack.

Figure 5 :
Figure 5: Flowchart of the method for crack parameters identi cation.

Figure 6 :
Figure 6: e mean value of objective function varies with evolutionary generations obtained by (a) original NSGA-III and (b) improved NSGA-III.

Figure 7 :
Figure 7: Comparison of Pareto optimal sets obtained by two methods: (a) crack position; (b) crack depth; (c) crack angle.

Figure 8 :
Figure 8: e objective function values corresponding to crack parameters in Pareto optimal solution.

Table 1 :
Parameters of the cracked rotor.
with di erent crack location as the crack depth is 0.25 cm.It shows that the 1X, 2X, and 3X amplitudes generated by transverse crack are greater than those by slant crack and the 1X, 2X, and 3X amplitudes generated by transverse crack are close to those by slant crack as the crack is near the bearing.

Table 2 ,
and it is worth noting that results are very accurate except for case 3, as the crack depth and crack position in case 3 are

Table 2 :
Crack parameters identification results using 300 sets of samples.