# Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings

David Noel, Mathieu Ritou, Benoît Furet, Sébastien Le Loch

► To cite this version:

David Noel, Mathieu Ritou, Benoît Furet, Sébastien Le Loch. Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings. Journal of Tribology, 2013, 135 (4), pp.041101. 10.1115/1.4024109. hal-00862772

HAL Id: hal-00862772

https://hal.science/hal-00862772

Submitted on 8 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings

Angular contact ball bearings are predominantly used for guiding high speed rotors such as machining spindles. For an accurate modeling, dynamic effects have to be considered, most notably in the bearings model. The paper is based on a dynamic model of angular contact ball bearings. Different kinematic hypotheses are discussed. A new method is proposed for the computation of the stiffness matrix: a complete analytical expression including dynamic effects is presented in order to ensure accuracy at high shaft speed. It is demonstrated that the new method leads to the exact solution, contrary to the previous ones. Besides, the computational cost is similar. The new method is then used to investigate the consequence of the kinematic hypotheses on bearing stiffness values. Last, the relevance of this work is illustrated through the computation of the dynamic behavior of a high speed milling spindle. The impact of this new computation method on the accuracy of a finite element spindle model is quantified.

**Keywords:** angular contact ball bearing, high speed rotor, stiffness matrix, bearing kinematics, dynamics, spindle

David Noel

Mathieu Ritou

Benoit Furet

Institut de Recherche en Communications et en Cybernetique de Nantes (IRCCyN),
UMR CNRS 6597,
1 rue de la Noe,
44321 Nantes, France;
LUNAM Universite,
Universite de Nantes,
1 quai de Tourville,
44035 Nantes, France

**Sebastien Le Loch**
LUNAM Universite,
Universite de Nantes,
1 quai de Tourville,
44035 Nantes, France

## 1 Introduction

The dynamic behavior of a high speed rotor needs to be studied during its design. The dynamic stability is one of the most important criteria. Experimental tests can determine the instability, however, their major drawbacks are their cost, thus limiting the possible number of tests. With a numerical modeling approach, the stability can be analyzed in detail. As an example, numerical models contribute to choosing the optimum set of operating conditions in high speed machining (HSM) in order to maximize the productivity while ensuring workpiece quality and spindle health [1]. Experimental steps are still necessary to update and validate the model. Updating would increase the accuracy, thus providing the values of parameters such as damping [2].

One main part of the model relies on the guidance of the rotor. The paper deals with the most widespread type: the angular contact ball bearing guided shaft. Its model, associated with a finite element model, enables the building of a global model of the rotor [3,4]. In this process, the accuracy of the bearing stiffness values is crucial. Indeed, the frequency response function and, in particular, the critical eigenfrequencies are directly linked to the bearings stiffness [5].

Bearings manufacturers have developed hybrid bearings for high speed rotors. The balls are made of ceramic instead of steel, which decreases the dynamic effects on the balls and the subsequent load on the outer race at high speed due to centrifugal forces. Thus, these high precision components can reach more than $2.510^6 d_m N$ (ball orbital diameter in millimeters times the shaft speed in rpm), which is characteristic of very critical applications [1].

The mechanical model of the bearing aims at obtaining not only the relation between the global loads **f** and global displacement **d**, but also the stiffness values. The models are characterized by the number of degrees of freedom (DOF). Since no explicit constitutive relation can be obtained for multi-DOF models, local variables (normal loads ($Q_i, Q_o$) and normal contact deformations ($\delta_i, \delta_o$)) need to be considered for each ball. Historically, two different approaches are possible.

The analytical approach is based on a hypothesis on the local load distribution. Sjoväll [6] first proposed the use of the load distribution factor “$\hat{e}$” and the Sjoväll integrals. Palmgren [7] popularized this method with a 2-DOF model. Houpert [8] generalized the method with a 5-DOF model. Hernot et al. [9] computed the stiffness matrix. Unfortunately, the dynamics effects on the balls cannot be taken into account with this approach.

The numerical approach is based on a rigid body displacement hypothesis of the inner ring (see Fig. 1). The approach has been introduced by Jones [10] with a 5-DOF model. The main steps of the model’s elaboration are exposed in Refs. [4,10,11] and the stiffness matrix computation is detailed in Refs. [4,10,12]. With the numerical approach, the dynamic effects on the balls are taken into account during local and global equilibriums. Two types of dynamic effects are considered (see Fig. 2). The centrifugal force $F_c$ is induced by the ball orbiting at the diameter $d_m$. A gyroscopic moment $M_g$ results from the shifting direction of the ball’s spin axis.

Previous works proved that dynamic effects play a significant role in spindle behavior [2,4]: the dynamic effects on the shaft and on the balls of the bearings. Indeed the stiffness varies according to the shaft speed. Therefore, the numerical approach has been chosen in the paper to express the relation between global loads and displacements.

After presenting a 5-DOF bearing model and the different kinematic assumptions, the numerical computation of the bearing model is explained. A new complete analytical expression of the stiffness matrix, taking into account the dynamic effects on balls, is detailed. In a second section, the improvement of the accuracy and the

<figure><img src="image_2.png" /><figcaption>**Fig. 1 Principle of the bearing model**</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2 Dynamic effects on the balls</figcaption></figure>

related computational cost is then studied. The impact of the kinematic hypotheses is discussed. Finally, the relevance of the method is shown by studying the first modes of a high speed spindle.

## 2 Dynamic Model of Bearing

**2.1 Elaboration of the Model.** This subsection explains the main steps of the model elaboration schematized in Fig. 1.

**2.1.1 Rigid Body Displacement and Geometric Relations.** The geometry at the contact is shown in Fig. 3, corresponding to an azimuth location $\psi$ of a ball defined in Fig. 4. The relative displacement at the ball's location is established from a rigid body displacement hypothesis of the inner ring

$$ \begin{cases} A_1 = BD \sin \alpha + \delta_x - \theta_z \Re_i \cos \psi + \theta_y \Re_i \sin \psi \\ A_2 = BD \cos \alpha + \delta_y \cos \psi + \delta_z \sin \psi \end{cases} \quad (1) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3 Position of the ball center and raceway groove curvature centers at angular position $\psi$ with and without applied load</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Ball equilibrium</figcaption></figure>

According to Fig. 3, the Pythagorean theorem gives two geometric equations

$$ \begin{cases} (A_1 - X_1)^2 + (A_2 - X_2)^2 - [(f_i - 0.5)D + \delta_i]^2 = 0 \\ X_1^2 + X_2^2 - [(f_o - 0.5)D + \delta_o]^2 = 0 \end{cases} \quad (2) $$

In this formulation, the global deformations of the rings are not taken into account. Raceways are assumed circular and without radial expansion. Ring deformations can be computed with a finite element model or with a continuum mechanics model. The present method can also be refined with radial ring expansions [13].

**2.1.2 Local Equilibrium.** For each ball, Newton's second law of motion leads to (see Fig. 4)

$$ \begin{cases} Q_i \sin \alpha_i - Q_o \sin \alpha_o + \frac{M_g}{D} (\lambda_i \cos \alpha_i - \lambda_o \cos \alpha_o) = 0 \\ Q_i \cos \alpha_i - Q_o \cos \alpha_o - \frac{M_g}{D} (\lambda_i \sin \alpha_i - \lambda_o \sin \alpha_o) + F_c = 0 \end{cases} \quad (3) $$

Equation (3) and Fig. 4 suppose that contact surfaces can provide sufficient reaction forces to the gyroscopic moment $M_g$, i.e., $\lambda_o M_g / D \le \mu Q_o$ and $\lambda_i M_g / D \le \mu Q_i$ where $\mu$ is the friction coefficient at the contacts. This assumption is reasonable for spindles since a high preload is applied on the bearings. The coefficients $\lambda_i$ and $\lambda_o$ express the gyroscopic moment distribution (see Sec. 2.2.3). The local variables of Eqs. (2) and (3) are $\mathbf{x} = (X_1, X_2, \delta_o, \delta_i)$. The contact angles ($\alpha_o, \alpha_i$) and local loads ($Q_o, Q_i$) are expressed as a function of $\mathbf{x}$, respectively, with the geometric relations of Fig. 3 and with the Hertz constitutive relation [14]: $\delta = KQ^{2/3}$.

**2.1.3 Global Equilibrium.** Once the local equilibrium has been solved, all of the local loads and contact angles are known. The global loads $\mathbf{f}$ are obtained by summing each ball's quantities, according to Newton's second law applied on the inner ring

$$ \begin{cases} F_x = \sum_z Q_i \sin \alpha_i + \lambda_i \frac{M_g}{D} \cos \alpha_i \\ F_y = \sum_z \left( Q_i \cos \alpha_i - \lambda_i \frac{M_g}{D} \sin \alpha_i \right) \cos \psi \\ F_z = \sum_z \left( Q_i \cos \alpha_i - \lambda_i \frac{M_g}{D} \sin \alpha_i \right) \sin \psi \\ M_y = \sum_z \left[ \Re_i \left( Q_i \sin \alpha_i + \lambda_i \frac{M_g}{D} \cos \alpha_i \right) - \lambda_i f_i M_g \right] \sin \psi \\ M_z = \sum_z \left[ -\Re_i \left( Q_i \sin \alpha_i + \lambda_i \frac{M_g}{D} \cos \alpha_i \right) + \lambda_i f_i M_g \right] \cos \psi \end{cases} \quad (4) $$
**Geometrical relation.** A simple expression of the pitch angle $\beta$ is found geometrically by considering the spin axis of the ball as the bisecting line of the two normal directions to the contact surfaces

$$ \beta = \frac{\alpha_i + \alpha_o}{2} \qquad (8) $$

This hypothesis corresponds to the assumption that the spin-to-roll ratios at each contact area are equal

$$ \left( \frac{\omega_s}{\omega_{\text{roll}}} \right)_i = \left( \frac{\omega_s}{\omega_{\text{roll}}} \right)_o = \tan\left(\frac{\alpha_i - \alpha_o}{2}\right) \times \left[ 1 + \gamma \sin\left(\frac{\alpha_i + \alpha_o}{2}\right) \sin\left(\frac{\alpha_i - \alpha_o}{2}\right) \right] \qquad (9) $$

The hypothesis leads to simpler equations for $\omega_m$ and $\omega_R$. Even if the hypothesis seems, at first, limited compared to the other ones, it is retained in the paper for further investigation. In the following, this theory will be referred to as the “geometric” theory.

### 2.2.3 Contact Hypothesis: Distribution of the Gyroscopic Moment.
The contact hypothesis concerns the gyroscopic moment distribution between the inner and outer raceways. The distribution is taken into account through the coefficients $\lambda_i$ and $\lambda_o$ (see Fig. 4).

In previous works, tangential loads on the ball, corresponding to the reaction to the gyroscopic moment, are assumed equal on both raceways [4,10,17]. Others associate a distribution according to the kinematic hypothesis [11,15,16,18]. The associated distribution parameters $\lambda_i$ and $\lambda_o$ are given in Table 2.

In this work, tangential loads on balls are considered equal on both raceways $\lambda_i = \lambda_o = 1$ in order to enable comparisons of the stiffness computations with Ref. [4].

### 2.3 Model Solving Methods.
The model has been established in the previous subsection, depending on the application case. Either the load **f** is known and the displacement **d** is required and vice versa. Thus, two types of global solving methods are needed and exposed in this subsection. For both, local solutions, i.e., for each ball, are obtained first. It is numerically solved.

#### 2.3.1 Local Solution.
In this first step, the global displacement **d** is known. The distance between the groove curvature centers is computed with Eq. (1). Each ball equilibrium needs to be solved. The four equations (see Eqs. (2) and (3)) expressed with the variables **x** are simultaneously solved. Since there is no explicit solution for the multi-DOF model, a numerical algorithm is required. The Newton-Raphson method is commonly chosen due to its relative simplicity. Here, $\epsilon_n$ is the residue vector of the four equations at iteration number $n$

$$ \mathbf{x}_{n+1} = \mathbf{x}_n - \mathbf{A}_n^{-1} \epsilon_n \quad \text{with} \quad \mathbf{x} = (X_1, X_2, \delta_o, \delta_i) \qquad (10) $$

The matrix **A** is built from the partial derivatives of the four equations with respect to the local variables **x**. The matrix **A** is discussed in detail in Refs. [4,10].

It should be noted that the convergence of the algorithm depends on the starting point. In the present case, all functions are not defined on the real field. To avoid any problem, the solution of the static uniaxial case analytically obtained should be chosen as a starting point.

<table><caption>**Table 2 Gyroscopic moment distributions**</caption><thead><tr><th></th><th>Equal distribution</th><th>Outer-race control</th><th>Inner-race control</th><th>Hybrid theory</th></tr></thead><tbody><tr><td>λ<sub>i</sub></td><td>1</td><td>0</td><td>2</td><td>2C/(1+C)</td></tr><tr><td>λ<sub>o</sub></td><td>1</td><td>2</td><td>0</td><td>2/(1+C)</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 7 Global displacement determination**</figcaption></figure>

#### 2.3.2 Global Load Determination.
The load determination consists of computing the load **f** that needs to be applied on the inner ring to obtain a given displacement **d**. Once each ball equilibrium is solved (see Sec. 2.3.1), the local loads ($Q_i, Q_o$) are summed according to Eq. (4) in order to obtain **f**.

#### 2.3.3 Global Displacement Determination.
The displacement determination is the inverse of the load determination. It consists of computing the displacement **d** resulting from the application of a given load **f**. It is required to study the behavior of the spindle bearings. A new Newton-Raphson algorithm is proposed. It is based on the load determination method (see Fig. 7)

$$ \mathbf{d}_{n+1} = \mathbf{d}_n - \mathbf{K}_n^{-1}(\mathbf{f}_n - \mathbf{f}) \qquad (11) $$

where $\mathbf{K}_n$ is the stiffness matrix and $\mathbf{f}_n$ the result of the load determination at iteration number $n$. The implementation of the algorithm is simple once the stiffness matrix is made available. Its computational cost is discussed in Sec. 3.1.

### 2.4 Computation of the Stiffness Matrix.
The stiffness matrix **K** is needed for the method of global displacement determination. Moreover, it is required in order to investigate the dynamic behavior of rotors, as discussed in Sec. 3.4. The **K** matrix is a linearized relation between the loads $\Delta\mathbf{f}$ and displacements $\Delta\mathbf{d}$, for a given loading state. This is a 5 × 5 Jacobian matrix built from the partial derivatives of loads with respect to the displacements: $\mathbf{K} = [\partial\mathbf{f}/\partial\mathbf{d}]$.

#### 2.4.1 Numerical Computation by Finite Differences.
The computation by finite differences is the easiest method to implement because no analytic expression is needed. The stiffness matrix is calculated from a numerical gradient obtained after five additional load determinations (for the five DOF) close to the imposed displacement configuration. The method does not work for a displacement determination and is time-consuming (see Sec. 3.1). However, the numerical computation has the benefit of providing reference values for the validation of analytical computations of the stiffness matrix.

#### 2.4.2 Analytical Computation.
This paragraph presents a new method for the analytical computation of the stiffness matrix. The method is based on Refs. [4,10] and is extended to accurately take into account the dynamic effects. In this way, the complete analytical expression of the stiffness matrix is achieved. This is particularly suitable for the computation of critical bearing applications (with high $d_mN$), such as high speed spindles.

The stiffness matrix is $\mathbf{K} = [\partial\mathbf{f}/\partial\mathbf{d}]$. Since the load **f** results from the sum of each ball load (see Eq. (4)), the bearing stiffness results from the stiffness contribution of each ball. Thus, the matrix is built from a sum of terms corresponding to each ball. Therefore, each entry of **K** is expressed by a linear combination of elementary partial derivatives of the dynamic variables
x_d = (X_1, X_2, \delta_o, \delta_i, F_c, M_g) and the trivial variables (A_1, A_2) with respect to displacements d

$$K_{ij} = \frac{\partial f^i}{\partial d^j} = \sum_z \left[ \sum_k \left( \frac{\partial f^i}{\partial x_d^k} \frac{\partial x_d^k}{\partial d^j} \right) + \frac{\partial f^i}{\partial A_1} \frac{\partial A_1}{\partial d^j} + \frac{\partial f^i}{\partial A_2} \frac{\partial A_2}{\partial d^j} \right] \quad (12)$$

Here, i, j, k, and z are, respectively, the indices of the global loads, global displacements, dynamic local variables, and balls. The $\partial f^i / \partial x_d^k$, $\partial f^i / \partial A_1$, and $\partial f^i / \partial A_2$ terms are analytically obtained by differentiating the equations of the global equilibrium (see Eq. (4)). Trigonometric manipulations related to Fig. 3 are necessary. The $\partial A_1 / \partial d^j$ and $\partial A_2 / \partial d^j$ terms are easily obtained from Eq. (1).

The difficulty consists of obtaining the remaining terms $\partial x_d^k / \partial d^j$. For this purpose, the six other Eqs. (2), (3), and (5) are differentiated with respect to the global displacements $d^j$. For instance the equation of the centrifugal force $F_c$ gives

$$\frac{\partial F_c}{\partial d^j} = \frac{\partial F_c}{\partial X_1} \frac{\partial X_1}{\partial d^j} + \dots + \frac{\partial F_c}{\partial \delta_i} \frac{\partial \delta_i}{\partial d^j} + \frac{\partial F_c}{\partial A_1} \frac{\partial A_1}{\partial d^j} + \frac{\partial F_c}{\partial A_2} \frac{\partial A_2}{\partial d^j} \quad (13)$$

The six differentiated equations lead to a linear combination of $\partial x_d / \partial d^j$. This can be formulated as a matrix product $B \cdot \partial x_d / \partial d^j = s^j$

$$\left[ \begin{array}{ccccccc} & & & & & & \\ & B_{1:4,1:4} & & & & & \\ & & 0 & 0 & & & \\ & & 0 & 0 & & & \\ & & & -1 & B_{36} & & \\ & & & & 0 & B_{46} & \\ & \frac{\partial F_c}{\partial X_1} & \frac{\partial F_c}{\partial X_2} & \frac{\partial F_c}{\partial \delta_o} & \frac{\partial F_c}{\partial \delta_i} & -1 & 0 \\ & \frac{\partial M_g}{\partial X_1} & \frac{\partial M_g}{\partial X_2} & \frac{\partial M_g}{\partial \delta_o} & \frac{\partial M_g}{\partial \delta_i} & 0 & -1 \end{array} \right] \begin{pmatrix} \frac{\partial X_1}{\partial d^j} \\ \frac{\partial X_2}{\partial d^j} \\ \frac{\partial \delta_o}{\partial d^j} \\ \frac{\partial \delta_i}{\partial d^j} \\ \frac{\partial F_c}{\partial d^j} \\ \frac{\partial M_g}{\partial d^j} \end{pmatrix} = s^j \quad (14)$$

with

$$\left\{ \begin{aligned} B_{36} &= \frac{1}{D} \left( \frac{\lambda_o X_1}{(f_o - 0.5)D + \delta_o} - \frac{\lambda_i (A_1 - X_1)}{(f_i - 0.5)D + \delta_i} \right) \\ B_{46} &= \frac{1}{D} \left( \frac{-\lambda_o X_2}{(f_o - 0.5)D + \delta_o} - \frac{\lambda_i (A_2 - X_2)}{(f_i - 0.5)D + \delta_i} \right) \end{aligned} \right. \quad (15)$$

The upper left terms $B_{ij}$ with $(i, j)$ from 1 to 4 correspond to the matrix A evoked in Eq. (10) and detailed in Refs. [4,10]. The vector $s^j$ contains known terms concerning $A_1$ and $A_2$. For instance

$$s_5^j = - \frac{\partial F_c}{\partial A_1} \frac{\partial A_1}{\partial d^j} - \frac{\partial F_c}{\partial A_2} \frac{\partial A_2}{\partial d^j} \quad (16)$$

Here, **B** is a 6 × 6 matrix that is, in fact, independent of the DOF j. The elementary partial derivatives $\partial x_d / \partial d^j$ are simultaneously computed by inverting the matrix (see Eq. (14)). The numerical computation is carried out for each ball and for each DOF j.

For each kinematic hypothesis seen in Sec. 2.2.2, the analytical expressions corresponding to the dynamic effects are different (the two last rows of **B** and $s^j$). A symbolic software is useful in order to establish these expressions. For example, $B_{64}$ for the outer-race control is given in the Appendix.

This method differs from the literature. The most accomplished method [4] is equivalent to neglecting the terms $B_{35}$, $B_{36}$, and $B_{46}$. In that work, $\partial x_d^k / \partial d^j$ for $k \in \{1, 2, 3, 4\}$ are computed apart and $\partial x_d^k / \partial d^j$ for $k \in \{5, 6\}$ are computed afterwards. This simplification decreases the complexity of the implementation because $F_c$ and $M_g$ do not need to be analytically expressed in the function of

<table><caption>**Table 3 Solving methods**</caption><thead><tr><th></th><th>Load determination</th><th>Displacement determination</th><th>Analytical computation of K (previous method)</th><th>Analytical computation of K (new method)</th><th>Numerical computation of K</th></tr></thead><tbody><tr><td>1</td><td>✓</td><td></td><td>✓</td><td></td><td></td></tr><tr><td>2</td><td></td><td>✓</td><td>✓</td><td></td><td></td></tr><tr><td>3</td><td>✓</td><td></td><td></td><td>✓</td><td></td></tr><tr><td>4</td><td></td><td>✓</td><td></td><td>✓</td><td></td></tr><tr><td>5</td><td>✓</td><td></td><td></td><td></td><td>✓</td></tr></tbody></table>

<table><caption>**Table 4 Bearing characteristics**</caption><tr><td>α</td><td>25 deg</td><td>D</td><td>9.52 mm</td></tr><tr><td>d<sub>m</sub></td><td>90 mm</td><td>f<sub>i</sub></td><td>0.56</td></tr><tr><td>z</td><td>25</td><td>f<sub>o</sub></td><td>0.54</td></tr><tr><td>E<sub>b</sub></td><td>315 GPa</td><td>E<sub>r</sub></td><td>210 GPa</td></tr><tr><td>ν<sub>b</sub></td><td>0.26</td><td>ν<sub>r</sub></td><td>0.3</td></tr><tr><td>ρ<sub>b</sub></td><td>3190 kg/m<sup>3</sup></td><td></td><td></td></tr></table>

<table><caption>**Table 5 Computational costs depending on the solving method**</caption><tr><td>Method</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td></tr><tr><td>Cost (%)</td><td>100</td><td>601</td><td>105</td><td>634</td><td>573</td></tr></table>

$x = (X_1, X_2, \delta_o, \delta_i)$. The relevance of the new method, based on the complete analytical expression of the stiffness matrix, is investigated in the following section.

### 3 Results and Comparison

The five computational methods studied in this section are presented in Table 3. Each method will be referred by its number. All configurations are covered.

The following high precision bearing has been retained for the study: SNFA VEX70/NS9CE3. It corresponds to a hybrid angular contact ball bearing that can be found in high speed and high power spindles. Table 4 presents its characteristics.

**3.1 Computational Cost.** In this subsection, computation costs are compared for the methods presented in Table 3. The computational costs for each method are given in Table 5 in percentages, relative to the reference case number 1.

The displacement determination is more costly because of the additional loop of the algorithm (see Fig. 7). In the example, the algorithm expressed by Eq. (11) converges in six steps. The differences between the analytical methods are minor for both the load and displacement determination methods (respectively, methods (1) to (3) and (2) to (4)). In fact, only some additional elementary operations are added. The new simultaneous computation of the partial derivatives of the dynamic effects, presented in Sec. 2.4.2, does not significantly increase the computational cost.

The finite difference method 5 is more costly than the analytical computations 1 and 3. Indeed, six load determinations have to be carried out to compute the **K** matrix. The implementation of the analytical stiffness matrix is, therefore, beneficial because the computational cost is almost divided by 6 (the difference between methods 5 to 1) and also because it enables us to use the displacement determination computation.

**3.2 Stiffness Matrix Accuracy.** The accuracy of the stiffness matrix is studied through the axial stiffness $K_{11}$ evaluated in the uniaxial case and considering the outer-race control hypothesis. The results from the load and displacement determination are identical because the stop criterion “tol” of Fig. 7 is small enough.
<figure><img src="image_7.png" /><figcaption>Fig. 8 Relative error on the axial stiffness ΔK₁₁ considering the uni-axial case under different loads Fₓ and speed ω</figcaption></figure>

Figure 8 presents the relative error on the axial stiffness of the previous analytical method 1, compared with the new one 3

$$ \Delta K_{11} = 100 \left| \frac{K_{11}^{(1)} - K_{11}^{(3)}}{K_{11}^{(3)}} \right| \quad (17) $$

The relative error increases with the shaft speed ω and decreases with the axial load Fₓ. This maximum stiffness error corresponds to the configuration for which dynamics effects are significant compared to local loads. The error reaches 16% in the example. The stiffness error is significantly higher for steel ball bearings due to the ball density (the magnitude of dynamic effects is 2.5 higher).

The simplification adopted in previous work [4] was correct because it dealt with moderate dₘN applications. Nevertheless, here, the stiffness error is too important for accurately modeling the dynamic behavior of bearings in the case of critical applications (high dₘN). Indeed, Fig. 8 attests that the new stiffness computation is particularly interesting beyond 15,000 rpm.

In the example, the maximum error between methods 3 and 5 is less than 0.01%, regardless of the loads and speeds. Besides, the finite difference method is known to provide the reference results after a proper convergence study. Thus, we can conclude that the new analytical expression is complete because it leads to the exact results, contrarily to the previous ones.

### 3.3 Kinematic Hypotheses.
The kinematic hypothesis concerns the computation of the pitch angle β (see Sec. 2.2.2). Its impact on the stiffness is presented here. The axial stiffness K₁₁ is evaluated with the new analytical method for a medium preload of 1220 N on the bearing (see Fig. 9).

The loss of stiffness at high speed is due to the dynamic effects on the balls. As seen in Fig. 9, the stiffness drop is significant, which is characteristic of elastically preloaded bearings. Note that the kinematic hypothesis has a notable impact on the estimated stiffness at high speed. The estimated stiffness with the hybrid hypothesis is located between the values obtained with the inner and outer-race control, which is logical because the control theory supposes extreme situations concerning ball spinning motions. The hybrid theory is closer to the outer-race control at high speed, which is consistent with the boundaries presented in Fig. 5. The stiffness obtained with the geometric theory is almost in the middle of the values obtained with the inner and outer-race control because of the identical spin-to-roll ratios on both raceways.

In the example, the drop in stiffness in relation to the shaft speed is 75.7% and 78.3%, respectively, for the hybrid theory and outer-race control. The drops would be smaller for a more loaded

<figure><img src="image_8.png" /><figcaption>Fig. 9 Axial stiffness versus shaft speed in relation to the kinematic hypotheses</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 10 Pitch angle β versus shaft speed in relation to the kinematic hypotheses</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 11 Gyroscopic moment Mg versus shaft speed for different kinematic hypotheses</figcaption></figure>

bearing. The relative error for the axial stiffness K₁₁ between the outer-race control and hybrid theory is 10.4% at 30,000 rpm.

The main difference in the model comes from the estimation of the pitch angle β in Fig. 10. This angle barely changes the values of the speeds ωR and ωm or centrifugal forces Fc. However, the impact on the gyroscopic moment Mg is important (see Fig. 11).
<figure><img src="image_11.png" /><figcaption>Fig. 12 CAD model of the spindle rotor</figcaption></figure>

<table><caption>Table 6 Influence of the stiffness error on the eigenfrequencies values</caption><thead><tr><th>Kinematic<br>hypothesis</th><th>Bending<br>modes</th><th>0 rpm<br>(Hz)</th><th>30,000 rpm<br>method 1 (Hz)</th><th>30,000 rpm<br>method 3 (Hz)</th><th><sup>f<sup>(3)</sup></sup>&frasl;<sub>f<sup>(3)</sup></sub> (%)</th></tr></thead><tbody><tr><td rowspan="2">Outer-race control</td><td>Mode 1</td><td>1073</td><td>384</td><td>362</td><td>6.16</td></tr><tr><td>Mode 2</td><td>1296</td><td>486</td><td>459</td><td>5.74</td></tr><tr><td rowspan="2">Hybrid theory</td><td>Mode 1</td><td>1073</td><td>406</td><td>390</td><td>3.95</td></tr><tr><td>Mode 2</td><td>1296</td><td>512</td><td>493</td><td>3.94</td></tr></tbody></table>

At high speed, the difference between the hybrid and control theory reaches 30% for $M_g$. Hence, with a neglected gyroscopic moment, the impact of the kinematic hypothesis is quasi null.

Let us compare kinematic theories at high speed. Harris [11] determined that the outer-race control theory is valid at high speed, but only with light loads on the bearing. In that case, stiffness values obtained from the hybrid theory tend towards the outer-race theory values (even more than in Fig. 9). Indeed, $\omega_{so}$ then tends towards zero; therefore, the hybrid theory is suitable. In the same way, results from the hybrid theory tend towards the geometric theory ones with large loads, due to an increase of $\omega_{so}$, which is also observed in Harris's book. Thus, the hybrid theory is also adapted with large loads.

In conclusion, the hybrid theory seems to be suitable for modeling the bearing kinematic because the behavior varies in relation to the bearing load. However, for a simplified implementation of the analytical stiffness matrix, the geometric hypothesis could provide a good estimation of the stiffness values for a significantly loaded bearing.

### 3.4 Impact on Spindle Dynamic Behavior.
In Sec. 3.2, a study of the stiffness accuracy is presented. It highlights the major error made at high speeds if the analytical expressions of the stiffness are not properly derived. The present subsection establishes the impact of this error on the computation of the eigenfrequencies of the spindle rotor. The impact of the kinematic hypothesis choice is also investigated.

The rotor volumetric model is presented in Fig. 12. The bearings are taken into account as linearized springs applied to the center of the inner ring (corresponding to each bearing DOF). The spring stiffness values are the diagonal terms of the stiffness matrix **K**. They are computed with the new analytical method proposed in the paper. The off diagonal terms of **K** have been neglected for the purposes of this paper. The effect of the off diagonal terms is discussed in Ref. [19].

The four bearings are identical (their characteristics are shown in Table 4). They are in a back-to-back arrangement of tandems and elastically preloaded by springs. Thus, the axial load $F_x$ of each bearing corresponds to half of the preload and is considered constant. For the example, the bearing preload is 1220 N, corresponding to the stiff preload recommended by the bearing manufacturer. The bearing inner rings, spacers, and motor parts are modeled as cylinder elements with mass but without stiffness. For the ease of comparison, the dynamic effects on the shaft are not considered in this simplified model of the rotor.

Catia® finite element simulations have been carried out for the kinematic hypothesis of the outer-race control and hybrid theory. The simulations have also been made with results from methods 1

and 3 in order to observe the impact of the new computation method. The eigenfrequencies of the first two bending modes are presented in Table 6.

As for a 1-DOF mass spring-system, a drop of stiffness means lower eigenfrequencies, which explains the significant frequency evolution between 0 and 30,000 rpm.

The uncertainty on the eigenfrequencies resulting from the kinematic hypotheses (outer-race control or hybrid theory) is about 7% in this example. Additional work focusing on refining the bearing kinematics would probably enable us to reduce this uncertainty.

The new stiffness computation of method 3 avoids errors of, respectively, 6% and 4% for the outer-race control and the hybrid theory compared with the previous analytical method. The bending mode frequency errors vary from 16 Hz to 27 Hz between the different methods. This is crucial in order to predict stable conditions for high speed milling. The new complete analytical method is, therefore, of great interest in the elaboration of high speed spindle models.

## 4 Conclusion

After presenting a 5-DOF model for angular contact ball bearings, a new method was proposed to compute a stiffness matrix that properly takes into account the dynamic effects on the balls. Complete analytical expressions of stiffness matrix have been obtained. The numerical results have revealed an enhanced accuracy for the stiffness values. The new method is particularly relevant for lightly loaded bearings for critical applications (both the high values of shaft speed and ball orbital diameter). The computational cost has been analyzed for several solving methods and has not revealed any significant cost increase for the new stiffness computation method. Kinematic hypotheses concerning the ball spinning motion have been addressed. The hybrid theory seems more suitable because it takes into account loads using a friction ratio. The impact of the new method has been illustrated through a high speed spindle model. The eigenfrequencies have confirmed that the new method is crucial in order to accurately predict the dynamic behavior of the spindle.

## Nomenclature

### Capital Letters

$A_1, A_2$ = distances between groove curvature centers (loaded)

$A$ = Jacobian matrix for the computation of the local variables **x**

$B = f_i + f_o - 1$ gap between groove curvature centers (unloaded)
**B** = Jacobian matrix for the computation of the stiffness matrix

**D** = ball diameter

**E** = modulus of elasticity

**F**<sub>c</sub> = centrifugal force on ball

**J** = mass moment of inertia of ball

**K** = contact constant according to Hertz's theory

**K** = bearing stiffness matrix

**L** = complete elliptic integral of the second kind calculated in accordance to the Hertz theory

**M**<sub>g</sub> = ball gyroscopic moment

**Q** = ball-raceway normal load

**R** = radius to locus of raceway groove curvature center
**R**<sub>i</sub> = 0.5d<sub>m</sub> + (f<sub>i</sub> - 0.5)D cos α<sub>o</sub>

**X**<sub>1</sub>, **X**<sub>2</sub> = distances between the inner groove curvature center and the ball center (loaded)

### Lowercase Letters

**d** = global displacement of the inner ring; d<sup>l</sup> corresponds to δ<sub>z</sub> for j = 3

**d** = (δ<sub>x</sub>, δ<sub>y</sub>, δ<sub>z</sub>, θ<sub>x</sub>, θ<sub>y</sub>)

**d**<sub>m</sub> = ball orbital diameter

**f** = r/D

**f** = global load on inner ring (f and d are expressed at O<sub>h</sub>, the center of the outer ring; see Fig. 4).

**f** = (F<sub>x</sub>, F<sub>y</sub>, F<sub>z</sub>, M<sub>y</sub>, M<sub>z</sub>)

**m** = ball mass

**r** = raceway groove curvature radius

**x** = local variables

**x** = (X<sub>1</sub>, X<sub>2</sub>, δ<sub>o</sub>, δ<sub>i</sub>)

**x**<sub>d</sub> = local dynamic variables

**x**<sub>d</sub> = (X<sub>1</sub>, X<sub>2</sub>, δ<sub>o</sub>, δ<sub>i</sub>, F<sub>c</sub>, M<sub>g</sub>)

**z** = number of balls

### Greek Symbols

**α** = contact angles

**β** = ball pitch angle

**γ** = D/d<sub>m</sub>

**δ** = ball-raceway normal displacement

**λ** = distribution parameter for the gyroscopic moment

**μ** = friction coefficient at the contact

**ν** = Poisson's ratio

**ρ** = mass density

**ψ** = ball angular position on the pitch circle, azimuth position

**ω** = shaft speed

**ω**<sub>m</sub> = orbital speed of ball

**ω**<sub>R</sub> = speed of ball about its own axis

**ω**<sub>roll</sub> = rolling motion

**ω**<sub>s</sub> = spinning motion

### Subscripts

**i** = inner ring

**o** = outer ring

## Appendix: Example of Analytical Expression in B Matrix

Considering the outer-race control hypothesis

$$B_{64} = \frac{\partial M_g}{\partial \delta_i} = \frac{J \omega^2 \tan \beta}{(f_i - 0.5)D + \delta_i} \times \left[ \gamma \left(\frac{\omega_m}{\omega}\right) \left(\frac{\omega_R}{\omega}\right)^2 \left( \frac{(X_1 - A_1) + \tan \beta(A_2 - X_2)) \sin \alpha_o}{\sin \alpha_i (1 + \gamma \cos \alpha_o)} - \frac{(1 + \gamma \tan \beta \sin \alpha_i)(A_2 - X_2)}{(1 - \gamma \cos \alpha_i)^2} + \frac{(2(X_1 - A_1) + \tan \beta(A_2 - X_2)) \tan \beta}{1 - \gamma \cos \alpha_i} \right) + \left( \frac{\omega_R}{\omega} \right) \left( (1 + \gamma \cos \alpha_o) \frac{(X_2 - A_2) \cos \alpha_o + (X_1 - A_1) \sin \alpha_o}{(1 + \cos(\alpha_o - \alpha_i))^2} + \frac{(\omega_m)(X_1 - A_1) + (A_2 - X_2) \tan \beta}{\sin \alpha_i} \right) \right]$$

## References

[1] Abele, E., Altintas, Y., and Brecher, C., 2010, “Machine Tool Spindle Units,” CIRP Ann., **59**(2), pp. 781–802.

[2] Gagnol, V., Bouzgarrou, B., Ray, P., and Barra, C., 2007, “Model-Based Chat- ter Stability Prediction for High-Speed Spindles,” Int. J. Mach. Tools Manuf., **47**(7–8), pp. 1176–1186.

[3] Nelson, H., 1980, “A Finite Rotating Shaft Element Using Timoshenko Beam Theory,” ASME J. Mech. Des., **102**(4), pp. 793–803.

[4] Cao, Y., and Altintas, Y., 2004, “A General Method for the Modeling of Spindle-Bearing Systems,” ASME J. Mech. Des., **126**(6), pp. 1089–1104.

[5] Swanson, E., Powell, C., and Weissman, S., 2005, “A Practical Review of Rotating Machinery Critical Speeds and Modes,” Sound Vib., **39**(5), pp. 10–17.

[6] Sjoväll, H., 1933, “The Load Distribution Within Ball and Roller Bearings Under Given External Radial and Axial Load,” Tek. Tidsk., pp. 97–102.

[7] Palmgren, A., 1959, *Ball and Roller Bearing Engineering*, 3rd ed., SKF industries, Philadelphia.

[8] Houpert, L., 1997, “A Uniform Analytical Approach for Ball and Roller Bearings Calculations,” ASME J. Tribol., **119**(4), pp. 851–859.

[9] Hernot, X., Sartor, M., and Guillot, J., 2000, “Calculation of the Stiffness Matrix of Angular Contact Ball Bearings by Using the Analytical Approach,” ASME J. Mech. Des., **122**(1), pp. 83–90.

[10] Jones, A., 1960, “A General Theory for Elastically Constrained Ball and Radial Roller Bearings Under Arbitrary Load and Speed Conditions,” ASME J. Basic Eng., **82**(2), pp. 309–320.

[11] Harris, T. A., and Kotzalas, M., 2007, *Advanced Concept of Bearing Technology*, 5th ed., CRC/Taylor & Francis, London.

[12] de Mul, J. M., Vree, J. M., and Maas, D. A., 1989, “Equilibrium and Associated Load Distribution in Ball and Roller Bearings Loaded in Five Degrees of Freedom While Neglecting Friction—Part I: General Theory and Application to Ball Bearings,” ASME J. Tribol., **111**(1), pp. 142–148.

[13] Noel, D., Ritou, M., Le Loch, S., and Furet, B., 2012, “Bearings Influence on the Dynamic Behavior of HSM Spindle,” 11th Biennial Conference on Engineering Systems Design and Analysis (ESDA), ASME.

[14] Hertz, H., 1881, “On the Contact of Elastic Solids,” J. Reine Angew. Math, **92**, pp. 156–171.

[15] Changan, D., Fuzhang, Z., Jun, Z., and Lei, Z., 2000, “Raceway Control Assumption and the Determination of Rolling Element Attitude Angle,” Proceedings of the Symposium of Mechanical Engineering, IEEE, New York, pp. 158–163.

[16] Lei, C., Rui, Z., Liu, J., Feng, R., and Zhao, J., 2010, “A New Method for Computing Contact Angle of High Speed Ball Bearing,” 2010 Third International Joint Conference on Computational Science and Optimization, Vol. 1, IEEE, New York, pp. 331–334.

[17] Jedrzejewski, J., and Kwasny, W., 2010, “Modelling of Angular Contact Ball Bearings and Axial Displacements for High-Speed Spindles,” CIRP Ann., **59**(1), pp. 377–382.

[18] Antoine, J. F., Abba, G., and Molinari, A., 2006, “A New Proposal for Explicit Angle Calculation in Angular Contact Ball Bearing,” ASME J. Mech. Des., **128**(2), pp. 468–478.

[19] Li, H., and Shin, Y., 2004, “Analysis of Bearing Configuration Effects on High Speed Spindles Using an Integrated Dynamic Thermo-Mechanical Spindle Model,” Int. J. Mach. Tools Manuf., **44**(4), pp. 347–364.