# Study on Simplified Model and Numerical Solution of High-Speed Angular Contact Ball Bearing

Authors: Guang Zeng; Chunjiang Zhao; Xiaokai Yu; Biao Sun; Zhigang Xiao; Qiang Bian; Bingyang Liu
Journal: Shock and Vibration
 (2020)

---

## Abstract

For the calculation model of high-speed angular contact bearing has many variables, the large root difference exists, and the Newton iterative method solving the convergence depends on the initial value problems; thus, the simplified calculation model is proposed and the algorithm is improved. Firstly, based on the nonlinear equations of variables recurrence method of the high-speed angular contact ball bearing calculation model, it is proved that the ultimate fundamental variables of calculation model are the actual inner and outer contact angles, the axial and radial deformations. According to this reason, the nonlinear equations are deformed and deduced, and the number of equations is reduced from 4Z + 2 to 2Z + 2 (Z represents the number of rolling bodies); a simplified calculation model is formed. Secondly, according to the small dependence of the artificial bee colony algorithm on the initial value, an improved artificial bee colony algorithm is proposed for the large root difference characteristics of high-speed ball bearings. The validity of the improved algorithm is verified by standard test function. The algorithm is used to solve the high-speed angular contact ball bearing calculation model. Finally, the deformations of high-speed angular contact ball bearings are compared and verified by experiments, and the results of improved algorithm show good agreement with the experiments results.


---

## Body

## 1. Introduction

Although the angular contact ball bearing has a simple structure, the internal force and motion relationship become complicated due to the gyroscopic moment and centrifugal force at high speed [1]. In order to study the dynamic characteristics of the ball bearing, Jones [2] first established a quasistatic model based on the raceway control theory and Hertz contact theory. Subsequently, Harris et al. [3] further improved the model, which is referred to as the Jones-Harris (J-H) model. Because the J-H model is effective and accurate for the calculation of internal load distribution, contact angle variation, fatigue life, and stiffness [4–11] of high-speed ball bearings, the model is still widely applied.The J-H model consists of a series of deformed geometric equations, rolling element equilibrium equations, and inner ring equilibrium equations, which form a high-dimensional nonlinear equation [12]. The system consists of a total of 4Z + 2 nonlinear equations [1]. In order to make J-H model solution more convenient, many scholars have conducted in-depth research. Antoine [13] proposed a new method for calculating contact angle of high-speed ball bearings under axial load. Zhang et al. [14] solved the dynamic characteristics of high-speed angular contact bearing based on various contact angles and hybrid theory. Kurvinen et al. [15] studied the method of solving the dynamic characteristics of angular contact bearings without considering gyroscopic moment and centrifugal force. Yu et al. [16] studied the particle swarm optimization (PSO) algorithm for solving the nonlinear equations of angular contact bearing. Wang et al. [17] studied the modeling and dynamic characteristics of angular contact bearings under axial, radial, and moment loads simultaneously.Although these methods effectively improve the convergence of the original J-H method, they do not fundamentally reduce the computational complexity of the model. Fang et al. [18] simplified the calculation model by deriving the internal balance equation of the bearing, but the solution still uses the Newton-Raphson method. Although these algorithms have proved to be effective in solving nonlinear systems of equations [1, 19], there are still some shortcomings such as high dependence on initial values, slow convergence speed, and being easy to fall into premature convergence and early convergence when solving nonlinear equations with large root differences [20]. Nowadays, advanced algorithms such as Shuffled Frog Leading algorithm [21], firefly algorithm [22], and artificial bee colony algorithm [23] have emerged with the rapid development of computer algorithm technology. These algorithms can transform the problem of solving nonlinear equations into optimization problems. Therefore, if they can reduce the equations number, meanwhile adopting the improved swarm intelligence algorithm for high-speed angular contact bearing model, they can effectively reduce the computational complexity and convergence and realize the fast and accurate solution of the high-speed ball bearing model.In order to solve the problem of too many nonlinear equations of high-speed angular contact bearing, the author determined the recursive relation among all variables in the model and deduced the calculation model according to the actual contact angle, axial deformation, and radical deformation. The number of solving equations is reduced from 4Z + 2 to 2Z + 2, which reduces the computational complexity by nearly half compared with the traditional algorithm. In addition, for the characteristics of large root difference of the angular contact ball bearing calculation, the artificial bee colony algorithm which initializes the food source in stages is adopted to improve the convergence precision and the solution speed of the solution. Finally, the deformation of angular contact ball bearings at high speed is compared and verified by experiments, and the improved algorithm results show good agreement with experiments results.


## 2. Establishing a Simplified Calculation Model

### 2.1. Traditional Computing Model

Solving the high-speed characteristics of angular contact ball bearings has to start from the contact deformation. Establishing the bearing deformation geometry equation, the rolling element balance equation, and the inner ring balance equation and then establishing the auxiliary equations including centrifugal force, gyroscopic moment, and contact angle are needed. Figure1 is a schematic diagram of the rolling elements in the different ball angular position ψj; Figure 2 is a schematic diagram of rolling body force; Figure 3 shows the relative position relationships between the centers of the inner and outer raceways groove curvature and the rolling elements centers before and after the bearing is loaded (according to the inner ring rotation and the outer ring fixation). According to the literature [12], combined with Hertz contact theory, bearing motion, and force relationship and raceway control theory, the calculation model of high-speed angular contact ball bearing is given.Figure 1
Position map of rolling elements at different azimuth anglesψj.Figure 2
Schematic diagram of the rolling element.Figure 3
Schematic diagram of the change of the center position of the groove before and after the load.


#### 2.1.1. Deformation Geometry Equation

Figure3 shows that the axial distance A1j and radial distance A2j between the curvature centers of inner and outer raceway grooves in the ball angular position ψj are as follows:(1)A1j=BDsinα0+δα+ℜiθcosψj,(2)A2j=BDcosα0+δrcosψj,where δa, δr, and θ, respectively, denote the relative axial deformation, relative radial deformation, and relative angular displacement of the inner and outer rings; ℜi denotes the radius of the groove curvature center trajectory of inner raceway; α0 is the initial contact angle of the bearing.The following geometric relationship can be determined from Figure3:(3)cosαoj=X2jfo−0.5D+δoj,(4)sinαoj=X1jfo−0.5D+δoj,(5)cosαij=A2j−X2jfi−0.5D+δij,(6)sinαij=A1j−X1jfi−0.5D+δij.In these equations,α is the contact angle of the bearing; f is the groove curvature of the bearing; the subscripts o and i are the outer ring and inner ring of the bearing, respectively; D is the ball diameter. According to the deformation geometric relationship between the components in Figure 3, combined with the Pythagorean theorem and (1) and (2), the position change relationship equation among the center of curvature of the rolling elements and the inner and outer grooves in the ball angular position ψj can be obtained:(7)A1j−X1j2+A2j−X2j2−fi−0.5D+δij2=0,(8)X1j2+X2j2−fo−0.5D+δoj2=0.


#### 2.1.2. Rolling Element Balance Equation

According to Hertz contact theory, the normal contact forcesQ of the rolling body and the inner and outer rings at the contact points are as follows:(9)Qij=Kijδij1.5,(10)Qoj=Kojδoj1.5.As shown in Figure2, in the ball angular position ψj, the forced balance equations of rolling elements in the radial and axial can be written as(11)Qijsinαij−Qojsinαoj−MgjDλijcosαij−λojcosαoj=0,(12)Qijcosαij−Qojcosαoj+MgjDλijsinαij−λojsinαoj+Foj=0,where Mg is the gyro moment of the ball. Substituting (3)–(6), (9), and (10) into (11) and (12) yields(13)Kijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij+λojMgjX2j/D−Kojδoj1.5X1jfo−0.5D+δoj=0,(14)Kijδij1.5A2j−X2j−λijMgj/DA1j−X1jfi−0.5D+δij−λojMgjX1j/D+Kojδoj1.5X2jfo−0.5D+δoj+Fcj=0.


#### 2.1.3. Inner Ring Balance Equation

According to the force balance relationship and the moment load balance relationship of the bearing inner ring in the axial and radial directions, the inner ring balance equation can be obtained:(15)Fa−∑j=1j=zQijsinαij−MgjDλijcosαij=0,Fr−∑j=1j=zQijcosαij+MgjDλijsinαijcosψj=0,M−∑j=1j=zQijsinαij−MgjDλijcosαijℜi+fiλijMgjcosψj=0.Substituting (3)–(6), (9), and (10) into (16)–(18) results in(16)Fa−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij=0,(17)Fr−∑j=1j=zKijδij1.5A2j−X2j+λijMgj/DA1j−X1jfi−0.5D+δijcosψj=0,(18)M−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δijℜi+fiλijMgjcosψj=0,where Fa, Fr, and M, respectively, are the axial, radial, and moment forces acting on the bearing. Given the initial iteration values of δa, δr, and θ, the iteration algorithm can be used to solve the nonlinear equations composed of formulas (7), (8), (13), and (14) in the ball angular position ψj and obtain X1j, X2j, δij, and δoj. By substituting these into the nonlinear equations (16)–(18) in every ball angular position ψj, the values of the main unknowns δa, δr, and θ can be obtained. Then, the iteration of X1j, X2j, δij, and δoj continued until the values of the main unknowns δa, δr, and θ meet the accuracy.


### 2.2. High-Speed Angular Contact Ball Bearings Calculation Model Analysis

Equations (7), (8), (13), (14), and (16)–(18) constitute a system of nonlinear equations, where there are 4Z + 2 equations. The calculation of contact stiffnesses Kij and Koj, centrifugal force Fcj, and gyroscopic moment Mgj of rolling elements and inner and outer rings at contact points all depends on the calculated values. And Kij, Koj, Fcj, and Mgj are functions of the actual internal and external contact angles αij and αoj; the actual internal angle αijand external contact angle αoj are functions of X1j and X2j. Therefore, the variables in the dynamic characteristic equations are highly coupled, which makes the iterative calculation process very complicated and the calculation amount is extremely large. The initial value is not chosen properly, and it is difficult to converge or even not converge.


#### 2.2.1. Variable Recursion Relation Analysis of Computational Model

(1) Variable Recursion Relation Analysis of Centrifugal ForceFcj. In ball angular position ψj, the centrifugal forces Fcj of the high-speed angular contact ball bearing rolling element can be expressed as follows: (19)Fcj=2.26×10−11D3dmnmj2.In the above formula, the formula for calculating the revolution speednmjof the rolling element is as follows:(20)nmj=nij1−γij1+γojΩi±Ωo1+γojcosαoj−βoj+1−γojcosαoj−βoj.According to the outer raceway control theory, the attitude angleβojof the rolling element is as follows:(21)βoj=arctancosαojsinαojcos2αoj+γoj.In the above formula, the formula for calculating the revolution speedγoj of the rolling element is as follows:(22)γoij=Dwcosαoijdm.As shown in Figure4, based on the above analysis, the following variable transfer relationship can be obtained, “where the left side of the arrow is the function to the right of the arrow.”Figure 4
Analysis of the variable transfer relationship of centrifugal forceFcj.(2) Variable Recursion Relation Analysis of Gyroscopic MomentMgj. In ball angular position ψj, the gyroscopic moment Mgj of the high-speed angular contact ball bearing can be expressed as (23)Mgj=4.5×10−11D5nbjnmjsinβj,where the calculation equations of the rotation speed nbjof the rolling element can be expressed:(24)nbj=dm1−γij1+γojΩi±ΩoDw1−γijcosαoj−β+Dw1+γojcosαi−β.As shown in Figure5, by substituting (20), (21), and (24) into (23), the variables recursive relation of gyroscopic moment Mgj can be obtained.Figure 5
Analysis of variable transfer relationship of gyro momentMgj.(3) Analysis of the Variables Recursive Relation of Contact StiffnessesKijandKojat Contact Points. In ball angular position ψj, the contact stiffnesses Kij and Koj of high-speed angular ball bearing at contact points between the rolling elements and the inner and outer rings can be expressed as follows: (25)Kij=2.15×105∑ρij−1/2nδ−3/2,(26)Koj=2.15×105∑ρoj−1/2nδ−3/2,where the function of principal curvature ∑ρ can be expressed:(27)∑ρ=ρ1Ι+ρ2Ι+ρ1Π+ρ2Π.For angular ball bearing:ρ1Ι=ρ1Π=2/Dw, ρ2Ι=−1/fjDw, and ρ2Π=∓2γ/Dw1±γ.Then, the equation can be obtained:(28)∑ρj=4Dw+2γjqDw1∓γjq−1fiDw,where minus was taken when j=i.nδ can be obtained by looking up the table based on the principal curvature difference function Fρ:(29)Fρ=ρ12−ρ11+ρ22−ρ21∑ρ.Based on the above analysis, the variable transfer relationship of contact stiffnessesKij and Koj are shown in Figure 6.Figure 6
Analysis of variable transfer relationship of contact point stiffnessesKij and Koj.(4) Analysis of the Variables Recursive Relation of Contact ForcesQijandQoj. In ball angular position ψj, the contact forces Qij and Qojof high-speed angular ball bearing at contact points between the rolling elements and the inner and outer rings can be expressed as follows: (30)Qij=Kijδij1.5,(31)Qoj=Kojδoj1.5.Through the relations among (25)–(31), the transfer relations of contact force and variables can be obtained, as shown in Figure 7.Figure 7
Analysis of variable transfer relationship of contact forcesQij and Qoj.(5) Overall Analysis of Variable Recursion Relations. Based on the analysis of the above variables, such as centrifugal force Fcj, gyro moment Mgj, contact stiffnesses Kij and Koj, and contact forces Qij and Qoj, the total transfer diagram of the variable transfer relationship of the dynamic characteristic equations can be obtained, as shown in Figure 8.Figure 8
Dynamic characteristic equations variable transfer relationship total analysis diagram.


#### 2.2.2. Derivation of Simplified Computational Models

Combining the relationship of variable transfer and dynamic characteristic equations, it is found that the eight basic variablesX1j, X2j, δij, δoj, Kij, Koj, Fcj, and Mgj are all functions of the actual internal and external contact angles αij and αoj. It can be concluded that the actual internal and external contact angles αij and αoj are the keys to solve and analyze the dynamic characteristic equations.The system of nonlinear equations consists of (7), (8), (13), (14), and (16)–(18). It transforms 4Z basic unknown variables about X1j, X2j, δij, and δoj in the equations into 2Z unknown variables about the actual internal αij and external contact angles αoj. Then, it makes 4Z + 2 unknown variables X1j, X2j, δij, δoj, δa, and δr of the original nonlinear equations reduced to 2Z + 2 unknown variables αij, αoj, δa, and δr. This reduces the dimension of the complex dynamic characteristic equations, makes the coupling hierarchy of the variables clear, reduces the iterative solution of the intermediate variables, and improves the convergence of the equations.


## 2.1. Traditional Computing Model

Solving the high-speed characteristics of angular contact ball bearings has to start from the contact deformation. Establishing the bearing deformation geometry equation, the rolling element balance equation, and the inner ring balance equation and then establishing the auxiliary equations including centrifugal force, gyroscopic moment, and contact angle are needed. Figure1 is a schematic diagram of the rolling elements in the different ball angular position ψj; Figure 2 is a schematic diagram of rolling body force; Figure 3 shows the relative position relationships between the centers of the inner and outer raceways groove curvature and the rolling elements centers before and after the bearing is loaded (according to the inner ring rotation and the outer ring fixation). According to the literature [12], combined with Hertz contact theory, bearing motion, and force relationship and raceway control theory, the calculation model of high-speed angular contact ball bearing is given.Figure 1
Position map of rolling elements at different azimuth anglesψj.Figure 2
Schematic diagram of the rolling element.Figure 3
Schematic diagram of the change of the center position of the groove before and after the load.


### 2.1.1. Deformation Geometry Equation

Figure3 shows that the axial distance A1j and radial distance A2j between the curvature centers of inner and outer raceway grooves in the ball angular position ψj are as follows:(1)A1j=BDsinα0+δα+ℜiθcosψj,(2)A2j=BDcosα0+δrcosψj,where δa, δr, and θ, respectively, denote the relative axial deformation, relative radial deformation, and relative angular displacement of the inner and outer rings; ℜi denotes the radius of the groove curvature center trajectory of inner raceway; α0 is the initial contact angle of the bearing.The following geometric relationship can be determined from Figure3:(3)cosαoj=X2jfo−0.5D+δoj,(4)sinαoj=X1jfo−0.5D+δoj,(5)cosαij=A2j−X2jfi−0.5D+δij,(6)sinαij=A1j−X1jfi−0.5D+δij.In these equations,α is the contact angle of the bearing; f is the groove curvature of the bearing; the subscripts o and i are the outer ring and inner ring of the bearing, respectively; D is the ball diameter. According to the deformation geometric relationship between the components in Figure 3, combined with the Pythagorean theorem and (1) and (2), the position change relationship equation among the center of curvature of the rolling elements and the inner and outer grooves in the ball angular position ψj can be obtained:(7)A1j−X1j2+A2j−X2j2−fi−0.5D+δij2=0,(8)X1j2+X2j2−fo−0.5D+δoj2=0.


### 2.1.2. Rolling Element Balance Equation

According to Hertz contact theory, the normal contact forcesQ of the rolling body and the inner and outer rings at the contact points are as follows:(9)Qij=Kijδij1.5,(10)Qoj=Kojδoj1.5.As shown in Figure2, in the ball angular position ψj, the forced balance equations of rolling elements in the radial and axial can be written as(11)Qijsinαij−Qojsinαoj−MgjDλijcosαij−λojcosαoj=0,(12)Qijcosαij−Qojcosαoj+MgjDλijsinαij−λojsinαoj+Foj=0,where Mg is the gyro moment of the ball. Substituting (3)–(6), (9), and (10) into (11) and (12) yields(13)Kijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij+λojMgjX2j/D−Kojδoj1.5X1jfo−0.5D+δoj=0,(14)Kijδij1.5A2j−X2j−λijMgj/DA1j−X1jfi−0.5D+δij−λojMgjX1j/D+Kojδoj1.5X2jfo−0.5D+δoj+Fcj=0.


### 2.1.3. Inner Ring Balance Equation

According to the force balance relationship and the moment load balance relationship of the bearing inner ring in the axial and radial directions, the inner ring balance equation can be obtained:(15)Fa−∑j=1j=zQijsinαij−MgjDλijcosαij=0,Fr−∑j=1j=zQijcosαij+MgjDλijsinαijcosψj=0,M−∑j=1j=zQijsinαij−MgjDλijcosαijℜi+fiλijMgjcosψj=0.Substituting (3)–(6), (9), and (10) into (16)–(18) results in(16)Fa−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij=0,(17)Fr−∑j=1j=zKijδij1.5A2j−X2j+λijMgj/DA1j−X1jfi−0.5D+δijcosψj=0,(18)M−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δijℜi+fiλijMgjcosψj=0,where Fa, Fr, and M, respectively, are the axial, radial, and moment forces acting on the bearing. Given the initial iteration values of δa, δr, and θ, the iteration algorithm can be used to solve the nonlinear equations composed of formulas (7), (8), (13), and (14) in the ball angular position ψj and obtain X1j, X2j, δij, and δoj. By substituting these into the nonlinear equations (16)–(18) in every ball angular position ψj, the values of the main unknowns δa, δr, and θ can be obtained. Then, the iteration of X1j, X2j, δij, and δoj continued until the values of the main unknowns δa, δr, and θ meet the accuracy.


## 2.1.1. Deformation Geometry Equation

Figure3 shows that the axial distance A1j and radial distance A2j between the curvature centers of inner and outer raceway grooves in the ball angular position ψj are as follows:(1)A1j=BDsinα0+δα+ℜiθcosψj,(2)A2j=BDcosα0+δrcosψj,where δa, δr, and θ, respectively, denote the relative axial deformation, relative radial deformation, and relative angular displacement of the inner and outer rings; ℜi denotes the radius of the groove curvature center trajectory of inner raceway; α0 is the initial contact angle of the bearing.The following geometric relationship can be determined from Figure3:(3)cosαoj=X2jfo−0.5D+δoj,(4)sinαoj=X1jfo−0.5D+δoj,(5)cosαij=A2j−X2jfi−0.5D+δij,(6)sinαij=A1j−X1jfi−0.5D+δij.In these equations,α is the contact angle of the bearing; f is the groove curvature of the bearing; the subscripts o and i are the outer ring and inner ring of the bearing, respectively; D is the ball diameter. According to the deformation geometric relationship between the components in Figure 3, combined with the Pythagorean theorem and (1) and (2), the position change relationship equation among the center of curvature of the rolling elements and the inner and outer grooves in the ball angular position ψj can be obtained:(7)A1j−X1j2+A2j−X2j2−fi−0.5D+δij2=0,(8)X1j2+X2j2−fo−0.5D+δoj2=0.


## 2.1.2. Rolling Element Balance Equation

According to Hertz contact theory, the normal contact forcesQ of the rolling body and the inner and outer rings at the contact points are as follows:(9)Qij=Kijδij1.5,(10)Qoj=Kojδoj1.5.As shown in Figure2, in the ball angular position ψj, the forced balance equations of rolling elements in the radial and axial can be written as(11)Qijsinαij−Qojsinαoj−MgjDλijcosαij−λojcosαoj=0,(12)Qijcosαij−Qojcosαoj+MgjDλijsinαij−λojsinαoj+Foj=0,where Mg is the gyro moment of the ball. Substituting (3)–(6), (9), and (10) into (11) and (12) yields(13)Kijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij+λojMgjX2j/D−Kojδoj1.5X1jfo−0.5D+δoj=0,(14)Kijδij1.5A2j−X2j−λijMgj/DA1j−X1jfi−0.5D+δij−λojMgjX1j/D+Kojδoj1.5X2jfo−0.5D+δoj+Fcj=0.


## 2.1.3. Inner Ring Balance Equation

According to the force balance relationship and the moment load balance relationship of the bearing inner ring in the axial and radial directions, the inner ring balance equation can be obtained:(15)Fa−∑j=1j=zQijsinαij−MgjDλijcosαij=0,Fr−∑j=1j=zQijcosαij+MgjDλijsinαijcosψj=0,M−∑j=1j=zQijsinαij−MgjDλijcosαijℜi+fiλijMgjcosψj=0.Substituting (3)–(6), (9), and (10) into (16)–(18) results in(16)Fa−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij=0,(17)Fr−∑j=1j=zKijδij1.5A2j−X2j+λijMgj/DA1j−X1jfi−0.5D+δijcosψj=0,(18)M−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δijℜi+fiλijMgjcosψj=0,where Fa, Fr, and M, respectively, are the axial, radial, and moment forces acting on the bearing. Given the initial iteration values of δa, δr, and θ, the iteration algorithm can be used to solve the nonlinear equations composed of formulas (7), (8), (13), and (14) in the ball angular position ψj and obtain X1j, X2j, δij, and δoj. By substituting these into the nonlinear equations (16)–(18) in every ball angular position ψj, the values of the main unknowns δa, δr, and θ can be obtained. Then, the iteration of X1j, X2j, δij, and δoj continued until the values of the main unknowns δa, δr, and θ meet the accuracy.


## 2.2. High-Speed Angular Contact Ball Bearings Calculation Model Analysis

Equations (7), (8), (13), (14), and (16)–(18) constitute a system of nonlinear equations, where there are 4Z + 2 equations. The calculation of contact stiffnesses Kij and Koj, centrifugal force Fcj, and gyroscopic moment Mgj of rolling elements and inner and outer rings at contact points all depends on the calculated values. And Kij, Koj, Fcj, and Mgj are functions of the actual internal and external contact angles αij and αoj; the actual internal angle αijand external contact angle αoj are functions of X1j and X2j. Therefore, the variables in the dynamic characteristic equations are highly coupled, which makes the iterative calculation process very complicated and the calculation amount is extremely large. The initial value is not chosen properly, and it is difficult to converge or even not converge.


### 2.2.1. Variable Recursion Relation Analysis of Computational Model

(1) Variable Recursion Relation Analysis of Centrifugal ForceFcj. In ball angular position ψj, the centrifugal forces Fcj of the high-speed angular contact ball bearing rolling element can be expressed as follows: (19)Fcj=2.26×10−11D3dmnmj2.In the above formula, the formula for calculating the revolution speednmjof the rolling element is as follows:(20)nmj=nij1−γij1+γojΩi±Ωo1+γojcosαoj−βoj+1−γojcosαoj−βoj.According to the outer raceway control theory, the attitude angleβojof the rolling element is as follows:(21)βoj=arctancosαojsinαojcos2αoj+γoj.In the above formula, the formula for calculating the revolution speedγoj of the rolling element is as follows:(22)γoij=Dwcosαoijdm.As shown in Figure4, based on the above analysis, the following variable transfer relationship can be obtained, “where the left side of the arrow is the function to the right of the arrow.”Figure 4
Analysis of the variable transfer relationship of centrifugal forceFcj.(2) Variable Recursion Relation Analysis of Gyroscopic MomentMgj. In ball angular position ψj, the gyroscopic moment Mgj of the high-speed angular contact ball bearing can be expressed as (23)Mgj=4.5×10−11D5nbjnmjsinβj,where the calculation equations of the rotation speed nbjof the rolling element can be expressed:(24)nbj=dm1−γij1+γojΩi±ΩoDw1−γijcosαoj−β+Dw1+γojcosαi−β.As shown in Figure5, by substituting (20), (21), and (24) into (23), the variables recursive relation of gyroscopic moment Mgj can be obtained.Figure 5
Analysis of variable transfer relationship of gyro momentMgj.(3) Analysis of the Variables Recursive Relation of Contact StiffnessesKijandKojat Contact Points. In ball angular position ψj, the contact stiffnesses Kij and Koj of high-speed angular ball bearing at contact points between the rolling elements and the inner and outer rings can be expressed as follows: (25)Kij=2.15×105∑ρij−1/2nδ−3/2,(26)Koj=2.15×105∑ρoj−1/2nδ−3/2,where the function of principal curvature ∑ρ can be expressed:(27)∑ρ=ρ1Ι+ρ2Ι+ρ1Π+ρ2Π.For angular ball bearing:ρ1Ι=ρ1Π=2/Dw, ρ2Ι=−1/fjDw, and ρ2Π=∓2γ/Dw1±γ.Then, the equation can be obtained:(28)∑ρj=4Dw+2γjqDw1∓γjq−1fiDw,where minus was taken when j=i.nδ can be obtained by looking up the table based on the principal curvature difference function Fρ:(29)Fρ=ρ12−ρ11+ρ22−ρ21∑ρ.Based on the above analysis, the variable transfer relationship of contact stiffnessesKij and Koj are shown in Figure 6.Figure 6
Analysis of variable transfer relationship of contact point stiffnessesKij and Koj.(4) Analysis of the Variables Recursive Relation of Contact ForcesQijandQoj. In ball angular position ψj, the contact forces Qij and Qojof high-speed angular ball bearing at contact points between the rolling elements and the inner and outer rings can be expressed as follows: (30)Qij=Kijδij1.5,(31)Qoj=Kojδoj1.5.Through the relations among (25)–(31), the transfer relations of contact force and variables can be obtained, as shown in Figure 7.Figure 7
Analysis of variable transfer relationship of contact forcesQij and Qoj.(5) Overall Analysis of Variable Recursion Relations. Based on the analysis of the above variables, such as centrifugal force Fcj, gyro moment Mgj, contact stiffnesses Kij and Koj, and contact forces Qij and Qoj, the total transfer diagram of the variable transfer relationship of the dynamic characteristic equations can be obtained, as shown in Figure 8.Figure 8
Dynamic characteristic equations variable transfer relationship total analysis diagram.


### 2.2.2. Derivation of Simplified Computational Models

Combining the relationship of variable transfer and dynamic characteristic equations, it is found that the eight basic variablesX1j, X2j, δij, δoj, Kij, Koj, Fcj, and Mgj are all functions of the actual internal and external contact angles αij and αoj. It can be concluded that the actual internal and external contact angles αij and αoj are the keys to solve and analyze the dynamic characteristic equations.The system of nonlinear equations consists of (7), (8), (13), (14), and (16)–(18). It transforms 4Z basic unknown variables about X1j, X2j, δij, and δoj in the equations into 2Z unknown variables about the actual internal αij and external contact angles αoj. Then, it makes 4Z + 2 unknown variables X1j, X2j, δij, δoj, δa, and δr of the original nonlinear equations reduced to 2Z + 2 unknown variables αij, αoj, δa, and δr. This reduces the dimension of the complex dynamic characteristic equations, makes the coupling hierarchy of the variables clear, reduces the iterative solution of the intermediate variables, and improves the convergence of the equations.


## 2.2.1. Variable Recursion Relation Analysis of Computational Model

(1) Variable Recursion Relation Analysis of Centrifugal ForceFcj. In ball angular position ψj, the centrifugal forces Fcj of the high-speed angular contact ball bearing rolling element can be expressed as follows: (19)Fcj=2.26×10−11D3dmnmj2.In the above formula, the formula for calculating the revolution speednmjof the rolling element is as follows:(20)nmj=nij1−γij1+γojΩi±Ωo1+γojcosαoj−βoj+1−γojcosαoj−βoj.According to the outer raceway control theory, the attitude angleβojof the rolling element is as follows:(21)βoj=arctancosαojsinαojcos2αoj+γoj.In the above formula, the formula for calculating the revolution speedγoj of the rolling element is as follows:(22)γoij=Dwcosαoijdm.As shown in Figure4, based on the above analysis, the following variable transfer relationship can be obtained, “where the left side of the arrow is the function to the right of the arrow.”Figure 4
Analysis of the variable transfer relationship of centrifugal forceFcj.(2) Variable Recursion Relation Analysis of Gyroscopic MomentMgj. In ball angular position ψj, the gyroscopic moment Mgj of the high-speed angular contact ball bearing can be expressed as (23)Mgj=4.5×10−11D5nbjnmjsinβj,where the calculation equations of the rotation speed nbjof the rolling element can be expressed:(24)nbj=dm1−γij1+γojΩi±ΩoDw1−γijcosαoj−β+Dw1+γojcosαi−β.As shown in Figure5, by substituting (20), (21), and (24) into (23), the variables recursive relation of gyroscopic moment Mgj can be obtained.Figure 5
Analysis of variable transfer relationship of gyro momentMgj.(3) Analysis of the Variables Recursive Relation of Contact StiffnessesKijandKojat Contact Points. In ball angular position ψj, the contact stiffnesses Kij and Koj of high-speed angular ball bearing at contact points between the rolling elements and the inner and outer rings can be expressed as follows: (25)Kij=2.15×105∑ρij−1/2nδ−3/2,(26)Koj=2.15×105∑ρoj−1/2nδ−3/2,where the function of principal curvature ∑ρ can be expressed:(27)∑ρ=ρ1Ι+ρ2Ι+ρ1Π+ρ2Π.For angular ball bearing:ρ1Ι=ρ1Π=2/Dw, ρ2Ι=−1/fjDw, and ρ2Π=∓2γ/Dw1±γ.Then, the equation can be obtained:(28)∑ρj=4Dw+2γjqDw1∓γjq−1fiDw,where minus was taken when j=i.nδ can be obtained by looking up the table based on the principal curvature difference function Fρ:(29)Fρ=ρ12−ρ11+ρ22−ρ21∑ρ.Based on the above analysis, the variable transfer relationship of contact stiffnessesKij and Koj are shown in Figure 6.Figure 6
Analysis of variable transfer relationship of contact point stiffnessesKij and Koj.(4) Analysis of the Variables Recursive Relation of Contact ForcesQijandQoj. In ball angular position ψj, the contact forces Qij and Qojof high-speed angular ball bearing at contact points between the rolling elements and the inner and outer rings can be expressed as follows: (30)Qij=Kijδij1.5,(31)Qoj=Kojδoj1.5.Through the relations among (25)–(31), the transfer relations of contact force and variables can be obtained, as shown in Figure 7.Figure 7
Analysis of variable transfer relationship of contact forcesQij and Qoj.(5) Overall Analysis of Variable Recursion Relations. Based on the analysis of the above variables, such as centrifugal force Fcj, gyro moment Mgj, contact stiffnesses Kij and Koj, and contact forces Qij and Qoj, the total transfer diagram of the variable transfer relationship of the dynamic characteristic equations can be obtained, as shown in Figure 8.Figure 8
Dynamic characteristic equations variable transfer relationship total analysis diagram.


## 2.2.2. Derivation of Simplified Computational Models

Combining the relationship of variable transfer and dynamic characteristic equations, it is found that the eight basic variablesX1j, X2j, δij, δoj, Kij, Koj, Fcj, and Mgj are all functions of the actual internal and external contact angles αij and αoj. It can be concluded that the actual internal and external contact angles αij and αoj are the keys to solve and analyze the dynamic characteristic equations.The system of nonlinear equations consists of (7), (8), (13), (14), and (16)–(18). It transforms 4Z basic unknown variables about X1j, X2j, δij, and δoj in the equations into 2Z unknown variables about the actual internal αij and external contact angles αoj. Then, it makes 4Z + 2 unknown variables X1j, X2j, δij, δoj, δa, and δr of the original nonlinear equations reduced to 2Z + 2 unknown variables αij, αoj, δa, and δr. This reduces the dimension of the complex dynamic characteristic equations, makes the coupling hierarchy of the variables clear, reduces the iterative solution of the intermediate variables, and improves the convergence of the equations.


## 3. Solution of Computational Model Based on Improved Artificial Bee Colony Algorithm

### 3.1. Artificial Bee Colony Algorithm and Its Algorithm Structure Improvement

When the standard artificial bee colony algorithm solves the high-dimensional optimization problem, there is a disadvantage that the convergence speed becomes slower as the dimension increases and it is easy to fall into the local optimum and premature convergence when solving the problem of large root difference. Therefore, the standard artificial bee colony algorithm needs to be improved.


#### 3.1.1. Improvements in Search Strategy

The standard artificial bee colony algorithm uses the same search formula for leading the bee phase and following the bee phase:(32)xm,i′=xm,i+rand−1,1∗xm,i−xk,i,where xk,i is randomly selected from the existing food source xm and k≠m, rand (−1, 1) represents a random number generated in the interval (−1, 1).The standard artificial bee colony algorithm uses the same search formula (32) in the leading bee stage and the following bee stage. Searching for the current food source xm, the neighborhood of the current solution xm, and in each cycle the global optimal solution is simply stored, but not in the search strategy, which plays the guiding role of the current global optimal solution to the real solution xbest. In order to improve the quality of the solution of the high-dimensional nonlinear problem and the global solution and convergence speed of the bee colony algorithm, the search formula (32) for the bee stage and follow-bee stages in the standard bee colony algorithm is improved to the following formula:(33)xm,i′=xbest,i+rand−1,1∗xbest,i−xm,i,(34)xm,i′=xk,i+rand−1,1∗xm,i−xk,i.After improving the standard artificial bee colony algorithm, the leader bee uses the search formula (33) to search. It is based on the global optimal solution xbest,i and plays its guiding role. Meanwhile, combining with the current solution xm,i, the leader bee carries out neighborhood search to obtain a better new solution xm,i′. Similarly, the follower bee uses search formula (34) to search. It is based on the random selection of other food source solutions xk,i and uses the current solution xm,i to search the neighborhood to obtain a better new solution xm,i′. Random selection of other food sources can enhance the global search ability of the algorithm and avoid single local search.


#### 3.1.2. Algorithm Structure Improvement

The geometric parameters of high-speed angular contact ball bearings belong to centimeter level and the loaded deformation parameters belong to micrometer level, so the results of calculation model have great root differences. However, the standard artificial bee colony algorithm does not consider individual diversity and difference, so it is necessary to make structural improvements to the standard artificial bee colony algorithm.Standard artificial bee colony algorithm generates N food sources directly in the initialization stage, that is, generating N D-dimensional initial solutionsxm=xm1,xm2,xmi,…,xmD, among them m=1,2,3,…,D, i=1,2,3,…,D. In the initialization stage, the cross-scale problem and priority relationship in each possible solution xm1,xm2,xmi,…,xmD are not considered step by step for the problem with large root differences. As a result, the algorithm converges to unreliable solutions prematurely and even falls into local dead cycles, which seriously affects the global convergence rate.The author uses the idea of segmentation to initialize the parameters of the food sourcexm: first generate parameters δa and δr and then generate parameters αij and αoj, which solves all the differences and cross-scale problems of the required parameters, so that the algorithm is reliable at the source and ensures the accuracy of the calculation results.Some codes after adjusting the structure of the algorithm are shown in Table1, which represents the parameters δa and δr; Table 2 represents the parameters αij and αoj to be determined.Table 1
Initialization phase algorithm structure adjustment part of the code.
Generateδa, δr(1)def calculate (self)(2) rnd = np.random.random (size = self.len-2∗self.z) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(3)fori, kin zip (range(0, self.len-2∗self.z), range (0, self.len-2∗ self.z)):(4)   self.chrom 3(i) = self. bound (0, 2∗self.z + i) + (self.bound (1, 2∗ self. z + i) − self.bound (0, 2 ∗ self. z + i))∗ rnd(k)(5)foriin range (0, self.len-2∗ self.z):(6)   self.X (2∗ self.z + i) = self.chrom3(i) # δa, δr#self.X (2∗self.z + 1) = self.chrom3(2∗self.z + 1)#self.X (4∗self.z + 2) = self.chrom(4∗self.z + 2)(7)  super().some_necessary_calculations()(8)r = 0(9)while r < self.z:(10)if self.A1 (r) < 0 or self.A2 (r) < 0:(11)       rnd = np.random.random (size = self.len-2∗ self.z) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(12)fori, kin zip (range (0, self.len-2 ∗ self.z), range (0, self.len - 2 ∗ self.z)):(13)        self.chrom3(i) = self.bound (0, 2 ∗ self.z + i) + (self.bound (1, 2 ∗ self.z + i) − self.bound (0, 2 ∗ self.z + i)) ∗ rnd (k) # xij=xmin, j + rand (0, 1) (xmax, j − xmin, j)(14)foriin range(0, self.len-2 ∗ self.z):(15)     self.X (2∗ self.z + i) = self.chrom3(i)#self.X (2∗self.z + 1) = self.chrom3(2∗self.z + 1)#self.X (4∗self.z + 2) = self.chrom(4∗self.z + 2)(16)     super().some_necessary_calculations()(17)r = 0Table 2
Initialization phase algorithm structure adjustment of another part of the code.
Generateαi, αo#αi(1)  rnd = np.random.random (size = self.z//2 + 1) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(2)  rand1 = sorted (rnd)#size = 8(3)foriin range (0,self.z//2 + 1):(4)   self.rand1(i) = rand1(i)(5)foriin range(self.z//2 + 1,self.z):(6)    self.rand1(i) = self.rand1(self.z − i) #size = 16(7)fori, k in zip (range(0,self.z//2 + 1), range(0, self.z//2 + 1)):(8)    self.chrom1(i) = self.bound (0, i) + (self.bound (1, i) − self.bound (0, i)) ∗ self.rand1(k) # Xij = Xmin, j + rand(9)foriin range (self.z//2 + 1,self.z):(10)    self.chrom1(i) = self.chrom1 (self.z − i)#αo(11)  rnd = np.random.random (size = self.z//2 + 1)# example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(12)  rand2 = sorted (rnd,reverse =True) #size = 8(13)foriin range(0,self.z//2 + 1):(14)   self.rand2(i) = rand 2(i)(15)foriin range (self.z//2 + 1,self.z):(16)   self.rand2(i) = self.rand2 (self.z − i) #size = 16(17)fori, k in zip (range (0, self.z//2 + 1), range (0, self.z//2 + 1)):(18)   self.chrom2(i) = self.boun ((0, i + self.z) + (self.bound (1, i + self.z) − self.bound (0, i + self.z)) ∗ self.rand2(k) # Xij = Xmin, j + rand (0, 1) (Xmax,j−Xmin, j) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(19)foriin range (self.z//2 + 1, self.z):(20)   self.chrom2(i) = self.chrom2 (self.z − i)


### 3.2. Testing and Analysis of Improved Algorithms

In order to verify the reliability and feasibility of the improved algorithm, four standard test functions are used to test the performance of the improved artificial bee colony algorithm. Among the four standard test functions,f1 is the Rastrigin function, f2 is the Ackley function, f3 is the Rosenbrock function, and f4 is the selfdef function, and its analytical expression is from (35) to (38).Rastrigin functionf1 is(35)f1x=10n+∑i=1nxi2−10cos2πxi.Ackley functionf2 is(36)f2x=−20exp0.21n∑i=1nxi2−exp1n∑i=1ncos2πxi+e+20.Rosenbrock functionf3 is(37)f3x=∑i=1n100xi2−xi+12+xi−12.Selfdef functionf4 is(38)f4x,y=x−99992+y2−52+x−9998y−3x1000+12.In Figure9, (a) is a three-dimensional diagram of a Rastrigin test function, (b) is a three-dimensional diagram of an Ackley test function, (c) is a three-dimensional diagram of a Rosenbrock test function, and (d) is a three-dimensional diagram of a Selfdef test function. f1 is a multimodal function. In the domain xi∈−5.12,5.12, this function has very many minima, and there is only one minimum, and the global minimum is minf1=f10,0,…,0=0. The f2 test function is used to verify the global convergence speed of the algorithm. The global minimum in the domain xi∈−5,5 is minf2=f20,0,…,0=0. f3 is a complex morbid function. Because it is difficult to find its global minimum, it is necessary to use this function to test the improved artificial bee colony algorithm. In the domain of xi∈−10,10, the global minimum is minf3=f30,0,…,0=0. f4 is used to test a function with a large root difference. This function is in the domain interval x∈−15000,15000, and the minimum value is minf4=f410000,2=0. From Figure 9, it can intuitively observe the extreme value of each function.Figure 9
Three-dimensional visualization of standard test functions.
(a)(b)(c)(d)Table3 shows the maximum number of iterations, test dimensions, algorithm control parameters, theoretical global optimal solutions, and search range for the four test functions where f1, f2, and f3 are used to test whether the algorithm has reliable global convergence speed and high accuracy in the face of high-dimensional, nonlinear, and multiextreme problems; f4 is used to test whether the improved algorithm has reliable search ability for large root difference problems.Table 3
Test functions and algorithm parameter settings.
FunctionMaximum number of iterationsFunctional dimensionAlgorithmic control parametersPopulation sizeTheoretical global optimal solutionSearch rangef1500050100300[−5.12, 5.12]f2500050100300[−5, 5]f3500050100300[−10, 10]f450002100300[−15000, 15000]The test results show that the standard artificial bee colony algorithm and the improved artificial bee colony algorithm are used to calculate the results of 10 tests for each test function, as shown in Table4. Meanwhile, the comparison curves of convergence rate of each test function are obtained, as shown in Figure 10.Table 4
Comparison of standard artificial bee colony algorithm and improved artificial bee colony algorithm.
Test functionAlgorithmOptimal solutionThe worst solutionMeanStandard deviationCalculating time (s)RastriginStandard1.95284E − 089.62104E − 063.24631E − 064.02246E − 0627.94Improved3.64366E − 102.30741E − 094.78643E − 107.56349E − 1023.43AckleyStandard4.14784E − 119.97926E − 112.38264E − 112.23417E − 1133.65Improved6.43285E − 137.15932E − 128.13264E − 132.78687E − 1228.26RosenbrockStandard1.41868E − 025.23661E − 014.23692E − 011.97058E − 0196.56Improved5.38452E − 041.35267E − 036.43285E − 043.17665E − 0473.24SelfdefStandard8.38282E − 072.05764E − 056.43657E − 067.05824E − 0623.68Improved1.57346E − 088.43697E − 083.85145E − 082.46283E − 0817.35Figure 10
Convergence rate contrast curve. (a) Rastrigin function. (b) Ackley function. (c) Rosenbrock function. (d) Selfdef function.
(a)(b)(c)(d)From Table4 and Figure 10, it can be seen that the improved algorithm improves in the order of magnitude compared with the standard algorithm in terms of the optimal solution accuracy, the worst solution accuracy, the mean value accuracy, and the variance accuracy. Meanwhile, the improved artificial bee colony algorithm reduces significantly the iteration calculation time compared with the standard artificial bee colony algorithm. These syntheses show that the improved artificial bee colony algorithm can solve the high-speed angular contact ball bearing calculation model with high dimension, nonlinearity, multiextremum, and large root difference.


### 3.3. Solution of Simplified Calculation Model

In order to adopt the improved artificial bee colony algorithm, the high-speed angular contact ball bearing calculation model needs to be simplified to the optimization calculation model. First, (7), (8), (13), (14), (16), and (17) are assembled into the following function groups:(39)F1j=A1j−X1j2+A2j−X2j2−fi−0.5D+δij2,F2j=X1j2+X2j2−fo−0.5D+δoj2,F3j=Kijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij+λojMgjX2j/D−Kojδoj1.5X1jfo−0.5D+δoj,F4j=Kijδij1.5A2j−X2j−λijMgj/DA1j−X1jfi−0.5D+δij−λojMgjX1j/D+Kojδoj1.5X2jfo−0.5D+δoj+Fcj,F5j=Fa−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij.F6j=Fr−∑j=1j=zKijδij1.5A2j−X2j−λijMgj/DA1j−X1jfi−0.5D+δijcosψj.Substituting (32)–(35) into function group (36), then transforming functions F1j, F2j, F3j, F4j, F5j, and F6j into objective functions, and finally using improved artificial bee colony algorithm to find Minimum value, the specific optimization model is as follows:Objective function:(40)F=∑j=1ZF1j2+∑j=1ZF2j2+∑j=1ZF3j2+∑j=1ZF4j2+F52+F62.The design variables are(41)x=αi1,αi2,αi3,…,αij,αo1,αo2,αo3,…,αoj,δa,δr.The basic constraints are(42)0≤αij≤α0,α0≤αij≤αrcsinA1j−X1jfi−0.5D+δij.Findxto minimize the Fx value; that is, find minFx.The design variables are 2Z + 2, which are nearly half less than the solution variables of the traditional numerical algorithm and are conducive for improving the global convergence rate. The specific flow of the high-speed angular contact ball bearing calculation model using the improved artificial bee colony algorithm is shown in Figure 11.Figure 11
Solution flow chart.


## 3.1. Artificial Bee Colony Algorithm and Its Algorithm Structure Improvement

When the standard artificial bee colony algorithm solves the high-dimensional optimization problem, there is a disadvantage that the convergence speed becomes slower as the dimension increases and it is easy to fall into the local optimum and premature convergence when solving the problem of large root difference. Therefore, the standard artificial bee colony algorithm needs to be improved.


### 3.1.1. Improvements in Search Strategy

The standard artificial bee colony algorithm uses the same search formula for leading the bee phase and following the bee phase:(32)xm,i′=xm,i+rand−1,1∗xm,i−xk,i,where xk,i is randomly selected from the existing food source xm and k≠m, rand (−1, 1) represents a random number generated in the interval (−1, 1).The standard artificial bee colony algorithm uses the same search formula (32) in the leading bee stage and the following bee stage. Searching for the current food source xm, the neighborhood of the current solution xm, and in each cycle the global optimal solution is simply stored, but not in the search strategy, which plays the guiding role of the current global optimal solution to the real solution xbest. In order to improve the quality of the solution of the high-dimensional nonlinear problem and the global solution and convergence speed of the bee colony algorithm, the search formula (32) for the bee stage and follow-bee stages in the standard bee colony algorithm is improved to the following formula:(33)xm,i′=xbest,i+rand−1,1∗xbest,i−xm,i,(34)xm,i′=xk,i+rand−1,1∗xm,i−xk,i.After improving the standard artificial bee colony algorithm, the leader bee uses the search formula (33) to search. It is based on the global optimal solution xbest,i and plays its guiding role. Meanwhile, combining with the current solution xm,i, the leader bee carries out neighborhood search to obtain a better new solution xm,i′. Similarly, the follower bee uses search formula (34) to search. It is based on the random selection of other food source solutions xk,i and uses the current solution xm,i to search the neighborhood to obtain a better new solution xm,i′. Random selection of other food sources can enhance the global search ability of the algorithm and avoid single local search.


### 3.1.2. Algorithm Structure Improvement

The geometric parameters of high-speed angular contact ball bearings belong to centimeter level and the loaded deformation parameters belong to micrometer level, so the results of calculation model have great root differences. However, the standard artificial bee colony algorithm does not consider individual diversity and difference, so it is necessary to make structural improvements to the standard artificial bee colony algorithm.Standard artificial bee colony algorithm generates N food sources directly in the initialization stage, that is, generating N D-dimensional initial solutionsxm=xm1,xm2,xmi,…,xmD, among them m=1,2,3,…,D, i=1,2,3,…,D. In the initialization stage, the cross-scale problem and priority relationship in each possible solution xm1,xm2,xmi,…,xmD are not considered step by step for the problem with large root differences. As a result, the algorithm converges to unreliable solutions prematurely and even falls into local dead cycles, which seriously affects the global convergence rate.The author uses the idea of segmentation to initialize the parameters of the food sourcexm: first generate parameters δa and δr and then generate parameters αij and αoj, which solves all the differences and cross-scale problems of the required parameters, so that the algorithm is reliable at the source and ensures the accuracy of the calculation results.Some codes after adjusting the structure of the algorithm are shown in Table1, which represents the parameters δa and δr; Table 2 represents the parameters αij and αoj to be determined.Table 1
Initialization phase algorithm structure adjustment part of the code.
Generateδa, δr(1)def calculate (self)(2) rnd = np.random.random (size = self.len-2∗self.z) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(3)fori, kin zip (range(0, self.len-2∗self.z), range (0, self.len-2∗ self.z)):(4)   self.chrom 3(i) = self. bound (0, 2∗self.z + i) + (self.bound (1, 2∗ self. z + i) − self.bound (0, 2 ∗ self. z + i))∗ rnd(k)(5)foriin range (0, self.len-2∗ self.z):(6)   self.X (2∗ self.z + i) = self.chrom3(i) # δa, δr#self.X (2∗self.z + 1) = self.chrom3(2∗self.z + 1)#self.X (4∗self.z + 2) = self.chrom(4∗self.z + 2)(7)  super().some_necessary_calculations()(8)r = 0(9)while r < self.z:(10)if self.A1 (r) < 0 or self.A2 (r) < 0:(11)       rnd = np.random.random (size = self.len-2∗ self.z) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(12)fori, kin zip (range (0, self.len-2 ∗ self.z), range (0, self.len - 2 ∗ self.z)):(13)        self.chrom3(i) = self.bound (0, 2 ∗ self.z + i) + (self.bound (1, 2 ∗ self.z + i) − self.bound (0, 2 ∗ self.z + i)) ∗ rnd (k) # xij=xmin, j + rand (0, 1) (xmax, j − xmin, j)(14)foriin range(0, self.len-2 ∗ self.z):(15)     self.X (2∗ self.z + i) = self.chrom3(i)#self.X (2∗self.z + 1) = self.chrom3(2∗self.z + 1)#self.X (4∗self.z + 2) = self.chrom(4∗self.z + 2)(16)     super().some_necessary_calculations()(17)r = 0Table 2
Initialization phase algorithm structure adjustment of another part of the code.
Generateαi, αo#αi(1)  rnd = np.random.random (size = self.z//2 + 1) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(2)  rand1 = sorted (rnd)#size = 8(3)foriin range (0,self.z//2 + 1):(4)   self.rand1(i) = rand1(i)(5)foriin range(self.z//2 + 1,self.z):(6)    self.rand1(i) = self.rand1(self.z − i) #size = 16(7)fori, k in zip (range(0,self.z//2 + 1), range(0, self.z//2 + 1)):(8)    self.chrom1(i) = self.bound (0, i) + (self.bound (1, i) − self.bound (0, i)) ∗ self.rand1(k) # Xij = Xmin, j + rand(9)foriin range (self.z//2 + 1,self.z):(10)    self.chrom1(i) = self.chrom1 (self.z − i)#αo(11)  rnd = np.random.random (size = self.z//2 + 1)# example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(12)  rand2 = sorted (rnd,reverse =True) #size = 8(13)foriin range(0,self.z//2 + 1):(14)   self.rand2(i) = rand 2(i)(15)foriin range (self.z//2 + 1,self.z):(16)   self.rand2(i) = self.rand2 (self.z − i) #size = 16(17)fori, k in zip (range (0, self.z//2 + 1), range (0, self.z//2 + 1)):(18)   self.chrom2(i) = self.boun ((0, i + self.z) + (self.bound (1, i + self.z) − self.bound (0, i + self.z)) ∗ self.rand2(k) # Xij = Xmin, j + rand (0, 1) (Xmax,j−Xmin, j) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(19)foriin range (self.z//2 + 1, self.z):(20)   self.chrom2(i) = self.chrom2 (self.z − i)


## 3.1.1. Improvements in Search Strategy

The standard artificial bee colony algorithm uses the same search formula for leading the bee phase and following the bee phase:(32)xm,i′=xm,i+rand−1,1∗xm,i−xk,i,where xk,i is randomly selected from the existing food source xm and k≠m, rand (−1, 1) represents a random number generated in the interval (−1, 1).The standard artificial bee colony algorithm uses the same search formula (32) in the leading bee stage and the following bee stage. Searching for the current food source xm, the neighborhood of the current solution xm, and in each cycle the global optimal solution is simply stored, but not in the search strategy, which plays the guiding role of the current global optimal solution to the real solution xbest. In order to improve the quality of the solution of the high-dimensional nonlinear problem and the global solution and convergence speed of the bee colony algorithm, the search formula (32) for the bee stage and follow-bee stages in the standard bee colony algorithm is improved to the following formula:(33)xm,i′=xbest,i+rand−1,1∗xbest,i−xm,i,(34)xm,i′=xk,i+rand−1,1∗xm,i−xk,i.After improving the standard artificial bee colony algorithm, the leader bee uses the search formula (33) to search. It is based on the global optimal solution xbest,i and plays its guiding role. Meanwhile, combining with the current solution xm,i, the leader bee carries out neighborhood search to obtain a better new solution xm,i′. Similarly, the follower bee uses search formula (34) to search. It is based on the random selection of other food source solutions xk,i and uses the current solution xm,i to search the neighborhood to obtain a better new solution xm,i′. Random selection of other food sources can enhance the global search ability of the algorithm and avoid single local search.


## 3.1.2. Algorithm Structure Improvement

The geometric parameters of high-speed angular contact ball bearings belong to centimeter level and the loaded deformation parameters belong to micrometer level, so the results of calculation model have great root differences. However, the standard artificial bee colony algorithm does not consider individual diversity and difference, so it is necessary to make structural improvements to the standard artificial bee colony algorithm.Standard artificial bee colony algorithm generates N food sources directly in the initialization stage, that is, generating N D-dimensional initial solutionsxm=xm1,xm2,xmi,…,xmD, among them m=1,2,3,…,D, i=1,2,3,…,D. In the initialization stage, the cross-scale problem and priority relationship in each possible solution xm1,xm2,xmi,…,xmD are not considered step by step for the problem with large root differences. As a result, the algorithm converges to unreliable solutions prematurely and even falls into local dead cycles, which seriously affects the global convergence rate.The author uses the idea of segmentation to initialize the parameters of the food sourcexm: first generate parameters δa and δr and then generate parameters αij and αoj, which solves all the differences and cross-scale problems of the required parameters, so that the algorithm is reliable at the source and ensures the accuracy of the calculation results.Some codes after adjusting the structure of the algorithm are shown in Table1, which represents the parameters δa and δr; Table 2 represents the parameters αij and αoj to be determined.Table 1
Initialization phase algorithm structure adjustment part of the code.
Generateδa, δr(1)def calculate (self)(2) rnd = np.random.random (size = self.len-2∗self.z) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(3)fori, kin zip (range(0, self.len-2∗self.z), range (0, self.len-2∗ self.z)):(4)   self.chrom 3(i) = self. bound (0, 2∗self.z + i) + (self.bound (1, 2∗ self. z + i) − self.bound (0, 2 ∗ self. z + i))∗ rnd(k)(5)foriin range (0, self.len-2∗ self.z):(6)   self.X (2∗ self.z + i) = self.chrom3(i) # δa, δr#self.X (2∗self.z + 1) = self.chrom3(2∗self.z + 1)#self.X (4∗self.z + 2) = self.chrom(4∗self.z + 2)(7)  super().some_necessary_calculations()(8)r = 0(9)while r < self.z:(10)if self.A1 (r) < 0 or self.A2 (r) < 0:(11)       rnd = np.random.random (size = self.len-2∗ self.z) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(12)fori, kin zip (range (0, self.len-2 ∗ self.z), range (0, self.len - 2 ∗ self.z)):(13)        self.chrom3(i) = self.bound (0, 2 ∗ self.z + i) + (self.bound (1, 2 ∗ self.z + i) − self.bound (0, 2 ∗ self.z + i)) ∗ rnd (k) # xij=xmin, j + rand (0, 1) (xmax, j − xmin, j)(14)foriin range(0, self.len-2 ∗ self.z):(15)     self.X (2∗ self.z + i) = self.chrom3(i)#self.X (2∗self.z + 1) = self.chrom3(2∗self.z + 1)#self.X (4∗self.z + 2) = self.chrom(4∗self.z + 2)(16)     super().some_necessary_calculations()(17)r = 0Table 2
Initialization phase algorithm structure adjustment of another part of the code.
Generateαi, αo#αi(1)  rnd = np.random.random (size = self.z//2 + 1) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(2)  rand1 = sorted (rnd)#size = 8(3)foriin range (0,self.z//2 + 1):(4)   self.rand1(i) = rand1(i)(5)foriin range(self.z//2 + 1,self.z):(6)    self.rand1(i) = self.rand1(self.z − i) #size = 16(7)fori, k in zip (range(0,self.z//2 + 1), range(0, self.z//2 + 1)):(8)    self.chrom1(i) = self.bound (0, i) + (self.bound (1, i) − self.bound (0, i)) ∗ self.rand1(k) # Xij = Xmin, j + rand(9)foriin range (self.z//2 + 1,self.z):(10)    self.chrom1(i) = self.chrom1 (self.z − i)#αo(11)  rnd = np.random.random (size = self.z//2 + 1)# example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(12)  rand2 = sorted (rnd,reverse =True) #size = 8(13)foriin range(0,self.z//2 + 1):(14)   self.rand2(i) = rand 2(i)(15)foriin range (self.z//2 + 1,self.z):(16)   self.rand2(i) = self.rand2 (self.z − i) #size = 16(17)fori, k in zip (range (0, self.z//2 + 1), range (0, self.z//2 + 1)):(18)   self.chrom2(i) = self.boun ((0, i + self.z) + (self.bound (1, i + self.z) − self.bound (0, i + self.z)) ∗ self.rand2(k) # Xij = Xmin, j + rand (0, 1) (Xmax,j−Xmin, j) # example: if self.len = 5, rnd = (0.29049083 0.96149427 0.75)(19)foriin range (self.z//2 + 1, self.z):(20)   self.chrom2(i) = self.chrom2 (self.z − i)


## 3.2. Testing and Analysis of Improved Algorithms

In order to verify the reliability and feasibility of the improved algorithm, four standard test functions are used to test the performance of the improved artificial bee colony algorithm. Among the four standard test functions,f1 is the Rastrigin function, f2 is the Ackley function, f3 is the Rosenbrock function, and f4 is the selfdef function, and its analytical expression is from (35) to (38).Rastrigin functionf1 is(35)f1x=10n+∑i=1nxi2−10cos2πxi.Ackley functionf2 is(36)f2x=−20exp0.21n∑i=1nxi2−exp1n∑i=1ncos2πxi+e+20.Rosenbrock functionf3 is(37)f3x=∑i=1n100xi2−xi+12+xi−12.Selfdef functionf4 is(38)f4x,y=x−99992+y2−52+x−9998y−3x1000+12.In Figure9, (a) is a three-dimensional diagram of a Rastrigin test function, (b) is a three-dimensional diagram of an Ackley test function, (c) is a three-dimensional diagram of a Rosenbrock test function, and (d) is a three-dimensional diagram of a Selfdef test function. f1 is a multimodal function. In the domain xi∈−5.12,5.12, this function has very many minima, and there is only one minimum, and the global minimum is minf1=f10,0,…,0=0. The f2 test function is used to verify the global convergence speed of the algorithm. The global minimum in the domain xi∈−5,5 is minf2=f20,0,…,0=0. f3 is a complex morbid function. Because it is difficult to find its global minimum, it is necessary to use this function to test the improved artificial bee colony algorithm. In the domain of xi∈−10,10, the global minimum is minf3=f30,0,…,0=0. f4 is used to test a function with a large root difference. This function is in the domain interval x∈−15000,15000, and the minimum value is minf4=f410000,2=0. From Figure 9, it can intuitively observe the extreme value of each function.Figure 9
Three-dimensional visualization of standard test functions.
(a)(b)(c)(d)Table3 shows the maximum number of iterations, test dimensions, algorithm control parameters, theoretical global optimal solutions, and search range for the four test functions where f1, f2, and f3 are used to test whether the algorithm has reliable global convergence speed and high accuracy in the face of high-dimensional, nonlinear, and multiextreme problems; f4 is used to test whether the improved algorithm has reliable search ability for large root difference problems.Table 3
Test functions and algorithm parameter settings.
FunctionMaximum number of iterationsFunctional dimensionAlgorithmic control parametersPopulation sizeTheoretical global optimal solutionSearch rangef1500050100300[−5.12, 5.12]f2500050100300[−5, 5]f3500050100300[−10, 10]f450002100300[−15000, 15000]The test results show that the standard artificial bee colony algorithm and the improved artificial bee colony algorithm are used to calculate the results of 10 tests for each test function, as shown in Table4. Meanwhile, the comparison curves of convergence rate of each test function are obtained, as shown in Figure 10.Table 4
Comparison of standard artificial bee colony algorithm and improved artificial bee colony algorithm.
Test functionAlgorithmOptimal solutionThe worst solutionMeanStandard deviationCalculating time (s)RastriginStandard1.95284E − 089.62104E − 063.24631E − 064.02246E − 0627.94Improved3.64366E − 102.30741E − 094.78643E − 107.56349E − 1023.43AckleyStandard4.14784E − 119.97926E − 112.38264E − 112.23417E − 1133.65Improved6.43285E − 137.15932E − 128.13264E − 132.78687E − 1228.26RosenbrockStandard1.41868E − 025.23661E − 014.23692E − 011.97058E − 0196.56Improved5.38452E − 041.35267E − 036.43285E − 043.17665E − 0473.24SelfdefStandard8.38282E − 072.05764E − 056.43657E − 067.05824E − 0623.68Improved1.57346E − 088.43697E − 083.85145E − 082.46283E − 0817.35Figure 10
Convergence rate contrast curve. (a) Rastrigin function. (b) Ackley function. (c) Rosenbrock function. (d) Selfdef function.
(a)(b)(c)(d)From Table4 and Figure 10, it can be seen that the improved algorithm improves in the order of magnitude compared with the standard algorithm in terms of the optimal solution accuracy, the worst solution accuracy, the mean value accuracy, and the variance accuracy. Meanwhile, the improved artificial bee colony algorithm reduces significantly the iteration calculation time compared with the standard artificial bee colony algorithm. These syntheses show that the improved artificial bee colony algorithm can solve the high-speed angular contact ball bearing calculation model with high dimension, nonlinearity, multiextremum, and large root difference.


## 3.3. Solution of Simplified Calculation Model

In order to adopt the improved artificial bee colony algorithm, the high-speed angular contact ball bearing calculation model needs to be simplified to the optimization calculation model. First, (7), (8), (13), (14), (16), and (17) are assembled into the following function groups:(39)F1j=A1j−X1j2+A2j−X2j2−fi−0.5D+δij2,F2j=X1j2+X2j2−fo−0.5D+δoj2,F3j=Kijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij+λojMgjX2j/D−Kojδoj1.5X1jfo−0.5D+δoj,F4j=Kijδij1.5A2j−X2j−λijMgj/DA1j−X1jfi−0.5D+δij−λojMgjX1j/D+Kojδoj1.5X2jfo−0.5D+δoj+Fcj,F5j=Fa−∑j=1j=zKijδij1.5A1j−X1j−λijMgj/DA2j−X2jfi−0.5D+δij.F6j=Fr−∑j=1j=zKijδij1.5A2j−X2j−λijMgj/DA1j−X1jfi−0.5D+δijcosψj.Substituting (32)–(35) into function group (36), then transforming functions F1j, F2j, F3j, F4j, F5j, and F6j into objective functions, and finally using improved artificial bee colony algorithm to find Minimum value, the specific optimization model is as follows:Objective function:(40)F=∑j=1ZF1j2+∑j=1ZF2j2+∑j=1ZF3j2+∑j=1ZF4j2+F52+F62.The design variables are(41)x=αi1,αi2,αi3,…,αij,αo1,αo2,αo3,…,αoj,δa,δr.The basic constraints are(42)0≤αij≤α0,α0≤αij≤αrcsinA1j−X1jfi−0.5D+δij.Findxto minimize the Fx value; that is, find minFx.The design variables are 2Z + 2, which are nearly half less than the solution variables of the traditional numerical algorithm and are conducive for improving the global convergence rate. The specific flow of the high-speed angular contact ball bearing calculation model using the improved artificial bee colony algorithm is shown in Figure 11.Figure 11
Solution flow chart.


## 4. Results and Discussion

Taking 218 high-speed angular contact ball bearings as an example, the author analyzes and calculates its dynamic characteristics. The specific geometric structure parameters and material parameters are shown in Table5, and the working condition parameters are given according to the specific analysis.Table 5
Structural parameters and material parameters of 218 angular contact ball bearings.
Bearing structural and material parametersNumerical valueInner raceway diameter (mm)102.79Outer raceway diameter (mm)147.73Rolling elements diameter (mm)22.23Number of rolling elements Z16The inner raceway curvature radius coefficientfi0.5232The outer raceway curvature radius coefficientfo0.5232Initial contact angelα0 (º)40Rolling elements and inner and outer ring elastic modulus (GPa)206Rolling elements and inner and outer ring Poisson’s ratioμ0.3Rolling elements and inner and outer ring densityρ (kg/m3)7850


### 4.1. Calculation Results Accuracy Verification

In order to verify the accuracy of the deformation calculation, the bearing stiffness test bench, as shown in Figure12, was used to measure the axial deformation of the high-speed angular contact bearing, under the action of both axial and radial forces.Figure 12
Bearing dynamic stiffness tester.Figure13 shows the axial deformation of the inner ferrule under different axial forces. Figure 14 shows the axial deformation of the inner ferrule when the bearing is subjected to both axial and radial forces. Meanwhile, Figures 13 and 14, respectively, compare the experimental results with the numerical results of the improved algorithm. It can be shown from Figure 13 that the axial deformation increases with the increase of axial force and velocity, respectively. When the axial load is low and the rotation speed is high, the axial deformation will be negative due to the centrifugal force of the rolling element. With the axial force increases, the axial deformation gradually changes from negative to positive. It can be seen from Figure 14 that the axial deformation also increases with the increase of the axial force and speed, respectively. This condition is consistent with the axial change trend of pure axial load, and it shows that the radial force has little effect on the axial deformation.Figure 13
Axial deformation only under axial load.Figure 14
Axial deformation under axial load and radial load.It can be seen from Figures13 and 14 that the numerical calculation results are basically consistent with the trend of the test results; there is a certain deviation between the two sets of data, which is due to the assumption of ring deformation in the numerical algorithm and the absence of consideration of bearing clearance and other factors. Furthermore, the results data of the change trend of the above axial deformation is completely consistent with the calculation results in the literature [17, 24, 25].


### 4.2. Analysis of High-Speed Characteristic Calculation Results of Angular Contact Ball Bearing

The improved artificial bee colony algorithm was used to calculate the contact force between the ball and the inner and outer channels under different speeds and axial preloads. The calculation results are shown in Figures15 and 16, respectively.Figure 15
The contact force between inner ring and rolling ball.Figure 16
The contact force between outer ring and rolling ball.It can be seen from Figure15 that the internal contact force increases linearly with the increase of axial force; at the same speed, the internal contact force increases substantially with the increase of axial force. Under the same axial force, the internal contact force decreases with the increase of rotating speed, and the decreasing range is smaller and smaller.It can be seen from Figure16 that when the rotation speed is low, the external contact force increases linearly with the increase of the axial force; thus the increase amplitude is small. That is because, at low speed, the centrifugal force is also small, and the external contact force is mainly provided by axial force. When the rotation speed is high, the external contact force decreases first and then increases with the increase of axial force. That is because, at high speed, the centrifugal force is large, and when the axial force is small, the external contact force is mainly provided by the centrifugal force. With the axial force increases, the revolution speed of the rolling element and the inner and outer rings and the rolling element changes, resulting in the decrease of the centrifugal force and the decrease of the external contact force. As the axial force continues to increase, the external contact force is mainly provided by the axial force, which in turn causes the external contact force to increase as the axial force increases.Figure17 shows the change of the actual contact angle with the ball angle position under the same axial load and the change of the radial force load. It can be seen that when the influence of rotational speed is not considered and the radial force is zero, the actual contact angle is greater than the initial contact angle and does not change with the change of the ball position. When the radial force is greater than zero, the contact angle changes with the change of the ball position, and the minimum contact angle appears at the top of the ball in the opposite direction of the maximum radial force. That is because the contact force is minimal when the ball is in this position. The maximum position of contact angle is located at the top of the maximum radial force positive direction; meanwhile, the contact load acting on the ball is the largest. In addition, when the axial force is fixed, with the increase of the radial force, the change range of the actual contact angle of the bearing increases.Figure 17
Change curve of actual contact angle with ball position under different radial loads.Figure18 shows the change of the actual contact angle with the ball angle position under the same radial load and the change of the axial force load. It can be seen that the axial force only affects the actual contact angle of the ball bearing and does not change the distribution of the actual contact angle. And the actual contact angle of a ball bearing increases with increasing axial force.Figure 18
Change curve of actual contact angle with ball position under different axial loads.Figure19 shows the influence of rotation speed on the contact force between the ball and the inner ring when the bearing subjects both axial and radial forces. It can be seen from the figure that, within the range of the ball angles I and IV zones, the internal contact force increases with the increase of the rotation speed and the increase of the IV zone becomes gentle with the increase of the rotation speed. In zone II of ball angle, the internal contact force decreases with the increase of rotation speed. In the range of zone III of ball angle, it is a mixed change zone. With the increase of rotation speed, the change of the internal contact force has no obvious regularity.Figure 19
Contact force between ball and inner ring under axial force and radial force.Figure20 shows the influence of rotation speed on contact force between ball and outer ring under the condition of bearing both axial force and radial force. It can be seen from the figure that, in the range of zone I and zone III of ball angle, the external contact force increases with the increase of rotating speed. With the increase of rotating speed, the range of change in zone III becomes smaller and smaller due to the influence of centrifugal force and gyroscopic moment. In the range of zone II of ball angle, it is a mixed change zone. With the increase of rotation speed, the change of the internal contact force has no obvious regularity.Figure 20
Contact force between ball and outer ring under axial force and radial force.Figure21 shows the influence of rotating speed on centrifugal force under the condition of bearing both axial force and radial force. As can be seen from the figure, the centrifugal force is particularly small at low speeds, and there is little difference in centrifugal force at different azimuths. As the speed increases, the centrifugal force increases sharply, and the gap between centrifugal forces at different azimuth angles also increases significantly.Figure 21
The variation of centrifugal force with different positions of rolling elements at different rotational speeds.Figure22 shows the effect of the rotational speed on the gyro torque under the condition that the bearing is subjected to both axial and radial forces. It can be seen from the figure that, under the same working conditions, the gyro torque of the rolling body increases with the increase of the rotation speed. When the speed n>5000r/min, the gyro torque increases sharply. The gyro torque will cause the gyro to rotate, and the gyro will cause friction and heat, which will affect the performance of the bearing. Moreover, the variation range of the gyro torque of the rolling element increases with the increase of the rotating speed, and the gyro torque in the bearing area is larger than that in the nonload bearing area.Figure 22
The variation of gyroscopic moment with different positions of rolling elements at different speeds.Meanwhile, the change trend results data of the above the actual contact angle, contact force, centrifugal force, and gyroscopic moment were completely consistent with the calculation result in the literatures [14, 24–27]. This verifies the correctness of the simplified model and the improved numerical algorithm to some extent.


## 4.1. Calculation Results Accuracy Verification

In order to verify the accuracy of the deformation calculation, the bearing stiffness test bench, as shown in Figure12, was used to measure the axial deformation of the high-speed angular contact bearing, under the action of both axial and radial forces.Figure 12
Bearing dynamic stiffness tester.Figure13 shows the axial deformation of the inner ferrule under different axial forces. Figure 14 shows the axial deformation of the inner ferrule when the bearing is subjected to both axial and radial forces. Meanwhile, Figures 13 and 14, respectively, compare the experimental results with the numerical results of the improved algorithm. It can be shown from Figure 13 that the axial deformation increases with the increase of axial force and velocity, respectively. When the axial load is low and the rotation speed is high, the axial deformation will be negative due to the centrifugal force of the rolling element. With the axial force increases, the axial deformation gradually changes from negative to positive. It can be seen from Figure 14 that the axial deformation also increases with the increase of the axial force and speed, respectively. This condition is consistent with the axial change trend of pure axial load, and it shows that the radial force has little effect on the axial deformation.Figure 13
Axial deformation only under axial load.Figure 14
Axial deformation under axial load and radial load.It can be seen from Figures13 and 14 that the numerical calculation results are basically consistent with the trend of the test results; there is a certain deviation between the two sets of data, which is due to the assumption of ring deformation in the numerical algorithm and the absence of consideration of bearing clearance and other factors. Furthermore, the results data of the change trend of the above axial deformation is completely consistent with the calculation results in the literature [17, 24, 25].


## 4.2. Analysis of High-Speed Characteristic Calculation Results of Angular Contact Ball Bearing

The improved artificial bee colony algorithm was used to calculate the contact force between the ball and the inner and outer channels under different speeds and axial preloads. The calculation results are shown in Figures15 and 16, respectively.Figure 15
The contact force between inner ring and rolling ball.Figure 16
The contact force between outer ring and rolling ball.It can be seen from Figure15 that the internal contact force increases linearly with the increase of axial force; at the same speed, the internal contact force increases substantially with the increase of axial force. Under the same axial force, the internal contact force decreases with the increase of rotating speed, and the decreasing range is smaller and smaller.It can be seen from Figure16 that when the rotation speed is low, the external contact force increases linearly with the increase of the axial force; thus the increase amplitude is small. That is because, at low speed, the centrifugal force is also small, and the external contact force is mainly provided by axial force. When the rotation speed is high, the external contact force decreases first and then increases with the increase of axial force. That is because, at high speed, the centrifugal force is large, and when the axial force is small, the external contact force is mainly provided by the centrifugal force. With the axial force increases, the revolution speed of the rolling element and the inner and outer rings and the rolling element changes, resulting in the decrease of the centrifugal force and the decrease of the external contact force. As the axial force continues to increase, the external contact force is mainly provided by the axial force, which in turn causes the external contact force to increase as the axial force increases.Figure17 shows the change of the actual contact angle with the ball angle position under the same axial load and the change of the radial force load. It can be seen that when the influence of rotational speed is not considered and the radial force is zero, the actual contact angle is greater than the initial contact angle and does not change with the change of the ball position. When the radial force is greater than zero, the contact angle changes with the change of the ball position, and the minimum contact angle appears at the top of the ball in the opposite direction of the maximum radial force. That is because the contact force is minimal when the ball is in this position. The maximum position of contact angle is located at the top of the maximum radial force positive direction; meanwhile, the contact load acting on the ball is the largest. In addition, when the axial force is fixed, with the increase of the radial force, the change range of the actual contact angle of the bearing increases.Figure 17
Change curve of actual contact angle with ball position under different radial loads.Figure18 shows the change of the actual contact angle with the ball angle position under the same radial load and the change of the axial force load. It can be seen that the axial force only affects the actual contact angle of the ball bearing and does not change the distribution of the actual contact angle. And the actual contact angle of a ball bearing increases with increasing axial force.Figure 18
Change curve of actual contact angle with ball position under different axial loads.Figure19 shows the influence of rotation speed on the contact force between the ball and the inner ring when the bearing subjects both axial and radial forces. It can be seen from the figure that, within the range of the ball angles I and IV zones, the internal contact force increases with the increase of the rotation speed and the increase of the IV zone becomes gentle with the increase of the rotation speed. In zone II of ball angle, the internal contact force decreases with the increase of rotation speed. In the range of zone III of ball angle, it is a mixed change zone. With the increase of rotation speed, the change of the internal contact force has no obvious regularity.Figure 19
Contact force between ball and inner ring under axial force and radial force.Figure20 shows the influence of rotation speed on contact force between ball and outer ring under the condition of bearing both axial force and radial force. It can be seen from the figure that, in the range of zone I and zone III of ball angle, the external contact force increases with the increase of rotating speed. With the increase of rotating speed, the range of change in zone III becomes smaller and smaller due to the influence of centrifugal force and gyroscopic moment. In the range of zone II of ball angle, it is a mixed change zone. With the increase of rotation speed, the change of the internal contact force has no obvious regularity.Figure 20
Contact force between ball and outer ring under axial force and radial force.Figure21 shows the influence of rotating speed on centrifugal force under the condition of bearing both axial force and radial force. As can be seen from the figure, the centrifugal force is particularly small at low speeds, and there is little difference in centrifugal force at different azimuths. As the speed increases, the centrifugal force increases sharply, and the gap between centrifugal forces at different azimuth angles also increases significantly.Figure 21
The variation of centrifugal force with different positions of rolling elements at different rotational speeds.Figure22 shows the effect of the rotational speed on the gyro torque under the condition that the bearing is subjected to both axial and radial forces. It can be seen from the figure that, under the same working conditions, the gyro torque of the rolling body increases with the increase of the rotation speed. When the speed n>5000r/min, the gyro torque increases sharply. The gyro torque will cause the gyro to rotate, and the gyro will cause friction and heat, which will affect the performance of the bearing. Moreover, the variation range of the gyro torque of the rolling element increases with the increase of the rotating speed, and the gyro torque in the bearing area is larger than that in the nonload bearing area.Figure 22
The variation of gyroscopic moment with different positions of rolling elements at different speeds.Meanwhile, the change trend results data of the above the actual contact angle, contact force, centrifugal force, and gyroscopic moment were completely consistent with the calculation result in the literatures [14, 24–27]. This verifies the correctness of the simplified model and the improved numerical algorithm to some extent.


## 5. Conclusion

According to the dynamic changes of various parameters of the angular contact ball bearing during high-speed operation, the equations for solving the high-speed dynamic characteristics are systematically analyzed.(1)
The transfer relationship of variables in the high-speed angular contact bearing model is explored, and a set of nonlinear equations with actual internal (inner) contact angle, actual external (outer) contact angle, radial deformation, and axial deformation as basic iterative variables are obtained. The number of nonlinear equations to be solved is greatly reduced.(2)
An artificial bee colony algorithm is proposed to solve the big root difference through algorithm improvement. The effectiveness of the improved algorithm is verified by the standard test function. The improved algorithm is applied to solve the dynamic characteristics of high-speed angular contact ball bearings and is verified by a practical project.(3)
The model is used to solve the influence of the rotational speed on the actual contact angle, contact force, centrifugal force, gyroscopic moment, and attitude angle, which provides a reference for the use, research, and design of the high-speed angular contact bearing.



---

*Source: 8843524-2020-12-08.xml*