Analytical Determination and Inﬂuence Analysis of Stiffness Matrix of Ball Bearing under Different Load Conditions

: Bearing stiffness, as one of the most important service characteristics for ball bearing, plays a crucial role in the bearing design and rotor dynamic analysis. To rapidly and accurately calculate the stiffness matrix of ball bearing under the arbitrary load conditions, a 5-DOF analytical model for bearing stiffness matrix analysis has been established by the ball–raceway contact analysis, implicit/explicit differential method, and matrix operations. The model has been validated comparing with the previous methods and experimental results. Based on this, the model has been used to investigate the inﬂuences of the load and operation conditions, the structural parameters variation on stiffness of ball bearing. The results show that property increasing axial preload can inhibit the attenuation of speed-varying stiffness, and the contact states between balls and raceways also have signiﬁcant inﬂuence on the change in the stiffness of ball bearings. Besides, a larger curvature coefﬁcient of inner raceway and a small curvature coefﬁcient of outer raceway can effectively improve the stiffness of ball bearing at high speed. Therefore, the proposed method can be a useful tool in bearing optimize design and performance analysis of ball and rotor system under various load conditions.


Introduction
As the key supporting parts in high-speed rotor system, rolling bearings are widely used across numerous rotating mechanical systems, including aerospace, automobile, machine tool, etc. With the rapid development of the modern industry, there is demand that the service performance of ball bearing must be improved to meet the higher operation requirement of mechanical system. To guarantee the optimum performance, the design and use parameters of rolling bearing should be adjusted by the designers and users according to the practical operation conditions, and an accurate and efficient mathematical model of bearing can be used to reduce the cost and cycle.
During the past decades, a large number of algorithms and methods have been presented in bearing modeling and analysis. In the earlier studies [1][2][3][4], the inertia forces of balls have been neglected, which may cause a considerable error in bearing analysis and characteristic prediction, especially for the high-speed operation condition. To improve the calculation accuracy of bearing analysis, Jones [5] proposed a classical mathematic model with the "Raceway Control Hypothesis" for ball bearing subjected to the combined radial, axial, and moment loads. In this model, the effect of the inertia forces of the balls are considered. Without the Raceway Control Hypothesis, Harris [6] developed a comprehensive model based on the lubrication analysis of ball-raceway contact area, and this model can be used in the bearing skidding behavior prediction and friction power loss calculation. In addition, the experimental results show that the partial results obtained from the "Raceway Control Hypothesis (RCH)" were not precise enough, thus the "Raceway Control Hypothesis" may be only suitable for ball bearing under some specific conditions. To overcome the above shortcomings, the more accurate mechanics and kinematics of local ball have been presented without the RCH in the following research [7][8][9]. In addition, a serial of dynamic models [10][11][12][13][14][15] of ball bearing have been presented in recently years to obtain and investigate the instantaneous motion and mechanics states of ball bearing under running condition, and the most representative study was proposed by Gupta [11,12]. In his model, the relative action relations among the bearing internal components and lubricant influences were fully taken in account. However, the model is time-consuming.
Bearing stiffness, as an important characteristic of ball bearing, has significant effects on the vibration, the carrying capacity, and the rotation accuracy of bearing and rotor system. An accuracy mathematical model can not only be used to guide the design of ball bearing [16] but can also be used in the dynamic properties prediction of rotor-bearing system [17][18][19]. Therefore, a great deal of models for bearing stiffness calculation have been proposed. In the earlier studies, researchers mainly focused on the radial and axial stiffness of rolling bearing. Gargiulo [20] and Wardle [21] separately gave a set of empirical formulas to predict the radial and axial stiffness of rolling bearing, and these formulas are still used in the estimation of bearing stiffness due to its simple and ease of use. To obtain the complete expressions of bearing stiffness matrix, a mathematic model of ball bearing has been presented by Houpert [22] and Hernot [23], but the dynamic forces of balls were ignored. Liu [24] gave a convenient method, called the Newton-Raphson method, for the analytical calculation of the stiffness of ball bearing stiffness under combined load but did not consider the influence of speed. However, the inertia forces of balls play a key role in the stiffness for ball bearing at the high-speed condition. In addition, the representative studies conducted by Lim and Singh [25][26][27] showed that only considering the axial and radial stiffness is not sufficient to explain the flexural and out-of-plane motions of ball bearing, so a complete stiffness matrix solution, including the tilting stiffness and off-diagonal cross-coupling terms, is necessary to describe the influences of ball bearing on vibration transmission. In order to quickly obtain bearing stiffness matrix, the finite difference method has been used to calculate the numerical differentiation of stiffness matrix firstly [28,29]. On this basis, Yang [30,31] calculated and analyzed the stiffness of the single bearing and complex rotor system. However, the calculation accuracy of the above method depends on the step and is difficult to guarantee. Therefore, the full-analytical and semi-analytical model of bearing stiffness solution have been proposed and have been widely used in practical engineering. Noel [32] gave a comprehensive algorithm for the analytical method of bearing stiffness matrix, and the inertia forces of balls have been fully considered to lead an exact solution. Fang [33] gave the complete analytical expression of the bearing stiffness matrix based on an improved model and classified the stiffness elements in the bearing stiffness matrix. However, there is no comparison and verification. Xia [34] introduced the implicit differential theory for the bearing stiffness calculation, but the detailed derivation process was not given. Gunduz and Singh [35] calculated and studied the stiffness matrix of the double-row ball bearings with different arrangements, but this model has neglected the dynamic effect of ball. Tsuha N [36] calculated the stiffness of the cylindrical roller bearing more accurately through an innovative EHL force model, but the model could not be directly applied to ball bearings.
Besides, based on the FE method, Guo and Parker [37] calculated and analyzed the stiffness of rolling bearing with different structural parameters, respectively, and the results have a good consistency with the analytical model and experiments. Based on the above stiffness model, some research [38][39][40] extended to study the time-varying stiffness of ball bearing caused by orbital motion and defects.
Besides the theoretical research, several experiments for bearing stiffness measuring were also conducted. Karus et al. [41] built a transmission test rig to measure the stiffness and damping of ball bearing under different preload. Walford and Stone [42] measured the radial and axial stiffness for ball bearings in a pair under oscillating condition on a test rig. Matti Rantatalo [43] and J.Jedrzejewski [44] also conducted a series of experiments on high spindle-bearing system, and the obvious stiffness attenuation with the rotating speed has been presented in both experimental results.
Above all, though the solution method of bearing stiffness matrix was widely studied, there still exist several problems that the influences factors on bearing stiffness have not fully considered [18,[20][21][22][23][24]32,34] and the stiffness calculation has numerical stability problems [28][29][30][31]. In addition, considering the various operation conditions of ball bearings, most of the previous models of bearing stiffness matrix calculation are only suitable for ball bearing under normal load conditions. Therefore, it is still needed to present a comprehensive analytical model for stiffness calculation of ball bearing under arbitrary load conditions. Besides, in order to guide bearing design, the details of effect of the speeding speed, structure parameters on the stiffness, and stiffness variations of the ball bearing under different conditions need to be studied.
This paper presents a quasi-static model without Raceway Control Hypothesis of ball bearing under arbitrary operation conditions, and a comprehensive analytical model for the 5-DOF stiffness matrix calculation of ball bearing is proposed. Based on the contact analysis of balls and raceways for ball bearing at operation conditions, the proposed algorithm can be applied to calculate and analyze the stiffness matrix for ball bearing under the arbitrary load conditions and rotating speeds. Besides, on this basis of proposed model, the influences of the external loads (including axial preload), rotating speeds, and structural parameters on bearing stiffness are analyzed in detail. The research of this method and the discussion of these influence laws provide a useful tool for the optimal design of bearings and the performance analysis of ball and rotor systems under various load conditions.

Quasi-Static Model of Ball Bearing
It can be seen from Figure 1 that ball bearing is usually composed of a fixed outer ring, a rotating inner ring, a cage and several balls, and the balls are uniformly distributed along the circle by the cage pockets. The angular position of kth ball can be written as: It is assumed that the bearing is acted by a combined external load vector F = {Fx, Fy, Fz, My, Mz}, and the relative displacement vector between the inner and outer rings is d = { , , , , }. Assuming that the rotating speed of inner ring is . It can be seen from Figure 2a that the position of ball center moves upward due to the ball centrifugal force. Consequently, the contact angle between the inner raceway and ball increases, while the contact angle between the outer raceway and ball decreases. It should be noticed that the It is assumed that the bearing is acted by a combined external load vector F = {F x , F y , F z , M y , M z }, and the relative displacement vector between the inner and outer rings is Machines 2022, 10, 238 4 of 21 d = {δ x , δ y , δ z , θ y , θ z }. Assuming that the rotating speed of inner ring is ω i . It can be seen from Figure 2a that the position of ball center moves upward due to the ball centrifugal force. Consequently, the contact angle between the inner raceway and ball increases, while the contact angle between the outer raceway and ball decreases. It should be noticed that the above analysis is only applicable to the ball-raceway contact situation, but when ball bearing is acted by the combined loads with a small axial force, parts of balls may no longer be loaded, and the bearing is divided into two parts: the loaded zone and the unloaded zone [45]. It can be seen from Figure 2b that parts of the balls are separate from the inner raceway by the action of centrifugal force. However, this phenomenon is usually ignored in modeling and stiffness calculation of ball bearing under the action of combined loads. It is assumed that the bearing is acted by a combined external load vector F = {Fx, Fy, Fz, My, Mz}, and the relative displacement vector between the inner and outer rings is d = { , , , , }. Assuming that the rotating speed of inner ring is . It can be seen from Figure 2a that the position of ball center moves upward due to the ball centrifugal force. Consequently, the contact angle between the inner raceway and ball increases, while the contact angle between the outer raceway and ball decreases. It should be noticed that the above analysis is only applicable to the ball-raceway contact situation, but when ball bearing is acted by the combined loads with a small axial force, parts of balls may no longer be loaded, and the bearing is divided into two parts: the loaded zone and the unloaded zone [45]. It can be seen from Figure 2b that parts of the balls are separate from the inner raceway by the action of centrifugal force. However, this phenomenon is usually ignored in modeling and stiffness calculation of ball bearing under the action of combined loads.  It also can be seen from Figure 2a that the expressions for the sine and cosine functions of the ball-raceways contact angles of ball in loaded zone are written as: where r i , r o denote the inner/outer raceway curvature radii, respectively. The curvature centers distances of the inner and outer raceways at angular position ψ k can be written: where D and d m are the dimeter of ball and the pitch radius of ball bearing, respectively. Besides, the geometry equilibrium equations can also be obtained according to Figure 2a: It can be seen from Figure 2b, when the ball-inner raceway separation occurs, that the center distance between the outer raceway curvature and ball is expressed [46]: Machines 2022, 10, 238

of 21
In addition, since the ball is not restrained by the inner ring, the following inequality about the center's distance between the inner ring curvature and ball needs to be satisfied [46]: where the distances PO d and PO i can be written as: when the ball is located in the loaded zone, the mechanical equilibrium equations of ball are expressed as: where Q ik and Q ok are the ball-inner/outer raceway contact loads, respectively. T ik and T ok are the ball-inner/outer frictions, respectively.
In the above equations, the K ik and K ok are the load-deformation coefficients in ballraceways contacts, the values of K ik and K ok are determined by the basic structural parameters and the contact angles, and the analytical calculation expressions K ik and K ok are given in Appendix A. δ ik and δ ok are the elastic deformations of ball-raceway contacts. In addition, according to the Raceway Control Theory of Jones, λ ik = 0 and λ ok = 2 are used in "Outer Raceway Control Hypothesis" and λ ik = 1, λ ok = 1 are used in "Inner Raceway Control Hypothesis". However, a lot of research show that the Inner and Outer Raceway Control Hypothesis shows a speed-discontinuity in bearing analysis and the "Raceway Control Hypothesis" is only suitable for the analysis of ball bearing under some specific conditions [6]. In order to overcome the above shortcomings, λ ik = 2 Q ik Q ik +Q ok and λ ok = 2 Q ok Q ik +Q ok can be selected according to equal friction coefficient assumption [8]. Besides, the inertia forces, including the gyroscopic moment M gk and the centrifugal force F ck of ball, are expressed: where m and J denote the mass and the mass moment of ball. ω mk and ω bk are the revolution and spin angular speeds of ball, respectively [9], ω mk and ω bk can be written: where γ ik = Dcosα ik d m and γ ok = Dcosα ok d m , and the pitch angle β k are given based on the d'Alembert principle [9]: with Machines 2022, 10, 238 where a ik and a ok the semi-major axes of the ball-raceways contacts. L ik and L ok are the second kind elliptic integrals of the ball-raceways contacts [9].
Besides, it can be seen from Figure 3b, if the ball separates from the inner raceway, the force equilibrium equation of ball is given as: where and the semi-major axes of the ball-raceways contacts. and are the second kind elliptic integrals of the ball-raceways contacts [9].
Besides, it can be seen from Figure 3b, if the ball separates from the inner raceway, the force equilibrium equation of ball is given as: Through the above analysis, it can be found that the ball in the loaded zone has 4 unknown variables = { , , , } in Equations (2) and (3). The Newton-Raphson algorithm has been used to calculate these variables for a given global displacement value vector d = { , , , , }. In addition, it can be seen from Equation (14) that only one variable is used to describe the force state of ball in the unloaded zone. However, it needs to be noted that the inequality (5) should be calculated first in the calculation of ball local equations. If the inequality (5) holds, Equation (14) is solved, or Equations (2) and (3) are solved. Through the above analysis, it can be found that the ball in the loaded zone has 4 unknown variables ε k = {X 1k , X 2k , δ ik , δ ok } in Equations (2) and (3). The Newton-Raphson algorithm has been used to calculate these variables for a given global displacement value vector d = {δ x , δ y , δ z , θ y , θ z }. In addition, it can be seen from Equation (14) that only one variable δ ok is used to describe the force state of ball in the unloaded zone. However, it needs to be noted that the inequality (5) should be calculated first in the calculation of ball local equations. If the inequality (5) holds, Equation (14) is solved, or Equations (2) and (3) are solved.
Based on the above analysis of local balls, the equilibrium relations between the external load vector F and the displacement vector d can be determined, and the mechanical equilibrium equations for rotating inner ring are given: where G[•] is the Heavisde function, it can be expressed as: In order to calculate the global displacement for a given external load, the nested-loop Newton-Raphson iteration algorithm is used as shown in Figure 4. Firstly, input the initial iteration values of inner displacement vector d and local variables ε k = {X 1k , X 2k , δ ik , δ ok }; Secondly, solve the local equilibrium equations update the values of local variables and ballraceway contact condition; Then, substitute new local variables into the global equations of inner ring and update the values of global displacements; Repeat the above iteration until the calculation errors are less than the set threshold. Finally, output the variable values and stiffness matrix. In order to complete the above iteration process, the expressions of J m in inner iteration and stiffness matrix K in outer loop iteration need to be determined, and the matrix J m denotes the Jacobian matrix of the ball local equations to ball local variables ε k = {X 1k , X 2k , δ ik , δ ok }. The detailed expressions are given as follows: Machines 2022, 10, x FOR PEER REVIEW 8 of 22

Computation of the Stiffness Matrix
The stiffness matrix of ball bearing can be calculated as Equation (19) through the partial derivatives between the load vector F to the displacement vector d: To obtain the complete expressions of the stiffness matrix, the relations between the external load vector F and global displacement vector d should be determined, the variable transfer graph is given in Ref [45]. Previous studies [45][46][47][48] show that the calculation process of the stiffness matrix should be divided into the partial derivatives of the load In addition, the stiffness matrix K of ball bearing is illustrated in the subsequent chapters.

Computation of the Stiffness Matrix
The stiffness matrix of ball bearing can be calculated as Equation (19) through the partial derivatives between the load vector F to the displacement vector d: To obtain the complete expressions of the stiffness matrix, the relations between the external load vector F and global displacement vector d should be determined, the variable transfer graph is given in Ref [45]. Previous studies [45][46][47][48] show that the calculation process of the stiffness matrix should be divided into the partial derivatives of the load vector F to the intermediated variables υ k and the intermediated variables υ k to global displacement vector d. The difference is that the selection of intermediate variables υ k is not the same, and the most accomplished methods [18,23] have chosen υ k = {X 1k , X 2k , δ ik , δ ok , M gk , F ck } as the intermediate variables. However, this selection strategy is very complex in the calculation of the partial derivatives between υ k and d, thus different simplification schemes have been adopted that may cause a certain error in the stiffness matrix calculation. To simplify the partial derivatives of υ k with respect to d, the intermediate variables υ k = {X 1k , X 2k , δ ik , δ ok , A 1k , A 2k } have been used in the proposed method. Therefore, the stiffness matrix is written as: where the partial derivatives between the loads vector F to the intermediated variables υ k can be written as: In the previous literatures [18,19], in order to simplify the calculation process, the contributions of the load-deformation coefficients K ik and K ok to bearing stiffness calculation have been ignored. However, due to the apparent change of the contact angles, it may lead to a larger error when ball bearing operates at the high-speed range. Therefore, the detailed calculation expression of ∂K ik ∂υ k and ∂K ok ∂υ k is given in Appendix A. In addition, the partial derivatives between the intermediated variables υ k = {X 1k , X 2k , δ ik , δ ok , A 1k , A 2k } to global displacement vector d is written as: where the partial derivatives ∂A 1k ∂d and ∂A 2k ∂d terms are calculated by the explicit differential of Equations (3) and (4). However, the relations among the ε k , A 1k , and A 2k are determined by Equations (2) and (3) in implicit form, so ∂A 1k ∂d and ∂A 2k ∂d terms can be obtained by the implicit differential method. Therefore, the detailed calculation process is given as follows: where the Jacobi determinant J can be written as: Above all, the detailed solving process and key steps of analytical calculation of bearing stiffness matrix are given, and all the influence factors on stiffness for bearing quasistatic analysis have been considered. In addition, since the above solution is implemented by matrix operation, the calculation process is more clear and easier to program.

Results
On this basis, the stiffness matrix for ball bearing under arbitrary speeds and loads conditions can be calculated. In view of the actual needs for the engineering and the project, only the radial and axial stiffness are chosen as the main research objects in this paper.

Comparison against Experimental and Theoretical Results
To show the effectiveness of the presented model, the calculation results of the highprecision ball bearings 7008C were compared to the results obtained by the most representative methods [18] at first. Then the static stiffness results for high-precision ball bearings 7008C and 7014C were compared against the experimental values in bearing user manual. At last, the dynamic stiffness results of bearing EEB3-2Z at 1000 rpm were compared against the published experimental and FEM results. As shown in Table 1, the key structural parameters of the above ball bearings are given. It can be seen from Figure 5 that the contrast results of the bearing stiffness vary with the speed calculated by the proposed model and the Cao-Jones model given. It can be found that the stiffness values obtained by the proposed model are located between the stiffness values calculated by the Cao-Jones model with the "Inner/Outer Raceway Control Hypothesis". As the "Raceway Control Theory" is based on extreme conditions, the results obtained by the proposed model are more reasonable and logical. To show the effectiveness of the presented model, the calculation results of the highprecision ball bearings 7008C were compared to the results obtained by the most representative methods [18] at first. Then the static stiffness results for high-precision ball bearings 7008C and 7014C were compared against the experimental values in bearing user manual. At last, the dynamic stiffness results of bearing EEB3-2Z at 1000 rpm were compared against the published experimental and FEM results. As shown in Table 1, the key structural parameters of the above ball bearings are given.
It can be seen from Figure 5 that the contrast results of the bearing stiffness vary with the speed calculated by the proposed model and the Cao-Jones model given. It can be found that the stiffness values obtained by the proposed model are located between the stiffness values calculated by the Cao-Jones model with the "Inner/Outer Raceway Control Hypothesis". As the "Raceway Control Theory" is based on extreme conditions, the results obtained by the proposed model are more reasonable and logical.  Then, as shown in Table 2, the comparison results of static stiffness between the proposed model and user manual values are given. One can find that the values of stiffness obtained by the proposed model are slightly more than the experimental values in bearing user manual, and it may be caused by the improvement of material properties as the special technology and heat-treatment are not considered in the model calculation.  Then, as shown in Table 2, the comparison results of static stiffness between the proposed model and user manual values are given. One can find that the values of stiffness obtained by the proposed model are slightly more than the experimental values in bearing user manual, and it may be caused by the improvement of material properties as the special technology and heat-treatment are not considered in the model calculation.  Figure 6 gives the comparison results of the dynamic stiffness for bearing EEB3-2Z obtained by different theoretical methods and experimental analysis. It can be found that the results given by the proposed method are more consistent with the experimental results.
Axial stiffness (N/µm)  Figure 6 gives the comparison results of the dynamic stiffness for bearing EEB3-2Z obtained by different theoretical methods and experimental analysis. It can be found that the results given by the proposed method are more consistent with the experimental results.

Stiffness Variation with Rotating Speed and Axial Preload
For the precision and high-speed ball bearing, a proper axial preload can be used to improve the rotating accuracy and supporting rigidity of ball bearing. According to Figure  7, when ball bearing is acted by a certain axial preload, the stiffness change with the rotating speed can be divided into three stages: decline slowly-drop rapidly-trend to gentle, and the same change rules have been found in previous theoretical analyses and experimental studies [34,43,44]. Besides, the decreasing trend of bearing stiffness with rotating speed has been eased. Figure 7 also gives the results of the stiffness change with the preload of ball at different speeds, and the drastic improvement of the stiffness due to the increase in the preload is presented.

Stiffness Variation with Different Load Conditions
As mentioned earlier, the ball-raceway contact analysis has been added in the bearing modeling to realize the stiffness calculation of ball bearing under arbitrary load. Therefore, to show the generality of the proposed model, three typical operation conditions of ball bearing SKF 6205 have been given to investigate the effects of the rotating speeds, the

Stiffness Variation with Rotating Speed and Axial Preload
For the precision and high-speed ball bearing, a proper axial preload can be used to improve the rotating accuracy and supporting rigidity of ball bearing. According to Figure 7, when ball bearing is acted by a certain axial preload, the stiffness change with the rotating speed can be divided into three stages: decline slowly-drop rapidly-trend to gentle, and the same change rules have been found in previous theoretical analyses and experimental studies [34,43,44]. Besides, the decreasing trend of bearing stiffness with rotating speed has been eased. Figure 7 also gives the results of the stiffness change with the preload of ball at different speeds, and the drastic improvement of the stiffness due to the increase in the preload is presented.  Figure 6 gives the comparison results of the dynamic stiffness for bearing EEB3-2Z obtained by different theoretical methods and experimental analysis. It can be found that the results given by the proposed method are more consistent with the experimental results.

Stiffness Variation with Rotating Speed and Axial Preload
For the precision and high-speed ball bearing, a proper axial preload can be used to improve the rotating accuracy and supporting rigidity of ball bearing. According to Figure  7, when ball bearing is acted by a certain axial preload, the stiffness change with the rotating speed can be divided into three stages: decline slowly-drop rapidly-trend to gentle, and the same change rules have been found in previous theoretical analyses and experimental studies [34,43,44]. Besides, the decreasing trend of bearing stiffness with rotating speed has been eased. Figure 7 also gives the results of the stiffness change with the preload of ball at different speeds, and the drastic improvement of the stiffness due to the increase in the preload is presented.

Stiffness Variation with Different Load Conditions
As mentioned earlier, the ball-raceway contact analysis has been added in the bearing modeling to realize the stiffness calculation of ball bearing under arbitrary load. Therefore, to show the generality of the proposed model, three typical operation conditions of ball bearing SKF 6205 have been given to investigate the effects of the rotating speeds, the

Stiffness Variation with Different Load Conditions
As mentioned earlier, the ball-raceway contact analysis has been added in the bearing modeling to realize the stiffness calculation of ball bearing under arbitrary load. Therefore, to show the generality of the proposed model, three typical operation conditions of ball bearing SKF 6205 have been given to investigate the effects of the rotating speeds, the external load conditions, and the internal contact states of balls and raceways on bearing stiffness. Figure 8 gives the stiffness change curves of ball bearing under one-direction radial load. One can find that both static and dynamic stiffness increase with the radial load. Besides, different from ball bearing acted by the axial load (Preload), the stiffness of ball bearing under a certain radial load (500 N) increases with the rotating speed. Figure 8 gives the stiffness change curves of ball bearing under one-direction radial load. One can find that both static and dynamic stiffness increase with the radial load. Besides, different from ball bearing acted by the axial load (Preload), the stiffness of ball bearing under a certain radial load (500 N) increases with the rotating speed.

The Combined Large Axial Load and Small Radial Load Condition
The above section has introduced the change law of ball bearing under the single direction load. On this basis, the stiffness change curves of bearing acted by the combined large axial load and small radial load at 2000 rpm are shown in Figure 9. According to the actual loading condition of ball bearing under operation condition, the angle displacement and and are set to zero. It can be seen from Figure 9 that both axial and radial stiffness decrease with the radial load.

The Combined Small Axial Load and Large Radial Load Condition
For contrast, the change curves of bearing stiffness acted by the combined small axial load and large radial load at 2000 rpm are shown in Figure 10. One can find that the axial stiffness reduces sharply with the radial load for a given axial load (50 N), and the radial stiffness firstly drops and then rises with the radial load. Besides, both the radial and axial stiffness increase with the axial load for ball bearing under a certain radial load (100 N). Different from the uniform increasing trend of axial stiffness with the axial load, the change curve of the radial stiffness contains two main pieces: the gradually increasing stage and the rapidly increasing stage. Compared with the change rule of bearing stiffness in Figure 9, the stiffness various in Figure 10 presents the diversity and complexity. In

The Combined Large Axial Load and Small Radial Load Condition
The above section has introduced the change law of ball bearing under the single direction load. On this basis, the stiffness change curves of bearing acted by the combined large axial load and small radial load at 2000 rpm are shown in Figure 9. According to the actual loading condition of ball bearing under operation condition, the angle displacement and θ y and θ z are set to zero. It can be seen from Figure 9 that both axial and radial stiffness decrease with the radial load. Figure 8 gives the stiffness change curves of ball bearing under one-direction radial load. One can find that both static and dynamic stiffness increase with the radial load. Besides, different from ball bearing acted by the axial load (Preload), the stiffness of ball bearing under a certain radial load (500 N) increases with the rotating speed.

The Combined Large Axial Load and Small Radial Load Condition
The above section has introduced the change law of ball bearing under the single direction load. On this basis, the stiffness change curves of bearing acted by the combined large axial load and small radial load at 2000 rpm are shown in Figure 9. According to the actual loading condition of ball bearing under operation condition, the angle displacement and and are set to zero. It can be seen from Figure 9 that both axial and radial stiffness decrease with the radial load.

The Combined Small Axial Load and Large Radial Load Condition
For contrast, the change curves of bearing stiffness acted by the combined small axial load and large radial load at 2000 rpm are shown in Figure 10. One can find that the axial stiffness reduces sharply with the radial load for a given axial load (50 N), and the radial stiffness firstly drops and then rises with the radial load. Besides, both the radial and axial stiffness increase with the axial load for ball bearing under a certain radial load (100 N). Different from the uniform increasing trend of axial stiffness with the axial load, the change curve of the radial stiffness contains two main pieces: the gradually increasing stage and the rapidly increasing stage. Compared with the change rule of bearing stiffness in Figure 9, the stiffness various in Figure 10 presents the diversity and complexity. In

The Combined Small Axial Load and Large Radial Load Condition
For contrast, the change curves of bearing stiffness acted by the combined small axial load and large radial load at 2000 rpm are shown in Figure 10. One can find that the axial stiffness reduces sharply with the radial load for a given axial load (50 N), and the radial stiffness firstly drops and then rises with the radial load. Besides, both the radial and axial stiffness increase with the axial load for ball bearing under a certain radial load (100 N). Different from the uniform increasing trend of axial stiffness with the axial load, the change curve of the radial stiffness contains two main pieces: the gradually increasing stage and the rapidly increasing stage. Compared with the change rule of bearing stiffness in Figure 9, the stiffness various in Figure 10 presents the diversity and complexity. In order to further reveal internal formation mechanism, the corresponding contact state between balls and raceways for ball bearing acted by the combined small axial load and large radial load are presented. As described in Figure 11a, the balls number in the loaded zone gradually decreases with the radial load. It should be noticed that, as the ball-raceway contact state change from full balls contact to partial balls contact, the change trend of bearing stiffness also changes. The axial stiffness changes from slow decline to rapid decline and the radial stiffness changes from decline to rising. In addition, as seen from Figure 11b, for the given radial load, the number of balls located at the loaded zone gradually increases with the axial load. Looking back at Figure 10, it can be found that the increasing trend of ball radial stiffness with axial load changes from slow rising to rapid rising when the ball-raceway contact state changes from partial balls contact to full balls contact. Above all, combining Figures 10 and 11, one can find that the changing trend of the stiffness for ball bearing in full balls contact condition is the same as the results in Figure 9, and the changing trend of bearing stiffness depends on the change of the internal contact state between balls and raceways. Therefore, one can conclude that the ball-raceway contact state has significant influence on the variation of bearing stiffness. way contact state change from full balls contact to partial balls contact, the change trend of bearing stiffness also changes. The axial stiffness changes from slow decline to rapid decline and the radial stiffness changes from decline to rising. In addition, as seen from Figure 11b, for the given radial load, the number of balls located at the loaded zone gradually increases with the axial load. Looking back at Figure 10, it can be found that the increasing trend of ball radial stiffness with axial load changes from slow rising to rapid rising when the ball-raceway contact state changes from partial balls contact to full balls contact. Above all, combining Figures 10 and 11, one can find that the changing trend of the stiffness for ball bearing in full balls contact condition is the same as the results in Figure 9, and the changing trend of bearing stiffness depends on the change of the internal contact state between balls and raceways. Therefore, one can conclude that the ball-raceway contact state has significant influence on the variation of bearing stiffness. . The contact states between balls and raceways for ball bearing acted by the small axial load and large radial load: (a) the contact states between balls and raceways varies with axial force; (b) the contact states between balls and raceways varies with radial force.

Stiffness Change with Internal Structural Parameters
In the above sections, the influence of operation conditions (speed and loads) on bearing stiffness has been described in detail. Due to the bearing internal structural parameter also having importance on bearing stiffness, in the following chapter, the high-precision way contact state change from full balls contact to partial balls contact, the change trend of bearing stiffness also changes. The axial stiffness changes from slow decline to rapid decline and the radial stiffness changes from decline to rising. In addition, as seen from Figure 11b, for the given radial load, the number of balls located at the loaded zone gradually increases with the axial load. Looking back at Figure 10, it can be found that the increasing trend of ball radial stiffness with axial load changes from slow rising to rapid rising when the ball-raceway contact state changes from partial balls contact to full balls contact. Above all, combining Figures 10 and 11, one can find that the changing trend of the stiffness for ball bearing in full balls contact condition is the same as the results in Figure 9, and the changing trend of bearing stiffness depends on the change of the internal contact state between balls and raceways. Therefore, one can conclude that the ball-raceway contact state has significant influence on the variation of bearing stiffness. . The contact states between balls and raceways for ball bearing acted by the small axial load and large radial load: (a) the contact states between balls and raceways varies with axial force; (b) the contact states between balls and raceways varies with radial force.

Stiffness Change with Internal Structural Parameters
In the above sections, the influence of operation conditions (speed and loads) on bearing stiffness has been described in detail. Due to the bearing internal structural parameter also having importance on bearing stiffness, in the following chapter, the high-precision Figure 11. The contact states between balls and raceways for ball bearing acted by the small axial load and large radial load: (a) the contact states between balls and raceways varies with axial force; (b) the contact states between balls and raceways varies with radial force.

Stiffness Change with Internal Structural Parameters
In the above sections, the influence of operation conditions (speed and loads) on bearing stiffness has been described in detail. Due to the bearing internal structural parameter also having importance on bearing stiffness, in the following chapter, the highprecision ball bearing 7008C is a case to discuss the influences of the various structural parameters on the stiffness.
As an important design parameter, the initial contact angle has great effects on the dynamic characteristic for ball bearings at high-speed range. The size of the initial contact angle can be adjusted as Equation (27) where P d denotes the internal radial clearance. It can be seen from Figure 12 that the influence rule of internal structural parameters variations on bearing initial contact angle is given. One can find that the initial contact angle decreases with the curvature center's distance between inner and outer raceways (both r i (f i ) and r o (f o ) increase). Besides, the initial contact angle increases with the radial clearance (d i decreases, d o increases, and D decreases).
including ri (fi), ro (fo), D, di, and do: where denotes the internal radial clearance. It can be seen from Figure 12 that the influence rule of internal structural parameters variations on bearing initial contact angle is given. One can find that the initial contact angle decreases with the curvature center's distance between inner and outer raceways (both ri (fi) and ro (fo) increase). Besides, the initial contact angle increases with the radial clearance (di decreases, do increases, and D decreases).

Single Structural Parameter Variation
In order to reveal the influence rules of single structural parameter variation on bearing stiffness under different speeds, the precision angular contact ball bearing 7008C under medium preload (290 N) has been chosen as the calculation example, and the key structural parameters are shown in Table 1. On this basis, the detailed calculation results are given in Figures 13-15. Figure 13 illustrates the influence of groove curvature coefficient (radius) variation on bearing stiffness under different speeds. It can be seen from Figure 13a,b that, for ball bearing at a static state or low speed, the axial stiffness decreases and radial stiffness remains basically unchanged with the rise of the curvature coefficient (radius) of inner raceway. Besides, when bearing under high-speed range, both the radial and axial stiffness increase with the curvature coefficient of inner raceway. Besides, it can be seen from Figure 13c,d that both the radial and axial stiffness decrease with the curvature radius of the outer raceway for the bearing at any speeds.

Single Structural Parameter Variation
In order to reveal the influence rules of single structural parameter variation on bearing stiffness under different speeds, the precision angular contact ball bearing 7008C under medium preload (290 N) has been chosen as the calculation example, and the key structural parameters are shown in Table 1. On this basis, the detailed calculation results are given in Figures 13-15. According to Figure 12, the size of the initial contact angle increases with the decrease of ri (fi) or ro (fo). Generally, for bearing at static state, the axial stiffness increases and radial stiffness decreases against the initial contact angle. However, for the raceway cur- with the increase of ball diameter, axial stiffness decreases and radial stiffness increases with the decrease of initial contact angle. When ball bearing is under the high-speed range, the axial stiffness does not decrease and slightly increases the ball diameter, it may be caused by the effect of the inertia forces of balls on the contact forces between outer raceways and balls.  with the decrease of initial contact angle. When ball bearing is under the high-speed range, the axial stiffness does not decrease and slightly increases the ball diameter, it may be caused by the effect of the inertia forces of balls on the contact forces between outer raceways and balls.   Figure 13 illustrates the influence of groove curvature coefficient (radius) variation on bearing stiffness under different speeds. It can be seen from Figure 13a,b that, for ball bearing at a static state or low speed, the axial stiffness decreases and radial stiffness remains basically unchanged with the rise of the curvature coefficient (radius) of inner raceway. Besides, when bearing under high-speed range, both the radial and axial stiffness increase with the curvature coefficient of inner raceway. Besides, it can be seen from Figure 13c,d that both the radial and axial stiffness decrease with the curvature radius of the outer raceway for the bearing at any speeds.
According to Figure 12, the size of the initial contact angle increases with the decrease of r i (f i ) or r o (f o ). Generally, for bearing at static state, the axial stiffness increases and radial stiffness decreases against the initial contact angle. However, for the raceway curvature coefficients variation, the change rule of bearing stiffness and contact angle is slightly different. Combining Figures 12 and 13, one can find that the radial stiffness basically remains unchanged with the changes of initial contact angle and raceway curvature coefficient. It is caused from, with the decrease of raceway curvature coefficient, the contact surfaces between ball and raceways becoming more fit and contact stiffness increasing, thus suppressing the decrease of bearing radial stiffness. However, for bearing at high speed, with the decrease of inner raceway curvature coefficient, radial stiffness decreases rapidly. With the decrease of outer raceway curvature coefficient, both radial and axial stiffness increase rapidly. It is caused by that the contact forces among balls and outer raceways increase rapidly as the rotating speed rise, and the sizes of bearing stiffness mainly depend on the contact angle and contact stiffness of ball-outer raceway that increases rapidly with the decrease of the curvature coefficient of outer raceway. Therefore, when ball bearing runs at high-speed range, when the outer raceway curvature coefficient decreases, both axial and radial stiffness increase rapidly.
Besides, according to Equation (27) and Figure 12, one can find that the raceway contact diameters also have great influence on the initial contact angle and radial clearance of ball bearing. Figure 14 illustrates the influences of the raceway contact diameters variations on bearing stiffness under different speeds is given. With the increase of the contact diameter of inner raceway, the axial stiffness decreases and radial stiffness increases, and axial stiffness tends to be constant for the bearing at high speed. The change rule of bearing stiffness with outer raceway contact diameter variation is exactly the opposite.
In addition, the influence of the ball diameter variation on bearing stiffness under different speeds is also given in Figure 15. When the bearing is at low and medium speeds, with the increase of ball diameter, axial stiffness decreases and radial stiffness increases with the decrease of initial contact angle. When ball bearing is under the high-speed range, the axial stiffness does not decrease and slightly increases the ball diameter, it may be caused by the effect of the inertia forces of balls on the contact forces between outer raceways and balls.

Constant Initial Contact Angle
In the last chapter, the influences of single structural parameter variation on bearing stiffness under different speeds were described in detail. One can find that, as the internal structural parameters change, the initial contact angle is changed. However, in the actual bearing design, the initial contact angle should be a specific value, for example, the initial contact angle is about 15 • for the 7008C series bearing. Therefore, the structural parameters of ball bearing are changed simultaneously to ensure a constant contact angle of the 7008C in this article research.
According to the working condition, the bearing clearance is divided into multi-levels, and the size of bearing clearance can be adjusted by changing the sizes of ball and raceway contact diameters. Without loss of generality, three different clearances, 30 µm, 50 µm, and 70 µm, are chosen, and the raceway curvature coefficients of inner and outer rings have been simultaneously changed to ensure the constant initial contact angle. It can be seen from Figure 16, for bearing at static state, with the increase of inner raceway curvature coefficient, the bearing stiffness first decreases and then increases, while for bearing at high speed, the bearing stiffness continuously increases with the curvature coefficient of the inner raceway. It means that the bearing designed with a smaller outer-raceway curvature coefficient and a larger inner-raceway curvature coefficient can improve the stiffness of ball bearing at high-speed range. Besides, with the increase of bearing clearance, both axial and radial stiffness is reduced, which means that a smaller clearance is also useful to improve the bearing stiffness. However, too small a clearance is more likely to lead to a thermal failure for ball bearing at the high-speed range.
bearing stiffness gradually decreases as the rotating speed rises, and the rates of decline depend on the raceway curvature coefficient. The projection of Figure 17a,b illustrates the relation between the radial stiffness and rotating speeds of ball bearing with different raceway curvature coefficients and the radial stiffness and raceway curvature coefficients of ball bearing under different speeds, respectively. One can find that, with the increase of inner raceway curvature coefficient, the downtrend of bearing stiffness with rotating speed tends to be gentle.  In order to further reveal the relations among the rotating speed, raceway curvature coefficient, and bearing stiffness, the three-dimension curved of the rotating speed, inner raceway curvature coefficient, and radial stiffness is given in Figure 17a. One can find that bearing stiffness gradually decreases as the rotating speed rises, and the rates of decline depend on the raceway curvature coefficient. The projection of Figure 17a,b illustrates the relation between the radial stiffness and rotating speeds of ball bearing with different raceway curvature coefficients and the radial stiffness and raceway curvature coefficients of ball bearing under different speeds, respectively. One can find that, with the increase of inner raceway curvature coefficient, the downtrend of bearing stiffness with rotating speed tends to be gentle.

Discussion and Conclusions
Based on the ball-raceway contact analysis, implicit/explicit differential method, and matrix operations, a comprehensive 5-DOF analytical model for stiffness matrix calculation of ball bearing under arbitrary loads condition has been proposed. Through comparison to the representative experimental and theoretical results in other literatures, the proposed model has been verified. Based on these, the influences of various operation condi- Figure 17. The influences of raceway curvature coefficient and rotating speed on bearing stiffness with constant initial contact angle and radial clearance: (a) the three-dimension curved of the rotating speed, inner raceway curvature coefficient, and radial stiffness; (b) the two-dimensional projection curves of three-dimensional surface (a). Q ik , Q ok contact loads of the kth ball-inner/outer raceway r i , r o inner/outer raceway curvature radius R xik , R xok equivalent curvature radius of the inner and outer ball-raceways contact T ik , T ok frictions of the kth ball-inner/outer raceway Z number of rolling balls As described in the classical Hertz contact theory, the load-deflection coefficients K ik and K ok can be expressed by using the least square linear regression algorithm: The effective contact radii of the ball-inner and ball-outer raceways are written as: