# Elliptical Crack Identification in a Nonrotating Shaft

Authors: B. Muñoz-Abella; L. Rubio; P. Rubio; L. Montero
Journal: Shock and Vibration
 (2018)

---

## Abstract

It is known that fatigue cracks are one of the most important problems of the mechanical components, since their propagation can cause severe loss, both personal and economic. So, it is essential to know deeply the behavior of the cracked element to have tools that allow predicting the breakage before it happens. The shafts are elements that are specially affected by the described problem, because they are subjected to alternative compression and tension stresses. This work presents, firstly, an analytical expression that allows determining the first four natural frequencies of bending vibration of a nonrotating cracked shaft, assumed as an Euler–Bernoulli beam, with circular cross section under pinned-pinned conditions, taking into account the elliptical shape of the crack. Second, once the direct problem is known, the inverse problem is approached. Genetic Algorithm technique has been used to estimate the crack parameters assuming known the natural frequencies of the cracked shaft.


---

## Body

## 1. Introduction

Shafts are a very common mechanical components that are usually subjected to very hard working conditions, specifically bending loads that originate periodic stresses, changing from tensile to compression stresses for each cycle of the shaft, often accompanied by high temperatures and aggressive environments. All of that frequently causes incipient fatigue cracks which can quickly lead to catastrophic failures [1–3]. So, it is very important to use diagnosis methodologies that allow detecting possible cracks and also predicting both the shape and size of cracks in order to replace the part before reaching the total breakdown. When a crack appears in a mechanical component it produces a change in the modal properties (natural frequencies and modes shapes) that is accompanied by an increase in its flexibility [1, 4]. Over the last four decades, from the 70s to the present day, many researches have focused on developing nondestructive maintenance methods founded on the vibration measurement characteristics, clearly the mode shapes and natural frequencies, for instance, some of them are presented in [3, 5–10], besides, some extensive reviews can be found in [11–14]. Some of these works are based on the fact that it is possible to estimate the crack parameters by measuring the changes in the natural frequency (inverse problem), since the natural frequencies can be measured accurately and most easily in comparison with other dynamic characteristics of a structure and the experimental errors do not significantly change their values [15, 16]. Before considering the crack identification inverse problem, the direct problem should be studied, that is, the determination of natural frequencies of the cracked element, assuming the crack parameters are known.First of all, to establish the relationship between crack parameters and natural frequencies, assuming that, during the test period the shaft does not rotate, countless authors consider shafts as beams with a circular cross section, among them [17–20]. According to Heydari et al. [21], one can find three approaches for the vibration modelling of cracked beams: continuous models with a continuous model for the crack, discrete models with a local flexibility model for the crack, and continuous models with a local flexibility model for the crack. With regard to the third approach, the continuous beam models with a local flexibility crack model many times consider two healthy segments that are linked by a massless rotational spring in order to take into account the flexibility due to the crack [5], so that the spring constant depends on the crack parameters by the fracture mechanics theory. This type of model has been successfully applied to Euler–Bernoulli cracked beams under different boundary conditions and rectangular cross sections [8, 22, 23], or circular cross sections [18, 24–27], but also concerning circular section beams or shafts, it is easy to find works that approach the problem by discrete methods [28–31].It is known that the increasing flexibility caused by the crack depends on the geometry of the crack, so it is important to estimate in the best possible way the actual crack shape. In the case of cracked shafts, frequently, cracks grow following transversal planes to their longitudinal axis acquiring elliptical shape fronts [32–37]. However, at present, in the knowledge of the authors, one cannot find in the literature any expression that allows calculating the relationship between crack parameters and natural frequencies in the case of elliptical cracks. In this paper, natural frequencies of bending vibrations of Euler–Bernoulli cracked beam with circular cross section, taking account the elliptical shape of the crack, have been obtained. The natural frequencies of the damaged beams have been calculated by the perturbation analysis, a method that was developed to calculate the eigenvalue and eigenvector of the damaged structures for decreasing the computation expense [38]. The perturbation techniques have been successfully used before by other authors to calculate the variation of natural frequencies induced by cracks in Euler–Bernoulli beams [8, 39–41].Once the direct problem is known, the inverse problem must be approached. It is usually considered as an optimization problem where the crack parameters are regarded as design variables, so that they are estimated by minimizing the differences between the reference values and the ones predicted by the model. Several optimization methods have been used to solve the inverse problem for crack identification, for example, classical optimization technique of minimizing a least square criterion [8, 42], the minimization of the mean-square and min-max formulations [43], or evolutionary algorithms [44]. Among them, genetic algorithm (GA) optimization techniques have been successfully included in fault detection methods applied to mechanical systems, for instance, in [45–50].In this paper, the perturbation theory has been used to obtain the analytical expressions of the natural frequencies of bending vibrations of nonrotating cracked shafts with circular cross section, taking into account location, size, and shape of the elliptical crack. The crack is represented by a rotational spring connecting the two parts of the beam produced by the crack, whose constant is evaluated from the expression of flexibility coefficient found in [35]. To solve the inverse problem of detecting and identifying the crack, the formulation previously obtained has been used in an optimization genetic algorithm which has been developed within the MatLab environment. Finally, in order to verify the GA, input data for the developed algorithm have been obtained from a modal 3D finite element (FE) model using Commercial Code ABAQUS [51].


## 2. Direct Problem

### 2.1. Model of the Cracked Shaft Beam

A cracked Euler–Bernoulli shaft beam that can vibrate in theX–Y plane is considered (Figure 1). It has length L and longitudinal uniform round cross section with diameter D. The beam contains a transversal crack with elliptical front characterized by semiaxes a and b, shown in Figure 1. It remains always open during vibration, and it is situated in an arbitrary location of abscissa xc. The crack characteristic parameters are as follows:(i)
α=a/D characteristic depth of the crack(ii)
γ=a/b shape factor of the crack (γ=1 corresponds to a semicircular crack and γ=0 corresponds to a straight crack)(iii)
δ=xc/L dimensionless location of the crackFigure 1
Cracked shaft beam.In order to take into account the presence of the crack, the shaft can be modelled with two beams connected by a massless torsional spring. The model leads to a discontinuity in the slope of the beam,Δθ, which is proportional to the bending moment transmitted by the crack section, Mxc:(1)Δθα,γ,xc=Λmα,γMxc,Λmα,γ is the bending flexibility coefficient that depends on the crack characteristics, both depth and shape front. It is given by(2)Λmα,γ=6421−ν2πED3λmα,γ,where E and v are Young’s modulus and Poisson’s ratio of the material and λmα,γ is the nondimensional flexibility coefficient for bending, which can be found in [35].The two beam segments can be treated separately. The vertical displacement of each part of the cracked beam connected by the spring is(3)yjx,t=ujxsinωct,j=1,2,where subscripts j=1 and j=2 mean left and right parts of the cracked beam, respectively, ujx is the transverse deflection, and ωc is the frequency of vibration of the cracked shaft. Both can be calculated solving the following eigenvalue problem:(4)EJd4ujdx4−ωc2ρAuj=0,j=1,2,where J is the moment of inertia, A is the area of the cross section of the beam, and ρ is the mass density of the material.The general solution for Equation (4) leads to a pair of equations one for each part of the cracked beam:(5)u1x=A1sinξcx+B1cosξcx+C1sinhξcx+D1coshξcx,0≤x≤xc,(6)u2x=A2sinξcx+B2cosξcx+C2sinhξcx+D2coshξcx,xc≤x≤L,where the coordinate x is measured from the left support of the beam and(7)ξc2=ωcAρEJ.The coefficientsAk, Bk, Ck, and Dk k=1,2 can be found using the boundary conditions. In this work, pinned-pinned end condition has been considered. The mathematical relations describing the boundary restraints are as follows:(8)u10=0,d2u10dx2=0,u2L=0,d2u2Ldx2=0.Also, four conditions corresponding to the cracked section are used:(i)
Transverse deflection continuity:(9)u1xc=u2xc.(ii)
Slope discontinuity:(10)Δθ=du2xcdx−du1xcdx=EJΛmd2u2xcdx2.(iii)
Bending moment continuity:(11)d2u1xcdx2=d2u2xcdx2.(iv)
Shear force continuity:(12)d3u1xcdx3=d3u2xcdx3.In consonance with other authors [8, 39, 40], the problem has been solved using the perturbation method. According to this method, the solution for the cracked beam can be expressed from the uncracked beam solution as(13)u1x=uox+εϕx+Oε2,(14)u2x=uox+εφx+Oε2,(15)ωc2=ωnc2+εμ+Oε2,where ωnc is the frequency of the harmonic vibration of the uncracked shaft and uox is its transverse deflection, ε is a positive parameter of the same order as Λm and ϕ, and φ and μ are the variables of the problem that have to be obtained as a part of the solution.Substituting (13), (14), and (15) into (4) and keeping the first-order terms only, the natural frequencies of the cracked shaft beam can be derived from the natural frequencies of the uncracked one and the flexibility of the spring that models the crack:(16)ωc2=ωnc2−ΛmEJd2uoxxc/dx22∫0LρAuox2dx.The displacementuo takes different expressions depending on the support conditions of the shaft beam [52]. From (16), the relation between ωc and ωnc (frequency ratio) can be calculated in terms of D, L, and λm and two parameters called g and f that depend on the support restrains:(17)ω∗=ωcωnc=1−641−ν2DλmLg4βf,where β is a value that also depends on the end conditions of the beam [52]. Table 1 shows β values for the first four natural frequencies of the shaft beam for pinned-pinned boundary condition.Table 1
Values ofβ.
Pinned-pinned
β
1
π
β
2
2
π
β
3
3
π
β
4
4
πThe expressions ofg and f until the 4th natural frequency are shown below (subscript i = 1, 2, 3, and 4 means the order of the natural frequency):(18)gi=sinβiδ2,fi=2βi−sin2βi.


### 2.2. Results of the Forward Problem

In this section, the expressions calculated before are applied to a nonrotating shaft beam whose characteristic parameters areL=2m, D=0.05m, E=72000MPa, ν=0.33, and ρ=2800kg/m3.


#### 2.2.1. Effects of the Crack Location

As a representative example of all crack shapes and sizes, for a cracked shaft with relative depthα=0.3 and shape factor γ=0.2, the variations of the first four normalized natural frequencies of the component with different crack locations (δ) are shown in Figure 2.Figure 2
Variations of the first four normalized frequencies in terms of crack position (δ) for α=0.3 and γ=0.2.
(a)
(b)
(c)
(d)As would be expected, the normalized natural frequencies depend on the position of the crack along the beam. If the crack is placed in a mode shape node, the natural frequency of the beam is not affected by the crack, that is, to say,ω∗=1, and also, the normalized natural frequencies decrease with the distance from the crack to the nearby nodes.


#### 2.2.2. Effects of the Crack Depth and Crack Shape

Figure3 shows the variations of the first four normalized natural frequencies of the component in terms of crack depth (α) and position (δ) for γ=0.5, and Figure 4 shows the variations of the first four normalized natural frequencies in terms of crack front shape (γ) and location (δ) for α=0.25; the conclusions are also valid for other values of γ and α, respectively. Both in Figures 3 and 4, axes have been oriented in order to achieve the best view of the represented surfaces.Figure 3
Variations of the first four normalized frequencies in terms of crack depth (α) and position (δ) for γ=0.5.
(a)
(b)
(c)
(d)Figure 4
Variations of the first four normalized frequencies in terms of crack front shape (γ) and position (δ) for α=0.25.
(a)
(b)
(c)
(d)With regard to the crack size influence in the frequency values, for all cases, the natural frequency of the cracked shaft decreases as crack is deeper. If one analyzes the variations of the frequency values depending on the crack shape front, it can be seen that, for all given frequencies, the natural frequencies of the element increase as crack shape is less straight, the largest values appear in the case of circular cracks (γ=1).


## 2.1. Model of the Cracked Shaft Beam

A cracked Euler–Bernoulli shaft beam that can vibrate in theX–Y plane is considered (Figure 1). It has length L and longitudinal uniform round cross section with diameter D. The beam contains a transversal crack with elliptical front characterized by semiaxes a and b, shown in Figure 1. It remains always open during vibration, and it is situated in an arbitrary location of abscissa xc. The crack characteristic parameters are as follows:(i)
α=a/D characteristic depth of the crack(ii)
γ=a/b shape factor of the crack (γ=1 corresponds to a semicircular crack and γ=0 corresponds to a straight crack)(iii)
δ=xc/L dimensionless location of the crackFigure 1
Cracked shaft beam.In order to take into account the presence of the crack, the shaft can be modelled with two beams connected by a massless torsional spring. The model leads to a discontinuity in the slope of the beam,Δθ, which is proportional to the bending moment transmitted by the crack section, Mxc:(1)Δθα,γ,xc=Λmα,γMxc,Λmα,γ is the bending flexibility coefficient that depends on the crack characteristics, both depth and shape front. It is given by(2)Λmα,γ=6421−ν2πED3λmα,γ,where E and v are Young’s modulus and Poisson’s ratio of the material and λmα,γ is the nondimensional flexibility coefficient for bending, which can be found in [35].The two beam segments can be treated separately. The vertical displacement of each part of the cracked beam connected by the spring is(3)yjx,t=ujxsinωct,j=1,2,where subscripts j=1 and j=2 mean left and right parts of the cracked beam, respectively, ujx is the transverse deflection, and ωc is the frequency of vibration of the cracked shaft. Both can be calculated solving the following eigenvalue problem:(4)EJd4ujdx4−ωc2ρAuj=0,j=1,2,where J is the moment of inertia, A is the area of the cross section of the beam, and ρ is the mass density of the material.The general solution for Equation (4) leads to a pair of equations one for each part of the cracked beam:(5)u1x=A1sinξcx+B1cosξcx+C1sinhξcx+D1coshξcx,0≤x≤xc,(6)u2x=A2sinξcx+B2cosξcx+C2sinhξcx+D2coshξcx,xc≤x≤L,where the coordinate x is measured from the left support of the beam and(7)ξc2=ωcAρEJ.The coefficientsAk, Bk, Ck, and Dk k=1,2 can be found using the boundary conditions. In this work, pinned-pinned end condition has been considered. The mathematical relations describing the boundary restraints are as follows:(8)u10=0,d2u10dx2=0,u2L=0,d2u2Ldx2=0.Also, four conditions corresponding to the cracked section are used:(i)
Transverse deflection continuity:(9)u1xc=u2xc.(ii)
Slope discontinuity:(10)Δθ=du2xcdx−du1xcdx=EJΛmd2u2xcdx2.(iii)
Bending moment continuity:(11)d2u1xcdx2=d2u2xcdx2.(iv)
Shear force continuity:(12)d3u1xcdx3=d3u2xcdx3.In consonance with other authors [8, 39, 40], the problem has been solved using the perturbation method. According to this method, the solution for the cracked beam can be expressed from the uncracked beam solution as(13)u1x=uox+εϕx+Oε2,(14)u2x=uox+εφx+Oε2,(15)ωc2=ωnc2+εμ+Oε2,where ωnc is the frequency of the harmonic vibration of the uncracked shaft and uox is its transverse deflection, ε is a positive parameter of the same order as Λm and ϕ, and φ and μ are the variables of the problem that have to be obtained as a part of the solution.Substituting (13), (14), and (15) into (4) and keeping the first-order terms only, the natural frequencies of the cracked shaft beam can be derived from the natural frequencies of the uncracked one and the flexibility of the spring that models the crack:(16)ωc2=ωnc2−ΛmEJd2uoxxc/dx22∫0LρAuox2dx.The displacementuo takes different expressions depending on the support conditions of the shaft beam [52]. From (16), the relation between ωc and ωnc (frequency ratio) can be calculated in terms of D, L, and λm and two parameters called g and f that depend on the support restrains:(17)ω∗=ωcωnc=1−641−ν2DλmLg4βf,where β is a value that also depends on the end conditions of the beam [52]. Table 1 shows β values for the first four natural frequencies of the shaft beam for pinned-pinned boundary condition.Table 1
Values ofβ.
Pinned-pinned
β
1
π
β
2
2
π
β
3
3
π
β
4
4
πThe expressions ofg and f until the 4th natural frequency are shown below (subscript i = 1, 2, 3, and 4 means the order of the natural frequency):(18)gi=sinβiδ2,fi=2βi−sin2βi.


## 2.2. Results of the Forward Problem

In this section, the expressions calculated before are applied to a nonrotating shaft beam whose characteristic parameters areL=2m, D=0.05m, E=72000MPa, ν=0.33, and ρ=2800kg/m3.


### 2.2.1. Effects of the Crack Location

As a representative example of all crack shapes and sizes, for a cracked shaft with relative depthα=0.3 and shape factor γ=0.2, the variations of the first four normalized natural frequencies of the component with different crack locations (δ) are shown in Figure 2.Figure 2
Variations of the first four normalized frequencies in terms of crack position (δ) for α=0.3 and γ=0.2.
(a)
(b)
(c)
(d)As would be expected, the normalized natural frequencies depend on the position of the crack along the beam. If the crack is placed in a mode shape node, the natural frequency of the beam is not affected by the crack, that is, to say,ω∗=1, and also, the normalized natural frequencies decrease with the distance from the crack to the nearby nodes.


### 2.2.2. Effects of the Crack Depth and Crack Shape

Figure3 shows the variations of the first four normalized natural frequencies of the component in terms of crack depth (α) and position (δ) for γ=0.5, and Figure 4 shows the variations of the first four normalized natural frequencies in terms of crack front shape (γ) and location (δ) for α=0.25; the conclusions are also valid for other values of γ and α, respectively. Both in Figures 3 and 4, axes have been oriented in order to achieve the best view of the represented surfaces.Figure 3
Variations of the first four normalized frequencies in terms of crack depth (α) and position (δ) for γ=0.5.
(a)
(b)
(c)
(d)Figure 4
Variations of the first four normalized frequencies in terms of crack front shape (γ) and position (δ) for α=0.25.
(a)
(b)
(c)
(d)With regard to the crack size influence in the frequency values, for all cases, the natural frequency of the cracked shaft decreases as crack is deeper. If one analyzes the variations of the frequency values depending on the crack shape front, it can be seen that, for all given frequencies, the natural frequencies of the element increase as crack shape is less straight, the largest values appear in the case of circular cracks (γ=1).


## 2.2.1. Effects of the Crack Location

As a representative example of all crack shapes and sizes, for a cracked shaft with relative depthα=0.3 and shape factor γ=0.2, the variations of the first four normalized natural frequencies of the component with different crack locations (δ) are shown in Figure 2.Figure 2
Variations of the first four normalized frequencies in terms of crack position (δ) for α=0.3 and γ=0.2.
(a)
(b)
(c)
(d)As would be expected, the normalized natural frequencies depend on the position of the crack along the beam. If the crack is placed in a mode shape node, the natural frequency of the beam is not affected by the crack, that is, to say,ω∗=1, and also, the normalized natural frequencies decrease with the distance from the crack to the nearby nodes.


## 2.2.2. Effects of the Crack Depth and Crack Shape

Figure3 shows the variations of the first four normalized natural frequencies of the component in terms of crack depth (α) and position (δ) for γ=0.5, and Figure 4 shows the variations of the first four normalized natural frequencies in terms of crack front shape (γ) and location (δ) for α=0.25; the conclusions are also valid for other values of γ and α, respectively. Both in Figures 3 and 4, axes have been oriented in order to achieve the best view of the represented surfaces.Figure 3
Variations of the first four normalized frequencies in terms of crack depth (α) and position (δ) for γ=0.5.
(a)
(b)
(c)
(d)Figure 4
Variations of the first four normalized frequencies in terms of crack front shape (γ) and position (δ) for α=0.25.
(a)
(b)
(c)
(d)With regard to the crack size influence in the frequency values, for all cases, the natural frequency of the cracked shaft decreases as crack is deeper. If one analyzes the variations of the frequency values depending on the crack shape front, it can be seen that, for all given frequencies, the natural frequencies of the element increase as crack shape is less straight, the largest values appear in the case of circular cracks (γ=1).


## 3. Inverse Problem

### 3.1. Genetic Algorithm

Once the direct problem is known, the crack detection and identification can be considered as an inverse problem. The direct problem, in this case, has consisted in obtaining the first four natural frequencies from the knowledge of the location and the crack properties. The inverse problem will be then, obtaining the location and the flexibility of the cracked section,δ and λmα,γ, knowing the natural frequencies. In the proposed methodology, the genetic algorithm (GA) optimization technique is applied to solve the inverse problem.The proposed GA is developed under the MatLab environment. To begin the GA, an initial population of 500 individuals are randomly generated. Every individual consists of 2 genes,λm and δ, encoded as a binary string (chromosome). Among the different selection techniques, the roulette wheel method has been chosen, that gives every member a chance of being selected according to each individual’s relative fitness. In relation to the crossover technique, the two-point operator has been employed; in this case, two points of the chromosome are randomly chosen in both individuals, and then the alleles falling in-between the two points are swapped to give two new offspring. On the contrary, the mutation method is used, which randomly creates adaptive rules with respect to the last successful or unsuccessful generation. At last, 85 is the percentage of the population to replace each generation.The fitness function is the most important element of a optimization genetic algorithm. According to the GA toolbox of MatLab, the fitness function,f, should be minimized. It is computed as the root mean-square (RMS) of the differences between the actual first four dimensionless natural frequencies and the computed dimensionless frequencies as follows:(19)fλm,δ=∑i=14ωi,a∗−ωi,c∗21/2,where i denotes the order of the natural frequency and ωi,a∗ and ωi,c∗ represent the actual and the computed values of ratio frequencies, respectively.To decide if the possible optimal values have been calculated, the best fitness value is checked in every step; if the fitness reaches a prefixed level or if the maximum number of iterations is achieved, the algorithm is stopped. As the stopping level, a value of10−6 has been chosen, and to guarantee that the algorithm converges to the optimal solution, each case has been run five times during 250 generations.


### 3.2. Numerical Results

In order to test the proposed algorithm, the same shaft beam previously described has been considered, with cracks in different positions and with different depths and front shapes, that is to say, different flexibilities. The study has been carried out for three couples ofα and γ and three crack locations δ. The three couples α,γ have been chosen considering the fact that the front of incipient cracks is practically circular and that when a crack grows its front becomes straight [53] (Figure 5). In summary, 9 cases have been considered. Table 2 shows the description of the considered cases.Figure 5
Considered cases in numerical results.Table 2
Input data: case description.
α
=
0.15
,
γ
=
0.73
α
=
0.28
,
γ
=
0.55
α
=
0.49
,
γ
=
0.12
(λm=7.04⋅10−4)
(λm=5.58⋅10−3)
(λm=3.74⋅10−2)
δ
=
0.25
Case 1
Case 2
Case 3
δ
=
0.46
Case 4
Case 5
Case 6
δ
=
0.82
Case 7
Case 8
Case 9Input data for the genetic algorithm have been obtained from a modal 3D finite element (FE) model using Commercial Code ABAQUS [51]. The mesh of the three-dimensional model is made by employing 8-node linear brick elements called C3D8, and in order to determine the level of mesh refinement, convergence analysis was carried out. The FE model finally included about 45000 elements and 50000 nodes. In Figure 6, as an example, in the case of healthy beam, both the first four bending vibration modes related to the first four natural frequencies obtained and undeformed beams obtained from the finite element model can be seen.Figure 6
Both first four bending vibration mode shapes and undeformed beams, for a healthy beam, obtained by the finite element model.In Table3, the values of the reference frequency ratios for 9 cases calculated by the FE model can be seen.Table 3
Values of the reference frequency ratios.
Pinned-pinned
ω
1
∗
ω
2
∗
ω
3
∗
ω
4
∗
Case 1
0.9995
0.9989
0.9995
1
Case 2
0.9956
0.9915
0.9959
1
Case 3
0.9717
0.9485
0.9767
0.9999
Case 4
0.9990
0.9999
0.9991
0.9998
Case 5
0.9915
0.9995
0.9927
0.9981
Case 6
0.9468
0.9967
0.9570
0.9889
Case 7
0.9997
0.9991
0.9990
0.9994
Case 8
0.9975
0.9930
0.9918
0.9952
Case 9
0.9835
0.9556
0.9529
0.9746It is important to take into account that it is possible to find that cracked shaft beams with different depths and front shape crack present identical vibration behavior as long as theλm value is the same in all cases. In Figure 7, it can be seen the surface that λm represents values in function of α and γ. In addition, Figure 8 shows several curves that have been calculated as the intersection of the aforementioned function and different planes characterized by constant λm values. As shown, there are infinite values of couples α,γ with which the same value of λm can be obtained.Figure 7
λ
m values in function of α and γ.Figure 8
λ
m values in function of α and γ for constant λm values.After applying the proposed GA, Table4 shows the obtained results, where the estimated values of λm and δ have been calculated as the average of the five repetitions of the same case, and the standard deviation values are given as(20)σx=∑i=1nxi−x¯2n−1,where x can be λm or δ and n is the number of repetitions; in this case, n=5.Table 4
Location and flexibility errors.
σ
λ
m
λ
¯
m
ε
λ
m
%
σ
δ
δ
¯
ε
δ
%
Case 1
5.14
⋅
10
−
7
7.37
⋅
10
−
4
4.7
2.15
⋅
10
−
5
0.249
0.6
Case 2
1.58
⋅
10
−
6
5.96
⋅
10
−
3
6.7
5.58
⋅
10
−
6
0.252
0.8
Case 3
2.90
⋅
10
−
6
3.58
⋅
10
−
2
4.3
6.29
⋅
10
−
6
0.256
2.5
Case 4
6.85
⋅
10
−
6
7.24
⋅
10
−
4
3.2
2.16
⋅
10
−
4
0.461
0.2
Case 5
1.80
⋅
10
−
6
5.95
⋅
10
−
3
6.5
1.24
⋅
10
−
5
0.46
0.0
Case 6
2.26
⋅
10
−
6
3.57
⋅
10
−
2
4.6
2.50
⋅
10
−
6
0.459
0.1
Case 7
1.09
⋅
10
−
5
7.28
⋅
10
−
4
3.4
2.62
⋅
10
−
3
0.819
0.1
Case 8
2.52
⋅
10
−
6
5.90
⋅
10
−
3
5.6
2.64
⋅
10
−
5
0.817
0.3
Case 9
3.05
⋅
10
−
6
3.48
⋅
10
−
2
7.1
3.72
⋅
10
−
6
0.811
1.1In all cases the error is calculated as(21)error%=|estimated value−true value|true value×100.Firstly, in the case of the crack location (δ), the estimation can be considered very good, since the error always keeps under 2.5% and the average error, taking into account all the studied cases, is less than 1%. Secondly, regarding the cracked section flexibility (λm), one can find errors ranging from 3.2% to 7.1%. With respect to flexibility, as it is explained in previous sections, for a given λm value, infinite couples α,γ can be found, so every λm value can be graphed as a curve composed by the corresponding infinite α,γ points. Figure 9 shows, for the studied cases that can be found in Table 2, that given a shape front γ, for every λm constant value, there is a difference between the real dimensionless depth of the crack obtained from λm real and the estimated crack size calculated from λm predicted. In other words, one can see the flexibility curves, both real (λm real) and predicted (λm predicted), as a function of their related infinite couples α,γ. Besides, in Figure 9, the specific points corresponding to the cases in Table 2 have also been highlighted.Figure 9
Flexibility curves. Real and predicted values in function ofα and γ for the studied cases. (a) Cases 1/4/7. (b) Cases 2/5/8. (c) Cases 3/6/9.
(a)
(b)
(c)In Table5, one can see the calculated errors (εα%), according to Equation (20), for the depth crack α corresponding to the previous analyzed cases and for the γ limit values, and these are straight front (γ=0) and circular front (γ=1).Table 5
Errors forα (εα%).
Cases
γ
=
0.73
γ
=
0.55
γ
=
0.12
γ
=
0
γ
=
1
1/4/7
0.33
—
—
2.09
1.54
2/5/8
—
5.46
—
6.04
5.04
3/6/9
—
—
0.67
0.00
0.00According to the results shown in Table5, it can be said that the errors for α estimated values, in the limit cases of straight and circular fronts, are less than the errors associated with flexibility λm for the same shape fronts. Ultimately, the proposed algorithm allows estimating well enough the size of the crack for all the studied boundary conditions since in the worst situations, the errors are under 7%.


## 3.1. Genetic Algorithm

Once the direct problem is known, the crack detection and identification can be considered as an inverse problem. The direct problem, in this case, has consisted in obtaining the first four natural frequencies from the knowledge of the location and the crack properties. The inverse problem will be then, obtaining the location and the flexibility of the cracked section,δ and λmα,γ, knowing the natural frequencies. In the proposed methodology, the genetic algorithm (GA) optimization technique is applied to solve the inverse problem.The proposed GA is developed under the MatLab environment. To begin the GA, an initial population of 500 individuals are randomly generated. Every individual consists of 2 genes,λm and δ, encoded as a binary string (chromosome). Among the different selection techniques, the roulette wheel method has been chosen, that gives every member a chance of being selected according to each individual’s relative fitness. In relation to the crossover technique, the two-point operator has been employed; in this case, two points of the chromosome are randomly chosen in both individuals, and then the alleles falling in-between the two points are swapped to give two new offspring. On the contrary, the mutation method is used, which randomly creates adaptive rules with respect to the last successful or unsuccessful generation. At last, 85 is the percentage of the population to replace each generation.The fitness function is the most important element of a optimization genetic algorithm. According to the GA toolbox of MatLab, the fitness function,f, should be minimized. It is computed as the root mean-square (RMS) of the differences between the actual first four dimensionless natural frequencies and the computed dimensionless frequencies as follows:(19)fλm,δ=∑i=14ωi,a∗−ωi,c∗21/2,where i denotes the order of the natural frequency and ωi,a∗ and ωi,c∗ represent the actual and the computed values of ratio frequencies, respectively.To decide if the possible optimal values have been calculated, the best fitness value is checked in every step; if the fitness reaches a prefixed level or if the maximum number of iterations is achieved, the algorithm is stopped. As the stopping level, a value of10−6 has been chosen, and to guarantee that the algorithm converges to the optimal solution, each case has been run five times during 250 generations.


## 3.2. Numerical Results

In order to test the proposed algorithm, the same shaft beam previously described has been considered, with cracks in different positions and with different depths and front shapes, that is to say, different flexibilities. The study has been carried out for three couples ofα and γ and three crack locations δ. The three couples α,γ have been chosen considering the fact that the front of incipient cracks is practically circular and that when a crack grows its front becomes straight [53] (Figure 5). In summary, 9 cases have been considered. Table 2 shows the description of the considered cases.Figure 5
Considered cases in numerical results.Table 2
Input data: case description.
α
=
0.15
,
γ
=
0.73
α
=
0.28
,
γ
=
0.55
α
=
0.49
,
γ
=
0.12
(λm=7.04⋅10−4)
(λm=5.58⋅10−3)
(λm=3.74⋅10−2)
δ
=
0.25
Case 1
Case 2
Case 3
δ
=
0.46
Case 4
Case 5
Case 6
δ
=
0.82
Case 7
Case 8
Case 9Input data for the genetic algorithm have been obtained from a modal 3D finite element (FE) model using Commercial Code ABAQUS [51]. The mesh of the three-dimensional model is made by employing 8-node linear brick elements called C3D8, and in order to determine the level of mesh refinement, convergence analysis was carried out. The FE model finally included about 45000 elements and 50000 nodes. In Figure 6, as an example, in the case of healthy beam, both the first four bending vibration modes related to the first four natural frequencies obtained and undeformed beams obtained from the finite element model can be seen.Figure 6
Both first four bending vibration mode shapes and undeformed beams, for a healthy beam, obtained by the finite element model.In Table3, the values of the reference frequency ratios for 9 cases calculated by the FE model can be seen.Table 3
Values of the reference frequency ratios.
Pinned-pinned
ω
1
∗
ω
2
∗
ω
3
∗
ω
4
∗
Case 1
0.9995
0.9989
0.9995
1
Case 2
0.9956
0.9915
0.9959
1
Case 3
0.9717
0.9485
0.9767
0.9999
Case 4
0.9990
0.9999
0.9991
0.9998
Case 5
0.9915
0.9995
0.9927
0.9981
Case 6
0.9468
0.9967
0.9570
0.9889
Case 7
0.9997
0.9991
0.9990
0.9994
Case 8
0.9975
0.9930
0.9918
0.9952
Case 9
0.9835
0.9556
0.9529
0.9746It is important to take into account that it is possible to find that cracked shaft beams with different depths and front shape crack present identical vibration behavior as long as theλm value is the same in all cases. In Figure 7, it can be seen the surface that λm represents values in function of α and γ. In addition, Figure 8 shows several curves that have been calculated as the intersection of the aforementioned function and different planes characterized by constant λm values. As shown, there are infinite values of couples α,γ with which the same value of λm can be obtained.Figure 7
λ
m values in function of α and γ.Figure 8
λ
m values in function of α and γ for constant λm values.After applying the proposed GA, Table4 shows the obtained results, where the estimated values of λm and δ have been calculated as the average of the five repetitions of the same case, and the standard deviation values are given as(20)σx=∑i=1nxi−x¯2n−1,where x can be λm or δ and n is the number of repetitions; in this case, n=5.Table 4
Location and flexibility errors.
σ
λ
m
λ
¯
m
ε
λ
m
%
σ
δ
δ
¯
ε
δ
%
Case 1
5.14
⋅
10
−
7
7.37
⋅
10
−
4
4.7
2.15
⋅
10
−
5
0.249
0.6
Case 2
1.58
⋅
10
−
6
5.96
⋅
10
−
3
6.7
5.58
⋅
10
−
6
0.252
0.8
Case 3
2.90
⋅
10
−
6
3.58
⋅
10
−
2
4.3
6.29
⋅
10
−
6
0.256
2.5
Case 4
6.85
⋅
10
−
6
7.24
⋅
10
−
4
3.2
2.16
⋅
10
−
4
0.461
0.2
Case 5
1.80
⋅
10
−
6
5.95
⋅
10
−
3
6.5
1.24
⋅
10
−
5
0.46
0.0
Case 6
2.26
⋅
10
−
6
3.57
⋅
10
−
2
4.6
2.50
⋅
10
−
6
0.459
0.1
Case 7
1.09
⋅
10
−
5
7.28
⋅
10
−
4
3.4
2.62
⋅
10
−
3
0.819
0.1
Case 8
2.52
⋅
10
−
6
5.90
⋅
10
−
3
5.6
2.64
⋅
10
−
5
0.817
0.3
Case 9
3.05
⋅
10
−
6
3.48
⋅
10
−
2
7.1
3.72
⋅
10
−
6
0.811
1.1In all cases the error is calculated as(21)error%=|estimated value−true value|true value×100.Firstly, in the case of the crack location (δ), the estimation can be considered very good, since the error always keeps under 2.5% and the average error, taking into account all the studied cases, is less than 1%. Secondly, regarding the cracked section flexibility (λm), one can find errors ranging from 3.2% to 7.1%. With respect to flexibility, as it is explained in previous sections, for a given λm value, infinite couples α,γ can be found, so every λm value can be graphed as a curve composed by the corresponding infinite α,γ points. Figure 9 shows, for the studied cases that can be found in Table 2, that given a shape front γ, for every λm constant value, there is a difference between the real dimensionless depth of the crack obtained from λm real and the estimated crack size calculated from λm predicted. In other words, one can see the flexibility curves, both real (λm real) and predicted (λm predicted), as a function of their related infinite couples α,γ. Besides, in Figure 9, the specific points corresponding to the cases in Table 2 have also been highlighted.Figure 9
Flexibility curves. Real and predicted values in function ofα and γ for the studied cases. (a) Cases 1/4/7. (b) Cases 2/5/8. (c) Cases 3/6/9.
(a)
(b)
(c)In Table5, one can see the calculated errors (εα%), according to Equation (20), for the depth crack α corresponding to the previous analyzed cases and for the γ limit values, and these are straight front (γ=0) and circular front (γ=1).Table 5
Errors forα (εα%).
Cases
γ
=
0.73
γ
=
0.55
γ
=
0.12
γ
=
0
γ
=
1
1/4/7
0.33
—
—
2.09
1.54
2/5/8
—
5.46
—
6.04
5.04
3/6/9
—
—
0.67
0.00
0.00According to the results shown in Table5, it can be said that the errors for α estimated values, in the limit cases of straight and circular fronts, are less than the errors associated with flexibility λm for the same shape fronts. Ultimately, the proposed algorithm allows estimating well enough the size of the crack for all the studied boundary conditions since in the worst situations, the errors are under 7%.


## 4. Conclusions

Analytical novel expressions available to calculate the first four natural frequencies of an Euler–Bernoulli cracked shaft beam have been presented in this work. These expressions consider that, most of the times, the cracks presented in a shaft have elliptical fronts, so the natural frequencies values of the cracked element depend not only on the location and the size of the crack but also on the shape of the crack front. The method employed is based on linking the two undamaged parts of the shaft by a rotational spring, located in the cracked section, which constantly depends on the crack parameters, that is to say, on the flexibility due to the crack; meanwhile, the formulated problem has been solved by the perturbation technique. Results obtained from the developed expressions show the natural frequency of the cracked shaft decreases as crack is deeper and increases as crack shape is less straight. Once the analytical expressions have been formulated, the inverse problem has been approached. The genetic algorithm technique has been used as an optimization method in order to provide an algorithm that allows to estimate the crack parameters; these are the location and the flexibility due to the crack, which depends on the size and shape of the defect. A finite element model has also been created in order to use its results to check the validity of the proposed algorithm. In regards to the dimensionless crack section flexibility, which depends on the size and shape of the crack, the calculated errors are in a range of values between from3.2% to 7.1%. However, because that it is possible to find that damaged shaft beams with different depths and front shape crack present identical vibration behavior as long as the λm value is the same in all cases, the proposed algorithm allows to estimate well enough the size of the crack since in the worst situations, the errors are under 6.5%, and the average error is 2.3%. Last, in the case of the crack location, the estimation can be considered very good, because the error always keeps under 2.5%, and the average error, taking into account all the studied cases, is less than 1%.



---

*Source: 4623035-2018-10-24.xml*