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

Authors: Danyang Wang; Chunrong Hua; Dawei Dong; Biao He; Zhiwen Lu
Journal: Shock and Vibration
 (2018)

---

## Abstract

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. This work proposes a new method 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.


---

## Body

## 1. 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. The 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–6], the methods can be split into two groups depending on whether the rotors are rotating or not. The 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 start-up 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 rotor-bearing 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–21], Markert et al. [22], Pennacchi and et al. [23]. The 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–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. Then, 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 Section2, 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.


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

### 2.1. Equations of Motion of a Cracked Two-Disc Rotor-Bearing System

As shown in Figure1, a two-disc rotor-bearing system is established by the finite element method, and the physical parameters are given in Table 1. The 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. The ball bearings are considered as linear springs and dampers.Figure 1
The two-disc rotor-bearing system with a breathing crack.Table 1
Parameters of the cracked rotor.
Parameters
Value (units)
Shaft length (l0)
0.6 m
Shaft diameter (D)
0.01 m
Position of bearing 1 (l4)
0.015 m
Position of bearing 2 (l2)
0.575 m
Bearing stiffness
9.6 × 105 N/m
Bearing damping
500 Ns/m
Position of disc 1 (l3)
0.095 m
Position of disc 2 (l2)
0.505 m
Disc diameter
0.074 m
Disc thickness
0.025 m
Disc eccentricity (e)
2 × 10−5 N/m
Unbalance orientation angle
0
Density of steel
7.8 × 103 kg/m3
Young’s modulus
3.11 × 1011 Pa
Poisson’s ratio
0.3
Gravitational acceleration
9.8 m/s2
Rayleigh damping coefficient (a)
0.6842
Rayleigh damping coefficient (b)
2.8012 × 10−5
First critical speed
2520 r/min
Rotating speed
840 r/min
Measuring point (lA)
0.145 m
Measuring point (lB)
0.245 m
Measuring point (lC)
0.345 m
Measuring point (lD)
0.445 mBy 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(1)Ms+Mrq¨+Cb+Cs+ΩGs+Grq˙+Kb+Ks+Kctq=Fu+Fg+Fex,where Ms and Mr are the mass matrices of the shaft and the two discs, respectively, Cs=αMs+βKs is the shaft Rayleigh damping matrix, Cb is the damping matrix of the bearings, Ω is the rotational speed of the shaft, Gs is the shaft gyroscopic matrix, Gr is the gyroscopic matrix of the two discs, q is the displacement matrix of the node, Ks is the stiffness matrix of the shaft, Kb is the stiffness matrix of the bearing, Kct is the stiffness matrix of the crack element, Fu is the excitation due to static unbalance of discs, Fg is the gravitational force, and Fex 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(2)Fui=Fuxi,Fuyi,Fuzi,Muxi,Muyi,MuziT.When rotors are rotating at a constant speedΩ, the elements in Fu can be expressed as [33](3)Fuxi=0,Muyi=0,Muzi=0,Fuyi=meΩ+θ˙xi2cosΩt+θxi+β−θ¨xisinΩt+θxi+β,Fuzi=meΩ+θ˙xi2sinΩt+θxi+β+θ¨xicosΩt+θxi+β,Muxi=mey¨sinΩt+θxi+β−z¨+gcosΩt+θxi+β,where θxt is the torsional angle, β is the unbalance orientation angle of the disc, and e is the disc eccentricity.


### 2.2. Model of the Shaft Element with Breathing Crack

Figure2(a) shows a cracked shaft element of length L and radius R, a is the crack depth, 2b is the crack width, P1−P12 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. xL 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. The 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°.Figure 2
Schematic diagram of cracked shaft element. (a) Cracked shaft element. (b) Crack cross section.
(a)
(b)The flexibility matrixG0 of the uncracked element and the additional flexibility matrix Gc 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(4)Gce=G0+Gc.The stiffness matrix of the cracked elementKce and uncracked element Kuce can be represented as(5)Kce=TGce−1TT,Kuce=TG0−1TT,where T is the transformation matrix written as(6)TT=100000−1000000100000−1000l00100000−10−l0000100000−1000000100000−1000000100000−1.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). The crack front is subdivided into m parts, the opening mode stress intensity factor KI is calculated by Equation (7), and a positive KI corresponds to the open crack state and a negative KI to the closed state. Thus, the location of the CCL is determined, and the stiffness matrix of the crack element can be obtained as(7)KI=∑i=16KIi.The Newmark method [35] is utilized to solve the equations numerically, setting the Newmark parameters δ = 0.5 and αn = 0.25 to ensure numerical stability.


### 2.3. Superharmonic Components Analysis of the Cracked Rotor System

The superharmonic components in subcritical speed region can be used as index to detect the cracks in the rotor-bearing systems [36, 37]. The 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. The 1X, 2X, and 3X superharmonic components amplitudes of the vibration response at the measuring point A are obtained through fast Fourier transform (FFT).Figure3 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 3
The amplitude curves of superharmonic components with different crack depths. (a) 1X amplitude; (b) 2X amplitude; (c) 3X amplitude.
(a)
(b)
(c)Figure4 depicts the amplitudes curves of superharmonic components of the vertical vibration response with different 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.Figure 4
The superharmonic components amplitudes with different crack locations. (a) 1X amplitude; (b) 2X amplitude; (c) 3X amplitude.
(a)
(b)
(c)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 identification, and the superharmonic features of the four measuring points are utilized to avoid the multisolution problem in the process of crack parameters identification.


## 2.1. Equations of Motion of a Cracked Two-Disc Rotor-Bearing System

As shown in Figure1, a two-disc rotor-bearing system is established by the finite element method, and the physical parameters are given in Table 1. The 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. The ball bearings are considered as linear springs and dampers.Figure 1
The two-disc rotor-bearing system with a breathing crack.Table 1
Parameters of the cracked rotor.
Parameters
Value (units)
Shaft length (l0)
0.6 m
Shaft diameter (D)
0.01 m
Position of bearing 1 (l4)
0.015 m
Position of bearing 2 (l2)
0.575 m
Bearing stiffness
9.6 × 105 N/m
Bearing damping
500 Ns/m
Position of disc 1 (l3)
0.095 m
Position of disc 2 (l2)
0.505 m
Disc diameter
0.074 m
Disc thickness
0.025 m
Disc eccentricity (e)
2 × 10−5 N/m
Unbalance orientation angle
0
Density of steel
7.8 × 103 kg/m3
Young’s modulus
3.11 × 1011 Pa
Poisson’s ratio
0.3
Gravitational acceleration
9.8 m/s2
Rayleigh damping coefficient (a)
0.6842
Rayleigh damping coefficient (b)
2.8012 × 10−5
First critical speed
2520 r/min
Rotating speed
840 r/min
Measuring point (lA)
0.145 m
Measuring point (lB)
0.245 m
Measuring point (lC)
0.345 m
Measuring point (lD)
0.445 mBy 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(1)Ms+Mrq¨+Cb+Cs+ΩGs+Grq˙+Kb+Ks+Kctq=Fu+Fg+Fex,where Ms and Mr are the mass matrices of the shaft and the two discs, respectively, Cs=αMs+βKs is the shaft Rayleigh damping matrix, Cb is the damping matrix of the bearings, Ω is the rotational speed of the shaft, Gs is the shaft gyroscopic matrix, Gr is the gyroscopic matrix of the two discs, q is the displacement matrix of the node, Ks is the stiffness matrix of the shaft, Kb is the stiffness matrix of the bearing, Kct is the stiffness matrix of the crack element, Fu is the excitation due to static unbalance of discs, Fg is the gravitational force, and Fex 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(2)Fui=Fuxi,Fuyi,Fuzi,Muxi,Muyi,MuziT.When rotors are rotating at a constant speedΩ, the elements in Fu can be expressed as [33](3)Fuxi=0,Muyi=0,Muzi=0,Fuyi=meΩ+θ˙xi2cosΩt+θxi+β−θ¨xisinΩt+θxi+β,Fuzi=meΩ+θ˙xi2sinΩt+θxi+β+θ¨xicosΩt+θxi+β,Muxi=mey¨sinΩt+θxi+β−z¨+gcosΩt+θxi+β,where θxt is the torsional angle, β is the unbalance orientation angle of the disc, and e is the disc eccentricity.


## 2.2. Model of the Shaft Element with Breathing Crack

Figure2(a) shows a cracked shaft element of length L and radius R, a is the crack depth, 2b is the crack width, P1−P12 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. xL 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. The 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°.Figure 2
Schematic diagram of cracked shaft element. (a) Cracked shaft element. (b) Crack cross section.
(a)
(b)The flexibility matrixG0 of the uncracked element and the additional flexibility matrix Gc 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(4)Gce=G0+Gc.The stiffness matrix of the cracked elementKce and uncracked element Kuce can be represented as(5)Kce=TGce−1TT,Kuce=TG0−1TT,where T is the transformation matrix written as(6)TT=100000−1000000100000−1000l00100000−10−l0000100000−1000000100000−1000000100000−1.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). The crack front is subdivided into m parts, the opening mode stress intensity factor KI is calculated by Equation (7), and a positive KI corresponds to the open crack state and a negative KI to the closed state. Thus, the location of the CCL is determined, and the stiffness matrix of the crack element can be obtained as(7)KI=∑i=16KIi.The Newmark method [35] is utilized to solve the equations numerically, setting the Newmark parameters δ = 0.5 and αn = 0.25 to ensure numerical stability.


## 2.3. Superharmonic Components Analysis of the Cracked Rotor System

The superharmonic components in subcritical speed region can be used as index to detect the cracks in the rotor-bearing systems [36, 37]. The 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. The 1X, 2X, and 3X superharmonic components amplitudes of the vibration response at the measuring point A are obtained through fast Fourier transform (FFT).Figure3 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 3
The amplitude curves of superharmonic components with different crack depths. (a) 1X amplitude; (b) 2X amplitude; (c) 3X amplitude.
(a)
(b)
(c)Figure4 depicts the amplitudes curves of superharmonic components of the vertical vibration response with different 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.Figure 4
The superharmonic components amplitudes with different crack locations. (a) 1X amplitude; (b) 2X amplitude; (c) 3X amplitude.
(a)
(b)
(c)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 identification, and the superharmonic features of the four measuring points are utilized to avoid the multisolution problem in the process of crack parameters identification.


## 3. Construction of the Kriging Surrogate Model

The 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. Thus, it is necessary to establish an uncomplicated relationship between the crack parameters and the superharmonic features for crack identification.The 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(8)X=x1x2x3⋮xn=l1d1β1l2d2β2l2d2β2⋮⋮⋮lndnβn,where xi=li,di,θi is the ith component of sampling crack parameters (li is the crack position, di 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 θ∗∈45,90.(9)Y=y1y2y3⋮yn=amp11Xx1amp21Xx1amp31Xx1amp41Xx1amp12Xx1amp22Xx1amp32Xx1amp42Xx1amp11Xx2amp21Xx2amp31Xx2amp41Xx2amp12Xx2amp22Xx2amp32Xx2amp42Xx2amp11Xx3amp21Xx3amp31Xx3amp41Xx3amp12Xx3amp22Xx3amp32Xx3amp42Xx3⋮⋮⋮⋮⋮⋮⋮⋮amp11Xxnamp21Xxnamp31Xxnamp41Xxnamp12Xxnamp22Xxnamp32Xxnamp42Xxn,where yi is the ith set of the output Y vector, ampnmXxi 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(10)ylxi=fxiTβ+εlxi,i=1,2,...,n,l=1,2,…,q,where fxi is a vector of a linear combination of selected functions, β is the matrix of regression coefficients, εlxi is the stochastic process with a variance of σl2 and a mean of zero, and q is the dimension of the Y vector. The covariance matrix is expressed by(11)Covεlxi,εlxj=σl2Rθ,xi,xj,i,j=1,2,...,n,l=1,2,…,q.R is the Gauss correlation function matrix and θ is the correlation parameter, given by(12)Ri,jθ,xi,xj=exp−θxi−xj,i,j=1,2,…,n.When existing an untested crack parameterx∗, the corresponding superharmonic component amplitudes can be estimated as the linear predictor, which is shown as(13)y^l∗x∗=cΤyl,l=1,2,…,q,where c is a n×1 coefficient vector and yl is the lth column of the matrix Y.The correlation vectorrx can be written as(14)rx=Ri,jθ,x1,x∗,Ri,jθ,x2,x∗,…,Ri,jθ,xn,x∗T,where xi is the initial crack samples and x∗ is the new sample.The predicted error is(15)y^l∗x∗−ylx∗=cTyl−ylx∗=cTZl−zlx∗+FTc−fx∗Tβl,in which(16)Zl=z1x1,z1x2,…,z1xn,F=f1x1⋯fpx1⋮⋱⋮f1xn⋯fpxn.To keep the predictor unbiased, it is required that(17)FTc−fx∗=0.Under this condition, the mean squared error of Equation (13) is(18)φlx∗=Eyl∗x∗−ylx∗2=EcTzl−zlx∗2=σl21+cTRc−2cTr.The Lagrangian function for the problem of minimizing the mean squared error with the constraint of Equation (17) is(19)Lc,λ=σl21+cTRc−2cTr−λTFTc−fx∗.The gradient of Equation (19) with respect to c is(20)L′cc,λ=2σl2Rc−r−Fλ.According to the first order necessary conditions for optimality(21)c=R−1r−FFTR−1F−1FTR−1r−fx∗.Substituting the above equation into Equation (10) gives(22)yl∗x∗=R−Tr−FFTR−1F−1FTR−1r−fx∗Tyl.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.(23)SC=1−∑j=1qy^l−yl2∑j=1qyl−y^2,EISE=1q∑l=1qy^l−yl2,where y^l is the lth set of the Y vector of the Kriging surrogate model, yl 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. This paper adopts the multipoint adding strategy to update the Kriging surrogate model if the predicted superharmonic components amplitudes cannot meet SC and EISE [41].


## 4. 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:(24)findx∗=l∗,d∗,θ∗,(25)minObject1=amp11Xx∗−amp11XxTarget,Object2=amp21Xx∗−amp21XxTarget,Object3=amp31Xx∗−amp31XxTarget,Object4=amp41Xx∗−amp41XxTarget,Object5=amp12Xx∗−amp12XxTarget,Object6=amp22Xx∗−amp22XxTarget,Object7=amp32Xx∗−amp32XxTarget,Object8=amp42Xx∗−amp42XxTarget,with(26)l∗∈5,55,l∗∈Z,d∗∈0.05,0.5,θ∗∈45,90,where ampnmXx∗ and ampnmXxTarget are the m× superharmonic component amplitudes at the nth measurement point. ampnmXx∗ is calculated via the Kriging surrogate model and ampnmXxTarget 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.The crack angles considered are 45° or 90°, and the crack is located in thel∗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.Therefore, 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:(27)lt∗=Roundlt,where lt∗ is the crack position obtained by preprocessing of lt at the tth generation and Round is the rounding function.(2)
Preprocessing of crack angle:(28)θt∗=Roundθt−45°45°+1×45°,where θt∗ is the crack angle obtained after preprocessing of θt at the tth generation. The preprocessing of the crack angle can make the crack angle close to 45° or 90° change to 45° or 90°.Replacelt and θt with lt∗ and θt∗ to form a new population.There are many Pareto optimal solutions by the improved NSGA-III. As it is very difficult 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 definingMobw as the selection index of the crack parameters identification results, Mobw can be obtained by weighing each objective function value:(29)Mobw=∑i=18Objecti×ki,where ki is the weight coefficient, and k1=k2=k3=k4=0.15 and k5=k6=k7=k8=0.1 in this paper.The corresponding identification result of the minimum value ofMobw is taken as the final crack parameters.


## 5. Numerical Investigation

### 5.1. Identification Results Based on a Kriging Surrogate Model

Numerical experiments are carried out to verify the effectiveness of the method. The detailed parameters of the rotor-bearing system are given in Table1, 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. Then, the Kriging surrogate model between crack parameters and corresponding 1X and 2X amplitudes is established.The 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.The flowchart of the method for identifying crack parameters based on the improved NSGA-III and Kriging surrogate model is shown in Figure5.Figure 5
Flowchart of the method for crack parameters identification.This paper considers eight cases (Table2) to evaluate the performance of the proposed method and selects case 1 to analyze the identification process of crack parameters for rotors in detail.Table 2
Crack parameters identification results using 300 sets of samples.
Case
Actual crack parameters
Predicted crack parameters
Percentage error
Position,l (element)
Depth,D (cm)
Angle,θ (degrees)
Position,l∗ (element)
Depth,D∗ (cm)
Angle,θ∗ (degrees)
Error (%)
Error (%)
Error (%)
1
20
0.375
90
20
0.37507
90
0.00
0.01
0.00
2
31
0.25
45
31
0.24988
45
0.00
0.05
0.00
3
50
0.1
45
49
0.15224
45
2.00
52.4
0.00
4
25
0.125
45
25
0.12145
45
0.00
2.84
0.00
5
21
0.225
90
21
0.20255
90
0.00
9.98
0.00
6
19
0.45
90
19
0.44813
90
0.00
0.42
0.00
7
49
0.225
45
49
0.22714
45
0.00
0.95
0.00
8
17
0.3
45
17
0.29411
45
0.00
1.96
0.00Figure6 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.Figure 6
The mean value of objective function varies with evolutionary generations obtained by (a) original NSGA-III and (b) improved NSGA-III.
(a)
(b)The crack parameters results identified with the improved NSGA-III are compared with results with the original NSGA-III (Figure7), 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.Figure 7
Comparison of Pareto optimal sets obtained by two methods: (a) crack position; (b) crack depth; (c) crack angle.
(a)
(b)
(c)Figures6 and 7 indicate that the improved NSGA-III can obtain Pareto optimal solutions accurately and efficiently.The values of the objective function corresponding to each group of crack parameters in Pareto optimal solution are shown in Figure8, and it is necessary to select a best set of solutions based on the index Mobw in Equation (29).Figure 8
The objective function values corresponding to crack parameters in Pareto optimal solution.As shown in Figure9, the selection index Mobw achieves the minimum value in the 90th population member, and the corresponding crack parameter is xi∗=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.Figure 9
IndexMobw selection of crack identification results.The identified crack parameters are listed in Table2, 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 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. The accuracy of crack parameter identification results is improved by increasing the number of samples in the boundary area.


### 5.2. 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.The configuration of ANN used in this paper is backpropagation algorithm (BPA), and the input and target output for training are expressed, respectively, asP=P1P2⋯Pn and T=T1T2⋯Tn, where P=YT, Y is the superharmonic components matrix in Equation (9) and Τ=XT, 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.The 300 sets of samples same as that in Section5.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.Table 3
Crack parameters identification results by ANN using 300 sets of samples.
Case
Actual crack parameters
Predicted crack parameters by ANN
Percentage error
Position,l (element)
Depth,D (cm)
Angle,θ (degrees)
Position,l∗ (element)
Depth,D∗ (cm)
Angle,θ∗ (degrees)
Error (%)
Error (%)
Error (%)
1
20
0.375
90
20
0.39386
90
0.00
5.03
0.00
2
31
0.25
45
31
0.25569
45
0.00
2.28
0.00
3
50
0.1
45
48
0.10314
45
4.00
3.14
0.00
4
25
0.125
45
28
0.11122
45
12.00
11.02
0.00
5
21
0.225
90
20
0.22909
45
4.76
1.81
50.00
6
19
0.45
90
23
0.49150
90
21.01
9.22
0.00
7
49
0.225
45
49
0.22645
45
0.00
0.64
0.00
8
17
0.3
45
17
0.27013
45
0.00
9.96
0.00In addition, the crack parameters identification results obtained by the training neural network with 600 sets of samples are shown in Table4, and it can be found that increasing the samples number trained by neural network can improve the identification accuracy of crack parameters.Table 4
Crack parameters identification results by ANN using 600 sets of samples.
Case
Actual crack parameters
Predicted crack parameters by ANN
Percentage error
Position,l (element)
Depth,D (cm)
Angle,θ (degrees)
Position,l∗ (element)
Depth,D∗ (cm)
Angle,θ∗ (degrees)
Error (%)
Error (%)
Error (%)
1
20
0.375
90
20
0.39146
90
0.00
4.39
0.00
2
31
0.25
45
31
0.23780
45
0.00
4.88
0.00
3
50
0.1
45
49
0.12001
45
2.00
20.01
0.00
4
25
0.125
45
27
0.11725
45
8.00
6.20
0.00
5
21
0.225
90
22
0.22614
90
4.76
0.51
0.00
6
19
0.45
90
22
0.49150
90
15.79
9.22
0.00
7
49
0.225
45
49
0.22591
45
0.00
0.40
0.00
8
17
0.3
45
17
0.29535
45
0.00
1.55
0.00However, 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 (Table2) in this paper.Based on the comparison result, it is clear to draw that the proposed method based on Kriging surrogate model and an improved NSGA-III can identify the crack parameters of operating rotors more accurately with fewer samples than ANN.


## 5.1. Identification Results Based on a Kriging Surrogate Model

Numerical experiments are carried out to verify the effectiveness of the method. The detailed parameters of the rotor-bearing system are given in Table1, 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. Then, the Kriging surrogate model between crack parameters and corresponding 1X and 2X amplitudes is established.The 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.The flowchart of the method for identifying crack parameters based on the improved NSGA-III and Kriging surrogate model is shown in Figure5.Figure 5
Flowchart of the method for crack parameters identification.This paper considers eight cases (Table2) to evaluate the performance of the proposed method and selects case 1 to analyze the identification process of crack parameters for rotors in detail.Table 2
Crack parameters identification results using 300 sets of samples.
Case
Actual crack parameters
Predicted crack parameters
Percentage error
Position,l (element)
Depth,D (cm)
Angle,θ (degrees)
Position,l∗ (element)
Depth,D∗ (cm)
Angle,θ∗ (degrees)
Error (%)
Error (%)
Error (%)
1
20
0.375
90
20
0.37507
90
0.00
0.01
0.00
2
31
0.25
45
31
0.24988
45
0.00
0.05
0.00
3
50
0.1
45
49
0.15224
45
2.00
52.4
0.00
4
25
0.125
45
25
0.12145
45
0.00
2.84
0.00
5
21
0.225
90
21
0.20255
90
0.00
9.98
0.00
6
19
0.45
90
19
0.44813
90
0.00
0.42
0.00
7
49
0.225
45
49
0.22714
45
0.00
0.95
0.00
8
17
0.3
45
17
0.29411
45
0.00
1.96
0.00Figure6 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.Figure 6
The mean value of objective function varies with evolutionary generations obtained by (a) original NSGA-III and (b) improved NSGA-III.
(a)
(b)The crack parameters results identified with the improved NSGA-III are compared with results with the original NSGA-III (Figure7), 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.Figure 7
Comparison of Pareto optimal sets obtained by two methods: (a) crack position; (b) crack depth; (c) crack angle.
(a)
(b)
(c)Figures6 and 7 indicate that the improved NSGA-III can obtain Pareto optimal solutions accurately and efficiently.The values of the objective function corresponding to each group of crack parameters in Pareto optimal solution are shown in Figure8, and it is necessary to select a best set of solutions based on the index Mobw in Equation (29).Figure 8
The objective function values corresponding to crack parameters in Pareto optimal solution.As shown in Figure9, the selection index Mobw achieves the minimum value in the 90th population member, and the corresponding crack parameter is xi∗=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.Figure 9
IndexMobw selection of crack identification results.The identified crack parameters are listed in Table2, 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 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. The accuracy of crack parameter identification results is improved by increasing the number of samples in the boundary area.


## 5.2. 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.The configuration of ANN used in this paper is backpropagation algorithm (BPA), and the input and target output for training are expressed, respectively, asP=P1P2⋯Pn and T=T1T2⋯Tn, where P=YT, Y is the superharmonic components matrix in Equation (9) and Τ=XT, 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.The 300 sets of samples same as that in Section5.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.Table 3
Crack parameters identification results by ANN using 300 sets of samples.
Case
Actual crack parameters
Predicted crack parameters by ANN
Percentage error
Position,l (element)
Depth,D (cm)
Angle,θ (degrees)
Position,l∗ (element)
Depth,D∗ (cm)
Angle,θ∗ (degrees)
Error (%)
Error (%)
Error (%)
1
20
0.375
90
20
0.39386
90
0.00
5.03
0.00
2
31
0.25
45
31
0.25569
45
0.00
2.28
0.00
3
50
0.1
45
48
0.10314
45
4.00
3.14
0.00
4
25
0.125
45
28
0.11122
45
12.00
11.02
0.00
5
21
0.225
90
20
0.22909
45
4.76
1.81
50.00
6
19
0.45
90
23
0.49150
90
21.01
9.22
0.00
7
49
0.225
45
49
0.22645
45
0.00
0.64
0.00
8
17
0.3
45
17
0.27013
45
0.00
9.96
0.00In addition, the crack parameters identification results obtained by the training neural network with 600 sets of samples are shown in Table4, and it can be found that increasing the samples number trained by neural network can improve the identification accuracy of crack parameters.Table 4
Crack parameters identification results by ANN using 600 sets of samples.
Case
Actual crack parameters
Predicted crack parameters by ANN
Percentage error
Position,l (element)
Depth,D (cm)
Angle,θ (degrees)
Position,l∗ (element)
Depth,D∗ (cm)
Angle,θ∗ (degrees)
Error (%)
Error (%)
Error (%)
1
20
0.375
90
20
0.39146
90
0.00
4.39
0.00
2
31
0.25
45
31
0.23780
45
0.00
4.88
0.00
3
50
0.1
45
49
0.12001
45
2.00
20.01
0.00
4
25
0.125
45
27
0.11725
45
8.00
6.20
0.00
5
21
0.225
90
22
0.22614
90
4.76
0.51
0.00
6
19
0.45
90
22
0.49150
90
15.79
9.22
0.00
7
49
0.225
45
49
0.22591
45
0.00
0.40
0.00
8
17
0.3
45
17
0.29535
45
0.00
1.55
0.00However, 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 (Table2) in this paper.Based on the comparison result, it is clear to draw that the proposed method based on Kriging surrogate model and an improved NSGA-III can identify the crack parameters of operating rotors more accurately with fewer samples than ANN.


## 6. Conclusions

A new crack parameters identification method based on a Kriging surrogate model and an improved NSGA-III is presented in the paper. The numerical results clearly show that the proposed approach can identify the crack location, depth, and angle for operating rotors accurately and efficiently.(1)
The 1X and 2X superharmonic components amplitudes of a rotor-bearing system in 1/3 subcritical speed region can be selected as the index vectors for crack parameters identification.(2)
Crack parameters and their corresponding 1X and 2X amplitudes of vibration response for rotors are used to establish the Kriging surrogate model, which avoids complicated finite element calculation, demanding a large number of iterative steps.(3)
NSGA-III is improved by adding the preprocessing of population members to avoid the disadvantages of the original NSGA-III for crack parameters identification, which can improve the accuracy of crack parameters identification and speed up the convergence of the multiobjective function values.(4)
Compared with ANN, the proposed method can identify crack parameters for operating rotors more accurately with fewer samples.It should be noted that only one crack is considered in this paper. There may be multiple cracks in rotors, and identification of multiple crack parameters is more complicated and requires further study.



---

*Source: 9274526-2018-10-08.xml*