A fast and reliable numerical method for analyzing loaded rolling element bearing displacements and stiffness

1. Introduction

The gear pair assembly has been considered as one of the major noise and vibration sources in the rotating machineries typically seen in automotive, aerospace and industrial applications [1-4]. The primary excitation force of gear pair vibration is the dynamic mesh force of engaged gear teeth caused by the transmission error (due to tooth profile, spacing error, and elastic deformation of the gear pair). The resultant vibration can be subsequently transmitted through the shaft-bearing system and excite the vibration of housing that radiates annoying noises. In order to design a reliable and quiet power transmission system, and/or trouble-shoot the noise and vibration issues, it is highly desirable to perform an understanding of the behavior of bearings and their interactions with the internal and housing components. In modeling the geared rotor systems, the bearing effect can be easily incorporated into the system model by introducing the relevant bearing stiffness set [5-7]. However, most of the rolling element bearings are precision elements with very complex components that inherently have nonlinear static/dynamic characteristics. Moreover, when a poorly designed bearing is considered, the load-displacement relation might become extremely complex and would cause vibrations. Hence, a fast and reliable bearing stiffness estimation method is necessary to facilitate the static or dynamic analysis of the rotating mechanical systems.

There have been numerous research efforts driven to determine the stiffness matrix of the rolling element bearings. Some of the earlier studies of rolling element bearings were performed by Jones [8], Harris [9], and Palmgern [10]. They investigated the radial and axial load-deflection relation using a nonlinear stiffness coefficient. Later, Gargiulo [11] provides an empirical formulae for radial and axial load-stiffness and deflection-stiffness relations by assuming rigid bearing races. The simplified bearing stiffness matrix obtained through these early studies is either based ideal boundary conditions assumption or neglecting certain degree of freedom (DOF). It was reported by Lim and Singh [12-14] that these early bearing stiffness formulations could not well represent the real bearing characteristics, especially for the coupled translational and rotational vibration. They proposed a more general bearing stiffness matrix with complete 5 DOF terms for ball and roller bearing elements. In their model, a discrete summation approach was adopted to obtain the total bearing forces and moments of all the loaded rolling elements. Hence, a set of nonlinear algebraic equations were formulated for the bearing stiffness, which was numerically solved by using Newton-Raphson method. Then, Hernot et al. [15] derived the stiffness matrix of a five-DOF (degree-of-freedom) angular contact ball bearing by using an analytical approach in which the load summation over ball elements is replaced by an integration. Similarly, the classic Newton-Raphson approach was applied to solve the whole matrix equations. Instead of using the analytical approaches, very recently, Guo and Parker [16] developed a finite element/contact mechanics model to obtain the bearing stiffness matrix for a wide range of bearing types and parameters. The accuracy of the results depends on contact control parameters and step size selected for the finite difference formulation. However, this method is very time-consuming. There are also several experimental techniques [17, 18] hat have been recently developed to determine the bearing stiffness matrix. One major concern of these FEM methods and experimental approaches is the efficiency especially from system design point of view. On the other hand, the general approach proposed by Lim and Singh [12] has demonstrated its efficiency and accuracy. Their theory has been widely adopted in general geared rotor dynamic analysis, such as spur and hypoid geared rotor systems [7, 19-22]. Very recently, Liew and Lim [23] extended the prior study to establish the time-varying stiffness formulation considering the orbital motion of the rolling elements. The time-varying bearing stiffness model has been implemented into both linear parallel and nonlinear non-parallel geared rotor system dynamics analysis [21].

In spite of these successes of obtaining stiffness matrix for bearing elements, one of the major concerns is the computational efficiency and robustness of these numerical models. The coefficients of stiffness matrix can be derived from the partial derivatives of the load expressions with respect to the displacements including translational and rotational coordinates of the inner ring of rolling element bearing relative to the housing support. In this formulation, the coefficients of stiffness matrix can be directly computed given a set of bearing displacement vectors. Yet, in practice, only the external forces applied to bearing system are known. In such cases, the bearing displacement vectors due to known bearing forces can be obtained by solving the system of nonlinear algebraic equations describing the load-displacement relation of the bearing system. For the solution to be reliable, an appropriate numerical method must be chosen to solve these nonlinear equations iteratively. However, the commonly used numerical methods, such as the Newton-Raphson and Powell’s hybrid methods [24], require careful application and can be cumbersome. Also, the accuracy of these algorithms typically relies on the trial and error of different initial estimates since the number of the numerical solutions is not known in advance. The classic Newton-Raphson method is a local minimal approach and the iterative solution tends to trap into the local point. Obviously, if the initial solution was chosen far away from the exact solution point, the whole numerical scheme will be very time-consuming.

In this paper, a reliable, fast and efficient energy method is developed to better quantify the relations of bearing displacements and applied forces. The proposed energy method is based on the principle of minimum potential energy derived from classical mechanics, which is a global searching method. The exact displacements can be found by searching for the displacements that yield the minimum total potential energy of the bearing system. This algorithm overcomes the deficiencies seen in classical iterative method requiring the trial-and-error of different initial estimates. Also, the trap of local minimum solution can be avoided especially for the existence of multiple numerical solutions typically seen in system of nonlinear algebraic equations. In addition, the proposed method can be used to guide the design of the bearings, especially for minimizing the vibrations of a poorly designed bearing. The minimized potential energy condition would give the desired system parameters to minimize the vibrations.

2. Bearing load-displacement relations

The schematic diagram for the relations between bearing forces and displacements for ball and roller bearings are shown in Fig. 1 and Fig. 2, respectively, detailed load-displacement relations were proposed by Lim and Singh [12, 14]. To recap, these relations are given briefly in this paper. In the proposed formulation, the load-displacement relations can be derived by considering the relations of (i) the displacements of inner ring and the deformation of outer raceway-ball/roller element-inner raceway, (ii) the load and the deformation for outer raceway-ball/roller element-inner raceway, and (iii) the normal loads on all ball/roller elements and bearing forces and moments.

Fig. 1Ball element bearing kinematics and co-ordinate system [12, 14]

a)

b)

c)

As seen from Fig. 1(a), for the ball bearing subjected to forces ($F_{bx}$, $F_{by}$, $F_{bz}$) and moments ($M_{bx}$, $M_{by}$), the resultant translational and rotational displacements generated in the bearing inner are (${\delta }_{bx}$, ${\delta }_{by}$, ${\delta }_{bz}$) and (${\beta }_{bx}$, ${\beta }_{by}$), respectively. Also, the ball and raceways will be displaced. Detailed schematic diagram of the deflection for the $j$th ball element and raceway located at angle ${\psi }_j$ from the $x$-axis is shown in Fig. 1(b). Here, $a_o$ is the position of the outer groove curvature centers, $a_i$ and $a_i^{\text{'}}$ are the initial and final locations of the inner raceway groove curvature center before and after the deflection of the $j$th ball element and raceways, respectively. The normal force $Q_j$ on ball element and the deflection of ball and raceways is given in detail in Fig. 1(c), where the bearing structures including ball element, inner raceway and outer raceway are ignored for purpose of simplification.

Based on Fig. 1(b)-(c), the total deformation of the $j$th ball element and raceway can be described in the following expressions:

where ${\psi }_j$ is the $j$th rolling element azimuth, ${\psi }_0$ is the angle between the first rolling element and $x$-axis, $N$ is the total number of ball elements for the ball bearing or roller elements for roller bearing, $A_0$ and $A_j$ are the unloaded and loaded relative distance between the inner and outer raceway groove curvature centers, and ${\alpha }_0$ is the unloaded contact angle for ball bearing or roller bearing. The $j$th ball element and raceways deflection in the axial ${\left(\delta \right)}_{zj}$ and the radial ${\left(\delta \right)}_{rj}$ directions are given as follow according to the displacements of inner ring:

where $r_j$ is the radial distance of the inner raceway groove curvature center for the ball bearing, and $r_L$ is the bearing radial clearance. Note ${\delta }_{bx}$, ${\delta }_{by}$ and ${\delta }_{bz}$ are the translational displacements of the inner ring along $x$, $y$, $z$ axis, respectively; ${\beta }_{bx}$, ${\beta }_{by}$ are angular displacements of the inner ring along $x$, $y$ axial directions.

The load-deformation relation for outer raceway-ball element-inner raceway can be defined by the Hertzian contact theory [25]:

where $Q_j$ is the resultant normal load on the ball element, and $K_n$ is the effective stiffness constant for the inner race-ball element-outer race contacts and is a function of the bearing geometry and material properties. The exponent $n$ is equal to 3/2 for ball type with elliptical contacts. The effects of centrifugal forces and gyroscopic moments of rolling element on ball bearing and roller bearing are ignored as these effects are considered only at extremely high rotational speeds.

The normal force $Q_j$ on the $j$th ball element can be divided into two components $Q_j\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_j$ and $Q_j\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_j$ along the negative direction of $r$ axis and $z$ axis, respectively, as shown in Fig. 1(c). Due to the $j$th rolling element azimuth ${\psi }_j$, the component $Q_j\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_j$ can be divided into $Q_j\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_j\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j$ (along the negative direction of $x$ axis) and $Q_j\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_j\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_j$ (along the negative direction of $y$ axis). The orthogonal load components of normal forces on all ball elements, including $Q_j\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_j\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j$ ($x$ axis), $Q_j\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_j\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_j$ ($y$ axis) and $Q_j\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_j$ ($z$ axis) and bearing forces and moments lead to the following forces and moments equilibrium relations:

where ${\alpha }_j$ is the loaded contact angle, and the two relations, $\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_j=\left(A_0\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_0+{\left(\delta \right)}_{zj}\right)/A_j$, $\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_j=\left(A_0\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_0+{\left(\delta \right)}_{rj}\right)/A_j$ can be obtained based on Fig. 1(c).

The load-displacement relation of ball bearing can be obtained by substituting the Eqs. (1)-(3) into the Eq. (4). This relation is very complicated if the normal force $Q_j$ is expressed explicitly with the displacements of inner ring. Here, the following simplified expression by Lim and Singh [12] is still used by substituting Eq. (1a), Eq. (1c) and Eq. (3) into Eq. (4):

Similarly, the load-displacement relations can be obtained for the roller bearing by following the derivation steps mentioned above. Figs. 2(a)-(c) show the undeformed (unloaded) and deformed (loaded) roller element. The point $o$ locating at the center of the effective roller length is the origins of $z^{\mathrm{\text{'}}}$ and $\xi$ axis, where $z^{\mathrm{\text{'}}}$ is the local rolling element axis coordinate and $\xi$ is the dimensionless local coordinate, $\xi =z^{\mathrm{\text{'}}}/L$, $L$ is the effective length of the roller. Note $z^{\mathrm{\text{'}}}$ varies from –$L/$2 to $L/$2, and $\xi$ from –0.5 to 0.5. The total roller and raceways elastic deformation at the origin $o$ for the $j$th roller is:

where $r_j$, $r_L$ and $r_c$ are the pitch bearing radius, bearing radial clearance and crown drop, respectively. Note the term ${\delta }_{bx}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_0\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j+{\delta }_{by}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_0\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_j+{\delta }_{bz}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_0$ is the deformation generated by the translational displacement of the inner ring, and ${\beta }_{bx}{r}_j\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_0\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_j-{\beta }_{by}{r}_j\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_0\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_j$ is the deformation due to the inner ring angular misalignment (${\beta }_{bx}$,${\beta }_{by}$). Both ${\beta }_{bx}$ and ${\beta }_{by}$ cause the rotation of the $j$th roller, and hence also contributing to the extra roller raceway deformation $∆V$ varying along the roller length as:

where:

where $\omega \left({\psi }_j\right)$ stands for the rotation angle of the $j$th roller due to inner ring angular misalignments (${\beta }_{bx}$,${\mathrm{}\beta }_{by}$). The elastic deformation of the $j$th roller effective length could be characterized by a dimensionless parameter $\xi$ and its azimuth angle ${\psi }_j$:

Fig. 2Roller element bearing kinematics and co-ordinate system

a)

b)

c)

The normal pressure along the roller length is usually different due to different deflections at each contact location, as shown in Fig. 2(c). The sum of normal pressure at dimensionless coordinate $\xi$ on the elemental length $d\xi$ can be calculated based on Hertzian contact theory [25]:

and the resultant normal force applied on the roller can be obtained by integration on $\xi$ as:

where $n$ is equal to 10/9 corresponding to line contact.

Fig. 3The rotation of roller element with respect to outer raceway

a)

b)

The subset (${\xi }_1,{\xi }_2$) can be bounded based on Fig. 3. For $L\omega \left({\psi }_j\right)>0$, it implies the deformation mode will be like Fig. 3(a) based on the Eq. (6c). In this case, the upper limit of integral ${\xi }_2$ is 0.5. The lower limit of integral ${\xi }_1$ should satisfy ${\delta }_{Rj{\xi }_1}=V\left({\psi }_j\right)+{\xi }_{1}L\omega \left({\psi }_j\right)\ge 0$ (The deflection ${\delta }_{Rj{\xi }_1}$ at the dimensionless coordinate ${\xi }_1$ is equal to or larger than zero), therefore, it could be obtained by choosing the larger one of both $-V\left({\psi }_j\right)/\left(L\omega \left({\psi }_j\right)\right)$ and –0.5. For $L\omega \left({\psi }_j\right)<0$, it refers to the Fig. 3(b). The lower limit of integral ${\xi }_1$ should be –0.5. The upper limit of integral ${\xi }_2$ should satisfy ${\delta }_{Rj{\xi }_2}=V\left({\psi }_j\right)+{\xi }_{2}L\omega \left({\psi }_j\right)\ge 0$, therefore, could be obtained by choosing the smaller one of both $-V\left({\psi }_j\right)/\left(L\omega \left({\psi }_j\right)\right)$ and 0.5. For $L\omega \left({\psi }_j\right)=0$, the roller does not rotate with respect to outer raceway and the lower and upper integral limits are –0.5 and 0.5, respectively. The normal force on the roller does not exist for $L\omega \left({\psi }_j\right)=0$ when the total elastic deformation of the roller raceways at the origin $o$ for the $j$th roller is less than zero. Based on the above discussion, the subset (${\xi }_1,{\xi }_2$) can be given as:

In order to obtain the equilibrium relations between normal pressures of all the rollers and bearing forces and moments, the normal pressure along the roller length needs to be transformed into point force as shown in Fig. 2(c). The value of the point force is the resultant normal force $Q_j$ and its location (the load eccentricity) is determined by following equation:

where the load eccentricity $e_j$ is the distance between the point load vector line of action and roller mid-point.The point load can be divided into three orthogonal load components along $x$, $y$, $z$ axis as done for ball element and the equilibrium relations between load components from all the rollers and bearing forces and moments are obtained by applying vector sum to the bearing inner ring:

The load-displacement relations for roller bearing can be easily attained by substituting the Eqs. (6)-(9) into Eq. (10). Similarly, the following reduced load-displacement expression can be obtained by substituting Eq. (7b) into Eq. (10):

Here, the bearing displacements for the ball and roller bearings can be solved from Eqs. (5) and (11) using classical iterative method. For example, the Newton-Raphson and Powell’s hybrid methods [24] have been applied to solve these nonlinear equations, which typically require a careful application and can be cumbersome due to the local minimal searching feature. In the following section, a reliable, fast and efficient energy method is applied to solve bearing displacements, which facilitates the quantification of the relations of bearing displacements and applied forces. The proposed energy method is a global searching method, which is based on the principle of minimum potential energy derived from classical mechanics.

3. Energy method

The energy method, based on the potential energy of the bearing system, is developed to solve the bearing displacements. The proposed energy method links the energy of the system to the displacements of the bearing, which can efficiently yield the correct displacement solution vectors. Further discussion of this proposed strategy is provided below.

Based on the theory of elasticity, the total potential energy of the rolling element bearing system is directly related to the displacements of the bearing system. The principle of minimum potential energy states that the potential energy corresponding to the correct solution satisfying all the differential equations and boundary condition is less than the potential energy corresponding to any other admissible displacement that satisfies geometrical equation and displacement boundary condition [26]. According to this principle, if the potential energy of the bearing system can be derived as a function of the displacements of inner ring relative to the housing support structure, the problem becomes finding the set of displacements that, satisfying the displacement boundary conditions and differential equations, make the potential energy minimum.

The potential energy of total bearing system consists of the elastic strain energy and the potential energy of bearing forces. Generally, the deformation is considered only at the “outer raceway-roller element-inner raceway” contacts when loads are applied to bearing. Therefore, the elastic strain energy is caused by the contact deflection and can be calculated by the following equation:

where $U_{strain}$ is the strain energy. Also, $\sigma ={\left({\sigma }_x{\sigma }_y{\sigma }_z{\tau }_{xy}{\tau }_{yz}{\tau }_{zx}\right)}^T$ are the six stress components, and $\epsilon ={\left({\epsilon }_x{\epsilon }_y{\epsilon }_z{\gamma }_{xy}{\gamma }_{yz}{\gamma }_{zx}\right)}^T$ are the six strain components corresponding to stress components. Note that ${\sigma }_i$ and ${\epsilon }_i$ ($i=x$, $y$, $z$) are the normal stress and strain components, and ${\tau }_j$ and ${\gamma }_j$ ($j=xy$, $yz$, $zx$) are the shear stress and strain components, respectively. Finally, $V$ is entire compressed volume, and $1/2{\sigma }^T\epsilon$ is the strain energy for each unit volume. It is almost impossible to determinate the elastic strain energy in applying Eq. (12) since it is very difficult to determine the stress and strain components accurately in the whole compressed volume. However, strain energy in the bearing can be calculated by the work due to forces causing contact deflection at the “outer raceway-roller element-inner raceway” contacts.

For the ball bearing, given the deformation for the “outer raceway-ball element-inner raceway” contacts as ${\delta }_j$ for the $j$th ball element, the corresponding normal load should be calculated in term of the Eq. (3) as:

where the elemental work due to $F_{jc}$ on the elemental displacement along the direction of ${\delta }_j$ can be describedas follow:

and the total work of $F_{jc}$ can be obtained by integrating the above expression:

where ${\delta }_{Bj}$ is the total deflection of the $j$th ball element and can be calculated by Eq. (1), (2) as a function of the displacements of inner ring of the bearing relative to the housing support. Therefore, the total strain energy for whole ball bearing can be obtained by adding $W_j$ for all ball elements:

The potential energy of bearing forces is subsequently calculated from:

Therefore, the total potential energy is the sum of elastic strain energy and potential energy of bearing forces given by:

Substituting Eqs. (16)-(17) into Eq. (18) yields the total potential energy described as a function of the displacements of inner ring of the bearing relative to the housing support structure.

For the roller bearing, provided the deformation at the dimensionless coordinate $\xi$, for the $j$th roller element is ${\delta }_{j\xi }$, then the contact load $F_{j\xi }$ on the elemental length $d\xi$ could be calculated as follow:

then, the elemental work due to $F_{j\xi }$ on the elemental displacement along the direction of ${\delta }_{j\xi }$ is calculated as:

and the total work for $j$th roller element can be attained by integrating with respect to ${\delta }_{j\xi }$ and $\xi$ as given below:

where ${\delta }_{Rj\xi }$ is the total deformation of the $j$th roller element at dimensionless coordinate $\xi$. Also, ${\delta }_{Rj\xi }$ can be calculated by the Eq. (6) and the limits of integration (${\xi }_1$, ${\xi }_2$) can be obtained from Eq. (8). Substituting Eq. (6d) into above Eq. (21) yields:

therefore, the total strain energy for the whole roller bearing could be attained by adding $W_j$ for all the roller elements:

The potential energy of bearing forces for roller bearing can be obtained by the same formula in Eq. (17) used for ball bearing:

The sum of the elastic strain energy and potential energy of the bearing forces contribute to the total potential energy as a function of bearing displacements can be written briefly as:

Given the above analysis, the energy method can be established employing the expression of potential energy in Eq. (18) and Eq. (26) as the objective functions for ball and roller bearing systems, respectively. Furthermore, the left expression of the given system of nonlinear algebraic Eq. (5) and Eq. (11) are actually the derivatives of the objective functions in Eq. (18) and Eq. (26), respectively. The global optimization method in MATLAB can be used to find the correct displacements that make total potential energy minimum and that also satisfy the differential equation and displacement boundary condition. The differential equation, in fact, does not need to be considered since bearing displacements here are the displacements of rigid mass point (i.e. the center of mass of rolling bearing) not a continuum body. The displacement boundary condition also does not exist since all the bearing displacements are unknown priori. It may be noted that since the energy method is based on the principle of minimum potential energy derived from the theory of elasticity, the approach for solving bearing displacements is assumed to be reliable. In addition, the energy method seeks the correct displacements from the energy minimization principle as opposed to earlier approaches [12] based on mathematical viewpoint since the nonlinear algebraic Eq. (5) and Eq. (11) are not solved directly. The feature of the proposed method will yield a fast and robust determination of the bearing stiffness, which can be employed in the vibration analysis of the geared rotor system modeling the bearing effect as supporting stiffness.

4. Computational study

In this section, numerical simulations were conducted to validate the efficiency of the proposed approach. The combination of loads for ball bearing and roller bearing has been given by Lim and Singh [12], which outlines the valid input combination of bearing loads. The translational and rotational displacements due to three different combinations of loads, that are (i) only constant axial force $F_{bz}$; (ii) constant radial force $F_{bx}$, axial force $F_{bz}$ and moment $M_{by}$; and (iii) constant axial force $F_{bz}$, and moments $M_{bx}$ and $M_{by}$ are calculated for ball bearing. And, the translational and rotational displacements for roller bearing subject to three different combinations of loads, that are (iv) only constant axial force $F_{bz}$; (v) constant radial force $F_{bx}$, axial force $F_{bz}$ and moment $M_{by}$; (vi) constant radial force $F_{bx}$, moments $M_{bx}$ and $M_{by}$ are also calculated in following computational cases. The design data for ball bearing and roller bearing is listed in Table 1 [12].

The energy method is used first to solve for the displacements of inner rings of ball bearing and roller bearing due to different combinations of loads. In the proposed computational process, the algorithm searches all the admissible displacements that make the potential energy of the bearing systems local minimum, and then finds the correct displacements that make the potential energy globally minimum from those admissible displacements found locally. In this manner, it is guaranteed that there is only one solution set for each one of these combinations of loads. The computational flowchart for solving bearing displacements using the energy method is shown in Fig. 4. Our new results are also compared to the calculations from the classical Newton method, Powell’s hybrid method [24] and Modified Powell’s method [27] applied to the systems of nonlinear algebraic Eq. (5) and Eq. (11). The comparison shows that all of the methods yield the same bearing displacements as these obtained from energy method. However, all of these classical approaches require multiple sets of initial guesses. As noted earlier in this paper, the classical iterative approaches necessitate the use of different initial guesses to find the proper displacements since the number of solutions of nonlinear equations is not known a priori. In contrast, these cumbersome procedures are not needed in the newly proposed energy method because all the admissible displacements solutions can be found even if multiple solutions exist for systems of nonlinear algebraic Eq. (5) and Eq. (11). Further discussions and comparisons are given below.

Fig. 4The computational flowchart of the energy method

For the case of only the axial displacement is non-zero when only the axial force is loaded onto the ball bearing, as shown in Fig. 5a, the results are as expected. The axial displacement decreases with greater unloaded contact angle ${\alpha }_0$, which implies the capacity to resist axial force improves as ${\alpha }_0$ increases. The same trend is shown when the combinations of loads, that are cases (ii) and (iii), are applied on the bearing system. For the combination of complex loads, denoted by case (ii), more non-zero displacements exist, as shown in Fig. 5b. The translational displacements in the $X$-direction increase with the unloaded contact angle, which is caused by the curvature of the raceway that provides lesser resistance to the radial loads. The rotational displacement ${\beta }_{by}$ increases in the beginning with increasing unloaded contact angle, and then decreases as the unloaded contact angle increases further. In Fig. 5c, the bearing inner ring generates rotational displacements about both the $X$-direction and $Y$-direction due to the existence of the moments $M_{bx}$ and $M_{by}$. The radial displacement here is obviously less than ${\beta }_{bx}$ for the combination of loads case (ii) because no radial forces are applied in third case.

Fig. 5aThe displacements of ball bearing given a constant axial force Fbz= 3000 N, as denoted by case (i)

Fig. 5bThe displacements of ball bearing given constant radial force Fbx= 1000 N, axial force Fbz= 3000 N and moment Mby= 5000 Nmm, as denoted by case (ii)

Fig. 5cThe displacements of ball bearing given constant axial force Fbz= 5000 N, moments Mbx= 3000 Nmm and Mby= 10000 Nmm, as denoted by case (iii)

For the roller bearing, it should be noted that ${\alpha }_0=$ 0° implies a cylindrical roller bearing type in which only the combinations of radial forces and transverse moments could be applied. This type of cylindrical roller bearing is not designed to carry axial forces without an axial flange. For this reason, the bearing structure is not expected to deform axially, and its potential energy is only a function of the radial and angular displacements. This is evident from the fact that the variation of axial displacement has no effect on the elastic deformation of roller element due to the vanishing term ${\delta }_{bz}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_0$ in the Eq. (6a). Therefore, axial displacement must vanish in the computational process. To deal with this special case, the lower and upper bounds of axial displacement are both set to zero, which guarantees that the axial displacement is kept constant at zero value and does not influence the calculation of the potential energy for each iterative step. This construction is necessary to attain a reasonable combination of radial and angular displacements with minimal potential energy. On the other hand, for thrust roller bearing (${\alpha }_0=$ 90°), the radial displacement should always be set to zero to avoid arbitrary radial translational motion and to ensure only the axial and angular displacements exist when subject to axial forces and moments. The ability to discern these special cases is necessary to facilitate the formulation of the proposed energy method to solve for the bearing displacements for any contact angles for roller bearing. Examples of computational studies of the roller bearing system are given next.

Fig. 6aThe displacements of roller bearing given a constant axial force Fbz= 10000 N, as denoted by case (iv)

Fig. 6bThe displacements of roller bearing given constant radial force Fbx= 3000 N, axial force Fbz= 10000 N, and moment Mby= 10000 Nmm, as denoted by case (v)

The bearing displacement results are shown in Fig. 6a for roller bearing due to constant axial force, as denoted by case (iv). Since the axial force carrying capacity provided by axial flange is not considered when contact angle approximates zero, only the displacements for ${\alpha }_0\ge$ 10° are calculated here. The axial displacement decreases with the increasing contact angle similar to the trend for axially loaded ball bearing. In this case, the axial displacement decreases very rapidly for ${\alpha }_0\le$30° and then slows down its decrease by approaching a horizontal asymptote for ${\alpha }_0\ge$ 30°. Note that all the other displacements except axial displacement are zero as expected. The displacements results in roller bearing due to the combination of loads, denoted by case (v), are obtained as shown in Fig. 6b only for contact angles varying from 10 to 30 degrees because of the practical limitation in the translational bearing motion for high contact angles. For example, the radial displacement along the $X$-direction reaches 2.2916 mm when contact angle is at 40 degrees, which is impractical. Also, as expected, there are more non-zero displacements for this case (v) than the pure axial loading case. Again, it is observed that the axial displacement decreases when the contact angle increases similar to case (iv). The radial displacement in the $X$-direction and the rotational angle along the $Y$-direction increase with the increasing contact angle, which implies that the roller bearings with greater angles provide less capacity to resist radial force, and the incremental radial displacement, contributes to the additional rotational angle.

Fig. 6cThe variation of radial displacements of cylindrical roller bearing (α0= 0°) given constant radial force Fbx= 5000 N, moment Mbx= 5000 Nmm, and Mby= 10000 Nmm, as denoted by case (vi)

Fig. 6dThe variation of rotational displacements of cylindrical roller bearing (α0= 0°) given constant radial force Fbx= 5000 N, moment Mbx= 5000 Nmm, and Mby= 10000 Nmm, as denoted by case (vi)

To study the effect of orbital position of roller elements on bearing displacements, time-varying bearing displacements for cylindrical roller bearing (${\alpha }_0=$ 0°) given constant radial force $F_{bx}$, moments $M_{bx}$ and $M_{by}$, as denoted by case (vi), are shown in Figs. 6c-6d. The radial and rotational displacements (axial displacements are zero) behave nearly like a sine wave with the normalized orbital position angle ($\psi /{\psi }_T$, ${\psi }_T=2\pi /N$ is the element-to-element angular distance).

5. Parametric studies

The stiffness for ball bearing and roller bearing can be derived as shown below applying Eq. (5) and Eq. (11):

The stiffness matrix of rolling element bearing systems includes radial, axial and rotational stiffness coefficients (diagonal), and coupling stiffness coefficients (off-diagonal). Explicit expressions of stiffness coefficients as the function of bearing displacements for ball bearing and roller bearing can be found in Lim and Singh’s [12, 14] papers, as shown in Appendix A1 and A2 (note that ${\left[K\right]}_b$ is symmetric, i.e., $k_{bij}=k_{bji}$).

For ball bearing, the axial and rotational stiffness coefficients given the constant axial bearing force denoted by case (i), as shown in Fig. 7, increase with the increase in unloaded contact angle ${\alpha }_0$, but the radial stiffness shows the inverse trend. The coupling stiffness coefficients remain almost constant when the unloaded contact angle is less than 30 degree, and decrease gradually when ${\alpha }_0$ exceeds 30 degrees. The results also show that all dominant stiffness coefficients mentioned above are very significant for deep groove ball bearing, but axial and rotational stiffness coefficients appear more prominent for angular contact ball bearing.

Fig. 7Dominant stiffness coefficients of ball bearing given a constant axial force Fbz= 3000 N, as denoted by case (i)

In the case when the ball bearing is subjected to the combinations of $F_{bx}$, $F_{bz}$ and $M_{by}$ denoted by case (ii), as shown in Fig. 8, there are more dominant stiffness terms as compared to the case with only axial preload. The radial, axial and rotational coefficients show the same trend as seen in Fig. 7. The coupling terms related to the axial displacement ${\delta }_z$ is insensitive to change in ${\alpha }_0$. When the ball bearing is loaded with the combinations of $F_{bz}$, $M_{bx}$ and $M_{by}$ denoted by case (iii), as shown in Fig. 9, the dominant stiffness coefficients are slightly different from the previous two cases. The coefficient $K_{bz\theta y}$ is also very significant for deep groove ball bearing except for the radial, axial and rotational stiffness terms. From middle to high value of ${\alpha }_0$, the axial and rotational stiffness coefficients are more dominant. The rotational terms, $K_{b\theta x\theta x}$ and $K_{b\theta y\theta y}$, are nearly the same even though the moment $M_{by}$ is much larger than $M_{bx}$, but $K_{bz\theta x}$ and $K_{bz\theta y}$ depict significant differences.

For roller bearing subjected to a constant axial force $F_{bz}$, as denoted by case (iv),the radial and rotational stiffness coefficients shown in Fig. 10 display the similar tendency as those seen for axially loaded ball bearing. The coupled stiffness coefficients are found to be significant in the middle contact angle range. The stiffness coefficients related to radial displacements vanish when contact angle is equal to 90° since this type of thrust bearings are not able to carry radial loads and hence no radial displacements are possible.

More dominant stiffness coefficients are shown in Fig. 11 for roller bearing given the combinations of $F_{bx}$, $F_{bz}$ and $M_{by}$, as denoted by case (v). The radial stiffness coefficients are dominant near the contact angle of 10 degrees but the rational and coupled stiffness coefficients except for those related to axial displacement are dominant for contact angles approaching 30 degrees.

Fig. 8Dominant stiffness coefficients of ball bearing given constant radial force Fbx= 1000 N, axial force Fbz= 3000 N and moment Mby= 5000 Nmm, as denoted by case (ii)

Fig. 9Dominant stiffness coefficients of ball bearing given constant axial force Fbz= 5000 N, moments Mbx= 3000 Nmm and Mby= 10000 Nmm, as denoted by case (iii)

Fig. 10Dominant stiffness coefficients of roller bearing given constant axial force Fbz= 10000 N, as denoted by case (iv)

Fig. 11Dominant stiffness coefficients of roller bearing given constant radial force Fbx= 3000 N, axial force Fbz= 10000 N and moment Mby= 10000 Nmm, as denoted by case (v)

The dominant time-varying stiffness coefficients for cylindrical roller bearing (${\alpha }_0=\mathrm{}$0°) due to combination of loads $F_{bx}$, $M_{bx}$ and $M_{by}$ denoted by case (vi) are shown in Fig. 12. The stiffness coefficients fluctuate with the orbital position of roller elements corresponding to the fluctuation of displacements observed in Figs. 6c-6d. The coupled stiffness coefficients vary more substantially than non-coupled stiffness coefficients. It is expected that incorporation of these time-varying stiffness coefficients can contribute to improvement in results of dynamic analysis of the bearing-supported geared rotor systems as compared to using only the averaged bearing stiffness values [13].

Fig. 12The variation of stiffness coefficients of cylindrical roller bearing (α0= 0°) given constant radial force Fbx= 5000 N, moments Mbx= 5000 Nmm and Mby= 10000 Nmm, as denoted by case (vi)

It is worth mentioning that the proposed energy method is highly efficient in calculating the bearing displacements and stiffness. In real time, all the cases showed that the bearing results only took several seconds to compute at specified orbital position of roller elements and unloaded contact angle using a typical desktop computer. The parameters studies show that most bearing stiffness terms vary significantly with unloaded contact angle except for those having zero values, and the time-varying characteristic of the stiffness terms may also be important.

6. Conclusions

A fast and reliable numerical-based energy method based on principle of minimum potential energy is developed and applied to compute for the displacements of ball and roller bearings under complex loads. In this energy method, a global optimal problem is analyzed in which the potential energy of the bearing system constitute the objective function. The global optimization algorithm is used to solve for the bearing displacements. The bearing displacements can be obtained by searching and finding the minimum total potential energy of the bearing system instead of solving the system of nonlinear algebraic equations directly. In addition, the proposed energy method derived from the well-established theory of elasticity runs quicker by avoiding the extra trial and error efforts in testing a range of initial estimates as needed in the classical iterative method. The computed bearing stiffness will be very important for the vibration analysis of the geared rotor system in order to design a quiet driveline. Also, the proposed energy method can be employed to analyze the sensitivity of bearing parameters on the potential energy, which can be applied to study the vibration problem of certain designed bearings.

The effect of unloaded contact angle on bearing stiffness of ball and roller bearings is analyzed. The analysis reveals the trends of the dominant stiffness coefficients. Time-varying characteristic of the stiffness coefficients are also studied using the proposed energy method. The analysis shows that the coupled stiffness coefficients are more sensitive to orbital position of roller elements as compared to non-coupled stiffness coefficients.