machines

Article
Analytical Determination and Influence Analysis of Stiffness
Matrix of Ball Bearing under Different Load Conditions
Qingbo Niu 1,2 , Yeteng Li 1 , Yongsheng Zhu 1, *, Shiyuan Pei 1 , Yanjing Yin 2 and Dongfeng Wang 2

1 Key Laboratory of Education Ministry for Modern Design and Rotor-Bearing System,
Xi’an Jiaotong University, Xi’an 710049, China; zys2008@126.com (Q.N.); liyeteng@stu.xjtu.edu.cn (Y.L.);
peishiyuan@xjtu.edu.cn (S.P.)
2 Luoyang Bearing Research Institute Co., Ltd., Luoyang 471039, China; zys_yyj@126.com (Y.Y.);
zyswdf@163.com (D.W.)
* Correspondence: yszhu@mail.xjtu.edu.cn; Tel.: +86-139-9114-9360

Abstract: 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, im-
plicit/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 influences 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
significant influence on the change in the stiffness of ball bearings. Besides, a larger curvature coeffi-
cient of inner raceway and a small curvature coefficient 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
Citation: Niu, Q.; Li, Y.; Zhu, Y.; Pei, optimize design and performance analysis of ball and rotor system under various load conditions.
S.; Yin, Y.; Wang, D. Analytical
Determination and Influence Keywords: stiffness matrix; ball bearing; analytical expressions; influence factors; optimization of
Analysis of Stiffness Matrix of Ball bearing design
Bearing under Different Load
Conditions. Machines 2022, 10, 238.
https://doi.org/10.3390/
machines10040238
1. Introduction
Academic Editor: Hui Ma As the key supporting parts in high-speed rotor system, rolling bearings are widely
Received: 16 February 2022
used across numerous rotating mechanical systems, including aerospace, automobile,
Accepted: 22 March 2022
machine tool, etc. With the rapid development of the modern industry, there is demand
Published: 28 March 2022
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
Publisher’s Note: MDPI stays neutral
use parameters of rolling bearing should be adjusted by the designers and users according
with regard to jurisdictional claims in
to the practical operation conditions, and an accurate and efficient mathematical model of
published maps and institutional affil-
bearing can be used to reduce the cost and cycle.
iations.
During the past decades, a large number of algorithms and methods have been
presented in bearing modeling and analysis. In the earlier studies [1–4], the inertia forces
of balls have been neglected, which may cause a considerable error in bearing analysis and
Copyright: © 2022 by the authors.
characteristic prediction, especially for the high-speed operation condition. To improve
Licensee MDPI, Basel, Switzerland. the calculation accuracy of bearing analysis, Jones [5] proposed a classical mathematic
This article is an open access article model with the “Raceway Control Hypothesis” for ball bearing subjected to the combined
distributed under the terms and radial, axial, and moment loads. In this model, the effect of the inertia forces of the
conditions of the Creative Commons balls are considered. Without the Raceway Control Hypothesis, Harris [6] developed a
Attribution (CC BY) license (https:// comprehensive model based on the lubrication analysis of ball–raceway contact area, and
creativecommons.org/licenses/by/ this model can be used in the bearing skidding behavior prediction and friction power loss
4.0/). calculation. In addition, the experimental results show that the partial results obtained from

Machines 2022, 10, 238. https://doi.org/10.3390/machines10040238 https://www.mdpi.com/journal/machines

Machines 2022, 10, 238 2 of 21

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–9]. In addition, a
serial of dynamic models [10–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–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–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–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.
 Machines 2022, 10, 238 3 of 21

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–24,32,34] and the stiffness calculation has numerical stability prob-
lems [28–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 com-
prehensive 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.

2. Theoretical Model
2.1. 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:
Machines 2022, 10, x FOR PEER REVIEW 4 of 22
2π
ψk = ( k − 1) (1)
Z

Figure
Figure 1. The
1. The coordinate
coordinate systems
systems for for
ballball bearing:
bearing: (a) Axial
(a) Axial section
section andand x axis
x axis coincided
coincided with
with thethe
rotation center line; (b) Traverse section and ball position angle (ψ = 2π ( k −
rotation center line; (b) Traverse section and ball position angle ( ψ k =k 2π (k − 1) Z ).
1 ) /Z).

It is assumed that the bearing is acted by a combined external load vector F = {Fx ,
FyIt
, Fisz , assumed
My , Mz }, that the relative
and the bearing displacement
is acted by a combined external
vector between theload vector
inner F = {Frings
and outer x, Fy,
is
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
 Figure 1. The coordinate systems for ball bearing: (a) Axial section and x axis coincided with the
rotation center line; (b) Traverse section and ball position angle ( ψ k = 2π (k − 1) Z ).
Machines 2022, 10, 238 4 of 21
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
d = {δx , δy , δz , θy , θz }. Assuming that the rotating speed of inner ring is ωi . It can be seen
Figure 2a that the position of ball center moves upward due to the ball centrifugal force.
from Figure 2a that the position of ball center moves upward due to the ball centrifugal
Consequently, the contactthe
force. Consequently, angle between
contact anglethe inner raceway
between the innerand ball increases,
raceway while thewhile
and ball increases,
contact
the contact angle between the outer raceway and ball decreases. It should bethat
angle between the outer raceway and ball decreases. It should be noticed the that
noticed
abovetheanalysis is only applicable to the ball–raceway contact situation, but when
above analysis is only applicable to the ball–raceway contact situation, but when ball ball bear-
ing isbearing
acted by the combined
is acted loads with
by the combined a small
loads with aaxial
smallforce,
axial parts
force,of balls
parts ofmay
ballsno
may longer
no longer
be loaded, and the bearing is divided into two parts: the loaded zone and
be loaded, and the bearing is divided into two parts: the loaded zone and the unloaded the unloaded
zone zone
[45]. It canItbe
[45]. seen
can from from
be seen Figure 2b that
Figure 2b parts of theofballs
that parts are separate
the balls from from
are separate the inner
the inner
raceway by the action of centrifugal force. However, this phenomenon
raceway by the action of centrifugal force. However, this phenomenon is usually is usually ignored
ignored
in modeling
in modelingand stiffness calculation
and stiffness of ball
calculation ofbearing under
ball bearing the action
under of combined
the action loads.
of combined loads.

FigureFigure The geometry

2. The2.geometry relations
relations among among the of
the centers centers of raceway
raceway curvaturecurvature
and ball: and ball:
(a) The (a)con-
ball The ball
tacts with inner
contacts raceway;
with (b) The ball
inner raceway; (b) separates from inner
The ball separates raceway.
from inner raceway.

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:
A1k − X1k A2k − X2k
(
sinαik = ri −0.5D +δik cosαik = ri −0.5D +δik
X1k X2k (2)
sinαok = ro −0.5D +δok cosαok = ro −0.5D +δok

where ri , ro 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:

A1k = BDsinα0 + 5dm θy sinψk − 0.5dm θz cosψk


(3)
A2k = BDcosα0 + δy cosψk + δz sinψk

where D and dm 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:

Eq1k = A1k − ∆ik sinαik − ∆ok sinαok = 0


(4)
Eq2k = A2k − ∆ik cosαik − ∆ok cosαok = 0

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]:

Oo Od0 = ( f o − 0.5) D + δok (5)

Machines 2022, 10, 238 5 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]:
r
 2  2
Oi0 Od0 = POd0 + POi0 ≤ ( f o − 0.5) D (6)

where the distances POd0 and POi0 can be written as:

(
POd0 = A2k − [( f o − 0.5) D + δok ]
(7)
POi0 = A1k

when the ball is located in the loaded zone, the mechanical equilibrium equations of ball
are expressed as:

Eq3k = Qik sinαik − Qok sinαok − Tik cosαik + Tok cosαok = 0
(8)
Eq4k = Qik cosαik − Qok cosαok − Tik sinαik + Tok sinαok + Fck = 0

where Qik and Qok are the ball–inner/outer raceway contact loads, respectively. Tik and Tok
are the ball–inner/outer frictions, respectively.

3 3
Qik = Kik δik2 Qok = Kok δok
2

M M (9)
T = λ gk gk
ik T =λ
ik D ok ok D

In the above equations, the Kik and Kok are the load-deformation coefficients in ball–
raceways contacts, the values of Kik and Kok are determined by the basic structural pa-
rameters and the contact angles, and the analytical calculation expressions Kik and Kok 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 Q+ikQ
ik ok
and λok = 2 Q Q+okQ can be selected according to equal friction coefficient assumption [8].
ik ok
Besides, the inertia forces, including the gyroscopic moment Mgk and the centrifugal force
Fck of ball, are expressed:
2

Fck = 0.5mdm ωmk
(10)
Mgk = Jωmk ωbk sinβ k
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:

 ω = ωi (1−γik )cos(αok − β k )
mk (1+γok )cos(αik − β k )+(1−γik )cos(αok − β k )
ωi (1−γik )(1+γok ) (11)
 ωbk = dm
D (1+γ )cos(α − β )+(1−γ )cos(α − β )
ok ik k ik ok k

where γik = Dcosα

dm
ik
and γok = Dcosαok
dm , and the pitch angle β k are given based on the
d’Alembert principle [9]:

M ( N + 1)sinαik + 2sinαok
tanβ k =   (12)
M ( N + 1)cosαik + 2 cosαok + dDm + G

with
 𝐺= 𝑀 𝑐𝑜𝑠(𝛼 − 𝛼 ) − 𝑁
𝑑 (13)
⎨ 𝐷
⎪ 1+
⎪𝑁 = 𝑑 𝑐𝑜𝑠 𝛼
⎩ 1 − 𝑐𝑜𝑠 𝛼
Machines 2022, 10, 238 where 𝑎 and 𝑎 the semi-major axes of the ball–raceways contacts. 𝐿 and 𝐿 6are
of 21
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:
 M = QQikaaik LLik
 𝐸𝑞 =D 𝑄 ok ok ok
− 𝐹(α =−0α ) − N]

G = dm M [cos (14)
ik ok (13)
 D
1+ d cosα
Through the above analysis, N it= can be found that the ball in the loaded zone has 4
m

 ok
1−cosαik
unknown variables 𝛆 = {𝑋 ,𝑋 ,𝛿 ,𝛿 } in Equations (2) and (3). The Newton–Raphson
where aikhas
algorithm aok the
andbeen usedsemi-major
to calculateaxes
theseofvariables
the ball–raceways
for a givencontacts. Lik and Lok value
global displacement are the
second kind elliptic integrals of the ball–raceways contacts [9].
vector d = {𝛿 ,𝛿 ,𝛿 ,𝜃 ,𝜃 }. In addition, it can be seen from Equation (14) that only one
Besides,
variable 𝛿 is itused
can to
be describe
seen fromthe Figure
force3b, if the
state of ball
ball separates from thezone.
in the unloaded innerHowever,
raceway, the
it
force equilibrium equation of ball is given as:
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)
Eq4k = Qok − Fck = 0 (14)
are solved.

Figure 3. Force analysis of the ball, including the inertia forces (Fck , Mgk ) and contact forces
(Qi/ok , Ti/ok ): (a) The ball contacts with inner ring; (b) The ball separates from inner ring.

Through the above analysis, it can be found that the ball in the loaded zone has
4 unknown variables εk = {X1k , X2k , δ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:
  
1 0
 
Fx 





 Fy  Z

 0 cosψ k
   

Q sinα − T cosα
   
 Fz  = ∑ G [δik ]  ik ik ik ik
   0 sinψ k

× (15)
Qik cosαik + Tik sinαik 


 My 

1
 dm
sinψ 0

  

  2 k  

Mz  − dm cosψ

0


2 k

where G[•] is the Heavisde function, it can be expressed as:

(
1 Oi0 Od0 > ro − 0.5D
G [δik ] = (16)
0 Oi0 Od0 ≤ ro − 0.5D

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 = {X1k , X2k , δik , δok };
Secondly, solve the local equilibrium equations update the values of local variables and ball–
raceway 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
 Machines 2022, 10, 238 7 of 21

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 Jm in
inner iteration and stiffness matrix K in outer loop iteration need to be determined, and
the matrix Jm denotes the Jacobian matrix of the ball local equations to ball local variables
εk = {X1k , X2k , δik , δok }. The detailed expressions are given as follows:
 ∂Eq1k   ∂Eq1k ∂Eq1k ∂Eq1k ∂Eq1k 
∂ε k ∂X1k ∂X2k ∂δik ∂δok
 ∂Eq2k   ∂Eq2k ∂Eq2k ∂Eq2k ∂Eq2k 
 ∂ε k
  ∂X1k ∂X2k ∂δik ∂δok

Machines 2022, 10, x FOR PEER REVIEW
Jm = 
 ∂Eq3k
=
  ∂Eq3k ∂Eq3k ∂Eq3k ∂Eq3k

 (17)
8 of 22
 ∂ε k   ∂X1k ∂X2k ∂δik ∂δok 
∂Eq4k ∂Eq4k ∂Eq4k ∂Eq4k ∂Eq4k
∂ε k ∂X1k ∂X2k ∂δik ∂δok

Figure 4. Calculation flow of proposed model, including the internal iteration of balls and the external
Figure 4. Calculation
iteration flow
of bearing of proposed model, including the internal iteration of balls and the exter-
ring.
nal iteration of bearing ring.
For the ball in the unloaded zone, and Oi0 Od0 ≤ ( f o − 0.5) D. Jm are given as follows:
2.2. Computation of the Stiffness Matrix
∂Eq4k
The stiffness matrix of ball bearing can
Jm =be calculated as Equation (19) through the(18)
∂δokdisplacement vector d:
partial derivatives between the load vector F to the
𝜕𝐹 𝜕𝐹 𝜕𝐹 𝜕𝐹 𝜕𝐹
In addition, the stiffness⎡matrix K of ball bearing is illustrated
⎤ in the subsequent chapters.
⎢ 𝜕𝛿 𝜕𝛿 𝜕𝛿 𝜕𝜃 𝜕𝜃 ⎥
⎢ 𝜕𝐹 𝜕𝐹 𝜕𝐹 𝜕𝐹 𝜕𝐹 ⎥
⎢ 𝜕𝛿 𝜕𝛿 𝜕𝛿 𝜕𝜃 𝜕𝜃 ⎥
⎢ 𝜕𝐹 𝜕𝐹 𝜕𝐹 𝜕𝐹 𝜕𝐹 ⎥
𝑲 × =⎢ ⎥ (19)
⎢ 𝜕𝛿 𝜕𝛿 𝜕𝛿 𝜕𝜃 𝜕𝜃 ⎥
⎢𝜕𝑀 𝜕𝑀 𝜕𝑀 𝜕𝑀 𝜕𝑀 ⎥
⎢ 𝜕𝛿 𝜕𝛿 𝜕𝛿 𝜕𝜃 𝜕𝜃 ⎥
⎢ ⎥
⎢ 𝜕𝑀 𝜕𝑀 𝜕𝑀 𝜕𝑀 𝜕𝑀 ⎥
⎣ 𝜕𝛿 𝜕𝛿 𝜕𝛿 𝜕𝜃 𝜕𝜃 ⎦
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 varia-
ble transfer graph is given in Ref [45]. Previous studies [45–48] show that the calculation
process of the stiffness matrix should be divided into the partial derivatives of the load
Machines 2022, 10, 238 8 of 21

2.2. 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:
 ∂F ∂F ∂F ∂F ∂F

x x x x x
∂δx ∂δy ∂δz ∂θy ∂θz
 ∂Fy ∂Fy ∂Fy ∂Fy ∂Fy 
 
 ∂δx ∂δy ∂δz ∂θy ∂θz 
 ∂Fz ∂Fz ∂Fz ∂Fz ∂Fz 
[ K ] 5×5 =
 ∂δx ∂δy ∂δz ∂θy ∂θz

 (19)
 ∂My ∂My ∂My ∂My ∂My 
 
 ∂δx ∂δy ∂δz ∂θy ∂θz 
∂Mz ∂Mz ∂Mz ∂Mz ∂Mz
∂δx ∂δy ∂δz ∂θy ∂θz

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–48] show that the cal-
culation 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 = {X1k , X2k , δik , δok , Mgk , Fck } as the intermediate variables. However, this selection strat-
egy 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 = {X1k , X2k , δik , δok , A1k , A2k } have been used in the proposed
method. Therefore, the stiffness matrix is written as:
 
 ∂Fxk  


 ∂υk 



  ∂F yk  

  
Z  ∂υk 
 h T T T T T
i 

[K]5×5 = ∑ G [δik ]  ∂υ 
 ∂Fzk  ∂υk ∂υk ∂υk ∂υk ∂υk
× ∂δx ∂δy ∂δz ∂θy ∂θz
(20)
1
  ∂Mk  6×5 
yk

   


  ∂υk  


 


 ∂M zk 

∂υk 5×6

where the partial derivatives between the loads vector F to the intermediated variables υk
can be written as:
 ∂Fxk 
 
∂υ 1 0
 ∂Fykk 

 ∂υk 
 
 0 cosψk  
 ∂Fzk  0 sinψk 
=  ×

 ∂υk 
 dm
 2 sinψk 0 
 ∂M  
yk 
 (21)
− d2m cosψk
 ∂υk 
∂Mzk 0 5×2
" ∂υk 5×6 #
∂Qik ∂Tik ∂sinαik
∂υk sinα ik − ∂υk cosα + Qik ∂υk ik − Tik ∂cosα
∂υ
ik
k
∂Qik ∂Tik ∂cosαik
∂υk cosαik + ∂υk sinαik + Qik ∂υk + Tik ∂sinα ik
∂υ k 2×6

with  3 h 1
i
∂Qik ∂Kik 2

∂υk = ∂υk δik + Kik 0 0 2
1.5δik 0 0 0
3 h 1
i (22)
∂Qok ∂Kok 2

∂υk = ∂υk δok + Kok 0 0 0 2
1.5δok 0 0
In the previous literatures [18,19], in order to simplify the calculation process, the
contributions of the load-deformation coefficients Kik and Kok to bearing stiffness calculation
have been ignored. However, due to the apparent change of the contact angles, it may lead
Machines 2022, 10, 238 9 of 21

to a larger error when ball bearing operates at the high-speed range. Therefore, the detailed
calculation expression of ∂K ik ∂Kok
∂υk and ∂υk is given in Appendix A.
In addition, the partial derivatives between the intermediated variables υk = {X1k , X2k ,
δik , δok , A1k , A2k } to global displacement vector d is written as:

∂υkT ∂υkT ∂υkT ∂υkT ∂υkT

h i
∂δx ∂δy ∂δz ∂θy ∂θz 6×5
 ∂X1k ∂X1k 
∂A1k ∂A2k
 ∂X2k ∂X2k 
 ∂A1k ∂A2k   
∂A1k ∂A1k ∂A1k ∂A1k ∂A1k

 ∂δik ∂δik

 ∂δx ∂δy ∂δz ∂θy ∂θz
(23)
= ∂A1k ∂A2k  × ∂A2k ∂A2k ∂A2k ∂A2k ∂A2k

 ∂δok ∂δok 
  ∂δx ∂δy ∂δz ∂θy ∂θz
 ∂A1k ∂A2k  2×5
 1 0 
0 1 6×2

where the partial derivatives ∂A 1k ∂A2k

∂d and ∂d terms are calculated by the explicit differential
of Equations (3) and (4). However, the relations among the εk , A1k , and A2k are determined
by Equations (2) and (3) in implicit form, so ∂A 1k ∂A2k
∂d and ∂d terms can be obtained by the
implicit differential method. Therefore, the detailed calculation process is given as follows:
( h i
1 0 0 d2m sinψk − d2m cosψk
∂A1k
∂d = (24)
∂A2k  
∂d = 0 cosψk sinψk 0 0

∂X1k
∂A1k = − 1J ∂(Eq 1k ,Eq2k ,Eq3k ,Eq4k )
∂( A ,X ,δ ,δ )
∂X1k
∂A2k = − 1J ∂(Eq 1k ,Eq2k ,Eq3k ,Eq4k )
∂( A ,X ,δ ,δ )
1k 2k ik ok 2k 2k ik ok
∂X2k
∂A1k = − 1J ∂(Eq 1k ,Eq2k ,Eq3k ,Eq4k )
∂( X ,A ,δ ,δ )
∂X2k
∂A2k = − 1J ∂(Eq 1k ,Eq2k ,Eq3k ,Eq4k )
∂( X ,A ,δ ,δ )
1k 1k ik ok 1k 2k ik ok
(25)
∂δik
∂A1k = − 1J ∂(∂Eq(X1k ,Eq 2k ,Eq3k ,Eq4k )
,X ,A ,δ )
∂δik
∂A2k = − 1J ∂(∂Eq(X1k ,Eq 2k ,Eq3k ,Eq4k )
,X ,A ,δ )
1k 2k 1k ok 1k 2k 2k ok
∂δok
∂A1k = − 1J ∂(Eq 1k ,Eq2k ,Eq3k ,Eq4k )
∂( X ,X ,δ ,A )
∂δok
∂A2k = − 1J ∂(Eq 1k ,Eq2k ,Eq3k ,Eq4k )
∂( X ,X ,δ ,A )
1k 2k ik 1k 1k 2k ik 2k

where the Jacobi determinant J can be written as:

 ∂Eq1k 
∂υk (1 : 4)
∂Eq2k
∂υk (1 : 4) 
 
∂( Eq1k , Eq2k , Eq3k , Eq4k ) 
J= = det ∂Eq3k
 (26)
∂( X1k , X2k , δik , δok )
∂υk (1 : 4) 

 
∂Eq4k
∂υk (1 : 4)

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 quasi-
static 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.

3. 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.

3.1. Comparison against Experimental and Theoretical Results

To show the effectiveness of the presented model, the calculation results of the high-
precision ball bearings 7008C were compared to the results obtained by the most rep-
resentative 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
Machines 2022, 10, 238 10 of 21

compared against the published experimental and FEM results. As shown in Table 1, the
key structural parameters of the above ball bearings are given.

Table 1. T Bearing parameters of NSK B7008C, B7014C, and SKF EEB3-2Z.

Machines 2022, 10, x FOR PEER REVIEW NSK NSK SKF SKF
11 of 22
Parameters
B7008C B7014C EEB3-2Z 6205
The curvature radius of inner raceway ri (mm) 4.000 5.330 2.064 4.108
The curvature radius of outer raceway ro (mm)
To show the effectiveness of the presented model, the calculation results of the4.108
3.790 5.050 2.064
high-
The contact diameter of inner raceway di (mm) 46.838 80.452 11.78 30.59
precision ball bearings 7008C were compared
The contact diameter of outer raceway do (mm) to the results
61.176 obtained
99.562 by the
19.72 most repre-
46.41
sentative methods
The number[18]ofatinternal
first. Then
ballsthe
Z static stiffness 19 results for
25 high-precision
7 ball bear-
9
ings 7008CTheanddiameter
7014C ofwereballscompared
D (mm) against the7.144 experimental
9.522 values3.969
in bearing user
7.900
The initial
manual. contact
At last, angle of ball
the dynamic α0 (◦ ) of bearing
bearingresults
stiffness 15.992 EEB3-2Z
15.942at 10006.429
rpm were 14.453
com-
paredThe pitch diameter
against of ball bearing
the published dm (mm)
experimental and FEM50.00results. As90.00
shown in15.75 38.50
Table 1, the key
structural parameters of the above ball bearings are given.
ItItcan
canbe
beseen
seenfrom
fromFigure
Figure55that
thatthe
thecontrast
contrastresults
resultsofofthe
thebearing
bearingstiffness
stiffnessvary
varywith
with
the
the speed
speed calculated
calculated by by the
the proposed
proposed model
model and
and the
the Cao–Jones
Cao–Jones model
model given.
given. ItItcan
can be
be
found
foundthat thatthe
thestiffness
stiffnessvalues
valuesobtained
obtainedbybythe
theproposed
proposedmodelmodelare
are located
locatedbetween
betweenthe the
stiffness
stiffness values
values calculated
calculated by by the
theCao–Jones
Cao–Jonesmodel
modelwithwiththe
the “Inner/OuterRaceway
“Inner/Outer Raceway Con-
Control
trol Hypothesis”.
Hypothesis”. As “Raceway
As the the “Raceway Control
Control Theory”
Theory” is based
is based on extreme
on extreme conditions,
conditions, the
the results
results obtained by the proposed model are more reasonable
obtained by the proposed model are more reasonable and logical. and logical.

(a) (b)
35 180
Radial stiffness (N/μm)

160 Cao-Jones
Axial stiffness (N/μm)

30 Cao-Jones model
(I-RCH) model
140
(I-RCH)
25 120
The proposed
The proposed
model 100
20 model
Cao-Jones model 80 Cao-Jones model
(O-RCH) (O-RCH)
15 60
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
Rotating speed (Krpm) Rotating speed (Krpm)

Figure
Figure 5. Comparison ofofbearing
5. Comparison bearingstiffness
stiffness varying
varying withwith the speed
the speed calculated
calculated by theby the proposed
proposed method
method with Cao-Jones method with the “Inner/Outer Raceway Control Hypothesis
with Cao-Jones method with the “Inner/Outer Raceway Control Hypothesis (I/O RCH)”: (a) axial (I/O RCH)”:
(a) axial stiffness; (b) radial stiffness.
stiffness; (b) radial stiffness.

Table Then,
1. T Bearing parameters
as shown of NSK
in Table B7008C,
2, the B7014C, and
comparison SKF EEB3-2Z.
results of static stiffness between the
proposed model and user manual values are given.NSK One can find
NSKthat theSKFvalues of SKF
stiffness
Parameters
obtained by the proposed model are slightly more than the experimental values in bearing
B7008C B7014C EEB3-2Z 6205
user manual, and it may be caused by the improvement of material properties as the special
The curvature radius of inner raceway ri (mm) 4.000 5.330 2.064 4.108
technology and heat-treatment are not considered in the model calculation.
The curvature radius of outer raceway ro (mm) 3.790 5.050 2.064 4.108
The contact diameter of inner raceway d i (mm) 46.838 80.452 11.78
Table 2. The comparison of bearing static stiffness (NSK B7008C, B7014C) with the experimental
30.59
The contact
values diameter
in bearing of outer
user manual raceway
(EL: extremelydo (mm)
light; L: 61.176 99.562 H: heavy).
light; M: medium; 19.72 46.41
The number of internal balls Z 19 25 7 9
Types The diameter7008C
of balls D (mm) 7.144 9.522 7014C3.969 7.900
Preload (N) The
60 initial
(EL) contact
120 (L)angle of ball(M)
10he290 bearing 590𝛼(H) (°) 145 15.992
(EL) 15.942
290 (L) 6.429
740 (M) 14.453
1470 (H)
Model The
37.5pitch diameter
48.7 of ball bearing
72.5 d m (mm)
99.3 50.00
65.1 90.00
83.9 15.75
124.4 38.50
170.1
Axial stiffness
Manual 39 51 77 110 68 88 135 190
(N/µm)
error Then, as shown
3.8% 4.5% in Table 2, the comparison
5.84% 9.72% results
4.26% of static
4.66%stiffness between10.47%
7.85% the pro-
posed 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 spe-
cial technology and heat-treatment are not considered in the model calculation.

Table 2. The comparison of bearing static stiffness (NSK B7008C, B7014C) with the experimental
 Axial stiff- Model 37.5 48.7 72.5 99.3 65.1 83.9 124.4 170.1
ness Manual 39 51 77 110 68 88 135 190
Axial stiff- Model
(N/µm) error 37.5
3.8% 48.7
4.5% 72.5
5.84% 99.3
9.72% 65.1
4.26% 83.9
4.66% 124.4
7.85% 170.1
10.47%
ness Manual 39 51 77 110 68 88 135 190
Machines 2022, 10, 238 (N/µm) error 3.8% 4.5% 5.84% 9.72% 4.26% 4.66% 7.85% 10.47% 11 of 21
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
Figure 6 gives the comparison results of the dynamic stiffness for bearing EEB3-2Z
the results given by the proposed method are more consistent with the experimental re-
obtained by different theoretical methods and experimental analysis. It can be found that
sults.Figure 6 gives the comparison results of the dynamic stiffness for bearing EEB3-2Z
the
obtained bygiven
results by the
different proposedmethods
theoretical methodand
are experimental
more consistent with the
analysis. experimental
It can re-
be found that
sults.
the results given by the proposed method are more consistent with the experimental results.

)
100

(N/μm
80
Radial stiffness

)
100

(N/μm
Stiffness 60
80
Radial stiffness
Stiffness
40
60
Axial/Radial

20
40
Axial stiffness
Axial/Radial

20050 100 150 200 250Axial300 350 400

stiffness
0 Axial load (N)
50 100 150 200 250 300 350 400
Axial load (N)
Figure 6. Comparison of the dynamic stiffness (1000 rpm) obtained by the proposed model (—□—
)Figure
to the 6.
Figure FEM
6. results (…△…)
Comparison
Comparison of the
of of Guo [37]
thedynamic
dynamic and the
stiffness
stiffness experiment
(1000
(1000 rpm) obtained
rpm) obtained (--○
resultsby
by --)proposed
the
the of Kraus model
proposed [41].
model ((——)
—□ —
to the FEM results ( . . . 4 . . . ) of Guo [37] and the experiment results (–#–) of Kraus [41].
) to the FEM results (…△…) of Guo [37] and the experiment results (--○--) of Kraus [41].
3.2. Stiffness Variation with Rotating Speed and Axial Preload
3.2. Stiffness Variation with Rotating Speed and Axial Preload
For the Variation
3.2. Stiffness precision andRotating
high-speed ball
andbearing, a proper axial preload can be used to
For the precisionwith and high-speed Speed ball Axial Preload
bearing, a proper axial preload can be used
improve the rotating accuracy and supporting rigidity of ball bearing. According to Figure
to improve
For the the rotating
precision and accuracy
high-speedand ball
supporting
bearing, rigidity
a properofaxial
ball preload
bearing.can According
be used to to
7, when ball bearing is acted by a certain axial preload, the stiffness change with the ro-
Figure 7,the
improve when ball bearing
rotating accuracyisand acted by a certain
supporting axialofpreload,
rigidity the stiffness
ball bearing. According change with
to Figure
tating speed can be divided into three stages: decline slowly—drop rapidly—trend to gen-
the
7, rotating
when ball speed
bearingcanis be divided
acted by a into three
certain stages:
axial decline
preload, the slowly—drop
stiffness change rapidly—trend
with the ro-
tle,
to and theand
gentle, samethechange
same rules have
change been
rules havefound
been infound
previous
in theoretical
previous analysesanalyses
theoretical and ex-
tating speed can be divided into three stages: decline slowly—drop rapidly—trend to gen-
perimental
andand studies
experimental [34,43,44].
studies rules Besides,
[34,43,44]. the decreasing trend of bearing stiffness with rotat-
tle, the same change have Besides,
been foundthe decreasing
in previoustrend of bearing
theoretical stiffness
analyses andwith
ex-
ing speedspeed
rotating has been eased.
has [34,43,44].
been Figure
eased. 7 also7 also
Figure givesgives
the results
thetrend of the
results of stiffness change withwith
the stiffness the
perimental studies Besides, the decreasing of bearing stiffnesschange
with rotat-
preload
the speed of ball
preload at different
of been
ball ateased. speeds,
different and and
speeds, the drastic
the the improvement
drastic improvement of the
of stiffness
thechange due due
stiffness to the
to
ing has Figure 7 also gives results of the stiffness with the
increase
the in
increase the
in preload
the is
preload presented.
is presented.
preload of ball at different speeds, and the drastic improvement of the stiffness due to the
increase in the preload is presented.
(a) (b)
50 50 250 250
(a) (b)
(N/μm)

10,000 rpm Fp=590 N 10,000 rpm

(N/μm)

40 Fp=590 N
50 50
40 250
200 250
200
(N/μm)

10,000 rpm
Static Fp=590 N Static 10,000 rpm
(N/μm)

30 Fp=590 N
Stiffness

40 40
30 200
150 200
150
Stiffness

Fp=290 N Static Static

Stiffness

30
20 30
20 150
100 Fp=290 N 150
100
Stiffness

Fp=290 N
RadialRadial
AxialAxial

20
10 20
10 20,000 rpm 100
50 Fp=290 N 100
50
Fp=120 N Fp=120 N 20,000 rpm
100 0 500 500
0Fp=120
10N 20 30 100 20020,000
400rpm 600 0 F =120
p
10 N 20 30 0 20020,000
400
rpm 600
0 Rotating speed (krpm) 0 Preload force (N) 0 Rotating speed (krpm) 0 Preload force (N)
0 10 20 30 0 200 400 600 0 10 20 30 0 200 400 600
Rotating speed (krpm)
Figure 7. Preload
Figure force (N)ofofthe
influences
The influences therotating Rotating
rotatingspeed
speed
andandspeed
axial
axial (krpm)
preload
preload Preload
on bearing
on bearing force
stiffness
stiffness (N)
(7008C):
(7008C): (a)
(a) bear-
bearing
ing axialaxial stiffness
stiffness varies
varies withwith speed
speed and
and preload;(b)
preload; (b)bearing
bearingradial
radialstiffness
stiffness varies
varies with
with speed
speed
Figure
and 7. The influences of the rotating speed and axial preload on bearing stiffness (7008C): (a)
and preload.
preload.
bearing axial stiffness varies with speed and preload; (b) bearing radial stiffness varies with speed
and
3.3. preload.
Stiffness Variation with Different Load Conditions
3.3. Stiffness Variation with Different Load Conditions
As mentioned
As mentioned earlier,
earlier, the
theball–raceway
ball–raceway contact
contactanalysis
analysishas
hasbeen
beenadded
added inin
the
thebearing
bear-
3.3. Stiffness Variation with Different Load Conditions
modeling to realize the stiffness calculation of ball bearing under arbitrary load.
ing modeling to realize the stiffness calculation of ball bearing under arbitrary load. There- Therefore,
As mentioned
to show earlier,
the generality of thethe ball–raceway
proposed model,contact analysis has been
three typical added in ofthe
ballbear-
fore, to show the generality of the proposed model, three operation conditions
typical operation conditions bear-
of
ing modeling
ing SKF to
6205SKF realize
have been the stiffness
given calculation
to given
investigate of ball
the effectsbearing under arbitrary load. There-
ball bearing 6205 have been to investigate theofeffects
the rotating speeds, the
of the rotating external
speeds, the
fore,
load to show theand
conditions, generality of thecontact
the internal proposed model,
states three
of balls andtypical operation
raceways conditions
on bearing of
stiffness.
ball bearing SKF 6205 have been given to investigate the effects of the rotating speeds, the
3.3.1. The Pure Radial Load Condition
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.
 3.3.1. The Pure Radial Load Condition
3.3.1. The Pure Radial Load Condition
Figure 8 gives the stiffness change curves of ball bearing under one-direction radial
load.Figure 8 gives
One can find the
thatstiffness change
both static and curves of ball
dynamic bearing
stiffness underwith
increase one-direction
the radialradial
load.
load. One
Besides, can findfrom
different thatball
both static and
bearing acteddynamic stiffness
by the axial load increase with
(Preload), the the radialofload.
stiffness ball
Machines 2022, 10, 238 Besides, different from ball bearing acted stiffness of12ball
of 21
bearing under a certain radial load (500 N)by the axialwith
increases loadthe
(Preload),
rotating the
speed.
bearing under a certain radial load (500 N) increases with the rotating speed.
(a) (b)
(a) (b)

(N/μm)

(N/μm)
120 85.5

(N/μm)

(N/μm)
120 85.5
100 20,000 rpm 85

Stiffness

Stiffness
100 20,000 rpm 85

Stiffness

Stiffness
80 84.5
80 Static 84.5
60 Static 84
Radial 84 Fy = 500 N

Radial
60
Fy = 500 N
Radial

Radial
40 83.5
40100 400 700 1000 83.5 0 5 10 15 20
100 Radial
400 load 700(N) 1000 0Rotating
5 speed
10 (Krpm)
15 20
Radial load (N) Rotating speed (Krpm)
Figure 8. The radial stiffness variation for ball bearing one-direction pure radial load: (a) the radial
Figure
Figure8.8.
stiffness The
Theradial
radial
changes stiffness
stiffness
against variation
variation
the radial for
for
force; ball
ball
(b) bearing
the radial one-direction
bearing one-direction
stiffness pure
pure
changes radial load:
radialthe
against load:(a)
(a)the radial
thespeed.
rotating radial
stiffness changes against the radial force; (b) the radial stiffness changes against the rotating speed.
stiffness changes against the radial force; (b) the radial stiffness changes against the rotating speed.
3.3.2.
3.3.2.The
TheCombined
CombinedLargeLargeAxial
AxialLoad
Loadand andSmall
SmallRadial
RadialLoad
LoadCondition
Condition
3.3.2. The Combined Large Axial Load and Small Radial Load Condition
The above
Theabove section
abovesection has
sectionhas introduced
hasintroduced
introducedthe the change
thechange law
changelaw of
lawof ball
ballbearing
ofball bearingunder
underthe
thesingle
single
The
direction load. On this basis, the stiffness change curves bearing
ofofbearing acted under
by the the single
combined
direction
direction load.
load. OnOn this basis,
this basis, the stiffness
the stiffness change curves
change curves bearing
of bearing acted by the
acted Accordingcombined
by the combined
large
largeaxial
axialload
loadand
andsmall
smallradial
radialload
loadat at2000
2000rpm
rpmare
areshown
shownin inFigure
Figure9. Accordingto
9.According the
large
actual axial load
loading and small
condition radial
of ball load at
bearing 2000 rpm
under are shown
operation in Figure
condition, 9.
the angle totothe
displace-
the
actual
actual loading
loading condition
condition of ball
of to
ball bearing
bearing under operation condition, the angle displacement
ment
and θand 𝜃 θ and
are 𝜃
set are set zero. canunder
It seenbe operation
seen condition,
from 9Figure 9 that the axial
both angleand
displace-
radial
ment y and
and 𝜃 zand 𝜃 to zero.
are set toItzero.
can be from
It can be Figure
seen that both
from Figure axial
9 that bothand radial
axial andstiffness
radial
stiffness
decreasedecrease
with thewith theload.
radial radial load.
stiffness decrease with the radial load.
(a) (b)
(a) (b)
23.7 116
(N/μm)
(N/μm)

23.7 116
115.5
23.6
(N/μm)
(N/μm)

23.6 115.5
115
23.5
23.5 115
114.5
Stiffness
Stiffness

23.4 114.5
Stiffness

114
Stiffness

23.4
23.3 114
23.3 113.5
23.2 113.5
113
Radial

23.2
Axial

23.1 113
Radial

112.5
Axial

23.1 Fx = 200 N Fx = 200 N

23 Fx =100
200 N 112.5
112 Fx =100
200 150
N 200
23 0 50 150 200 112 0 50
0 50 100
Radial load150
(N) 200 0 Radial
50 100load 150
(N) 200
Radial load (N) Radial load (N)
Thestiffness
Figure9.9.The
Figure stiffnessvariation
variationfor
forball
ballbearing
bearingacted
actedbybythe
thecombined
combinedlarge
largeaxial
axialload
loadand
andsmall
small
Figure
radial 9. The
radialload:
load: stiffness
(a)
(a)axial variation
axialstiffness
stiffness for with
varies
varies ball
withbearing
radial acted
radialload;
load; by
(b)
(b) the combined
radial
radialstiffness large with
stiffnessvaries
varies axial load
radialand
withradial small
load.
load.
radial load: (a) axial stiffness varies with radial load; (b) radial stiffness varies with radial load.
3.3.3.The
3.3.3. TheCombined
CombinedSmall SmallAxial
AxialLoad
Loadand andLarge
LargeRadial
RadialLoad
LoadCondition
Condition
3.3.3. The Combined
Forcontrast,
contrast,the Small
thechange Axial
changecurves Load and
curvesofofbearing Large Radial
bearingstiffness Load
stiffnessacted
actedby Condition
bythe
thecombined
combinedsmall smallaxial
axial
For
load For
and contrast,
large the
radial change
load at curves
2000 of
rpm bearing
are shownstiffness
in acted
Figure by
10. the
One
load and large radial load at 2000 rpm are shown in Figure 10. One can find that the axialcombined
can find small
that theaxial
axial
load andreduces
stiffness
stiffness large
reducesradial loadwith
sharply
sharply at 2000
with therpm
the areload
radial
radial shown
load forain
for Figure
agiven
given 10. load
axial
axial One
loadcan find
(50N),
(50 thatthe
N),and
and theradial
the axial
radial
stiffness reduces
stiffnessfirstly
stiffness sharply
firstlydrops
dropsand with
andthen the
thenrises radial
riseswith
withtheload for
theradial a given
radialload. axial
load.Besides, load
Besides,both (50
boththe N), and
theradial the
radialand radial
andaxial
axial
stiffness
stiffness firstly drops
increase and
with then
the axialrises
loadwith
for the
ballradial
bearingload. Besides,
under a both
certain
stiffness increase with the axial load for ball bearing under a certain radial load (100 N). the radial
radial loadand axial
(100 N).
stiffness
Differentincrease
from thewith the
uniform axial load
increasing for ball
trend ofbearing under
axial stiffness a certain
with radial
the axial
Different from the uniform increasing trend of axial stiffness with the axial load, the load
load, the(100 N).
change
curve of the
Different fromradial stiffness contains two mainofpieces: the gradually increasing stage the
and
change curve ofthe
theuniform increasing
radial stiffness trend
contains twoaxial
mainstiffness
pieces: thewithgradually
the axialincreasing
load,
the rapidly
change curve increasing
of the stage.stiffness
radial Compared with the
contains twochange
main rule of bearing
pieces: the stiffnessincreasing
gradually in Figure 9,
stage and the rapidly increasing stage. Compared with the change rule of bearing stiffness
the stiffness
stage and the various
rapidly in Figure 10stage.
increasing presents the diversity
Compared andchange
with the complexity. ofIn
ruleand order to
bearing further
stiffness
in Figure 9, the stiffness various in Figure 10 presents the diversity complexity. In
in Figure 9, the stiffness various in Figure 10 presents the diversity and complexity.and
reveal internal formation mechanism, the corresponding contact state between balls In
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
 way
way contact
contact state
state change
change fromfrom full
full balls
balls contact
contact to to partial
partial balls
balls contact,
contact, thethe change
change trend
trend
of bearing stiffness also changes. The axial stiffness changes from
of bearing stiffness also changes. The axial stiffness changes from slow decline to rapid slow decline to rapid
decline
decline andand the
the radial
radial stiffness
stiffness changes
changes fromfrom decline
decline to to rising.
rising. In
In addition,
addition, as as seen
seen from
from
Figure
Figure 11b, for the given radial load, the number of balls located at the loaded zone grad-
11b, for the given radial load, the number of balls located at the loaded zone grad-
Machines 2022, 10, 238 ually
ually increases
increases withwith thethe axial
axial load.
load. Looking
Looking back back at at Figure
Figure 10,
10, it
it can
can bebe found
found thatthat the
the
13 of 21
increasing
increasing trend
trend of of ball
ball radial
radial stiffness
stiffness with
with axial
axial load
load changes
changes from
from slow
slow rising
rising toto rapid
rapid
rising
rising when
when thethe ball–raceway
ball–raceway contact
contact state
state changes
changes fromfrom partial
partial balls
balls contact
contact toto full
full balls
balls
contact.
contact. Above all, combining Figures 10 and 11, one can find that the changing trend of
Figures 10 and 11, one can find that the changing trend of the stiffness for ball bearing of
Above all, combining Figures 10 and 11, one can find that the changing trend in
the
the stiffness for ball bearing in full balls contact condition is the same as the results in
full stiffness for ball
balls contact bearingis in
condition thefull
sameballs
as contact
the resultscondition
in Figureis the same
9, and theas the results
changing trendin
Figure
Figure 9,
9, and
and the
the changing
changing trend
trend of
of bearing
bearing stiffness depends on
on the change of
of the internal
of bearing stiffness depends on the change stiffness depends
of the internal contact the change
state between theballs
internal
and
contact
contact state
state between
between balls
balls and
and raceways.
raceways. Therefore,
Therefore, one
one can
can conclude
conclude that
that the
the ball–race-
ball–race-
raceways. Therefore, one can conclude that the ball–raceway contact state has significant
way
way contact
contact
influence
state
on state
has
has significant
the variationsignificant influence
influence
of bearing
on
on the
stiffness. the variation
variation of of bearing
bearing stiffness.
stiffness.

(a)
(a) (b)
(b)
14
14 12
12 85
85 65
65
F
Fxx == 50
50 N
N Fyy == 100
F 100 N
N F
Fxx == 50
50 N
N F
Fyy == 100
100 N
N

(N/μm)
Stiffness(N/μm)
(N/μm)

80
Stiffness(N/μm)

12
12 10
10 80
75 60
60
10 88 75
10

RadialStiffness
70
70
AxialStiffness

88 66 55
55
65
65
66 44
60
60 50

Radial
50
Axial

44 22 55
55
22 00 50
500 100 200 300 400 500 45
45
00 100
100 200
200 300
300 400
400 500
500 00 10
10 20
20 30
30 40
40 50
50 0 100 200 300 400 500 00 10
10 20
20 30
30 40
40 50
50
Radial
Radial load
load (N)
(N) Axial
Axial load
load (N)
(N) Radial
Radial load
load (N)
(N) Axial
Axial load
load (N)
(N)
Figure
Figure 10. The stiffness variation for ball bearing acted by small axial load and large radial load: (a)
Figure 10.
10. The
Thestiffness
stiffnessvariation
variationfor
forball
ballbearing
bearingacted byby
acted small axial
small load
axial andand
load large radial
large load:
radial (a)
load:
The
The axial
axial stiffness
stiffness varies
varies against
against radial
radial or
or axial
axial load;
load; (b)
(b) The
The radial
radial stiffness
stiffness varies
varies against
against radial
radial or
or
(a) The axial stiffness varies against radial or axial load; (b) The radial stiffness varies against radial
axial
axial load.
load.
or axial load.
(a)
(a) (b)
(b)
zone

zone

10 10
contactzone

contactzone

10 10
99 99
ballsinincontact

ballsinincontact

88 88
77 77
66 66
ofballs

ofballs

55 55
numberof

numberof

44 44
The number

The number

33 33
22 22 F
Fxx == 50
F 50 N
N Fyy == 100
100 N
N
11 11
The

The

00 100
100 200
200 300300 400
400 500
500 00 10
10 20
20 30
30 40
40 50
50
Radial
Radial load
load (N)
(N) Axial
Axial load
load (N)
(N)
Figure 11.
Figure
Figure 11. The contact
11. The
The contact states
contact states between
states between balls
between balls and
balls and raceways
and raceways for
raceways for ball
for ball bearing
ball bearing acted
bearing acted by
acted by the
by the small
the small axial
small axial
axial
load
load and
and large
large radial
radial load:
load: (a)
(a)the
the contact
contactstates
statesbetween
between balls
balls and
and raceways
raceways varies
varieswith
with axial
load and large radial load: (a) the contact states between balls and raceways varies with axial force;axialforce;
force;
(b)
(b) the
the contact
contact states
states between
between balls
balls and
and raceways
raceways varies
varies with
with radial
radial force.
force.

3.4. Stiffness Change

3.4. Change with Internal
Internal Structural Parameters
Parameters
3.4. Stiffness
Stiffness Change with
with Internal Structural
Structural Parameters
In theabove
In above sections, the influence of operation conditions (speedloads)
and loads) on
In the
the above sections,
sections, thethe influence
influence ofof operation
operation conditions
conditions (speed
(speed and
and loads) onon bear-
bear-
bearing
ing stiffness has been described in detail. Due to the bearing internal structural
ing stiffness
stiffness has
has been
been described
described in in detail.
detail. Due
Due to
to the
the bearing
bearing internal
internal structural
structural parameter
parameter
parameter
also also having importance on bearinginstiffness, in the chapter,
followingthechapter, the high-
also having importance on bearing stiffness, in the following chapter, the high-precision
having importance on bearing stiffness, the following high-precision
precision 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) through changing internal structural parameters,
including ri (fi ), ro (fo ), D, di , and do :

Pd
cosα0 = 1 − , P = do − di − 2D A = r o + ri − D (27)
2A d
where Pd 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
 including ri (fi), ro (fo), D, di, and do:
𝑃
𝑐𝑜𝑠 𝛼 = 1 − ，𝑃 = 𝑑 − 𝑑 − 2𝐷 𝐴 = 𝑟 + 𝑟 − 𝐷 (27)
2𝐴

Machines 2022, 10, 238 where 𝑃 denotes the internal radial clearance. 14 of 21
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
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
(both ri (fi ) and ro (fo ) increase). Besides, the initial contact angle increases with the radial
clearance
clearance (di(d
decreases,
decreases, do increases, and
d increases, D decreases).
and D decreases).
i o
𝑟 𝑟r
𝑓f == ri，𝑓
, f =
=
o (28)
(28)
i 𝐷 o 𝐷D
D

(a) (b) (c)

24 18 20
Initial contact angle (deg)

22 Inner Inner
17 18
20
18 16 16
16
15 14
14 Outer Outer
12 14 12
0.51 0.54 0.57 −10 0 10 −10 0 10
Raceway curvature Raceway contact Ball diameter
coefficient diameter (μm) (μm)

Figure
Figure TheThe
12.12. influences
influences of of structural
structural parameters
parameters onon bearing
bearing initial
initial contact
contact angle:
angle: (a)(a) inner/outer
inner/outer
raceway
raceway curvature
curvature coefficients fi/o;f(b)
coefficients i/o ; (b) raceway
raceway contact
contact diameters di/o;d(c)
diameters i/o ; ball
(c) ball D. D.
diameter
diameter
3.4.1. Single Structural Parameter Variation
3.4.1. Single Structural Parameter Variation
In order to reveal the influence rules of single structural parameter variation on bearing
In order
stiffness to reveal
under the influence
different rules
speeds, the of singleangular
precision structural parameter
contact variation
ball bearing on bear-
7008C under
ingmedium
stiffnesspreload
under different speeds, the precision angular contact ball bearing 7008C un-
(290 N) has been chosen as the calculation example, and the key structural
Machines 2022, 10, x FOR PEER REVIEW 16 of 22
derparameters
medium preload (290 N) has been chosen as the calculation example, and the key
are shown in Table 1. On this basis, the detailed calculation results are given in
structural
Figures parameters
13–15. are shown in Table 1. On this basis, the detailed calculation results
are given in Figures 13–15.
(a) Figure 13 illustrates the influence(b) of groove curvature coefficient (radius) variation
70 on bearing stiffness under different speeds.
180 It can be seen from Figure 13a,b that, for ball
Radial stiffness(N/μm)

bearing at a static state or low speed, the axial stiffness decreases and radial stiffness re-
Axial stiffness(N/μm)

60 150
mains basically unchanged with the rise of the curvature coefficient (radius) of inner race-
50 120
way. Besides, when bearing under high-speed range, both the radial and axial stiffness
40
increase with the curvature coefficient90of inner raceway. Besides, it can be seen from Fig-
30 ure 13c,d that both the radial and axial60stiffness decrease with the curvature radius of the
20 outer raceway for the bearing at any speeds.30
10 0
0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.51 0.52 0.53 0.54 0.55 0.56 0.57
Inner raceway curvature coefficient fi Inner raceway curvature coefficient fi
(c) (d)
50 180
Radial stiffness(N/μm)
Axial stiffness(N/μm)

160
40
140
30 120
100
20
80
10 60
0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.51 0.52 0.53 0.54 0.55 0.56 0.57
Outer raceway curvature coefficient fo Outer raceway curvature coefficient fo

Figure 13. The

Figure 13. influences of raceway
The influences curvature
of raceway coefficient
curvature on bearing
coefficient stiffness:
on bearing (a) bearing
stiffness: axial axial
(a) bearing
stiffness
stiffness variesvaries
with with
innerinner
racewayraceway curvature
curvature coefficient fi ; (b)fibearing
coefficient ; (b) bearing
radialradial stiffness
stiffness variesvaries
with with
inner raceway curvature coefficient f ; (c) bearing axial stiffness varies with outer
inner raceway curvature coefficient fi ; (c) bearing axial stiffness varies with outer raceway curvature
i raceway curvature
coefficient fo ; (d)fobearing
coefficient ; (d) bearing
radialradial stiffness
stiffness varies varies
with with
outerouter raceway
raceway curvature
curvature coefficient
coefficient fo (—○—
fo (—#—
rotating speedspeed
rotating = 0 rpm,
=0— ∗——
rpm, rotating speed speed
— rotating = 5000=rpm,
5000— — —
rpm, □— rotating
rotating speed =10,000 rpm, —◇
rpm, —3—
speed =10000
rotating speed = 20,000 rpm).
— rotating speed = 20,000 rpm).

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,
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
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 race-
caused by the effect of the inertia forces of balls on the contact forces between outer race-
ways and balls.
ways and balls.
Machines 2022, 10, 238 15 of 21

(a) (b)
(a) (b)
40 200
40 200

stiffness(N/μm)
stiffness(N/μm)

175

Radialstiffness(N/μm)
Axialstiffness(N/μm)

35 175
35
150
30 150
30 125
25 125
25 100
100
20

Radial
75
Axial

20 75
15 50
15−10 −5 0 5 10 −10
50 −5 0 5 10
−10 −5 0 −10 −5
Inner raceway contact diameter d5i (μm) 10
Inner raceway 0
contact diameter di5(μm) 10
(c) Inner raceway contact diameter di (μm) Inner raceway contact diameter di (μm)
(d)
(c)40 (d)
200
40 200

stiffness(N/μm)
stiffness(N/μm)

Radialstiffness(N/μm)
35
Axialstiffness(N/μm)

35
150
30 150
30
25
25 100
100
20
Radial
Axial

20
15 50
−10
15 −5 0 5 10 −10
50 −5 0 5 10
−10 −10 −5
Outer −5 0 diameter d5 (μm)
raceway contact o
10 Outer raceway 0
contact diameter do5(μm) 10
Outer raceway contact diameter do (μm) Outer raceway contact diameter do (μm)
Figure 14. The influences of the contact diameter of inner or outer raceway on bearing stiffness: (a)
Figure
Figure14.14. The
The influences of thethe contact
contactdiameter
diameterofofinner
innerororouter
outerraceway
raceway onon bearing
bearing stiffness:
stiffness: (a)
bearing axial stiffness varies with the contact diameter of inner raceway di; (b) bearing radial stiff-
bearing
(a) bearingaxial stiffness
axial varies
stiffness with
varies thethe
with contact
contactdiameter of inner
diameter of innerraceway
racewaydi; (b)
di ; bearing radial
(b) bearing stiff-
radial
ness varies with the contact diameter of inner raceway di; (c) bearing axial stiffness varies with the
ness varies
stiffness with
varies thethe
with contact diameter
contact diameter of inner raceway di; (c) bearing
di ; (c) axial stiffness varies with the
contact diameter of outer raceway do; (d) of inner
bearing raceway
radial stiffness bearing
varies axial
with thestiffness
contactvaries with
diameter of
contact
the diameter
contact diameterof outer raceway d o; (d) bearing radial stiffness varies with the contact diameter of
do ; (d)
outer raceway do; (—of○outer raceway
— rotating speed = 0bearing
rpm, —radial stiffnessspeed
— rotating varies= with
5000 the
rpm, contact
—□—diameter
rotating
ofouter racewaydod;o(;—
outerraceway ○— rotating
(—#— rotatingspeed
speed==00rpm, rpm,— —∗— rotating speed
— rotating speed == 5000
5000 rpm,
rpm, — —□— —rotating
rotating
speed =10,000 rpm, —◇— rotating speed = 20,000 rpm).
speed =10,000 rpm, — ◇ — rotating speed
speed =10,000 rpm, —3— rotating speed = 20,000 rpm). = 20,000 rpm).

(a) (b)
(a) (b)
45 200
45 200
stiffness(N/μm)
stiffness(N/μm)

40
Radialstiffness(N/μm)
Axialstiffness(N/μm)

40
35 150
35 150
30
30
25 100
25 100
Radial

20
Axial

20
15 50
−10
15 −5 0 5 10 −10
50 −5 0 5 10
−10 −5 Ball diameter
0 D (μm) 5 10 −5 Ball diameter −10
0 D (μm) 5 10
Ball diameter D (μm) Ball diameter D (μm)
Figure15.
15. The influences
influencesof
ofball
balldiameter
diameteron bearing stiffness: (a) (a)
bearing axial stiffness varies with
Figure 15. The
Figure The influences of ball diameter ononbearing
bearing stiffness:
stiffness: bearing
(a) bearing axial
axial stiffness
stiffness varies
varies with
ball diameter
with D; (b)D;
ball diameter bearing radial stiffness
(b) bearing varies with ball diameter D (—○—Drotating
(—#— speed =0
ball diameter D; (b) bearing radial radial stiffness
stiffness varies
varies with with
ball ball diameter
diameter D (—○— rotating rotating
speed =0
rpm,  rotating speed = 5000 rpm, □ rotating speed =10,000 rpm, ◇ rotating speed =
speed —
= 0 —
rpm, — ∗ — rotating speed = —
5000 — rpm, —  — rotating speed —
=10,000— rpm, —3— ro-
rpm, —— rotating speed = 5000 rpm, —□— rotating speed =10,000 rpm, —◇— rotating speed =
tating
20,000speed
rpm).= 20,000 rpm).
20,000 rpm).
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 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 curvature
Machines 2022, 10, 238 16 of 21

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 coef-
ficient. 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.

3.4.2. 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
Machines 2022, 10, 238 17 of 21
of inner raceway curvature coefficient, the downtrend of bearing stiffness with rotating
speed tends to be gentle.

(a) (b)
Pd = 30μm Pd = 50μm Pd = 70μm Pd = 30μm Pd = 50μm Pd = 70μm
50 40 40 250 200 200

Radial Stiffness (N/μm)

Axial Stiffness (N/μm)

40 200
30 30 150 150
30 150
20 20 100 100
20 100
10 10 50 50
10 50
0 0 0 0 0 0
0.51 0.545 0.52 0.57 0.52 0.6 0.51 0.545 0.52 0.57 0.52 0.6
Raceway curvature Raceway curvature Raceway curvature Raceway curvature Raceway curvature Raceway curvature
coefficient fi coefficient fi coefficient fi coefficient fi coefficient fi coefficient fi

Figure
Figure 16.
16. The
Theinfluences
influences of of raceway
racewaycurvature
curvature coefficient
coefficient on
on bearing
bearing stiffness
stiffness with
with different
different radial
radial
clearance
clearance and
and constant
constant initial
initial contact
contact angle:
angle: (a)(a) bearing
bearing axial
axial stiffness
stiffness varies
varies with
with inner
inner raceway
raceway
curvature
curvature coefficient fi; (b)
coefficient fi; (b) bearing
bearing radial
radial stiffness
stiffnessvaries
varieswith
withinner
innerraceway
racewaycurvature
curvaturecoefficient
coefficientfi
fi(——
(------ rotating
rotating speed
speed = 0=rpm,
0 rpm,
—— —— rotating
rotating speed
speed = 20,000
= 20,000 rpm).
rpm).

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
Machines 2022, 10, x FOR PEER REVIEW
inner raceway curvature coefficient, the downtrend of bearing stiffness with19 of 22
rotating speed
tends to be gentle.

Figure
Figure TheThe
17.17. influences
influences of raceway
of raceway curvature
curvature coefficient
coefficient and speed
and rotating rotating speed on
on bearing bearing stiffness
stiffness
with constant initial contact angle and radial clearance: (a) the three-dimension curved of the rotat-
with constant initial contact angle and radial clearance: (a) the three-dimension curved of the rotating
ing speed, inner raceway curvature coefficient, and radial stiffness; (b) the two-dimensional projec-
speed, inner raceway curvature coefficient, and radial stiffness; (b) the two-dimensional projection
tion curves of three-dimensional surface (a).
curves of three-dimensional surface (a).
4. 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 calcula-
tion of ball bearing under arbitrary loads condition has been proposed. Through compar-
ison to the representative experimental and theoretical results in other literatures, the pro-
posed model has been verified. Based on these, the influences of various operation condi-
Machines 2022, 10, 238 18 of 21

4. 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 conditions (ro-
tating speed and loads) and structural parameters on bearing stiffness have been addressed.
Some novel conclusions have been given as follows:
1. Property increasing the axial preload for ball bearing can effectively inhibit the stiffness
attenuation phenomenon with rotating speed.
2. The internal contact states of balls and raceways have significant effects on the stiff-
ness variation.
3. The influences of initial contact angle variation determined by different structure
parameters on stiffness are different.
4. A smaller raceway curvature coefficient of outer ring and a larger raceway curvature
coefficient of inner ring can be designed for the high-speed ball bearing range to
improve its stiffness property.
This paper gives a detailed analysis of the effect factors, including the speeding speed
and structure parameters on the stiffness and stiffness variations of the ball bearing under
different conditions, and the conclusions drawn above can be used for ball bearing design
guidance and operation and maintenance.

Author Contributions: Formal analysis, funding acquisition, writing—original draft, Q.N.; formal
analysis, writing—review & editing Y.L.; conceptualization, writing—original draft, Y.Z.; data cu-
ration, software, S.P.; validation, Y.Y. and D.W. All authors have read and agreed to the published
version of the manuscript.
Funding: This research was funded by National key research and development program, grant
number: 2018YFB2000205.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: No date need be reported.
Conflicts of Interest: The authors declare no conflict of interest.

Nomenclature

A distance between curvature centers of raceway groove

aik , aok semi-major axes of the inner and outer ball–raceways contact
D rolling ball diameter
di , do contact diameter of inner/outer raceway
dm pitch diameter
E equivalent elastic modulus
Fck centrifugal force of the kth ball
Fik , F ok The second kind of elliptic integral of the inner and outer ball–raceways contact
fi , fo ri /D; ro /D
J mass moment of the ball
J Jacobi determinant
Jm Jacobian matrix of the ball local equations to ball local variables
K stiffness matrix
Kik , Kok load-deformation coefficients in the kth ball–inner/outer raceway contact
Lik , Lok the second kind elliptic integral of the inner and outer ball–raceways contact
m mass of the ball
Mgk gyroscopic moment of the kth ball
Pd internal radial clearance
Machines 2022, 10, 238 19 of 21

Qik , Qok contact loads of the kth ball–inner/outer raceway

ri , ro inner/outer raceway curvature radius
R xik , R xok equivalent curvature radius of the inner and outer ball–raceways contact
Tik , Tok frictions of the kth ball–inner/outer raceway
Z number of rolling balls
α0 initial contact angle
αik , αok contact angle of the kth ball
βk pitch angle of the kth ball
θy , θz The angular displacements of ball bearing
δik , δok elastic deformations of the kth ball–inner/outer raceway contacts
δx , δy , δz elative displacement of the inner and outer rings
λik , λok the distribution factor of the contact load of the inner and outer ball–raceways contact
γik , γok Dcosαik /dm ; Dcosαok /dm
ωi rotating speed of inner ring
ωbk spin angular speed of the kth ball
ωmk revolution angular speed of the kth ball
ψk angular position of the kth rolling ball

Appendix A. The Detailed Solution Process of Bearing Static Model

As described in the classical Hertz contact theory, the load-deflection coefficients Kik
and Kok can be expressed by using the least square linear regression algorithm:
!1/2 !1/2
Rik Eik Rok Eok
Kik = πEkik 3
, Kok = πEkok 3
(A1)
4.5Fik 4.5F ok

where E denotes the equivalent elastic modulus, and k, E, F can be expressed as follows:
0.6360 0.6360
Ryik Ryok
 
kik = 1.0339 , kok = 1.0339 (A2)
R xik R xok

0.5986 0.5986
Eik = 1.0003 + , Eok = 1.0003 + (A3)
Ryik /R xik Ryok /R xok
   
Fik = 1.5277 + 0.6023 ln Ryik /R xik , F ok = 1.5277 + 0.6023 ln Ryok /R xok (A4)

R xik Ryik R xok Ryok

Rik = , Rok = (A5)
R xik + Ryik R xok + Ryok
The effective contact radii of the ball–inner and ball–outer raceways are written as:

D (dm − Dcosαik ) ri D
R xik = , Ryik = (A6)
2dm 2ri − D

D (dm + Dcosαok ) ro D
R xok = , Ryok = (A7)
2dm 2ro − D

References
1. Stribeck, R. Ball bearing for various loads. Trans. ASME 1907, 29, 420–463.
2. Sjovall, H. The load distribution within ball and roller bearings under given external radial and axial load. Tek. Tidskr. 1933, 19, 72–75.
3. Lundberg, G.; Palmgren, A. Dynamic Capacity of Roller Bearings. Acta Polytechnica, Generalstabens Litograf. Anst. Förl. 1952, 2, 96–127.
4. Palmgren, A.; Ruley, B. Ball and Roller Bearing Engineering; SKF Industries, Inc.: Philadelphia, PA, USA, 1959.
5. Jones, A.B. A general theory for elastically constrained ball and radial roller bearings under arbitrary load and speed conditions.
J. Basic Eng. 1960, 82, 309–320. [CrossRef]
6. Harris, T.A. Rolling Bearing Analysis, 3rd ed.; John Wiley & Sons: Hoboken, NJ, USA, 1991; pp. 256–296.
7. Foord, C.A. High-speed ball bearing analysis. Proc. Inst. Mech. Eng. Part G J. Aerosp. Eng. 2006, 220, 537–544. [CrossRef]
8. Wang, W.-Z.; Hu, L.; Zhang, S.-G.; Zhao, Z.-Q.; Ai, S. Modeling angular contact ball bearing without raceway control hypothesis.
Mech. Mach. Theory 2014, 82, 154–172. [CrossRef]
Machines 2022, 10, 238 20 of 21

9. Ding, C.A.; Zhou, F.; Zhu, J.; Zhan, L. Raceway control assumption and the determination of rolling element attitude angle. Chin.
J. Mech. Eng. En. 2001, 37, 58–61. [CrossRef]
10. Walters, C.T. Dynamics of ball bearings. J. Lubr. Technol. 1971, 93, 1–10. [CrossRef]
11. Gupta, P.K. Dynamics of rolling-element bearings—Part III: Ball-bearing analysis. J. Lubr. Technol. Trans. ASME 1979, 101, 312–318.
[CrossRef]
12. Gupta, P.K. Dynamics of rolling-element bearings—Part IV: Ball-bearing results. J. Lubr. Technol. Trans. ASME 1979, 101, 319–326.
[CrossRef]
13. Meeks, C.R.; Ng, K.O. The dynamics of ball separators in ball bearings—Part I: Analysis. ASLE Trans. 1985, 28, 277–287. [CrossRef]
14. Meeks, C.R.; Long, T. Ball bearing dynamic analysis using computer methods—Part I: Analysis. J. Tribol. Trans. ASME 1996, 118,
52–58. [CrossRef]
15. Wang, Y.; Wang, W.; Zhang, S.; Zhao, Z. Investigation of skidding in angular contact ball bearings under high speed. Tribol. Int.
2015, 92, 404–417. [CrossRef]
16. Panda, S.; Nanda, P.; Mishra, D. Comparative study on optimum design of rolling element bearing. Tribol. Int. 2015, 92, 595–604.
[CrossRef]
17. Li, H.; Shin, Y.C. Analysis of bearing configuration effects on high speed spindles using an integrated dynamic thermo-mechanical
spindle model. Int. J. Mach. Tool. Manuf. 2004, 44, 347–364. [CrossRef]
18. Cao, Y.; Altintas, Y. A general method for the modelling of spindle-bearing systems. J. Mech. Des. Trans. ASME 2004, 126,
1089–1104. [CrossRef]
19. Cao, H.; Holkup, T.; Altintas, Y. A comparative study on the dynamics of high speed spindles with respect to different preload
mechanisms. Int. J. Adv. Manuf. Tech. 2011, 57, 871–883. [CrossRef]
20. Gargiulo, E.P. A simple way to estimate bearing stiffness. Mach. Des. 1980, 52, 107–110.
21. Wardle, F.P.; Lacey, S.J.; Poon, S.Y. Dynamic and static characteristics of a wide speed range machine tool spindle. Precis. Eng.
1983, 5, 175–183. [CrossRef]
22. Houpert, L. A uniform analytical approach for ball and roller bearings calculations. J. Tribol. 1997, 119, 851–858. [CrossRef]
23. Hernot, X.; Sartor, M.; Guillot, J. Calculation of the Stiffness Matrix of Angular Contact Ball Bearings by Using the Analytical
Approach. J. Mech. Des. Trans ASME 2000, 122, 83–90. [CrossRef]
24. Liu, J.; Tang, C.; Wu, H.; Xu, Z.; Wang, L. An analytical calculation method of the load distribution and stiffness of an angular
contact ball bearing. Mech. Mach. Theory 2019, 142, 103597. [CrossRef]
25. Lim, T.C.; Singh, R. Vibration transmission through rolling element bearings, part I: Bearing stiffness formulation. J. Sound Vib.
1990, 139, 179–199. [CrossRef]
26. Lim, T.C.; Singh, R. Vibration transmission through rolling element bearings, part II: System studies. J. Sound Vib. 1990, 139,
201–225. [CrossRef]
27. Lim, T.C.; Singh, R. Vibration transmission through rolling element bearings, Part III: Geared rotor system studies. J. Sound Vib.
1991, 151, 31–54. [CrossRef]
28. Bollinger, J.G.; Geiger, G. Analysis of the static and dynamic behavior of lathe spindles. Int. J. Mach. Tool. Des. Res. 1964, 3,
193–209. [CrossRef]
29. While, M.F. Rolling Element Bearing Vibration Transfer Characteristics: Effect of Stiffness. J. Appl. Mech. 1979, 46, 677–684.
[CrossRef]
30. Yang, Z.; Li, B.; Yu, T. Influence of structural parameters and tolerance on stiffness of high-speed ball bearings. Int. J. Precis. Eng.
Man. 2016, 17, 1493–1501. [CrossRef]
31. Yang, Z.; Chen, H.; Yu, T. Effects of rolling bearing configuration on stiffness of machine tool spindle. Proc. Inst. Mech. Eng. Part C
J. Mech. Eng. Sci. 2017, 232, 775–785. [CrossRef]
32. Noël, D.; Ritou, M.; Furet, B.; Le Loch, S. Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings.
J. Tribol. Trans ASME 2013, 135, 1–8. [CrossRef]
33. Fang, B.; Yan, K.; Hong, J.; Zhang, J. A comprehensive study on the off-diagonal coupling elements in the stiffness matrix of the
angular contact ball bearing and their influence on the dynamic characteristics of the rotor system. Mech. Mach. Theory 2021, 158,
104251. [CrossRef]
34. Sheng, X.; Li, B.; Wu, Z.; Li, H. Calculation of ball bearing speed-varying stiffness. Mech. Mach. Theory 2014, 81, 166–180.
[CrossRef]
35. Gunduz, A.; Singh, R. Stiffness matrix formulation for double row angular contact ball bearings: Analytical development and
validation. J. Sound Vib. 2013, 332, 5898–5916. [CrossRef]
36. Tsuha, N.; Cavalca, K.L. Stiffness and damping of elastohydrodynamic line contact applied to cylindrical roller bearing dynamic
model. J. Sound Vib. 2020, 481, 115444. [CrossRef]
37. Guo, Y.; Parker, R.G. Stiffness matrix calculation of rolling element bearings using a finite element/contact mechanics model.
Mech. Mac. Theory 2012, 51, 32–45. [CrossRef]
38. Zhang, X.; Han, Q.; Peng, Z.; Chu, F. Stability analysis of a rotor–bearing system with time-varying bearing stiffness due to finite
number of balls and unbalanced force. J. Sound Vib. 2013, 332, 6768–6784. [CrossRef]
39. Petersen, D.; Howard, C.; Prime, Z. Varying stiffness and load distributions in defective ball bearings: Analytical formulation and
application to defect size estimation. J. Sound Vib. 2015, 337, 284–300. [CrossRef]
Machines 2022, 10, 238 21 of 21

40. Petersen, D.; Howard, C.; Sawalhi, N.; Ahmadi, A.M.; Singh, S. Analysis of bearing stiffness variations, contact forces and
vibrations in radially loaded double row rolling element bearings with raceway defects. Mech. Syst. Signal. Processing 2015, 50–51,
139–160. [CrossRef]
41. Kraus, J.; Blech, J.J.; Braun, S.G. In Situ Determination of Rolling Bearing Stiffness and Damping by Modal Analysis. J. Vib. Acoust.
1987, 109, 235. [CrossRef]
42. Stone, B.J.; Walford, T. The measurement of the radial stiffness of rolling element bearings under oscillating conditions. J. Mech.
Eng. Sci. 1980, 22, 175–181.
43. Matsubara, M. Computational modelling of precision spindles supported by ball bearings. Int. J. Mach. Tools. Manuf. 1998, 28,
429–442. [CrossRef]
44. Jedrzejewski, J.; Kwasny, W. Modelling of angular contact ball bearings and axial displacements for high-speed spindles. CIRP
Ann. Manuf. Technol. 2010, 59, 377–382. [CrossRef]
45. Zhang, J.; Fang, B.; Hong, J.; Zhu, Y. Effect of preload on ball-raceway contact state and fatigue life of angular contact ball bearing.
Tribol. Int. 2017, 114, 365–372. [CrossRef]
46. Zhang, J.; Fang, B.; Zhu, Y.; Hong, J. A comparative study and stiffness analysis of angular contact ball bearings under different
preload mechanisms. Mech. Mach. Theory 2017, 115, 1–17. [CrossRef]
47. Fang, B.; Zhang, J.; Yan, K.; Hong, J.; Wang, M.Y. A comprehensive study on the speed-varying stiffness of ball bearing under
different load conditions. Mech. Mach. Theory 2019, 136, 1–13. [CrossRef]
48. Zhang, Y.; Fang, B.; Kong, L.; Li, Y. Effect of the ring misalignment on the service characteristics of ball bearing and rotor system.
Mech. Mach. Theory 2020, 151, 103889. [CrossRef]