Study on the extraction method of transverse open crack’s information

1. Introduction

It is pretty common for large-unit rotor system such as air separation equipment to suffer from fatigue crack problem. Longtime running, over heat and additional moment caused by machining error could all lead to transverse crack of a running rotor system. Further expansion of the crack will cause paroxysmal damage of the rotor system and furthermore, unstable movement and shaft damage of the system. Analyzing the vibration character of the cracked rotor system and extracting the crack information from vibration signal can effectively estimate the position and depth of a crack. Study on the method of extracting crack information has become an important direction of the study on cracked rotor problem.

For years, scholars both home and abroad have made loads of related studies on the complicated dynamic phenomenon of cracked rotor system and made some important achievement. Dimarogonas, Papadopoulos [1], by studying the cracked rotor system, brought forward that the impact of the crack upon the vibration characteristics of rotor lies on the variation of rotor stiffness matrix, and used the crack flexibility function to express that variation. Gounaris, Dimarogonas [2], based on the methods of fracture mechanics, brought forward the crack unit of Euler-Bernoulli beam model. Jibing Zheng and Guang Meng [3] analyzed the dynamic behaviors of the cracked rotor under the impact of nonlinear whirl, especially the bifurcation of system response and chaos characteristics, and found many occasions where period 3 solution becomes other periodic solution as the initial value varies. Darpe, Gupta, etc [4], using the nonlinear respiratory crack model, studied the responses under the influence of the coupling effect of vertical vibration, transverse vibration and torsional vibration of the cracked shaft. This paper established the dynamic model of cracked rotor system by using the crack flexibility matrix and obtained the residual dynamic equation of the cracked rotor system by comparing the dynamic equations of healthy and cracked rotor system. Based on that equation, it deduced the method for determining the position and depth of the open crack, and verified its feasibility through simulation and experiment.

2. Dynamic model of fault rotor system

2.1. Crack model

Suppose that a shaft segment unit, given a certain stiffness, has a transverse crack whose depth is $a$, and force $P_1$-$P_{12}$ applied on either side of it. As shown by Fig. 1(a), $P_1$ and $P_7$ are axial forces; $P_2$, $P_3$, $P_8$, $P_9$ are shear forces; $P_4$, $P_5$, $P_{10}$, $P_{11}$ are bending moments circling $x$ and $y$ axis, and $P_6$ and $P_{12}$ are torques. The section’s semidiameter is $R$, the crack’s width is $\text{2}b$, and the meanings of $\alpha$, $h$ and $dx$ are shown as in Fig. 1(b).

Fig. 1Parameters of crack model

a)

b)

Based on Paris’ equation [5], we can get the additional displacement $u_i$ caused by the crack in $i$ direction:

1

$u_i=\frac{\partial }{\partial P_i}\left[{\int }_0^{a}J\left(\alpha \right)d\alpha \right],$

where $J\left(\alpha \right)$ – strain energy density function, find in reference [6-9] for detail. $P_i$ – relative load in $i$ direction.

The crack flexibility coefficient is (2$b$ in width):

2

$c_{ij}=\frac{\partial u_i}{\partial P_j}=\frac{{\partial }^2}{\partial P_i\partial P_j}\left[{\int }_0^{a}J\left(\alpha \right)d\alpha \right]=\frac{{\partial }^2}{\partial P_i\partial P_j}\left[{\int }_{-b}^b{\int }_0^{a}J\left(\alpha \right)d\alpha dx\right].$

The flexibility coefficients can constitute the crack flexibility matrix, whose order rests with the degree of freedom of the nodes on the shaft segment. When bending vibration is considered solely, the crack flexibility matrix is:

$C_c^e=\left[\begin{array}{cccc}{c}_{22}& 0& 0& 0\\ 0& c_{33}& 0& 0\\ 0& 0& c_{44}& c_{45}\\ 0& 0& c_{54}& c_{55}\end{array}\right].$

The relative flexibility coefficients in matrix $C_c^e$ are:

3

$\left\{\begin{array}{l}{c}_{22}=\frac{1-{\nu }^2}{\pi ER}\times 4{\int }_0^{\stackrel{-}{a}}{\int }_0^{\stackrel{-}{b}}\stackrel{-}{y}{F}_{III}^2\left(\stackrel{-}{h}\right)d\stackrel{-}{x}d\stackrel{-}{y},\\ c_{33}=\frac{1-{\nu }^2}{\pi ER}\times 4{\int }_0^{\stackrel{-}{a}}{\int }_0^{\stackrel{-}{b}}\stackrel{-}{y}{F}_{II}^2\left(\stackrel{-}{h}\right)d\stackrel{-}{x}d\stackrel{-}{y},\\ c_{44}=\frac{1-{\nu }^2}{\pi E{R}^3}\times 32{\int }_0^{\stackrel{-}{a}}{\int }_0^{\stackrel{-}{b}}{\stackrel{-}{x}}^2\stackrel{-}{y}{F}_I^2\left(\stackrel{-}{h}\right)d\stackrel{-}{x}d\stackrel{-}{y},\\ c_{45}=\frac{1-{\nu }^2}{\pi E{R}^3}\times 64{\int }_0^{\stackrel{-}{a}}{\int }_0^{\stackrel{-}{b}}\stackrel{-}{x}\stackrel{-}{y}\sqrt{1-{\stackrel{-}{x}}^2}{F}_1\left(\stackrel{-}{h}\right)F_2\left(\stackrel{-}{h}\right)d\stackrel{-}{x}d\stackrel{-}{y},\\ c_{55}=\frac{1-{\nu }^2}{\pi E{R}^3}\times 64{\int }_0^{\stackrel{-}{a}}{\int }_0^{\stackrel{-}{b}}\stackrel{-}{y}(1-{\stackrel{-}{x}}^2)F_2^2\left(\stackrel{-}{h}\right)d\stackrel{-}{x}d\stackrel{-}{y},\end{array}\right.$

where $v$ is material shearing coefficient, $E$ is the elastic modulus of the shaft segment material. The details of the other coefficients as follows:

$\stackrel{-}{a}=a/R,\stackrel{-}{b}=b/R,\stackrel{-}{y}=y/R,\stackrel{-}{x}=x/R,\stackrel{-}{h}=y/h,$

$F_{\text{1}}\left(\stackrel{-}{h}\right)={\left(\frac{\text{tan}\lambda }{\lambda }\right)}^{\frac{\text{1}}{\text{2}}}\frac{\left[\text{0}.\text{752}+\text{2}.\text{02}\stackrel{-}{h}+\text{0}.\text{37}{\left(\text{1}-\text{sin}\lambda \right)}^{\text{3}}\right]}{\text{cos}\lambda },$

$F_{\text{2}}\left(\stackrel{-}{h}\right)={\left(\frac{\text{tan}\lambda }{\lambda }\right)}^{\frac{\text{1}}{\text{2}}}\frac{\left[\text{0}.\text{923}+\text{0}.\text{199}{\left[\text{1}-\text{sin}\lambda \right]}^{\text{4}}\right]}{\text{cos}\lambda },$

$F_{II}\left(\stackrel{-}{h}\right)=\frac{\left(\text{1}.\text{122}-\text{0}.\text{561}\stackrel{-}{h}+\text{0}.\text{085}{\stackrel{-}{h}}^{\text{2}}+\text{0}.\text{18}{\stackrel{-}{h}}^{\text{3}}\right)}{{\left(\text{1}-\stackrel{-}{h}\right)}^{\frac{\text{1}}{\text{2}}}},$

$F_{III}\left(\stackrel{-}{h}\right)={\left(\frac{\text{tan}\lambda }{\lambda }\right)}^{\frac{\text{1}}{\text{2}}},$

$\lambda =\frac{\pi }{2}\stackrel{-}{h}.$

The influence upon the rotor system by the crack can be denoted by the variation of the stiffness of the cracked shaft segment. Suppose that the flexibility matrix of a non-cracked rotor shaft segment is:

4

${\mathbf{C}}_u^e=\left[\begin{array}{cccc}\frac{l^3}{3EI}& 0& 0& \frac{l^2}{2EI}\\ 0& \frac{l^3}{3EI}& -\frac{l^2}{2EI}& 0\\ 0& -\frac{l^2}{2EI}& \frac{l}{EI}& 0\\ \frac{l^2}{2EI}& 0& 0& \frac{l}{EI}\end{array}\right],$

where $I$ is the section inertia moment of the shaft segment, $E$ is the elastic modulus of the shaft material and $l$ is the length of the shaft segment.

So the whole flexibility matrix of the rotor shaft segment with crack can be expressed as:

5

${\mathbf{C}}^e={\mathbf{C}}_u^e+{\mathbf{C}}_c^e.$

Accordingly, the stiffness matrix of the cracked rotor shaft segment is:

6

${\mathbf{K}}_u^e={\mathbf{T}}_c{\left({\mathbf{C}}^e\right)}^{-1}{\mathbf{T}}_c^T={\mathbf{T}}_c{({\mathbf{C}}_u^e+{\mathbf{C}}_c^e)}^{-1}{\mathbf{T}}_c^T,$

where, ${\mathbf{T}}_c$ is the stiffness transition matrix of the cracked shaft, details of which can be found as follows:

${\mathbf{T}}_c=\left[\begin{array}{cccc}-1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& l& -1& 0\\ -l& 0& 0& -1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right].$

3. Determining the coefficients of the crack position

As the crack only exists at the $R$-th segment, $\mathrm{\Delta }\mathbf{K}$ in detail is as follows:

where ${\mathbf{O}}_{i\times j}$ is a $i$×$j$-dimensioned 0-matrix.

As only $\mathrm{\Delta }{\mathbf{K}}^e$ contains non-zero factors in Eq. (10), the left end of Eq. (9) can transform into:

where ${\mathbf{F}}_c^e$ is the crack-caused additional force at either end of the cracked shaft segment.

As the flexibility factors of the open crack do not alter in rotation, from Eq. (11) we know that ${\mathbf{F}}_c^e$ is in simple harmonic variation. The Eq. (9) can be expressed in complex number:

Suppose $\mathrm{\Delta }\mathbf{q}=\mathrm{\Delta }\mathbf{A}\mathrm{e}\mathrm{x}\mathrm{p}\left(j\omega t\right)$, from Eq. (12) we get:

Based on the stress condition of the two ends of the shaft segment, ${\mathbf{F}}_c^e$ can be expressed by the force at one side of the cracked shaft and moment component [11-15], so Eq. (13) can be changed into:

where $\mathbf{T}$ is shown as follows:

In light of the actual form of ${\mathbf{F}}_c$ in Eq. (11), only the column from 4$R-$3 to 4$R+$4 have practical sense. If the residual vibration quantity of nodes $i$ and $j$ of the shaft segment in $x$ and $y$ direction and node $k$ in either direction ($x$ or $y$) are known, from Eq. (14) we have:

Based on Eqs. (14) and (16), we can derive the expression of the residual vibration quantity of node $k$ in $x$ and $y$ direction ($y$ direction as an example):

where:

Thus the position of the crack can be estimated by the following expression:

where ${\mathrm{\Delta }y}_k^{\text{'}}$ is the response of node k obtained from simulation or experiment.

Based on Eq. (19), go through $R=$1, 2,…, $N$ and get the value of ${\lambda }_R$ of each shaft segment, and the one with the smallest value of ${\lambda }_R$ is where the crack happens.

4. Extraction of the depth of the crack

According to Eq. (3), the information about the depth of the crack $a$ is contained in the flexibility coefficient, thus it is essential to obtain the crack flexibility coefficient of the rotor system to solve the crack depth problem [16-18].

From Eq. (11) we get:

where ${\mathbf{q}}_c^e$ is the response vector at the two ends of the segment where the crack happens, and ${\mathbf{F}}_c^e$ can be obtained from Eq. (16).

From Eq. (20) we have:

Clean up and we get:

where, ${\left({\mathbf{T}}_c\right)}^{-1}={\left({\mathbf{T}}_c^{\text{T}}{\mathbf{T}}_c\right)}^{-1}{\mathbf{T}}_c^T$ so from Eq. (22) we have:

Assume $\left\{\begin{array}{cccc}{a}_1& a_2& a_3& a_4\end{array}\right\}={\left({\mathbf{T}}_c^{\text{T}}{\mathbf{T}}_c\right)}^{-1}{\mathbf{T}}_c^{\text{T}}{\mathbf{F}}_c^e$, $\left\{\begin{array}{cccc}{b}_1& b_2& b_3& b_4\end{array}\right\}={\mathbf{T}}_c^{\text{T}}{\mathbf{q}}_c^e$, so from Eq. (23) we get:

After the position R of the cracked segment is determined, the ${\mathbf{q}}_c^e$ is known, thus $a_1$, $a_2$, $a_3$, $a_4$, $b_1$, $b_2$, $b_3$, $b_4$ can all be derived. The four Eq. (24) have five unknown crack flexibility coefficients. In order that the crack flexibility coefficients can be solved, we take two different rotation frequency $\omega$, solve its corresponding ${\mathbf{q}}_c^e$ by Eq. (8). From Eq. (3) we know that the crack flexibility coefficients do not vary according to rotation frequency $\omega$, thus, given two different frequencies, we have eight equations to solve five unknown numbers. This is an incompatible equations and we can derive the crack flexibility coefficients $c_{22}$, $c_{33}$, $c_{44}$, $c_{45}$ and $c_{55}$ by the method of least squares.

Fig. 2 shows that the crack flexibility coefficient vary as the dimensionless crack depth. We can see from the figure that each crack flexibility coefficient corresponds with dimensionless crack depth $\stackrel{-}{a}$. Substitute the obtained $c_{22}$, $c_{33}$, $c_{44}$, $c_{45}$ and $c_{55}$ into Eq. (3), we have the corresponding dimensionless crack depth ${\stackrel{-}{a}}_{22}$, ${\stackrel{-}{a}}_{33}$, ${\stackrel{-}{a}}_{44}$, ${\stackrel{-}{a}}_{45}$ and ${\stackrel{-}{a}}_{55}$. Then we can calculate the mean of them and get dimensionless crack depth $\stackrel{-}{a}$, which is further used to derive the crack depth $a$.

Fig. 2Crack flexibility coefficient

5. Numeric simulation

Take the single span rotor with a single disk as an example, the cracked rotor system going through finite-element model of discretization is shown by Fig. 3, in which the dots are nodes and the number denotes the segment number. The disk’s diameter is of 60 mm at segment 3; the two bearings at each sides are at segment 1 and 13 and the two supporting bearings are simplified as spring damper.

Coefficient of each shaft segment unit of the rotor system is shown in Table 1.

The elastic modulus of the rotor system is 2.1×10¹¹ Pa, density being 7850 kg/m³. Proportional damping is chosen with its number being 2×10³ N·s/m; stiffness of the bearing is 2×10⁶ N/m; unbalancing value is 156×10^-6 kg·m.

Fig. 3Single span rotor with a single disk model

The open crack’s flexibility coefficient do not vary as the rotation angle does and there is only 1X frequency in its response, so we need Eq. (9), i.e. residual dynamic equation, to determine the position of the crack. Assume that the crack is at the center of the 6th segment and the dimensionless crack depth is 0.5. We first get the displacements of nodes 4 and 11 in $x$ and $y$ direction and that of the node 5 in y direction, and by Eq. (19) we can get the value of ${\lambda }_R$ of each segment. As shown in Fig. 4, the crack happens where the value of ${\lambda }_R$ is the smallest, which indicates that the crack happens at segment 6.

Fig. 4Determination of the position of the crack

After that, by obtaining the responses at the two ends of the cracked shaft segment of the rotor system at the rotating speed of 6000 rpm and 6600 rpm, we can get from Eq. (24) $c_{22}$, $c_{33}$, $c_{44}$, $c_{45}$ and $c_{55}$ and their dimensionless crack depth ${\stackrel{-}{a}}_{22}$, ${\stackrel{-}{a}}_{33}$, ${\stackrel{-}{a}}_{44}$, ${\stackrel{-}{a}}_{45}$ and ${\stackrel{-}{a}}_{55}$, as shown in Table 2.

By averaging ${\stackrel{-}{a}}_{22}$, ${\stackrel{-}{a}}_{33}$, ${\stackrel{-}{a}}_{44}$, ${\stackrel{-}{a}}_{45}$ and ${\stackrel{-}{a}}_{55}$ we get the dimensionless crack depth $\stackrel{-}{a}=$0.508, which differs slightly from the default value 0.5. Referring to the dimensionless form of the Eq. (3), with the radius of the shaft being 5 mm, we can determine that the actual crack depth is 2.54 mm.

6. Test of robustness

The structure of the rotor system is shown in Fig. 3. And the coefficient of each shaft segment unit of the rotor system is shown in Table 1. Assume that the crack was at the center of the 6th segment and the dimensionless crack depth was 0.5 equally. And due to the influence of the noise, assumed that higher harmonic response [19-21] of the measured signal followed normal distribution, its variance was 0.05. Thus, the 20 times simulations were done. And the result is shown in Fig. 5. Can be seen from the diagram, the crack location was diagnosed in the 6th segment by the method every time.

Fig. 5Result of robustness test

And then, by obtaining the responses at the two ends of the cracked shaft segment of the rotor system at the rotating speed of 6000 rpm and 6600 rpm, we can get the average value of $c_{22}$, $c_{33}$, $c_{44}$, $c_{45}$, $c_{55}$ and their dimensionless crack depth ${\stackrel{-}{a}}_{22}$, ${\stackrel{-}{a}}_{33}$, ${\stackrel{-}{a}}_{44}$, ${\stackrel{-}{a}}_{45}$ and ${\stackrel{-}{a}}_{55}$, as shown in Table 3.

So we get the dimensionless crack depth $\stackrel{-}{a}=$0.522. The set value is 0.5. The error is within 5 %. Error rate is small. Can be seen from the above analysis, the method has good robustness.

7. Experiment verification

This experiment performs on the Bently rotor test-bed as shown in Fig. 6. The shaft of the rotor system is 450 mm in length and the coefficients of each segment unit are shown in Table 1. There is a crack of 3 mm in depth at segment unit 5. The displacement sensor is installed at node 5 (135 mm from the initial end of the coupling) and node 11 (355 mm from the initial end of the coupling). Use B&K3560B acquisition instrument to collect, filter, magnify and A/D transfer, the sampling frequency being 3.2 kHz.

Fig. 6Bently rotor test-bed

Collect displacements of node 5 and 11 in $x$ and $y$ direction with rotating speed being 5000 rpm and 5500 rpm respectively. Fig. 7 shows the time-domain waveform, frequency spectrogram and Poincare sectional drawing at node 5 of the cracked rotor system.

Fig. 7Response of the cracked rotor system

a)

b)

c)

We can see that the rotor system performs cycle 1 motion in the Poincare section. And the frequency spectrogram shows that it is mainly working frequency when 2$X$ frequencies seldom appears. Thus we still use the residual vibration quantity equation as the method to extract the crack’s position and use Eq. (19) to obtain the value of ${\lambda }_R$ of each segment. As shown in Fig. 8, the crack appears at segment 5 where ${\lambda }_R$ takes on its smallest value.

Fig. 8Determination of the position of the crack

Then we measure the response of the two ends of segment 5 (node 5 and 6) when the rotating speed is at 5000 rpm and 5500 rpm respectively. By Eq. (24) we get $c_{22}=$0.4059, $c_{33}=$0.5050, $c_{44}=$0.7755, $c_{45}=$2.4056 and $c_{55}=$5.9346, as shown in Table 4.

Combining with Eq. (3) we get the dimensionless crack depth as 0.623, and then we get the actual depth of the crack as 3.115 mm. The extracted position of the crack and its depth (3 mm) are in line with the actual situation.

8. Discussion

The limit of the method:

1) This method needs to obtain the system frequency response function matrix $\mathbf{E}\left(\text{j}\omega \right)={[\mathbf{K}+j\omega (\mathbf{C}+\omega \mathbf{G})-{\omega }^{\text{2}}\mathbf{M}]}^{-\text{1}}$ in advance through the structural parameters or the experimental test.

2) If there is not the steady-state period response caused by strong nonlinear such as quasi-periodic, chaotic and so on in the rotor system, the method is ineffective.

The advantages of this method:

1) This method will be able to identify the location of the crack through the only two positions of vibration response. Due to the sensors can be only placed in the bearing displacement of the large rotating machinery, this method is applicable.

2) This method is not only effective for piecewise linear local nonlinear system, but also for other types of local nonlinear in the engineering. As long as it can be fitted to a polynomial form, this method can be used to identify the local nonlinear position.

9. Conclusion

1) Based on crack flexibility matrix, the dynamic model and residual dynamic model of the cracked rotor system were established in this paper. By analyzing the stress conditions of the two ends of the cracked shaft segment and combining the residual quantity dynamic equation, it provides a method by which information about the open crack’s position and its depth can be extracted from the vibration signals of the rotor system.

2) The information extracting method put forward by this paper can be used to diagnose open crack fault, as of rotor system, cantilever beam or truss structure. The method needs to collect the vibration signals of three nodes at two different rotating speeds, which is easy and feasible.

3) For a real-life industrial environment, this method remains to be tested. But in the theoretical point of view, the method is feasible.