**1. Introduction**

Inter-shaft bearings are widely used in the supporting drive systems of twin-rotor aero engines. They are one of the necessary parts of the areoengine-supporting drive system. Due to the harsh working conditions of aeroengine, faults often occur in inter-shaft bearings.

At present, the dynamic modeling method of rolling bearing with local defects is widely studied. Rolling elements are equivalent to nonlinear springs, and a 2-DOF model is mainly for the dynamic analysis of a transient state during the running of the rolling bearing [1]. At the same time, Gupta [2–9] published many models to describe motion states and force states of each part and to consider the speed changes of each part and the corresponding effect of inertial force. However, the model does not consider the influence of damping, and it only simplifies the collision contact according to the elastic contact. Based on the McFadden model [10], Su [11,12] established the dynamic model of bearing with defects and distribution defects under variable loads and used this model to reveal the frequency characteristics of the two defects. Hu [13] proposed a 5-DOF dynamic model of deep-groove ball bearings. The model theoretically formulated the elastic deformation and nonlinear contact forces of bearings coupling dual rotors. Patel [14] established a 6-DOF dynamic model for deep-groove ball bearings and calculated the vibration response of

**Citation:** Tian, J.; Ai, X.; Zhang, F.; Wang, Z.; Wang, C.; Chen, Y. Dynamic Modeling and Simulation Analysis of Inter-Shaft Bearings with Local Defects Considering Elasto-Hydrodynamic Lubrication. *Coatings* **2022**, *12*, 1735. https:// doi.org/10.3390/coatings12111735

Academic Editors: Andrew J. Pinkerton and Georgios Skordaris

Received: 30 September 2022 Accepted: 8 November 2022 Published: 13 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

single-point and multi-point faults under constant load. Based on the Hertz contact theory, Xi [15], Ma [16], Patel [17], Liu [18] and Cao [19] established bearing dynamic models with different degrees of freedom and analyzed the nonlinear response of rolling bearing with local defects.

Patil [20] considered the rolling element as a nonlinear contact spring and proposed a dynamic model to research the effects of the defect size on the vibration characteristics of the bearing. Based on the Patil model, Kulkarni [21] used cubic Hermite spline interpolation difference functions to simulate the pulse generated by bearing faults. Cui [22] established a nonlinear dynamic model, which used rectangular displacement excitation to simulate local defects. Khanam [23] proposed a dynamic model using a semi-sinusoidal curve to describe the displacement excitation. Wang [24] came up with a multi-body dynamic model to study the vibration response of cylindrical roller bearings with local defects and analyzed the influence of time-varying contact on bearing defects. The above research shows that the dynamic model of rolling bearing based on the Hertz contact theory, combined with the displacement excitation function, can accurately simulate the vibration signals generated by bearing surface faults.

Sassi [25] established a 3-DOF bearing vibration equation considering the influence of an oil film. The interaction between the rollers and the fault area were analyzed. Wang [26] developed a coupling model of the rotor bearings system, considering the oil film. The variation laws of oil films were analyzed based on this model. Yan [27] and Zhang [28] also considered the influence of the lubricating oil film of the bearing in the established dynamic model and proved that it can more accurately describe the real-time state of bearing operation when considering the influence of EHL on the bearing.

In addition to establishing models for rolling bearing with local defects, the researchers also proposed a variety of coupling modeling methods for fault bearing and rotor systems [29–32]. Niu [33] established a dynamic model considering the change in contact force direction of roller sliding and entering defects. Cao [34] proposed a dynamic model method of bearing, considering a bearing pedestal and rotor system. To sum up, research on the dynamic model of rolling bearing fault carried out by experts and scholars mostly focuses on the faults of conventional rolling bearing. There is not much research on the modeling of inter-shaft bearing fault dynamics. Rolling bearing and inter-shaft bearing have both similarities and differences. The modeling method of rolling bearing can be used as the basis of inter-shaft bearing modeling, but it can not completely describe the motion state of inter-shaft bearing. In this paper, the inter-shaft bearing is taken as the main research object. Based on the nonlinear Hertz contact theory, a local defect dynamic model of the inter-shaft bearing, considering the influence of time-varying displacement excitation (TVDE) and elasto-hydrodynamic lubrication (EHL), is proposed. The model is used to simulate and analyze the time–frequency distribution rules of fault signals and the evolution law of micro-faults.

Other parts of this paper are arranged as follow. The fault simulation experiment of co-rotating and counter rotating on the birotor bearing test rig is shown in Part 3. The fault characteristics of inter-shaft bearing with local defects are studied in Part 4, and the time–frequency characteristics of micro faults are also studied in this section. Finally, part 5 summarizes the conclusion of this paper.

#### **2. Establishment of Fault Dynamic Model**

#### *2.1. Subsection Simplification and Assumption of Inter-Shaft Bearings*

Inter-shaft bearing is an important component of aeroengine rotor systems and works between an HP (high-pressure) rotor and an LP (low-pressure) rotor. Unlike ordinary bearings, the inner ring (IR) and outer ring (OR) of bearings rotate at the same time. Depending on the rotation direction, the working speed of the bearing is the speed difference or sum of the two rotors speed. The supporting form of typical dual-rotor aeroengine inter-shaft bearing is shown in Figure 1.

**Figure 1.** Supporting form of inter-shaft bearing.

The IR and OR of the inter-shaft bearing rotate at the same time, so it is considered to generate vibration frequently during operation. Considering that both the IR and OR have vertical and horizontal vibrations, a 4-DOF dynamic model is established in this paper. The model assumes that the rollers do not slip. Referring to the Patil model [20], the contact between rollers and rings is simplified as a nonlinear spring-mass system. The model is shown in Figure 2.

**Figure 2.** Simplified model of inter-shaft bearing.

#### *2.2. Nonlinear Hertz Contact Forces*

On the basis of the nonlinear Hertz contact theory [35], Harris [36] derived the relationship between the nonlinear load and the displacement of the rolling bearing as shown in Equation (1).

$$
\zeta = \mathbb{K} \delta^n \tag{1}
$$

where *K* is the load deformation coefficient. *δ* is the radial displacement. Cylindrical roller bearing *n* = 10/9.

*F*

Since the rollers in the inter-shaft bearing has contact with both the IR and OR, the total contact stiffness is as follows.

$$K = \left[\frac{1}{\left(1/K\_i\right)^{1/n} + \left(1/K\_0\right)^{1/n}}\right]^n \tag{2}$$

where *Ki* is the contact stiffness of the IR. *K*0 is the contact stiffness of the OR. For cylindrical roller bearing, *K* = 8.06 × 10<sup>4</sup> × *l*8/9, where *l* is the roller length.

*δj*, the radial deformation of the *j* roller is

$$\delta\_j = (\mathbf{x}\_0 - \mathbf{x}\_i)\cos\theta\_j + (y\_0 - y\_i)\sin\theta\_j - \mathbf{C}\_r - H\_d \tag{3}$$

where *x*0 and *y*0 are the displacement of the center of the OR along the x-axis and y-axis. *xi* and *yi* is the displacement of the center of the IR along the x-axis and the y-axis. *θj* is the angle between the center of the j*th* roller and the x-axis. *Cr* is the initial clearance between the roller and the raceway. *Hd*is the TVDE when the roller passes through the defect.

$$
\theta\_j = \omega\_\text{c} t + \frac{2\pi(j-1)}{Z} + \theta\_{t0} \tag{4}
$$

$$
\omega\_{\varepsilon} = \frac{1}{2} \left[ \omega\_i (1 - \frac{d}{D\_m} \cos a) + \omega\_o (1 + \frac{d}{D\_m} \cos a) \right] \tag{5}
$$

where *ωi* is the angular velocity of IR of inter-shaft bearing, *ω*o is the angular velocity of OR of inter-shaft bearing, *ωc* is the angular velocity of inter-shaft bearing retainer, *θt*0 is the initial Angle of the first roller relative to the axis, *Z* is the number of roller of inter-shaft bearing, *d* is the roller diameter, *Dm* is the pitch diameter of inter-shaft bearing, and *α* is the bearing pressure angle.

Substitute Equations (2) and (3) into Equation (2) to obtain the nonlinear Hertz contact force of the *j* roller:

$$F\_{\hat{l}} = K \left[ (\mathbf{x}\_0 - \mathbf{x}\_i) \cos \theta\_{\hat{l}} + (y\_0 - y\_i) \sin \theta\_{\hat{l}} - G\_r - H\_d \right]^{10/9} \tag{6}$$

The total nonlinear Hertz contact force received by the inter-shaft bearing is the sum of the nonlinear Hertz contact forces of all the rollers. The components of the total contact force on the x and y axes are:

$$\begin{cases} F\chi = K \sum\_{j=1}^{Z} \left[ (\mathbf{x}\_0 - \mathbf{x}\_i) \cos \theta\_j + (y\_0 - y\_i) \sin \theta\_j - \mathbf{C}\_r - H\_d \right]^{10/9} \cos \theta\_j\\ F\_Y = K \sum\_{j=1}^{Z} \left[ (\mathbf{x}\_0 - \mathbf{x}\_i) \cos \theta\_j + (y\_0 - y\_i) \sin \theta\_j - \mathbf{C}\_r - H\_d \right]^{10/9} \sin \theta\_j \end{cases} \tag{7}$$

#### *2.3. Time-Varying Displacement Excitation*

This paper assumes that the defect of the inter-shaft bearing runs through the entire raceway. Therefore, the defect length is larger than the roller, but the width of the fault is far less than the diameter of the roller. When the roller passes through the defect of the raceway, the roller will not contact the bottom of the defect, and the falling displacement of the roller center will be less than the defect depth.

The relationship between and bearing movement is shown in Figure 3. In the figure, *Hd* is the time-varying displacement when the roller passes through the defect, *H* is the defect depth, *B* is width of the defect, and *He* is the maximum deviation of the roller center when the roller passes through the defect. It can be seen from Figure 4 that the roller can be divided into three stages when it passes through the defect. The first stage is that the roller starts to enter the defect and completely enters the defect. At this time, the time-varying displacement *Hd* increases gradually as the roller enters the defect; in the meantime, the roller always contacts the left side of the defect. In stage 2, the roller completely enters the defect. At this time, the roller contacts the left and right sides of the defect, and *Hd* reaches the peak value of *He*. In stage 3, the roller enters the defect completely and leaves the defect just after. At this time, the roller only contacts the right side of the defect, and *Hd* gradually decreases. To sum up, while the roller passes through the defect, *Hd* increases first and then decreases with motion process. Therefore, the semi-sinusoidal

function is used to describe the TVDE of the local defect to the roller. The maximum displacement excitation *He* of the roller is:

$$H\_{\mathfrak{k}} = d/2 - \sqrt{(d/2)^2 - (B/2)^2} \tag{8}$$

**Figure 3.** TVDE of the roller.

**Figure 4.** Roller passing defect state diagram.

#### 2.3.1. The TVDE of Defects in the OR

Figure 5 reveals the relationship between the rollers and the defect angle position when the IR and OR of the bearing rotation in the opposite direction. The angular velocity of the IR *wi* is clockwise, and the angular velocity of the OR *w*0 is counterclockwise. When *wi* is smaller than *w*0, the angle velocity of the cage *wc* is counterclockwise at this time, and *ϕf* is the defect angle. The counterclockwise direction is specified in the figure as the positive direction. *θb*0 is the angle of the first roller relative to the X-axis at the initial moment. At moment *t*, *θbi* is the rotational angle of the *ith* roller. *θ f* 0 is the angle relative to the X-axis at the initial moment of the defect. *θ f* is the rotational angle of the defect. The TVDE *Hd* generated when the roller passes through the defect is:

$$H\_d = \begin{cases} H\_c \sin\left(\frac{\pi}{\theta\_f^\*} (\text{mod}(\theta\_f, 2\pi) - \text{mod}(\theta\_{bi}, 2\pi))\right) \\ 0 \le \text{mod}(\theta\_f, 2\pi) - \text{mod}(\theta\_{bi}, 2\pi) \le \varphi\_f \\\\ 0 \qquad \qquad \qquad \qquad \end{cases} \tag{9}$$

where mod is the remainders. The expressions for *θbi* and *θ f* are:

$$
\theta\_{bi} = \frac{2\pi}{Z}(i - 1) + w\_c t + \theta\_{b0}(i = 1 \dots Z) \tag{10}
$$

$$
\theta f = w\_0 t + \theta\_{f0} \tag{11}
$$

**Figure 5.** Roller and defect angular position of OR fault.

#### 2.3.2. The TVDE with Defects in the IR

Figure 6 shows the relationship of rollers and defects when the IR and OR rotate in the opposite direction. When *wi* is smaller than *w*0, both *wc* and *w*0 are in the same direction. Counterclockwise is the positive direction. The expression of TVDE *Hd* is:

$$H\_d = \begin{cases} H\_\varepsilon \sin\left(\frac{\pi}{\theta\_f} \left( \text{mod}(\theta\_{bi}, 2\pi) - \text{mod}\left(\theta\_{f0} + 2\pi - \text{mod}(w\_i t, 2\pi), 2\pi\right) \right) \right) \\\quad 0 \le \text{mod}(\theta\_{bi}, 2\pi) - \text{mod}\left(\theta\_{f0} + 2\pi - \text{mod}(w\_i t, 2\pi), 2\pi\right) \le \rho\_f \\\\ 0 & \text{else} \end{cases} \tag{12}$$

$$\theta\_{b\bar{i}} = \frac{2\pi}{Z}(\bar{i} - 1) + w\_{\bar{c}}t + \theta\_{b0}(\bar{i} = 1 \dots Z) \tag{13}$$

**Figure 6.** Roller and defect angular position of IR fault.

2.3.3. TVDE of Rolling Element Fault

Figure 7 shows the relationship between rollers and the defect angle position when the IR and OR rotate in reverse. Set *wi* to be clockwise rotation and *w*0 to be counterclockwise rotation. When *w*0 is greater than *wi*, *wc* is counterclockwise, and the angle velocity of the roller *wb* is counterclockwise. The counterclockwise direction is the positive direction. Assume the roller 1 has a local defect fault. The defect angle is *ϕb f* . rack the motion of the first roller and establish another moving coordinate system on axis *om*. The counterclockwise direction is taken as the positive direction, and *θb f* 0 is the initial angle of right defect relative to the *om*-axis. At moment t, the angle of the relative moving *om*-axis is *θb f* , and the TVDE *Hd* of the bearing is:

$$H\_d = \begin{cases} \begin{array}{l} H\_\varepsilon \sin\left(\frac{\pi}{\theta\_{bf}} \text{mod}(2\pi + w\_b t - \theta\_{bf0}, 2\pi)\right) \\ 0 \le \text{mod}(2\pi + w\_b t - \theta\_{bf0}, 2\pi) \le \varphi\_{bf} \\\ H\_\varepsilon \sin\left(\frac{\pi}{\theta\_{bf}} (\text{mod}(2\pi + w\_b t - \theta\_{bf0}, 2\pi) - \pi)\right) \\ \pi \le \text{mod}(2\pi + w\_b t - \theta\_{bf0}, 2\pi) \le \pi + \varphi\_{bf} \end{array} \tag{14}$$
 
$$0 \qquad \qquad \qquad 0 \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad 0$$

**Figure 7.** Roller and defect angular position of roller fault.

#### *2.4. Dynamic Model for Inter-Shaft Bearings with Local Defects*

Based on the hypothesis theory of rigid rings and considering the eccentric load on the bearing, the inter-shaft bearing is established as a 4-DOF dynamic model. The contact stiffness, damping, nonlinear Hertz contact force, eccentric load, time-varying displacement and constant radial load are substituted into the dynamic equation to establish the dynamic model of an inter-shaft bearing fault.

$$\begin{cases} M\_0 \ddot{\mathbf{x}}\_0 + c \dot{\mathbf{x}}\_0 + K \sum\_{j=1}^Z \lambda \delta^{10/9} \cos \theta\_j = W + F\_{l0} \cos(\omega\_0 t) \\\ M\_0 \ddot{\mathbf{y}}\_0 + c \dot{\mathbf{y}}\_0 + K \sum\_{j=1}^Z \lambda \delta^{10/9} \sin \theta\_j = F\_{l0} \sin(\omega\_0 t) \\\ M\_i \ddot{\mathbf{x}}\_i + c \dot{\mathbf{x}}\_i - K \sum\_{j=1}^Z \lambda \delta^{10/9} \cos \theta\_j = W + F\_{li} \cos(\omega\_i t) \\\ M\_i \ddot{\mathbf{y}}\_i + c \dot{\mathbf{y}}\_i - K \sum\_{j=1}^Z \lambda \delta^{10/9} \sin \theta\_j = F\_{li} \sin(\omega\_i t) \end{cases} \tag{15}$$

where *Mi* and *M*0 is the quality of IR and OR; *Fli* is the eccentric load applied to the IR; *c* is the damping coefficient; *Fl*0 is the eccentric load applied to the OR. *W* is the radial load perpendicular to inner surface of IR; *λ* is the switch signal to indicate if the roller and the raceway are in contact, expressed as:

$$
\lambda = \begin{cases} 1 & \delta > 0 \\ 0 & \delta \le 0 \end{cases} \tag{16}
$$

The dynamic differential equation is solved by the Newmark-*β* method [37]. Taking *γ* = 1/2 and *β* = 1/4, this method has second-order accuracy and is an unconditionally stable method, which can accurately solve the fault dynamic model of inter-shaft bearing.

#### *2.5. Fault Dynamic Model Considering the Influence of EHL*

The inter-shaft bearing is running at a high speed between high-pressure rotors and low-pressure rotors. Therefore, the inter-shaft bearing needs real-time lubrication. There is an oil film between the rollers and the raceway, and it is always in the state of EHL under the dynamic pressure. The distribution of oil-film pressure in the inter-shaft bearing is shown in Figure 8. The oil-film thickness under the EHL state affects the tribological and dynamic characteristics of the friction pair. The oil pressure generated by the high-speed rotation of the bearing makes the oil film produce a "stiffening effect". The EHL pressure curve is similar to the Hertz pressure curve in distribution. Therefore, in the dynamic model considering the influence of EHL established in this paper, the stiffness of bearing is considered as the series stiffness of oil-film stiffness and contact stiffness [37,38].

**Figure 8.** Elasto-hydrodynamic pressure distribution of inter-shaft bearing.

The paper assumed that the lubricating oil is at a constant temperature and that there is no end leakage, and the contact points of rollers and the rings do not slip. The Dowson– Higginson line-contact film-thickness formula is used to calculate the oil-film thickness between the roller and the IR and OR [39] as follows,

$$h\_{\rm min} = \frac{2.56 a^{0.54} (\eta\_0 \ell I)^{0.7} R^{0.43} L^{0.13}}{E^{\prime 0.03} W^{0.13}} \tag{17}$$

where *α* is the viscosity coefficient of lubricating oil; *η*0 is the dynamic viscosity of lubricating oil; *U* is the average linear velocity of roller; *R* is the equivalent radius of curvature; *L* is the line contact length of roller; *W* is the radial load; *E* is the composite modulus of elasticity, *E* = *E* 1−*ν*<sup>2</sup> .

The inter-shaft bearing is characterized by the same rotation direction of the IR and OR. Assume the rotation speed of the IR is *ni*, the radius is *Ri*, the rotation speed of the OR is *no*, and the roller radius is *r*. By defining the constant *γ* = *dDm* cos *α*, the average linear velocity *U* between the IR and OR of the roller can be obtained according to Equation (5):

$$\Delta U = \frac{\pi}{120} D\_m [n\_i(1-\gamma) + n\_0(1+\gamma)] \tag{18}$$

Since the inter-shaft bearing pressure angle *α* = 0, then

$$\gamma = \frac{d}{D\_{\text{m}}} \cos \kappa = \frac{r}{R\_i + r} \tag{19}$$

The equivalent radius is:

$$R = \left(\frac{1}{r} \pm \frac{1}{R\_i + r \mp r}\right)^{-1} \tag{20}$$

where "−" is the contact equivalent radius of the IR and the roller; "+" is the contact equivalent radius of the OR and the roller.

Substitute Equation (19) into Equation (20), then:

$$R = r(1 \mp \gamma) \tag{21}$$

By substituting each parameter into Equation (17), the minimum oil-film thickness *hi* and *h*0 of the rollers and the IR and OR can be obtained.

$$\begin{cases} h\_i = 0.21a^{0.54} (\eta\_0 D\_m)^{0.7} [n\_i(1-\gamma) + n\_0(1+\gamma)]^{0.7} \\ \quad \ast r^{0.43} (1-\gamma)^{0.43} L^{0.13} E'^{-0.03} W^{-0.13} \\ h\_o = 0.21a^{0.54} (\eta\_0 D\_m)^{0.7} [n\_i(1-\gamma) + n\_0(1+\gamma)]^{0.7} \\ \quad \ast r^{0.43} (1+\gamma)^{0.43} L^{0.13} E'^{-0.03} W^{-0.13} \end{cases} \tag{22}$$

The total thickness of the oil film is:

$$h = h\_i + h\_0 = C\_i \mathcal{W}^{-0.13} + C\_0 \mathcal{W}^{-0.13} \tag{23}$$

where *Ci* is the coefficient in front of *W*−0.13 in *hi*; *C*0 is the coefficient in front of *W*−0.13 in *h*0.

The oil-film stiffness of the rollers and IR and OR of the inter-shaft bearing can be obtained by Equation (23):

$$K\_H = \frac{\text{d}W}{\text{d}h} = \left(0.13(\text{C}\_i + \text{C}\_0)\mathcal{W}^{-1.13}\right)^{-1} \tag{24}$$

According to the calculation formula of series stiffness, the total stiffness of the rolling bearing considering the influence of EHL is:

$$K' = \frac{1}{1/K + 1/K\_H} = \frac{KK\_H}{K + K\_H} \tag{25}$$

#### **3. Experimental Validation of the Numerical Model**
