Nonlinear responses analysis caused by slant crack in a rotor-bearing system

1. Introduction

In rotating machinery, the fatigue shaft crack may be appearing due to manufacturing flaws, corrosive or thermal loading, which can be identified as the main cause of many catastrophic failures. The local crack will introduce local flexibilities, which has a great influence on the vibration behaviors of the rotor system. A large amount of researches on dynamics of cracked rotor has been performed, and work in this area is still continuing. The early researches on crack can date back to the 1970s. Some review papers summarized the main research progress, such as Wauer[1], Gasch [2] and Dimarogonas [3]. Recent researches mainly focused on the crack breathing mechanism, different modeling approaches for the cracked rotor element, cracked shaft vibrations and crack identification methods in [4-6].

Generally, cracks propagate in surfaces with are roughly planar and perpendicular to the rotating axis of the shaft. However, the crack may also propagate along a helicoidal path under large torque loading and these cracks are referred in the literature as helicoidal or slant cracks [5]. Assuming that the slant crack opens and closes synchronously with torsional excitation frequency, Ichimonji et al. [7, 8] firstly analyzed the dynamic characteristics of a simple rotor by a qualitative study. Sekhar and Prasad [9] established a finite element model of a rotor-bearing system for flexural vibrations by including a shaft having a slant crack. In their model, a flexibility matrix for a slant crack and later the stiffness matrix of a slant cracked element were developed. Prabhakar et al.[10] studied vibration characteristics of a slant-cracked rotor passing through its flexural critical speed by using finite element method for flexural vibrations and analyzed the transient response of a cracked rotor by applying an unbalance force and a harmonically varying torque on the rotor. Through the modeling; the dynamic analysis; and detection and monitoring techniques, Sekar et al.[11] compared the vibration behavior of rotors with the transverse crack and slant crack. Darpe[12] presented a simple Jeffcott rotor with a slant crack. In his model, the flexibility matrix of the rotor with slant crack is developed and the stiffness coefficients based on the flexibility values are used in the equations of motion. In another paper [13], Darpe derived a new flexibility matrix for the slant crack that accounts for the additional stress intensity factors due to orientation of the crack compared to the transverse crack and compared the stiffness coefficients and coupled vibration response characteristics between rotor with the slant and transverse crack. Lin and Chu [14, 15] derived four motion equations of two transversal, one torsional and one longitudinal directions of a Jeffcott rotor system with a 45° slant crack on the shaft, and the main stiffnesses and coupling ones are derived based on the relationship between stress intensity factors and strain energy density function. Han et al. [16] established a finite element model of a geared rotor with slant crack and carried out the dynamic analysis of a geared rotor-bearing system with a breathing slant crack is performed. In another paper [17], by using the direct spectral method, Han et al. analyzed the forced responses of a geared rotor system with slant cracked shaft and time-varying mesh stiffness and discussed the effects of straight or slant crack and crack depth on the forced response of the system without and with torsional excitation. Liu et al. [18]compared the dynamic characteristics of rotor with bearing transverse and slant cracks by using finite element method.

A lot of researches on the rotor crack are based on rigid bearing supports, and there is less work reported on the cracked flexible rotor supported fluid film bearings [19-24]. However, the stability and dynamic characteristics of this kind of system are very important as dangerous cracks were mainly reported in rotors running on fluid film bearings. By using the eight rotordynamic coefficients to simulate the effects of journal bearing, Gomez-Mancilla et al. [19] developed an extended Jeffcott rotor on lubricated journal bearings having masses and imbalances at disc and bearings, and analyzed the orbit evolution and vibration patterns at the local resonances. Based on the same simulated method of journal bearing in [19], Sekhar [20] developed equations of motion for transient response and carried out dynamic analysis by considering the effects of fluid film bearings.

Considering the nonlinear oil-film force of journal bearing, many researchers[21-24] studied the dynamic characteristics of cracked rotor system. Meng and Gasch [21] investigated the stability and the stability degree of a flexible cracked rotor supported on different kinds of journal bearings and analyzed the influences of the crack stiffness ratio, the fixed Sommerfeld number, the gravity parameter, and the mass ratio on the system stability. Papadopoulos et al. [22]investigated the dynamic behavior of a rotor-bearing system based a continuous approach, where the shaft rotates on two journal bearings that are simulated as forces acting on the rotor. Yang and Suh[23] studied the various nonlinear responses and dynamic instabilities by a comprehensive rotor model incorporating translational and rotational inertia, bending stiffness, gyroscopic movements, and shear deformation and experiencing slow crack growth and nonlinear bearing film forces. Bachschmid and Dellupi [24]used a model-based identification procedure to identify the nonlinear forces of linearized and nonlinear oil films in two lobe journal bearings. Their presented method utilizes the linear model of the rotor and the nonlinear oil-film forces in the identification procedure.

Most of the previous work focused on the steady-state vibrations of the rotor-bearing system. However, the transient response during run-up and run-down is also important as steady-state response to detect cracks[20, 25-29]. Sekhar [20]analyzed the coast-down phenomenon caused by transverse rotor crack by considering the dissipation through the journal film and by evaluating the deceleration for each speed, and found characteristic sub-critical response peaks when the cracked rotor decelerates through its critical speed. Prabhakar et al. [28]analyzed the vibration characteristics of a slant-cracked rotor passing through its flexural critical speed by using finite element method for flexural vibrations. Sekhar [29] analyzed both start-up and run-down phenomena and compared for the same rotor for the same angular acceleration and deceleration together with the effects of crack, and extract the features of transverse cracks from the time-domain signals of rotor-bearing system by using the continuous wavelet transform (CWT).

Many literatures also showed that the unbalance has a great effect on the vibration response of the cracked rotor[19, 30-32]. Sinou [31] studied the influence of the crack-unbalance interactions and more particularly the relative orientation between the front crack and the unbalance on the vibration responses of the cracked rotor system. Cheng et al.[32] investigated the influence of nonlinear breathing of the crack and the imbalance orientation angle on the stability, critical speed and peak response of the rotor.

Based on the analysis of the above literatures, it is clear that the dynamics of an elastically supported cracked rotor is physically and fundamentally different from that supported using fluid-film bearings [23]. Although the vibration responses of the rotor system supported by fluid-film bearings have been investigated [33-41], the researches on nonlinear dynamic responses caused by crack in a rotor-bearing system supported by fluid bearings are not sufficient. On the basis of [12, 13, 33, 37], this study investigates the nonlinear dynamic characteristics induced by slant crack and linear/nonlinear bearing force in a flexible rotor system attached two disks are analyzed by using finite element method. Moreover, the study also discusses the effects of eccentric phase differences of two disks on the nonlinear responses of the rotor-bearing system under steady-state process (constant rotating speed) and run-up process.

2. Finite element model of a rotor system with a slant crack

2.1. FE model of a rotor segment with slant crack

In order to simplify the modeling process, some assumptions for the rotor crack are made as follows:

(1) A non-propagating slant crack is assumed, and the crack surfaces are assumed to be planar and smooth, and the crack thickness is ignorable, i.e., the gap between two crack surfaces is zero [42].

(2) The material properties of the cracked shaft are assumed to be linearly elastic.

A shaft segment containing a slant crack (see Fig. 1) is modeled by a Timoshenko beam element with six degrees of freedom per node. However, considering the presence of a slant crack, the stiffness matrix of the beam element needs to be modified, which is different with the usual beam element. The coupling phenomena including bending-longitudinal, bending-torsion and longitudinal-torsion have been taken into consideration to modify the stiffness matrix. Then the modified crack shaft element is fitted into the finite element assemblage of the whole rotor system.

A shaft element with a slant crack having a depth of $a$ oriented at an angle of ${\theta }_1$ relative to the axis of the shaft is shown in Fig. 1(b) and the crack center is situated at a distance $x$ from the left end of the element. The element is loaded with shear forces $P_2$, $P_3$ and $P_8$, $P_9$, bending moments $P_5$, $P_6$ and $P_{11}$, $P_{12}$, axial forces $P_1$ and $P_7$, and torsional moments $P_4$ and $P_{10}$ (see Fig. 1(a)). In Fig. 1(c), $R$ is the radius of the shaft; $b=\sqrt{R^2-(R-a{)}^2}/\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_1$ is a half of the crack width; $\alpha =\sqrt{R^2-(\beta \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_1{)}^2}-(R-a)$ is the depth of crack at any distance $\beta$ from the center along the crack edge; $h=2\sqrt{R^2-(\beta \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_1{)}^2}$ is the total height of the strip of width $d\beta$.

Based on fracture mechanics, the total strain energy $U$ is expressed as:

1

$U=U^0+U^c,$

where $U^0$ and $U^c$ are the strain energy of the uncracked shaft element and the additional strain energy due to crack, respectively. The strain energy $U^0$ can be expressed as follows:

2

$U^0=\frac{1}{2}\int \left[\frac{{\alpha }_s{V}_1^2}{GA}+\frac{{\alpha }_s{V}_2^2}{GA}+\frac{M_1^2}{EI}+\frac{M_2^2}{EI}+\frac{T^2}{G{I}_0}+\frac{F^2}{AE}\right]dx,$

where $V_1$ and $V_2$ are the shear forces, $M_1$ and $M_2$ are the bending moments, $T$ is the torsional moment, $F$ is the axial force acting at the crack cross-section, $G$ is the modulus of rigidity, $E$ is the Young’s Modulus, $I$ is the area moment of inertia of the cross-section, $I_0$ is the polar moment of inertia of the cross-section and ${\alpha }_s$ is the shear coefficient.

Fig. 1Shaft finite element: a) the crack element showing forces acting and coordinate system, b) top-view of the slant crack element, c) cross-section of slant crack

Based on the loaded forces and moments acting on the crack element (see Fig. 1), $V_1=P_2$, $V_2=P_3$, $T=P_4$, $F=P_1$, $M_1=P_{2}x-P_6$, and $M_2=P_{3}x+P_5$, the $U^0$ can be rewritten as follow:

3

$U^0=\frac{1}{2}\left[\frac{P_1^{2}l}{AE}+\frac{{\alpha }_s{P}_2^{2}l}{GA}+\frac{P_2^2{l}^3}{3EI}+\frac{{\alpha }_s{P}_3^{2}l}{GA}+\frac{P_3^2{l}^3}{3EI}+\frac{P_4^{2}l}{G{I}_0}+\frac{P_5^{2}l}{EI}+\frac{P_6^{2}l}{EI}-\frac{P_2{P}_6{l}^2}{EI}+\frac{P_3{P}_5{l}^2}{EI}\right].$

The additional strain energy due to the slant crack can be expressed as:

4

$U^c=\frac{1}{E\text{'}}{\int }_A\left[{\left({\sum }_{i=1}^6{K}_i^{\mathrm{{\rm I}}}\right)}^2+{\left({\sum }_{i=1}^6{K}_{{}_i}^{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right)}^2+(1+\upsilon ){\left({\sum }_{i=1}^6{K}_{{}_i}^{\mathrm{{\rm I}}\mathrm{{\rm I}}\mathrm{{\rm I}}}\right)}^2\right]dA,\left(i=\mathrm{1,2},\dots ,6\right),$

where, $E^{\text{'}}=E/1-\upsilon$, $\upsilon$ is the Poisson’s ratio. $K_i^{\mathrm{{\rm I}}}$, $K_i^{\mathrm{{\rm I}}\mathrm{I}}$ and $K_i^{\mathrm{I}\mathrm{I}\mathrm{I}}$represent the stress intensity factors corresponding to the opening, sliding and shearing modes of crack displacement, respectively. Additionally, the detailed expressions of these stress intensity factors can be found in [12].

Assumed that $u_i$ and $P_i$ are the displacement and force, the flexibility coefficients of the slant crack element are expressed as follows:

5

$g_{ij}=\frac{{\partial }^{2}U}{\partial P_i\partial P_j},\mathrm{}\mathrm{}\left(i,j=\mathrm{1,2},\dots ,6\right).$

The flexibility matrix of the slant crack element is given by:

6

$\mathbf{g}=\left[\begin{array}{llllll}{g}_{11}^{}& & & & & \\ g_{21}^{}& g_{22}^{}& & \mathrm{S}\mathrm{y}\mathrm{m}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{c}& & \\ g_{31}^{}& g_{32}^{}& g_{33}^{}& & & \\ g_{41}^{}& g_{42}^{}& g_{43}^{}& g_{44}^{}& & \\ g_{51}^{}& g_{52}^{}& g_{53}^{}& g_{54}^{}& g_{55}^{}& \\ g_{61}^{}& g_{62}^{}& g_{63}^{}& g_{64}^{}& g_{65}^{}& g_{66}^{}\end{array}\right].$

The stiffness matrix of cracked shaft element ${\mathbf{K}}^c$ can be written as:

7

${\mathbf{{\rm K}}}^c=\mathbf{T}{\mathbf{g}}^{-1}{\mathbf{T}}^T,$

where $\mathbf{T}$ is the transformation matrix is given by:

8

${\mathbf{T}}^T=\left[\begin{array}{llllllllllll}1& 0& 0& 0& 0& 0& -1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& -1& 0& 0& 0& l_{ce}\\ 0& 0& 1& 0& 0& 0& 0& 0& -1& 0& -l_{ce}& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& -1\end{array}\right],$

here, $l_{ce}$ is the length of the slant crack element.

2.2. Modeling of breathing behavior of the slant crack

During the operation of the cracked rotor, the crack is not keeping a constant state, but changing at each instant. In practice, part of the slant crack opens and closes in different rotating angles (see Fig. 2), which is called the breathing behavior of the crack. In order to study the flexibility variation with amount of crack opening and closing, Darpe et al. [43] proposed a concept of crack closure line (CCL), as is shown in Fig. 2(a). The crack closure line is an imaginary line perpendicular to the crack edge and it separates the open and closed parts of the crack. The crack edge is divided into 50 points in this study. When the rotor system stays in a rotating state, the position of the CCL keeps changing along the crack edge as the rotor rotates anticlockwise, namely, from 1 to 50 while opening from A to B and from 50 to 100 while closing from A to B. For example, at $\omega =$840 rev/min, the law of the crack breath is shown in Fig. 2b. The figure indicates that the crack breathes from fully closed state in the range of $\theta \in \mathrm{}$[308°, 360°]∪[0°, 42°], and opens gradually from fully closed to fully open state when $\theta \in [$42°, 116$°]$, and thereafter fully opens at $\theta \in$ [116°, 237°], then closes gradually from fully open to fully closed state at $\theta \in [$237°, 308°$]$.

The value of the stress intensity factors ($K^0$) at the crack front is applied to describe the crack breathing behavior, which is used for confirming critical line separating the open or closed part:

9

$K^0={\sum }_{i=1}^6{K}_i^{\mathrm{{\rm I}}}=K_1^{\mathrm{{\rm I}}}+K_2^{\mathrm{{\rm I}}}+K_3^{\mathrm{{\rm I}}}+K_4^{\mathrm{{\rm I}}}+K_5^{\mathrm{{\rm I}}}+K_6^{\mathrm{{\rm I}}}.$

The sign of $K^0$ is employed to indicate the state of crack along the crack edge. If the sign of $K^0$ at any point along the crack edge is negative, it indicates compressive stress at that point, and hence the crack is assumed closed there. If the sign of $K^0$ is positive, it indicates tensile stress and the crack stays in open state.

According to the sign of $K^0$, the crack closure line position can be ascertained. Furthermore, the integration limits responsible for flexible matrix coefficients can also be determined. Hence, the time-varying stiffness matrix ${\mathbf{K}}^r$ in rotating coordinates is obtained. Through using global transformation matrix ${\mathbf{T}}_g$, the ${\mathbf{K}}^r$ is transformed from the rotating coordinate to the ${\mathbf{K}}^s$ in stationary coordinate. ${\mathbf{T}}_g$ is assembled using the elemental coordinate transformation matrix ${\mathbf{T}}_e$, which is shown as:

10

${\mathbf{T}}_e=\left[\begin{array}{llllll}1& 0& 0& 0& 0& 0\\ 0& \mathrm{c}\mathrm{o}\mathrm{s}\theta & \mathrm{s}\mathrm{i}\mathrm{n}\theta & 0& 0& 0\\ 0& -\mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{c}\mathrm{o}\mathrm{s}\theta & 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& \mathrm{c}\mathrm{o}\mathrm{s}\theta & \mathrm{s}\mathrm{i}\mathrm{n}\theta \\ 0& 0& 0& 0& -\mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right],$

where $\theta$ is the angle of rotation:

11

${\mathbf{K}}^s={\mathbf{T}}_g^T{\mathbf{K}}^r{\mathbf{T}}_g.$

Fig. 2Variation of crack closure line position with angular position of rotor

2.3. FE model of the rotor-bearing system

A rotor-bearing system, attached with two identical disks and supported by two bearings, is shown in Fig. 3. In order to study the rotor-bearing system efficiently, the FE model of the rotor-bearing system is simplified according to the following assumptions:

(a) The shaft is divided into 10 Timoshenko beam elements and 11 nodes. Every node has six degrees of freedom as is shown in Fig. 4(a). The slant crack occurs at the fifth element.

(b) The rigid disks are simulated by lumped mass elements which are superimposed upon the corresponding shaft nodes (see Fig. 4(b)). In Fig. 4(b), ${\phi }_1$ is the orientation of unbalance of disk 1 relative to the crack direction ($\xi$) and ${\phi }_2$ is the eccentric phase difference between the unbalance of disk 1 ($m_{d1}e$) and the unbalance of disk 2 ($m_{d2}e$).

These elements are simulated by the mass $m_d$, the diametral and polar moments of inertia ($I_d$, and $I_{pd}$), meanwhile the gyroscopic effects of the disks are also considered.

(c) The left and right bearings are identical and simulated by linear bearing forces and nonlinear oil-film forces presented in [33].

The general displacement vector of a beam element for the shaft ${\mathbf{q}}^e$ can be expressed as:

12

${\mathbf{q}}^e=[x_A\mathrm{}\mathrm{}{y}_A\mathrm{}\mathrm{}{z}_A\mathrm{}{\theta }_{xA}{\theta }_{yA}{\theta }_{zA}{x}_B{y}_B{z}_B{\theta }_{xB}{\theta }_{yB}{\theta }_{zB}{]}^T,$

where the superscript $e$ stands for an element.

The general displacement vector of a rigid disk ${\mathbf{q}}_d^e$ is given as:

13

${\mathbf{q}}_d^e=[x_C\mathrm{}\mathrm{}{y}_C\mathrm{}\mathrm{}{z}_C\mathrm{}{\theta }_{xC}{\theta }_{yC}{\theta }_{zC}{]}^T.$

Fig. 3Schematic diagram and physical dimensions of the simplified rotor-bearing system

Fig. 4FE model schematics of a shaft element and rigid disk

The mass, stiffness and gyroscopic matrices of shaft and disk elements are denoted as ${\mathbf{M}}^e$, ${\mathbf{K}}^e$, ${\mathbf{G}}^e$ and ${\mathbf{M}}_d^e$, ${\mathbf{G}}_d^e$, respectively. The detailed expressions for these matrices can be found in [44].

Considering the effects of the bearing oil-film forces, unbalance exciting force, the rotor gravity, the dynamic equations of the rotor-bearing system with slant crack can be written as:

14

${\mathbf{M}}^{\mathrm{s}}\ddot{\mathbf{q}}+\left(\omega {\mathbf{G}}^{\mathrm{s}}+{\mathbf{C}}^{\mathrm{s}}\right)\dot{\mathbf{q}}+{\mathbf{K}}^s\mathbf{q}={\mathbf{F}}_u+{\mathbf{F}}_b-{\mathbf{F}}_g,$

where ${\mathbf{M}}^s$, ${\mathbf{K}}^s$, ${\mathbf{C}}^s$ and $\omega {\mathbf{G}}^s$ are the mass, stiffness, damping and gyroscopic matrices of the system in stationary coordinates. $\mathbf{q}$ denotes the displacement vector. ${\mathbf{F}}_u$ represents the excitation force vector due to the disk unbalance, here only the disk unbalances are applied at nodes 4 and 7. ${\mathbf{F}}_b$ denotes the oil-film force vectors of the sliding bearings at nodes 1 and 10 (see Fig. 2). ${\mathbf{F}}_g$ denotes the static gravitational force vector in $y$ direction.

The nodal force vectors ${\mathbf{F}}_u$ at nodes 4 (disk 1) and 7 (disk 2) are given as follows:

15

${\mathbf{F}}_u={\left[\begin{array}{c}⋮\\ f_{ux4}\\ f_{uy4}\\ f_{uz4}\\ T_{ux4}\\ M_{uy4}\\ M_{uz4}\\ ⋮\\ f_{ux7}\\ f_{uy7}\\ f_{uz7}\\ T_{ux7}\\ M_{uy7}\\ M_{uz7}\\ ⋮\end{array}\right]}_{60\times 1}={\left[\begin{array}{c}⋮\\ 0\\ m_{d1}e{\omega }^2\mathrm{c}\mathrm{o}\mathrm{s}\theta +m_{d1}e{a}_{ac}\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ m_{d1}e{\omega }^2\mathrm{s}\mathrm{i}\mathrm{n}\theta -m_{d1}e{a}_{ac}\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ 0\\ 0\\ 0\\ ⋮\\ 0\\ m_{d2}e{\omega }^2\mathrm{c}\mathrm{o}\mathrm{s}(\theta +{\phi }_2)+m_{d2}e{a}_{ac}\mathrm{s}\mathrm{i}\mathrm{n}(\theta +{\phi }_2)\\ m_{d2}e{\omega }^2\mathrm{s}\mathrm{i}\mathrm{n}(\theta +{\phi }_2)-m_{d2}e{a}_{ac}\mathrm{c}\mathrm{o}\mathrm{s}(\theta +{\phi }_2)\\ 0\\ 0\\ 0\\ ⋮\end{array}\right]}_{60\times 1},$
$\theta =\left\{\begin{array}{l}{\omega }_{0}t+0.5a_{cc}{t}^2,a_{ac}\ne 0,\\ \omega t,a_{ac}=0,\end{array}\right.$

where $f$, $T$ and $M$ denote the node force, torque and bending moment. The subscripts $ux$, $uy$, $uz$ denote the directions of forces, torque and bending moment and the number subscript denotes node. $\theta$ is the angle of rotation. ${\omega }_0$ is the initial angular velocity, in this paper ${\omega }_0=$10 rad/s.

It should be noted that only stiffness matrix is constantly updated, usually after every degree of rotation. The mass and damping matrices are assumed constant. Rayleigh damping form is applied and obtained by the following formula [37]:

16

$\mathbf{C}={\alpha }_1\mathbf{M}+{\beta }_1\mathbf{K},$

17

${\alpha }_1=\frac{\pi \left({\omega }_{n2}{\xi }_1-{{\omega }_n}_1{\xi }_2\right){{\omega }_n}_1{{\omega }_n}_2}{15({\omega }_{n2}^2-{\omega }_{n1}^2)},{\beta }_1=\frac{60\left({\omega }_{n2}{\xi }_2-{\omega }_{n1}{\xi }_1\right)}{\pi ({\omega }_{n2}^2-{\omega }_{n1}^2)},$

where ${\omega }_{n1}$ and ${\omega }_{n2}$ are the first and second natural frequencies (rev/min); ${\xi }_1$ and ${\xi }_2$ are the first and second modal damping ratios, respectively. In this paper, ${\xi }_1=$0.005 and ${\xi }_2=$0.01.

Two methods for modeling the supporting forces of bearing are proposed, one is linear stiffness-damping model, and the other is nonlinear oil-film force model based on short bearing theory. The specific demonstrations are as follows:

2.3.1. Linear bearing force model

Only the stiffness and damping coefficients terms in horizontal ($z$-coordinate) and vertical ($y$-coordinate) directions are considered, while the cross stiffness and damping terms are neglected. The horizontal and vertical stiffness terms $K_{b1y}\text{,}$$K_{blz}\text{,}$$K_{b2y}\text{,}$$K_{b2z}$ have equal values (2.5×10⁷ N/m), and the horizontal and vertical damping terms $C_{b1y}$, $C_{blz}$, $C_{b2y}$, $C_{b2z}$ have equal values (2.1×10³ N∙s/m). Based on the simulation data above, the first and second bending natural frequencies of the system (${\omega }_1$) are 28.38 and 102.83 Hz, and the first torsional frequency (${\omega }_e$) is 205.08 Hz. In order to verify the model validly, FE model in ANSYS is also established (see Fig. 5). In ANSYS software, the rotor shaft is simulated by using BEAM188 element; the rigid disks are modeled as concentrated masses by using MASS21 element; two identical bearings are modeled by using COMBI214 elements. The results in Fig. 5 indicate that the first three natural frequencies and mode shapes obtained by two models are almost exactly the same.

Fig. 5The comparison of mode shapes in bending and torsional directions between this paper and ANSYS software

2.3.2. Nonlinear bearing force (oil-film force) model

According to the short bearing theory, the nonlinear oil-film force vectors ${\mathbf{F}}_b$ at nodes 1 and 10 are expressed as follows:

18

${\mathit{F}}_b={\left[\begin{array}{l}\begin{array}{c}\begin{array}{c}{F}_{bx1}\\ F_{by1}\end{array}\\ \begin{array}{c}{F}_{bz1}\\ T_{bx1}\\ M_{by1}\\ M_{bz1}\\ ⋮\end{array}\end{array}\\ \begin{array}{l}\begin{array}{c}{F}_{bx10}\\ F_{by10}\end{array}\\ \begin{array}{c}{F}_{bz10}\\ T_{bx10}\\ M_{by10}\\ M_{bz10}\end{array}\end{array}\end{array}\right]}_{60\times 1}=\sigma {\left[\begin{array}{c}\begin{array}{l}\begin{array}{c}0\\ f_{by}\end{array}\\ \begin{array}{c}{f}_{bz}\\ 0\\ 0\\ 0\end{array}\end{array}\\ ⋮\\ \begin{array}{l}\begin{array}{c}0\\ f_{by}\end{array}\\ \begin{array}{c}{f}_{bz}\\ 0\\ 0\\ 0\end{array}\end{array}\end{array}\right]}_{60\times 1},$

where $f_{bz}$ and $f_{by}$ are dimensionless oil-film forces in horizontal and veridical directions, and $\sigma$ is as follow:

19

$\sigma =\eta \omega \frac{D}{2}L{\left(\frac{D}{2c}\right)}^2{\left(\frac{L}{D}\right)}^2,$

here, $\eta$, $L$, $D$ and $c$ are oil viscosity, bearing length, journal diameter and mean radial clearance, respectively. The detailed expressions of $f_{bz}$ and $f_{by}$ can be found in [33].

3. Vibration responses of the rotor system with slant crack under linear and nonlinear bearing forces

Model parameters used in numerical simulation are shown in Table 1. In this section, in order to indicate the effects of the simplified forms of bearing on the vibration responses of the rotor system with slant crack, two simplified bearing models: linear and nonlinear bearing forces are compared under steady-state ($a_{ac}=$0 rad/s²) and transient ($a_{ac}=$10 rad/s²) conditions. The detailed simulation cases are listed in Table 2. The detailed flow chart for vibration response calculation of a rotor-bearing system with slant crack is displayed in Fig. 6.

3.1. Vibration responses in horizontal and torsional directions under linear bearing forces

3.1.1. Case 1: ${\mathit{a}}_{\mathit{a}\mathit{c}}=\mathbf{}$0 rad/s² under linear bearing forces

With the increasing eccentric phase difference of two disks ${\phi }_2$, spectrum cascades of the right journal (node 10) in horizontal and torsional directions under $a_{ac}=$0 rad/s² and $\stackrel{-}{a}=$0.3 are shown in Figs. 7 and 8, which describe the steady-state unbalance response at rotating speed ranging from 360 rev/min to 3840 rev/min in steps of 120 rev/min. For the sake of clarity, some enlarged views from 360 rev/min to 1440 rev/min are also provided. These figures show the following dynamic phenomena:

1) The horizontal vibration responses in Fig. 7 indicate that the amplitude of the rotating frequency component ($f_r$) reduces with the increasing ${\phi }_2$. By the enlarged views, it is clear that the distinctive superharmonic resonance phenomena can be observed at the 1/2nd first bending critical speed $\text{(}{f}_1/2=$ 840 rev/min). Under the out-of-phase eccentricities of two disks (${\phi }_2=$180°), superharmonic resonance can also observed at 600 rev/min, which is close to the 1/3rd first bending critical speed ($f_1/3=$560 rev/min).

Fig. 6Flow chart for vibration response calculation of a rotor-bearing system with slant crack

2) The torsional vibration responses in Fig. 8 indicate that the distinctive superharmonic resonance phenomena appear near 2400 rev/min and 3200 rev/min, which is close to 1/5th and 1/4th first torsional natural frequencies ($f_3/5=$2460 rev/min and $f_3/4=$3075 rev/min). With the increasing ${\phi }_2$, superharmonic resonance at 1/7th and 1/6th first torsional natural frequencies can also observed, such as at ${\phi }_2=$120° and 180° (see Figs. 8(c) and 8(d)).

Table 2A description of simulation

Simplified bearing forms	Cases	Angular acceleration	Crack parameters	Unbalance parameters	Bearing parameters	Figures concerned
Linear bearing force	Case 1	$a_{ac}=$0 rad/s2	$\stackrel{-}{a}=$0.3, ${\phi }_1=$0°	${\phi }_2=$0°, 60°, 120°, 180º	$K_{b1y}=K_{blz}=K_{b2y}=K_{b2z}=$25 MN/m, $C_{b1y}=C_{blz}=C_{b2y}=C_{b2z}=$2.1 kN·s/m	Figs. 7, 8
	Case 2	$a_{ac}=$10 rad/s2	$\stackrel{-}{a}=\text{0.3}$, ${\phi }_1=$0°	${\phi }_2=$0°, 60°, 120°, 180º	$K_{b1y}=K_{blz}=K_{b2y}=K_{b2z}=$25 MN/m, $C_{b1y}=C_{blz}=C_{b2y}=C_{b2z}=$2.1 kN·s/m	Figs. 9-11
Nonlinear bearing force	Case 3	$a_{ac}=$0 rad/s2	$\stackrel{-}{a}=\text{0.3}$, ${\phi }_1=$0°	${\phi }_2=$0°, 60°, 120°, 180°	$\eta =$0.04 Pa·s, $L=$10 mm, $D=$25 mm, $c=$0.03 mm	Figs. 12,13
	Case 4	$a_{ac}=$10 rad/s2	$\stackrel{-}{a}=$0.3, ${\phi }_1=$0°	${\phi }_2=$0°, 60°, 120°, 180°	$\eta =$0.04 Pa·s, $L=$10 mm, $D=$25 mm, $c=$0.03 mm	Figs. 14-16

Fig. 7Spectrum cascades in horizontal direction under linear bearing condition (case 1: aac= 0 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

3.1.2. Case 2: ${\mathit{a}}_{\mathit{a}\mathit{c}}=$ 10 rad/s² under linear bearing forces

Spectrum cascades of the right journal (node 10) in horizontal and torsional directions under $a_{ac}=$10 rad/s² and $\stackrel{-}{a}=$0.3 are indicated in Figs. 9, 10 and 11, which describe the transient unbalance response at rotating speed ranging from 140 rev/min to 3830 rev/min in steps of 90 rev/min. For the sake of clarity, some enlarged views from 340 rev/min to 1490 rev/min are also provided. These figures show the following dynamic phenomena:

1) The horizontal vibration responses in Fig. 9(a) and 10 indicate that angular acceleration has some effect on the first bending critical speed in comparison with the steady-state process ($a_{ac}=$0 rad/s²) because of the tangential inertia force caused by run-up process, and the superharmonic resonance phenomena also appear at the 1/2nd first bending critical speed.

Fig. 8Spectrum cascades in torsional direction under linear bearing condition (case 1: aac= 0 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

Fig. 9Time responses through its critical speed under linear bearing condition (aac= 10 rad/s2, a-= 0.3): a) horizontal direction z, b) torsional direction θx. Note: the subgraphs from left to right are φ2= 0°, 60°, 120° and 180°, respectively

2) The torsional vibration responses in Fig. 9(b) and 11 indicate that the distinctive superharmonic resonance phenomena appear near the 1/5th and 1/4th first torsional natural frequencies. When the rotor vibration amplitude is small, the superharmonic resonances near the 1/6th and 1/7th first torsional natural frequencies are also obvious (see Fig. 11(d)). Moreover, the large torsional amplitude of 2$f_r$ is also a typical feature at ${\phi }_2=$0°, 60° and 120° during run-up process (see Fig. 11(a), 11(b) and 11(c)).

Fig. 10Spectrum cascades in horizontal direction under linear bearing condition (aac= 10 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

Fig. 11Spectrum cascades in torsional direction through its critical speed under linear bearing condition (aac= 10 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

3.2. Vibration responses in horizontal and torsional directions under nonlinear bearing forces

3.2.1. Case 3: ${\mathit{a}}_{\mathit{a}\mathit{c}}=\mathit{}$0 rad/s² under nonlinear bearing forces

Applying the nonlinear bearing force model (Eq. (18)), spectrum cascades of the right journal (node 10) in horizontal and torsional directions under $a_{ac}=$0 rad/s² and $\stackrel{-}{a}=$0.3 are shown in Figs. 12 and 13, which describe the steady-state unbalance response at rotating speed ranging from 360 rev/min to 3840 rev/min in steps of 120 rev/min. Some dynamic phenomena are described as follows:

1) The horizontal vibration responses in Fig. 12 show that the lateral vibration is mainly affected by nonlinear oil-film force. Typical oil-film instability features, such as half-speed whirl and subharmonic components of $f_r$/2 and 3$f_r$/2, can be observed in Fig. 12(a) and 12(b). Superharmonic resonance phenomena disappear compared with those under linear bearing force.

2) The torsional vibration responses in Fig. 13 display that the subharmonic components of $f_r$/2 and 3$f_r$/2 caused by nonlinear oil-film forces still exist. This also demonstrates the lateral-torsional coupling vibration phenomenon caused by slant crack. Moreover, the distinctive superharmonic resonance phenomena related to the first torsional natural frequency can also observed, which is similar to the laws under linear bearing force.

Fig. 12Spectrum cascades in horizontal direction under nonlinear bearing condition (aac= 0 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

3.2.2. Case 4: ${\mathit{a}}_{\mathit{a}\mathit{c}}=\mathit{}$10 rad/s² under nonlinear bearing forces

Spectrum cascades of the right journal (node 10) in horizontal and torsional directions under $a_{ac}=$10 rad/s² and $\stackrel{-}{a}=$0.3 are shown in Figs. 14, 15 and 16, which describe the transient unbalance response. In the figure, $f_{n1}$ denotes the first mode whip frequency. Figs. 14(a) and 15 show that the oil-film features are restrained during run-up due to the tangential inertia forces and multiple frequency components 2$f_r$, 3$f_r$ and 4$f_r$ become dominant. Fig. 14(b) and 16 display that the oil-film features appear and the superharmonic resonance phenomena related to the first torsional natural frequency can also observed.

Fig. 13Spectrum cascades in torsional direction under nonlinear bearing condition (aac= 0 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

Fig. 14Time response through its critical speed under nonlinear bearing condition (aac= 10 rad/s2, a-= 0.3): a) horizontal direction z, b) torsional direction θx. Note: the subgraphs from left to right are φ2= 0°, 60°, 120° and 180°, respectively

4. Conclusions

The nonlinear dynamic characteristics induced by slant crack and linear/nonlinear bearing force in a flexible rotor system attached two disks, are analyzed by using finite element method. The effects of eccentric phase differences of two disks on the nonlinear responses of the rotor-bearing system are analyzed under steady-state process (constant rotating speed) and run-up process. Some conclusions are summarized as follows:

1) The lateral vibration shows that superharmonic resonance phenomenon related to the first bending critical speed can be observed under linear bearing forces. However, it is almost unseeable under nonlinear bearing forces. In the analyzed frequency range, the superharmonic resonance near the 1/2nd first bending critical speed is obvious.

2) Under linear and nonlinear bearing forces, the torsional vibration all indicates that the superharmonic resonance phenomena related to the first torsional natural frequency appear. In the analyzed frequency range, the superharmonic resonances near the 1/5th and 1/4th first torsional natural frequencies are obvious, which can be identified as distinctive features to diagnose the slant crack.

Fig. 15Spectrum cascades in horizontal direction under nonlinear bearing condition (aac= 10 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

Fig. 16Spectrum cascades in torsional direction under nonlinear bearing condition (aac= 10 rad/s2, a-= 0.3): a) φ2= 0°, b) φ2= 60°, c) φ2= 120°, d) φ2= 180°

3) The eccentric phase difference of two disks ${\phi }_2$ and angular acceleration $a_{ac}$ during run-up have some effects on the vibration of the rotor system with slant crack. Large ${\phi }_2$ can decrease the rotor vibration and restrain the oil-film instability. The angular acceleration $a_{ac}$ can restrain the oil-film instability due to the tangential inertia force. Moreover, the large torsional amplitude of the second harmonic frequency (2$f_r$) is also a typical feature during run-up.

In this study, the vibration responses of cracked rotor systems are simulated by numerical simulation. In this practical engineering, the measured crack signals are usually contaminated by noise, in my future work, the effects of noise on the feature extraction of crack fault will be focused by removing the noise signal using wavelet analysis.