Abstract
With the continuous advancement in bearing speed, bearing assembly clearance (pocket, radial, and guide clearances) has a particularly significant impact on dynamic characteristics of bearing. In order to improve its performance, it is necessary to generate an optimized design for assembly clearance based on dynamic analysis. In this paper, a dynamic model of cylindrical roller bearing has been put forth based on the variable step fourth-order Runge–Kutta method, while the simulation results have been verified by a high-speed bearing cage motion testing machine. Considering pocket, radial, and guide clearances as independent variables and maximum impact force, whirl deviation ratio, and minimum power loss as the objectives, a multiobjective optimized model of cylindrical roller bearing has been developed using central composite experimental design (CCD) and response surface method. After optimization, the maximum impact force was reduced by 9.26%, the whirl deviation ratio was reduced by 2.87%, and power loss was reduced by 1.45%.
1. Introduction
The cylindrical roller bearing is widely used in aeroengines, high-speed railways, robotics, and other fields because of its small friction coefficient and excellent high-speed performance. According to failure case analysis and statistics, cage fatigue fracture is one of the serious failure forms of high-speed cylindrical roller bearings. The whirl of the cage and frequent collision and friction between rollers and cage are the main causes for cage failure [1, 2].
Bearing assembly clearance plays a vital role in the stability of bearings. Ghaisas et al. [3] studied the effects of cage pocket clearance, guide clearances, roller diameter, and inner ring inclination angle on cage whirl trajectory and cage speed deviation ratio and further used cage speed deviation to evaluate the cage stability. Tu et al. [4] proposed a nonlinear model to investigate the dynamic interactions between the rolling element and cage under rotational speed fluctuation conditions. In this regard, the effects of cage pocket clearance and amplitude and frequency of fluctuation on interaction forces between the rolling element and cage were analyzed. Yang et al. [5] presented a five degree-of-freedom (5-DOF) quasidynamic model to analyze the relationship between cage clearance and heating characteristics. The results showed the existence of a critical value for both the guide and pocket-hole clearances. Zhang et al. [6] developed differential equations for high-speed cylindrical roller bearings based on dynamic theory. It was found that larger or smaller guide land clearance or pocket circumference clearance raised the vibration of the cage. Therefore, assigning a clearance value within an optimal range for guide land and pocket circumference clearances will reduce the cage’s vibration effectively. Choe et al. [7] analyzed the effects of clearance and rotation speed on the cage stability and performance using the probability density function of cage whirling frequency. The cage wear loss increased with decreasing cage ball-pocket clearance because of the collisions between the balls and cage pocket. Kim [8] developed a simplified roller bearing dynamic model that reflected radial clearance. The study results showed that clearance of bearing components (pocket, radial, and guide clearances) had a significant impact on the dynamic characteristics of the cage and should be considered in designing cylindrical roller bearings to improve its dynamic performance. With the increase of bearing speed, bearing clearance will also cause nonlinear resonances, bifurcations, stabilities, and additional dynamic loads of the bearings in the rotor system. Yang et al. [9] established an analytical model, and the vibration responses of the system are calculated in a large speed range with the consideration of different ball numbers and different rotor eccentricities. Gao et al. [10, 11] developed a nonlinear dynamic model of intershaft bearing in a dual-rotor system, and then the dynamic load of the intershaft bearing is obtained by solving the dynamic equations of the system numerically. An exhaustive parametric analysis for temperatures and nonlinear thermal behaviors of the intershaft bearing affected by dynamic parameters (including the rotation speed ratio, unbalances of rotors, the radial clearance, the stiffness, and the roller number of the intershaft bearing) and thermal parameters (including the lubricant viscosity and the ambient temperature) was carried out.
The performance of bearings depends largely on the design. Choi and Yoon et al. [12] have optimized the design variables for double row angular contact ball bearings of an automobile wheel-bearing unit using a genetic algorithm (GA). Chakraborty et al. [13] applied optimized deep groove ball bearing (DGBB) using GA; however, certain conditions like assembly angle were kept constant for all pairs of solutions. Rao and Tiwari et al. [14] attempted a parametric study on the design variables specified by Chakraborty et al. [13] and applied bounds to the five constant variables involved in the constraints. Gupta et al. [15] applied nondominated sorting genetic algorithm-II (NSGA-II) to a mathematical model which was a multiobjective function, and a parametric study was carried out. Tiwari and Chandran et al. [16] took three primary objectives, namely, the dynamic capacity, the elastohydrodynamic lubrication (EHL), minimum film-thickness, and the maximum temperature which have been optimized with an evolutionary algorithm. Kalyan et al. [17] proposed a design methodology to optimize the dynamic capacity and the elastohydrodynamic minimum film thickness, and in the optimization problem, a total of seven design variables (which comprises bearing pitch diameter, roller mean diameter, effective roller length, number of rollers, and three other unknown constraint factors), two objectives, and 16 constraints are considered. Hingawe et al. [18] adopted response surface methodology based Box–Behnken design to study the effect of texture geometrical parameters, namely, position, depth, orientation, and slender ratio on tribological characteristics such as load-carrying capacity and friction coefficient; furthermore, a multiobjective optimization study is carried out using grey relational analysis. Ahmad et al. [19] took five internal geometrical parameters which are taken as design variables; moreover, five design constraint factors are also included. Thirty-six constraint equations are considered, which are mainly based on geometry and strength considerations. The objective functions were optimized individually (i.e., the single-objective optimization) and then simultaneously (i.e., the multiobjective optimization), and NSGA-II has been used as the optimization tool. Pareto-optimal fronts are obtained for one of the bearings. Of many points on the Pareto-optimal front, only the knee solutions have been presented.
In most literature studies, the optimization design of bearings mainly optimizes the structural parameters such as the diameter of the roller and the radius coefficient of ditch rate of inner and outer rings, and the optimization result is a certain value of bearing structural parameters. But in the actual production process, no matter how high precision the equipment is, the bearing elements always had errors, which will lead to a few individual bearings having excellent performance. According to the previous analysis, bearing clearance has a great impact on bearing performance. If the optimal bearing clearance range can be obtained through optimization, then, according to the principle of assembly dimension chain, the bearing clearance is used as a closed ring to determine the basic size and deviation of bearing elements. Then, the qualified product rate of the bearing will be improved greatly.
In view of these, the present paper expounds on a model of high-speed cylindrical roller bearing based on the dynamic theory. Consequently, considering maximum impact force, whirl deviation ratio, and minimum power loss as the objectives, CCD and RSM have been used to optimize the clearance of cylindrical roller bearing components. The research results can be used in the structural design of high-speed cylindrical roller bearing.
2. Dynamic Modeling and Experimental Verification
A dynamic model of cylindrical roller bearing based on the Gupta method [20] describes the three-dimensional and instantaneous motion of bearing elements. However, the interaction between elements is complex, and a large number of coordinate transformations and numerical operations of differential equations are required to delineate the relationship. In order to facilitate effective analysis and simplify the calculation, the study makes the following reasonable assumptions:(i)Ignoring the axial motion of all parts(ii)Ignoring the influence of temperature(iii)Centroid of each element of bearing coincides with the geometric center
2.1. Relationship between Bearing Elements
2.1.1. Interaction between Roller and Rings
Considering the interaction between the roller and outer ring as an example, the interaction between the roller and the inner ring is similarly adopted. Firstly, a series of coordinate systems such as inertial coordinate system OiXiYiZi, roller fixed body coordinate system ObXbYbZb, and ferrule fixed body coordinate system OrXrYrZr was created to describe the relative position and motion of bearing elements, as shown in Figure 1.
The inertial coordinate system was fixed in space whose origin was located at the curvature center of the outer ring, and the Xi axis was along the bearing axis. The roller fixed body coordinate system was fixed on the roller and moved along with it, whose origin was located at the center of mass of the roller, and Xb was along the axis of the roller. The ring fixed body coordinate system was fixed on the ring and moved along with it, whose origin was located at the ring centroid, and Xr was along the ring axis. In addition, the roller orientation coordinate system OaXaYaZa needed to be constructed to describe the track position of the roller centroid, whose origin was located in the roller center.
Under the initial conditions, Xa was parallel to Xi, Za passed through the roller center and intersected vertically with the axis where Xi was located, and Ya was determined by the right-hand rule. Other local coordinate systems were constructed according to modeling needs.
As shown in Figure 1, the position vectors and of the roller center and ring center relative to the inertial coordinate system were obtained according to the geometric relationship (superscript i represents the inertial coordinate system). The position vector from the center of the roller to the center of the ring (superscript r represents the ring fixed body coordinate system) in the ring fixed body coordinate system is expressed as follows:where Tir was the transformation matrix from inertial coordinate system to ring fixed body coordinate system.
During the operation of the roller, if the roller is skewed, that is, it rotates around the Yb axis, the contact length between the roller and the ring changes along the axis direction. Therefore, it is necessary to cut the roller with length l into m round flakes on average. The coordinate of the kth round flakes center relative to the roller center was ((−0.5+(k−0.5)/m) l,0,0). The position vector from the center of round flakes to the center of the ring is expressed as follows:where Tbr was the transformation matrix from roller body-fixed coordinate system to ring body-fixed coordinate system.
In order to describe the track position of round flakes center in bearing, it was necessary to establish the azimuth coordinate system OasXasYasZas of round flakes whose origin was located at its center, and other directions were consistent with the direction of the roller orientation coordinate system. According to the components of in Yr and Zr directions, the azimuth angle t of the position vector from the center of the round flakes to the ring center was determined using the following formula:
The transformation matrix Tras from ring coordinate system to round flakes azimuth coordinate system was determined by azimuth angle. Considering the point P to be a point on the round flakes close to the outer raceway, the position vector of point P relative to the ring center in the round flakes in the azimuth coordinate system was calculated as follows:where D was the diameter of the roller and ϕ was the angle at which point P turned around the roller axis.
To determine whether the round flakes were in contact with the outer raceway, the approach between point P and the outer raceway was calculated as follows:
If was positive, there exists a negative gap, and deformation occurs between roller and raceway. The contact load Fbnk between the kth round flakes and the raceway was calculated according to the Hertz line contact theory, as shown belowwhere Eeq was the equivalent elastic modulus, leq was the effective contact length of the roller, and dx was the thickness of the round flakes.
Given the contact between the roller round flakes and raceway at point P, it was necessary to obtain the relative sliding speed and relative position vector between the center of the raceway contact area and roller contact center to calculate the tangential traction force on the roller. The contact coordinate system OcXcYcZc was established at P, and the direction of each coordinate axis was consistent with that of the orientation of the round flakes coordinate system with the coordinate origin at point P. Assuming that contact deformation of roller and raceway was equal, the position vector from contact point P to roller center in roller fixed body coordinate system was obtained as
The translation speed of the ring in the inertial coordinate system was , rotation speed was , translation speed of the roller in the inertial coordinate system was , and rotation speed was . In the contact coordinate system, the velocity at the contact point P of the round flakes is calculated using the following:where Tib was the transformation matrix from inertial coordinate system to roller fixed body coordinate system, Tbc was the transformation matrix from roller fixed body coordinate system to contact coordinate system.
The linear velocity of contact point P on the ring in the contact coordinate system was calculated aswhere Tac was the transformation matrix from slice azimuth coordinate system to contact coordinate system, Tias was the transformation matrix from inertial coordinate system to slice azimuth coordinate system, Tbas was the transformation matrix from roller fixed body coordinate system to round flakes azimuth coordinate system. Consequently, the sliding speed of point P on ring relative to point P on round flakes was written as follows:
The relative sliding speed u was calculated by the following formula:
The relative sliding speed was substituted into (12) (Gupta lubrication traction model) to calculate the lubrication traction coefficient μ as follows:where A, B, C, and D are lubricant coefficients.
In the contact coordinate system, the force vector acting on the round flakes was written as follows:
In the roller fixed body coordinate system, the torque vector on the round flakes was written as follows:
Therefore, the resultant force and resultant torque vector of the raceway acting on the roller was written as
According to the principle of force and reaction force, the resultant force and resultant torque vector of the roller acting on the bearing raceway were obtained.
2.1.2. Interaction between Roller and Cage
There is a gap between the roller and cage pocket; however, they collide during movement at high speed. In order to analyze the interaction between roller and cage, cage fixed body coordinate system OcaXcaYcaZca and pocket coordinate system OcpXcpYcpZcp were formed based on the interaction formulae developed in the previous section. As shown in Figure 2, the origin Oca of the fixed body coordinate system of the cage coincided with the center of the cage.
In the initial state, the directions of the three coordinate axes were consistent with the inertial coordinate system wherein the pocket coordinate system was fixed on each pocket (origin Ocp at the center of the pocket). The Xcp axis was perpendicular to the pocket wall, while the Zcp axis was located along the radial direction of the cage, pointing outward.
Vectors and , as observed in Figure 2, were the position vectors of the roller center and cage center, respectively, in the inertial coordinate system. From this vector relationship, the position vector of the roller center relative to the cage center was obtained. The represented the position vector from the cage pocket center to the cage center in the cage fixed body coordinate system.
Considering the inclination of the roller (rotation around Zb axis) and partial contact between roller and front wall of the cage pocket, the slicing method has been used with the contact point as point . The solution flow was similar to the one used in the previous section. After a series of the vector sum and difference operations, the position vector (superscript cp represents the cage pocket coordinate system) of the contact point on round flakes relative to the cage pocket in the pocket coordinate system has been obtained.
Without considering the axial movement of the roller, it collides only with the front wall of the cage. Therefore, the approach between the roller and cage was calculated by the following formula:where lc was the pocket width of the cage.
When , that is, the deformation has occurred between roller and cage, the contact load can be obtained by using Hertz contact theory. Subsequently, interaction force and torque between roller and cage could be obtained according to the treatment method similar to the interaction between roller and raceway.
2.1.3. Interaction between Cage and Guide Ring
In this section, as outer and inner ring guidance modes are similar, the outer ring guidance mode has been discussed as an example. The interaction between the cage and the guide ring is shown in Figure 3.
As shown in Figure 3, in the ring fixed body coordinate system, was the position vector of cage center relative to guide ferrule center, while was the position vector of any point P on the edge of the cage relative to the cage center (superscript r represents the ring fixed body coordinate system). Then, in the ring fixed body coordinate system, was the position vector of any point P on the edge of the cage relative to the guide ring center and it is expressed as
Setting the guide radius of the ferrule as Rcg and the approach as between the edge of the cage, the guide ring was calculated as follows:
When , deformation occurs between the cage and the guide ring while the contact load was obtained using the Hertz contact theory. The interaction force and torque between the cage and the guide ring were obtained using the treatment method similar to the interaction between the roller and the raceway.
2.1.4. Solution and Verification of Dynamic Equations
The motion of bearing elements is generally classified into translation and rotation around the center of mass. For convenience, dynamic equations of different elements have been described using different coordinate systems with respect to their motion characteristics. Without considering the axial motion of the roller, the translational motion of the roller was described as follows:where mb is the mass of the roller, rb is the roller radius, Fr and Fθ are the radial and axial components of the resultant force on the rolling element in the inertial cylindrical coordinate system.
Ignoring the axial motion, the translational motion of the inner ring, outer ring, and cage were described aswhere Fy and Fz were the components of the resultant force on the inner ring, outer ring or cage in the Y and Z directions, respectively, m was the mass of the inner ring, outer ring or cage.
The rotation of any bearing elements was described by the Euler equation aswhere Ix, Iy, and Iz were the rotational inertia of the roller in three directions of the coordinate axis, M1, M2 and M3 were the components of the resultant torque on each coordinate axis of the roller fixed body coordinate system.
Dynamic equations of roller bearing system were numerically assessed using the step-changing fourth-order Runge–Kutta method. A flowchart of the dynamic simulation program used in the current study is shown in Figure 4.
Compared with the simplified model, the complex dynamic model based on the Gupta modeling method considered the issues of roller sliding and cage collision and simulated transient motion of the cage. The model verification was performed by measuring the speed of the cage. To validate the model, the following test was carried out:(i)A radial load of 1200 N was applied(ii)The rotating speed of the shaft was increased from 20000 r/min to 32000 r/min(iii)Group of tests for every 2000 r/min was conducted(iv)Actual speed of cage under 7 groups of different speeds was obtained(v)The results were compared with simulation values
The test equipment consisted of a high-speed bearing cage motion state test machine and a high-speed camera. The high-speed bearing cage motion state testing machine was composed of test main parts, hydraulic loading system, high-temperature lubrication system, and cooling system, as shown in Figure 5.
The high-speed camera was composed of an ACS camera, F-Port adapter, tripod quick mount board, standard communication cable, power adapter, control software MLink, etc., as shown in Figure 6.
The high-speed camera was used to capture the position of mark point 1 on the cage in real-time, and MoviasNeo postprocessing software was used to obtain the displacement curve of mark point 1 relative to the center point 0 of the main shaft. As shown in Figure 7, the time difference corresponding to the adjacent wave crest position of the displacement curve was the time required for the cage to rotate for one cycle so as to obtain the actual speed of the cage.
The structural parameters and working condition parameters of the object bearing in this paper are shown in Tables 1 and 2.
Further, the test results were compared with model simulation results to verify the model. As shown in Figure 8, the simulation results were generally close to the test results, and the trend was mostly the same. However, there was a large deviation in some places, which was due to the difference between dynamic model parameters and actual parameters. Although there may be errors in the actual measurement process of cage speed, the comparison of values proved that the bearing dynamic model showed high overall accuracy.
3. Design and Optimization
From the review of research discussed in Section 1, the present paper considered the following three factors that had the most impact on the dynamic performance of bearing as the optimization independent variables: radial clearance (A), guide clearance (B), and pocket clearance (C). The optimization objectives considered were maximum impact force (Y1), whirl deviation ratio (Y2), and power loss (Y3). RSM was the design approach used for modeling and analysis of the problem in the present study. Using this approach, the correlation between variable input parameters and output responses was developed. This approach was found to be effective and efficient [18].
In this research, modeling and analysis of clearances of cylindrical roller bearing components were carried out using the CCD design of RSM. As shown in Figure 9, the CCD center point (0,0,0) was used to determine the test error and reproducibility, while (−1,1) represented the high and low levels of variables. The range of axial points was (±α,0,0), (0,±α,0), (0,0,±α), and the value of axial points away from the center was related to the factor value k. From the above analysis, the bearing boundary design parameters of k = 3, α = 1.682 have been considered in this paper. Table 3 shows the codes and levels of test design factors.
For the three parameters and five levels, a parametric set of 20 runs (typical working conditions: rotation speed 24000 r/min, radial load 1300 N) were obtained using the dynamic model in Section 2, as shown in Table 4.
According to the data samples in Table 4, the regression models of maximum impact force, whirl deviation ratio, and power loss were obtained.
It can be observed from Tables 5–7 that the three regression models were highly significant (P-value＜0.01), whereas the mismatch of models was not significant ( value >0.05). Therefore, the fitting degree of these regression models was high. The three determination coefficient R2 and correction coefficient of the three models were both close to 1, indicating the maximum impact force, whirl deviation ratio, and power loss fitting regression model have high reliability.
Based on the analysis results of the regression model, the 3D response surface diagram of the interaction effects of various factors has been drawn. The response surface of radial and guide clearance to the maximum impact force is shown in Figure 10(a). During the increase of radial clearance from 0.015 mm to 0.025 mm, the maximum impact force increased with guide clearance. The response surface of radial and pocket clearance to the maximum impact force is shown in Figure 10(b). When radial clearance increased from 0.015 mm to 0.025 mm, the maximum impact force increased with pocket clearance.
(a)
(b)
As indicated in Figure 11(a), when radial clearance increased from 0.015 mm to 0.025 mm, the whirl deviation ratio decreased with the increase of guide clearance. As depicted in Figure 11(b), when the radial clearance increased from 0.015 mm to 0.025 mm, the whirl deviation ratio increased with pocket clearance.
(a)
(b)
As illustrated in Figure 12(a), when the radial clearance increased from 0.015 mm to 0.025 mm, the power loss decreased with the increase of guide clearance. As observed in Figure 12(b), when the radial clearance increased from 0.015 mm to 0.025 mm, the power loss decreased with the increase of pocket clearance.
(a)
(b)
Combining the regression models of maximum impact force, whirl deviation ratio, and power loss, a multiobjective optimization model was developed, as follows:
NSGA-II algorithm was used for the solution. As the algorithm was relatively mature, the solution process was no longer described. Performance comparison between combination of optimized parameters and original design parameters are shown in Table 8.
In practical work, the bearing clearance interval is often given according to the bearing working conditions, that is, a certain tolerance was given on the basis of the optimal clearance. For example, the value range of pocket clearance in Table 8 can be set as 0.280–0.290 mm. Then, as shown in Figure 13, the assembly dimension chain between roller and cage was established, in which pocket diameter lc and roller diameter D were the unit ring, pocket clearance δbc was the closed ring, which was solved by the extreme value method, pocket diameter lc was and roller diameter D was . The cage and roller manufactured according to the above two dimensions can ensure that the pocket clearance was in the allowable range. Compared with directly optimizing the bearing structure size and giving tolerance in other literature, the yield of this method was higher, which has been verified in engineering applications.
4. Discussion
The results show that after optimization, the maximum impact force reduced by 9.26%, the whirl deviation ratio decreased by 2.87%, and the power loss lowered by 1.45%. However, the above optimization process was conducted for typical service conditions. The validity of the optimization results in the case of change in parameters of working conditions is worthy of discussion. Therefore, maximum impact force reduction, whirl deviation ratio, and power loss corresponding to the original design value and optimal parameter combination under different working conditions have been compared.
The following inferences were derived from Figures 14(a) and 14(f):(i)With constant radial load (Fr = 1300 N), the maximum impact force and power loss corresponding to both design value and optimization value increased with the increase in speed, while the whirl deviation ratio first decreased and then increased(ii)With constant rotating speed and increasing load (N = 24000 r/min), maximum impact force decreased, power loss increased, and whirl deviation ratio first decreased and then increased
(a)
(b)
(c)
(d)
(e)
(f)
However, under all working condition parameter combinations, maximum impact force, whirl deviation ratio, and power loss corresponding to the optimal parameter combinations were lower than the design values, which showed that optimized bearing dynamic performance was better than the original design.
5. Conclusion
A dynamic model of cylindrical roller bearing was developed and solved by variable step fourth-order Runge–Kutta method. In order to verify the effectiveness of the dynamic model, the actual speed of the cage was obtained by using a high-speed bearing cage motion state testing machine and a high-speed camera. Then the central composite experimental design (CCD) method was used to design the simulation test while regression models of maximum impact force, whirl deviation ratio, and power loss were obtained by the response surface method. A multiobjective optimized model of cylindrical roller bearing has been developed based on the regression models, and the NSGA-II method was used to solve the optimized model. The conclusions of this research were as follows:(i)The actual speed of the cage was obtained by using a high-speed bearing cage motion state testing machine and high-speed camera to verify the effectiveness of the dynamic model.(ii)The optimal parameter combination of clearance (pocket, radial, and guide clearances) of cylindrical roller bearing components was obtained. After optimization, the maximum impact force decreased by 9.26%, the whirl deviation ratio reduced by 2.87%, and the power loss lowered by 1.45%.(iii)The maximum impact force reduction, whirl deviation ratio, and power loss corresponding to the original design value and optimal parameter combination under different working conditions have been compared. The optimized bearing is better than the original design value under various working conditions.
Data Availability
Relevant data have been presented in the paper. If necessary, the data are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The grant was obtained from the National Science Funds of China (no. 51805151), Key Scientific Research Projects of Colleges and Universities in Henan Province (no. 19A460019), and National Key Research and Development Project (no. 2018YFB2000203).