Thermally affected stiffness matrix of angular contact ball bearings in a high-speed spindle system

Abstract

Bearing stiffness directly affects the dynamic characteristics in a high-speed spindle system and plays an important role in terms of manufacturing quality. We developed a new approach for predicting the thermal behavior of a high-speed spindle, calculated the thermal expansion, and generated a bearing stiffness matrix for angular contact ball bearings. The heat convection of spindle housing in air, the balls in lubricant, the spindle shaft in quiescent air, and the bearing inner ring surfaces were determined. Heat sources such as bearing friction, and the heat contributed by the built-in motor, were simulated using an analysis systems (ANSYS) steady-state thermal model. The results were imported into a static ANSYS structural model. Ball thermal expansion was calculated based on changes in the coordinates of nodal points on the ball surface. Finally, a thermally affected bearing stiffness matrix was generated by applying the Newton–Raphson technique. Decreases in the bearing radial, axial, angular, and coupling stiffness values as rotational spindle speed increased were calculated. Also, the stiffness coefficients at a specific speed increased significantly caused by the thermal effects. Finally, for validation, the bearing stiffness was compared to that calculated using an earlier thermal network approach.

Introduction

Machining requires high-level precision. The dynamic spindle characteristics of a machining tool such as speed, torque, dynamic stiffness, and thermal expansion affect the throughput and quality of machining operations. It is difficult to eliminate heat generated by friction between the contacts of machine tool spindles. ¹ Therefore, for the computation of the spindle bearing dynamic stiffness, it is essential to analytically model thermally affected parameters.

Palmgren ² was the first to explore the relationship between bearing load and the corresponding displacement. Jones ³ determined inner and outer ring displacements according to the position of balls operating within a five degrees of freedom (DOF) support. The 5 DOF bearing loads were then developed using the model of De Mul et al., ⁴ who used a Jacobian matrix to solve the nonlinear contact force equations of bearings and defined ball contact stiffness using a five-by-five matrix. More recently, a 5 DOF, quasi-static Jones bearing model with different bearing configurations was used to study the stiffness characteristics of duplex angular contact ball bearings. ⁵ To obtain the five-by-five dynamic stiffness matrix, a finite element model ⁶ was constructed using LS-DYNA software to calculate bearing displacement and stress. However, thermal effects were not considered in these previous works, although heat generation within bearings may trigger failure and compromise the reliability, performance, and rotational speed of a machine tool spindle system.^7,8

Burton and Staph ⁷ performed one of the first studies on the thermal behavior of angular contact ball bearings. To obtain the amount of thermal expansion in the angular contact ball bearings, Burton and Staph drew a heat-transfer diagram of the thermal resistance associated with deformation caused by elastic contacts. Stein and Tu ⁹ developed a state-space model based on the heat transmission laws governing different bearing configurations under various lubrication conditions. Heat sources, heat transfers, and sinks were included in the closed-loop, finite-difference thermal spindle model of Bossmanns and Tu. ¹⁰ Takabi and Khonsari ¹¹ studied thermal behavior using a thermal resistance network to estimate the thermal expansion of bearing components. Recent studies^12,13 created comprehensive thermal grid models to explore the thermal performance of angular contact ball bearings. Ye et al. ¹⁴ analyzed the thermally dynamic characteristics of angular contact ball bearings. Overall, although the various models introduced above used imaginative approaches to investigate thermal-mechanical behavior, the effects of heat on stiffness were not considered.

Jorgensen and Shin ¹⁵ further developed the model of De Mul et al. by adding a heat-transfer network to the spindle bearing system. Jorgensen and Shin combined the centrifugal expansion and inner ring-shaft thermal expansion into the ball and groove center separation in order to compute the ball bearing stiffness matrix. In addition, a finite element approach was taken, assuming that the spindle geometry was axis-symmetric. ¹⁶ Assessments of thermal behavior based on thermal networks have certain limitations; the thermal resistance network and associated calculations are not precise because the geometric factors in thermal resistance formula change with thermal expansion, especially for a multi-bearing motorized spindle system with complicated shapes.

In this article, we describe a workflow using an ANSYS steady-state model. This model simplifies observations of heat transfer within a spindle system. The thermal conditions of the model are determined based on the volume or surface area of the parts and include coefficients reflecting internal heat generation and convection. We used a finite element approach to solve the complicated multi-bearing spindle system. Thus, the drawbacks of thermal resistance networks are overcome. The model predicts the thermal behavior of an experimental spindle system and provides data on the thermal expansion of balls, the inner ring-shaft, and the outer ring-housing. These parameters can then be used to calculate bearing stiffness.

Bearing stiffness matrix

Centrifugal expansion

As spindle rotational velocity increases, the orbital speed of all balls in angular contact ball bearings also increases. The centrifugal force on each ball is in accordance with Harris ¹⁷

$$F_c=\frac{1}{2}m{d}_m{\omega }^2{\left(\frac{{\omega }_{\mathrm{obt}}}{\omega }\right)}_j^2$$

(1)

where j denotes the jth bearing ball at a certain attitudinal angle and ${\omega }_{\mathrm{obt}}$ is the ball orbital velocity. Using the equation of Timoshenko and Goodier, the expansion of the bearing inner ring and spindle shaft caused by centrifugal force is given by Jorgensen and Shin ¹⁵

$${\epsilon }_c=\frac{\rho {\omega }^2}{32E}{d}_m[D_i^2(3+\upsilon )+d_m^2(1-\upsilon )]$$

(2)

Ball loads and displacement

Figure 1 shows the loads acting on the ball, including the centrifugal force $F_c$ , the bearing inner ring contact force $Q_i$ , and the bearing outer ring contact force $Q_e$ . Displacement of the ball and bearing rings is governed by these loads and thermal expansion. Figure 2 ¹⁵ shows the geometry of the ball and the bearing inner ring displacement caused by the loads and thermal expansion. On the contact line, the distances between the ball center and the bearing inner and outer ring groove centers are $l_{\mathrm{oi}}$ and $l_{\mathrm{oe}}$ , respectively. The motion of the ball center is $(v_r,v_z)$ and the motion of the inner ring groove center is $(u_r,u_z)$ . The contact deformations between the ball and the raceways are given by

$${\delta }_i=l_i-l_{\mathrm{oi}}-{\sigma }_i,{\delta }_e=l_e-l_{\mathrm{oe}}-{\sigma }_e$$

(3)

where ${\sigma }_i$ and ${\sigma }_e$ are the initial clearances between the ball and the bearing inner and outer rings, respectively. $l_i$ and $l_e$ represent the separation between the ball and the inner and outer groove centers, respectively, calculated based on the ball thermal expansion ${\psi }_b$ , the centrifugal expansion ${\epsilon }_c$ , and the initial geometric parameters

$$\begin{array}{c}{l}_i={\psi }_b+\hfill \\ {[{(l_{\mathrm{oi}}\cos {\alpha }_o+u_r-v_r+{\psi }_{\mathrm{ir}}+{\epsilon }_c)}^2+{(l_{\mathrm{oi}}\sin {\alpha }_o+u_z-v_z)}^2]}^{1/2}\hfill \\ l_e={\psi }_b+{[{(l_{\mathrm{oe}}\cos {\alpha }_o+v_r-{\psi }_{\mathrm{er}})}^2+{(l_{\mathrm{oe}}\sin {\alpha }_o+v_z)}^2]}^{1/2}\hfill \end{array}$$

(4)

where ${\psi }_b$ is the difference between the initial and final ball radius under static thermal conditions, here defined by fitting the nodal point coordinates to the ANSYS static structural model. The inner ring-spindle thermal expansion is in accordance with Timoshenko and Goodier ¹⁸

$${\psi }_{\mathrm{ir}}=\frac{{\varsigma }_s}{3}\mathrm{\Delta }{T}_{\mathrm{sm}}(1+{\upsilon }_s)r_s+\frac{{\varsigma }_i}{3}\mathrm{\Delta }{T}_i(1+{\upsilon }_i)r_i$$

(5)

where ${\varsigma }_s$ and ${\varsigma }_i$ are the thermal expansion coefficients (m/m/°C) of the materials from which the spindle shaft and the inner bearing ring, respectively, are fabricated. $\mathrm{\Delta }{T}_{\mathrm{sm}}$ and $\mathrm{\Delta }{T}_i$ are the temperature changes of the spindle shaft where the built-in motor located and the bearing inner ring, respectively. $r_s=1/2D_i$ is the radius of the spindle shaft at the point where the bearing is located.

Figure 1.

Loads on a ball.

Figure 2.

Thermal expansion and displacement.

The outer ring-housing thermal expansion is given by

$${\psi }_{\mathrm{er}}=\frac{{\varsigma }_e}{3}\frac{(1+{\upsilon }_e)r_e}{r_h+r_e}[\mathrm{\Delta }{T}_e(2r_e+r_h)+\mathrm{\Delta }{T}_h(2r_h+r_e)]$$

(6)

where $r_h$ is the radius of the housing in which the bearing is located. $\mathrm{\Delta }{T}_e$ is the temperature change of the bearing outer ring, and $\mathrm{\Delta }{T}_h$ is the change in housing temperature.

The Hertzian theory of spherical contact states that if the bearing ring contact forces $Q_i$ and $Q_e$ are functions of bearing deflection, then^15,17,19

$$Q_i=C_i{\delta }_i^{1.5},Q_e=C_e{\delta }_e^{1.5}$$

(7)

where $C_i$ and $C_e$ are the contact stiffness coefficients according to the material properties and bearing geometry, respectively. Harris ¹⁷ and Sopanen and Mikkola ²⁰ expressed $C_i$ and $C_e$ using elliptical integral parameters

$$C_i=\pi {\Im }_{i}E\prime \sqrt{\frac{R_i{\xi }_i}{4.5{\zeta }_i^3},}{C}_e=\pi {\Im }_{e}E\prime \sqrt{\frac{R_e{\xi }_e}{4.5{\zeta }_e^3}}$$

(8)

where $\Im$ , $\xi$ , and $\zeta$ are an ellipticity parameter and the first and second elliptic integrals, respectively, defined using the least squares method

$$\begin{array}{c}\Im =1.0339{\left(\frac{{\Re }_y}{{\Re }_x}\right)}^{0.636},\xi =1.0003+0.5968\left(\frac{{\Re }_x}{{\Re }_y}\right),\hfill \\ \zeta =1.5277+0.6023\ln \left(\frac{{\Re }_y}{{\Re }_x}\right)\hfill \end{array}$$

(9)

Parameters ${\Re }_x$ and ${\Re }_y$ are the effective radii of curvature and calculated as functions of the radii of curvature for ball-inner, ball-outer race contact.

$E\prime$ is obtained from the modulus of elasticity and the Poisson’s ratios of the balls and the inner and outer rings of the angular contact ball bearing

$$\frac{1}{E\prime }=\frac{1}{2}(\frac{1-{\upsilon }_a^2}{E_a}+\frac{1-{\upsilon }_b^2}{E_b})$$

(10)

where the subscripts a and b refer to the contact pair. In angular contact ball bearings, a contact pair refers to the contact between the ball and either the bearing inner or outer groove.

Local equilibrium equations

From Figure 1, the loads on the ball are given by

$$\left\{F_b\right\}=\left\{\begin{array}{c}{F}_r\hfill \\ \multicolumn{1}{c}{F_z}\end{array}\right\}=\left\{\begin{array}{c}{Q}_i\cos {\alpha }_i-Q_e\cos {\alpha }_e+F\hfill \\ \multicolumn{1}{c}{Q_i\sin {\alpha }_i-Q_e\sin {\alpha }_e}\end{array}\right\}=\left\{\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right\}$$

(11)

The Newton–Raphson technique is used to solve the nonlinear equilibrium equation (11) by considering changes $(\mathrm{\Delta }{v}_r,\mathrm{\Delta }{v}_z)$ in ball center coordinates. The changes $(\mathrm{\Delta }{v}_r,\mathrm{\Delta }{v}_z)$ after each iteration are obtained using local equilibrium equations

$$\left\{F_b\right\}+\left[J_b\right]·\mathrm{\Delta }v=\left\{\begin{array}{c}{F}_r\hfill \\ \multicolumn{1}{c}{F_z}\end{array}\right\}+\left[\begin{array}{cc}\frac{\partial F_r}{\partial v_r}\hfill & \frac{\partial F_r}{\partial v_z}\hfill \\ \multicolumn{1}{c}{\frac{\partial F_z}{\partial v_r}}& \frac{\partial F_z}{\partial v_z}\hfill \end{array}\right]·\left\{\begin{array}{c}\mathrm{\Delta }{v}_r\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }{v}_z}\end{array}\right\}=0$$

(12)

The Jacobian matrix of the local equilibrium equations is the derivative of the ball loads of equation (11) as the ball moves

$$\left[J_b\right]=\left[\frac{\partial \left\{F_b\right\}}{{\partial \left\{v\right\}}^T}\right]=\left[\begin{array}{cc}\frac{\partial F_r}{\partial v_r}\hfill & \frac{\partial F_r}{\partial v_z}\hfill \\ \multicolumn{1}{c}{\frac{\partial F_z}{\partial v_r}}& \frac{\partial F_z}{\partial v_z}\hfill \end{array}\right]$$

(13)

By deriving partial derivatives and applying the chain rule, the elements of the Jacobian matrix of equation (13) are functions of the contact stiffness coefficients, bearing deformations, inner and outer contact angles, and thermal expansion

$$\begin{array}{cc}\multicolumn{1}{c}{\frac{\partial F_r}{\partial v_r}=}& -\frac{3}{2}{C}_i{\delta }_i^{1/2}\frac{l_i{\cos }^2{\alpha }_i}{l_i-{\psi }_b{\sin }^2{\alpha }_i}-Q_i\frac{{\sin }^2{\alpha }_i}{l_i-{\psi }_b{\sin }^2{\alpha }_i}\hfill \\ \multicolumn{1}{c}{}& -\frac{3}{2}{C}_e{\delta }_e^{1/2}{\cos }^2{\alpha }_e-\frac{Q_e}{l_e-{\psi }_b}{\sin }^2{\alpha }_e\hfill \end{array}$$

(14)

$$\begin{array}{cc}\multicolumn{1}{c}{\frac{\partial F_z}{\partial v_z}=}& -\frac{3}{2}{C}_i{\delta }_i^{1/2}{\sin }^2{\alpha }_i-\frac{Q_i}{l_i-{\psi }_b}{\cos }^2{\alpha }_i\hfill \\ \multicolumn{1}{c}{}& -\frac{3}{2}{C}_e{\delta }_e^{1/2}{\sin }^2{\alpha }_e-\frac{Q_e}{l_e-{\psi }_b}{\cos }^2{\alpha }_e\hfill \end{array}$$

(15)

$$\begin{array}{cc}\multicolumn{1}{c}{\frac{\partial F_r}{\partial v_z}}& =\frac{\partial F_z}{\partial v_r}=-\frac{3}{2}{C}_i{\delta }_i^{1/2}\sin {\alpha }_i\cos {\alpha }_i\hfill \\ \multicolumn{1}{c}{}& +\frac{Q_i}{l_i-{\psi }_b}\sin {\alpha }_i\cos {\alpha }_i-\frac{3}{2}{C}_e{\delta }_e^{1/2}\sin {\alpha }_e\cos {\alpha }_e\hfill \\ \multicolumn{1}{c}{}& +\frac{Q_e}{l_e-{\psi }_b}\sin {\alpha }_e\cos {\alpha }_e\hfill \end{array}$$

(16)

Global equilibrium equations

Figure 3 shows the bearing configuration of our experimental multi-bearing spindle system, in which a global Cartesian system (x, y, z) is used with the z-axis running along the bearing rotational axis; the origin is at the bearing centroid. A cylindrical system $(r,\phi ,z)$ is used to track the inner ring groove center (the transformation reference). $\phi$ is the angle between the positive x-axis and the r-axis. Specific teams of the local and global systems can be described as follows

$$\begin{array}{c}\mathrm{Groove}\mathrm{center}\mathrm{displacement}\mathrm{of}\mathrm{a}\mathrm{bearing}:\hfill \\ \left\{u\right\}=[\begin{array}{ccc}{u}_r\hfill & u_z\hfill & \phi \hfill \end{array}{]}^T\hfill \end{array}$$

(17)

$$\begin{array}{c}\mathrm{Inner}\mathrm{groove}\mathrm{center}\mathrm{loading}\mathrm{of}\mathrm{a}\mathrm{bearing}:\hfill \\ \left\{Q\right\}=[\begin{array}{ccc}{Q}_r\hfill & Q_z\hfill & M\hfill \end{array}{]}^T\hfill \end{array}$$

(18)

$$\begin{array}{c}\mathrm{Transformation}\mathrm{matrix}:\hfill \\ \left[R\theta \right]={\left[\begin{array}{ccccc}\cos \theta \hfill & \sin \theta \hfill & 0\hfill & -z_p\sin \theta \hfill & z_p\cos \theta \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & r_p\sin \theta \hfill & -r_p\cos \theta \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -\sin \theta \hfill & \cos \theta \hfill \end{array}\right]}^T\hfill \end{array}$$

(19)

$$\begin{array}{c}\mathrm{Bearing}\mathrm{global}\mathrm{displacement}:\hfill \\ \{\delta \}={\left[\begin{array}{ccccc}{\delta }_x\hfill & {\delta }_y\hfill & {\delta }_z\hfill & {\gamma }_x\hfill & {\gamma }_y\hfill \end{array}\right]}^T={\left[R\theta \right]}^T\left\{u\right\}\hfill \end{array}$$

(20)

$$\begin{array}{c}\mathrm{Bearing}\mathrm{global}\mathrm{force}:\hfill \\ \{F\}={\left[\begin{array}{ccccc}{F}_x\hfill & F_y\hfill & F_z\hfill & M_x\hfill & M_y\hfill \end{array}\right]}^T={\left[R\theta \right]}^T\left\{Q\right\}\hfill \end{array}$$

(21)

$$\begin{array}{c}\mathrm{Global}\mathrm{force}\mathrm{equilibrium}\mathrm{equation}:\hfill \\ \left\{F\right\}+\sum _{j=1}^n{\left[R\theta \right]}_j^T{\left\{Q\right\}}_j=\left\{0\right\}\hfill \end{array}$$

(22)

$$\begin{array}{c}\mathrm{Five}-\mathrm{by}-\mathrm{five}\mathrm{Jacobian}\mathrm{matrix}\mathrm{of}\mathrm{a}\mathrm{ball}:\hfill \\ \left[\frac{\partial \{F\}}{\partial {\left\{\delta \right\}}^T}\right]=-\sum _{j=1}^n{\left[R\theta \right]}_j^T·\frac{\partial {\left\{Q\right\}}_j^T}{\partial {\left\{u\right\}}_j}·{\left[R\theta \right]}_j\hfill \end{array}$$

(23)

Figure 3.

Bearing configuration and coordinates.

Bearing stiffness is given by

$$K_b=\left[\begin{array}{ccccc}{K}_{\mathrm{xx}}\hfill & K_{\mathrm{xy}}\hfill & K_{\mathrm{xz}}\hfill & K_{x{\theta }_x}\hfill & K_{x{\theta }_y}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{yx}}}& K_{\mathrm{yy}}\hfill & K_{\mathrm{yz}}\hfill & K_{y{\theta }_x}\hfill & K_{y{\theta }_y}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{zx}}}& K_{\mathrm{zy}}\hfill & K_{\mathrm{zz}}\hfill & K_{z{\theta }_x}\hfill & K_{z{\theta }_y}\hfill \\ \multicolumn{1}{c}{K_{{\theta }_{x}x}}& K_{{\theta }_{x}y}\hfill & K_{{\theta }_{x}z}\hfill & K_{{\theta }_x{\theta }_x}\hfill & K_{{\theta }_x{\theta }_y}\hfill \\ \multicolumn{1}{c}{K_{{\theta }_{y}x}}& K_{{\theta }_{y}y}\hfill & K_{{\theta }_{y}z}\hfill & K_{{\theta }_y{\theta }_x}\hfill & K_{{\theta }_y{\theta }_y}\hfill \end{array}\right]$$

(24)

The physical meaning and symmetric characteristics of the dominant terms in this matrix may be found in Li and Shin. ¹⁹

Thermal behavior

Bearing heat generation

Total friction torque

As the spindle rotates, the balls and the inner ring come into contact; the contact of the balls with the outer ring creates frictional losses and heat is thus generated. The amount of heat generated depends on the friction torque, which varies with spindle geometry, rotational speed, the type of lubricant, and the internal and external loads.^15,17,19 The total friction torque is in accordance with Palmgren ²

$$M_t=M_l+M_v+M_f$$

(25)

where $M_l$ is the torque caused by internal load and $M_v$ is the torque attributable to lubricant viscosity. For ball bearings, the torque caused by flanges $(M_f)$ is not considered.

After considering the inner and outer ring factors, $M_t$ can be rewritten as ¹⁵

$$M_{\mathrm{ti}}=\frac{1}{n}\frac{D_a}{r_e}[\frac{1}{2}{f}_{\mathrm{om}}{\left(\frac{F_{\mathrm{si}}}{C_{\mathrm{si}}}\right)}^{1/3}{F}_{\beta i}{d}_m+0.225f_0{\left(\eta \omega \right)}^{2/3}{d}_m^3]$$

(26)

$$M_{\mathrm{te}}=\frac{1}{n}\frac{D_a}{r_i}[\frac{1}{2}{f}_{\mathrm{om}}{\left(\frac{F_{\mathrm{se}}}{C_{\mathrm{se}}}\right)}^{1/3}{F}_{\beta e}{d}_m+0.225f_0{\left(\eta \omega \right)}^{2/3}{d}_m^3]$$

(27)

where $f_0$ depends on the bearing type and $f_0=1$ depends on angular contact ball bearings. In this study, we use oil-mist lubricant, then $f_{\mathrm{om}}=0.001$ .

Using the $C_s$ , $F_{\beta }$ , and $F_s$ parameters of Harris, ¹⁷ where $C_s$ is the static load rating, we derive

$$C_{\mathrm{si}}={\phi }_{\mathrm{si}}\mathrm{kn}{D_a}^2\cos {\alpha }_o{C}_{\mathrm{se}}={\phi }_{\mathrm{se}}\mathrm{kn}{D_a}^2\cos {\alpha }_o$$

(28)

where k is the number of bearing rows and n is the number of balls in a bearing. For the angular contact ball bearings used here, k = 1. ${\phi }_s$ is related to the load rating.

For angular contact ball bearings, ${\phi }_s$ is defined using the dimensionless values of the semi-major and semi-minor axes of the contact ellipse ( $a^{*}$ and $b^{*}$ , respectively)

$${\phi }_{\mathrm{si}}=23.8\frac{{\left(a_i^{*}{b}_i^{*}\right)}^3}{{(4-\frac{D_a}{r_{\mathrm{ig}}}+\frac{2\lambda }{1-\lambda })}^2},{\phi }_{\mathrm{se}}=23.8\frac{{\left(a_e^{*}{b}_e^{*}\right)}^3}{{(4-\frac{D_a}{r_{\mathrm{eg}}}-\frac{2\lambda }{1+\lambda })}^2}$$

(29)

where $\lambda$ is

$$\lambda =\frac{D_a}{d_m}\cos {\alpha }_o$$

(30)

$a^{*}$ and $b^{*}$ can be derived from the ellipticity parameters, that is, the first elliptic integrals

$$a^{*}={\left(\frac{2{\Im }^2\xi }{\pi }\right)}^{1/3},b^{*}={\left(\frac{2\xi }{\pi \Im }\right)}^{1/3}$$

(31)

$F_{\beta }$ reflects the dependence of the applied load on load magnitude and load direction

$$F_{\beta }=0.9F_a\mathrm{ctn}{\alpha }_o-0.1F_r$$

(32)

where $\mathrm{ctn}{\alpha }_o$ is the cotangent function of the contact angle $\alpha$ . $F_a$ and $F_r$ are the axial and radial loads, respectively. $F_s$ is the static equivalent load

$$F_s=X_s{F}_r+Y_s{F}_a$$

(33)

$X_s$ and $Y_s$ depend on the type of bearing used and the initial contact angle. For angular contact ball bearings with an initial contact angle of 15°, $X_s=0.5$ and $Y_s=0.47$ .

Spin moment

Heat is also generated by the spin moment^15,17,19

$$M_{\mathrm{spi}}=\frac{3}{8}{\chi }_s{Q}_i{\zeta }_i{\vartheta }_i,M_{\mathrm{spe}}=\frac{3}{8}{\chi }_s{Q}_e{\zeta }_e{\vartheta }_e$$

(34)

where ${\chi }_s$ is the frictional coefficient and $\vartheta$ is the third elliptic integral. $\vartheta$ is estimated according to several factors, including the dimensionless semi-major axes of the contact ellipses between the ball and the inner and outer races.

The major axes of the contact ellipses $\zeta$ between each ball and the inner and outer raceways are

$${\zeta }_i=0.236a_i^{*}\left(\frac{Q_i}{\sum {\rho }_i}\right),{\zeta }_e=0.236a_e^{*}\left(\frac{Q_e}{\sum {\rho }_e}\right)$$

(35)

where $a_i^{*}$ and $a_e^{*}$ are calculated from equation (31); these correspond to the bearing raceways. $\sum {\rho }_i$ and $\sum {\rho }_e$ represent the sum of the curvature of contact deformations between the ball and the inner and outer rings, respectively. The sum of the curvature can be derived from the inner groove radius $r_{\mathrm{ig}}$ and the outer groove radius $r_{\mathrm{eg}}$

$$\begin{array}{c}\sum {\rho }_i=\frac{1}{D_a}(4-\frac{D_a}{r_{\mathrm{ig}}}+\frac{2\lambda }{1-\lambda }),\hfill \\ \sum {\rho }_e=\frac{1}{D_a}(4-\frac{D_a}{r_{\mathrm{eg}}}-\frac{2\lambda }{1+\lambda })\hfill \end{array}$$

(36)

Total heat generation

Heat generation in each contact area is given by Jorgensen and Shin ¹⁵

$$H_i={\omega }_{\mathrm{obt}}\frac{D_a}{d_m}{M}_{\mathrm{ti}}+{\omega }_{\mathrm{si}}{M}_{\mathrm{spi}},H_e={\omega }_{\mathrm{obt}}\frac{D_a}{d_m}{M}_{\mathrm{te}}+{\omega }_{\mathrm{se}}{M}_{\mathrm{spe}}$$

(37)

From equation (37), the total contact heat generation $H_t$ is theoretically the sum of the heat generated at the inner and outer rings

$$H_t=H_i+H_e$$

(38)

The orbital velocity ${\omega }_{\mathrm{obt}}$ and spin angular velocities ${\omega }_{\mathrm{si}}$ and ${\omega }_{\mathrm{se}}$ are in accordance with Hirano ²¹

$${\omega }_{\mathrm{obt}}=\omega {[1+\frac{(\sigma +\cos {\alpha }_e)(\cos \beta \cos \beta \prime \cos {\alpha }_i+\sin \beta \sin {\alpha }_i)}{(\sigma -\cos {\alpha }_i)(\cos \beta \cos \beta \prime \cos {\alpha }_e+\sin \beta \sin {\alpha }_e)}]}^{-1}$$

(39)

$$\begin{array}{c}{\omega }_{\mathrm{si}}=-{\omega }_b\sin (\beta -{\alpha }_i)+(\omega -{\omega }_{\mathrm{obt}})\sin {\alpha }_i\hfill \\ {\omega }_{\mathrm{se}}={\omega }_b\sin ({\alpha }_e-\beta )-{\omega }_{\mathrm{obt}}\sin {\alpha }_e\hfill \end{array}$$

(40)

where ${\omega }_b$ is the angular velocity of the ball

$${\omega }_b={\omega }_{\mathrm{obt}}\frac{\sigma +\cos {\alpha }_e}{\cos \beta \cos \beta \prime \cos {\alpha }_i+\sin \beta \sin {\alpha }_e}$$

(41)

The angles $\beta$ and $\beta \prime$ in Figure 4 are the ball rotational angles in the x–z and y–z planes, respectively. Hirano ²¹ and Jorgensen and Shin ²² indicated that although $\beta \prime$ varies markedly with rotational speed, the assumption $\beta \prime =0$ affords reasonable predictions of both $\beta$ and rotational velocity.

Figure 4.

Rotational velocities and angles.

Heat contribution of the built-in motor

The specifications of the three-phase alternating current (AC) induction motor are listed in Table 1. The total heat generated by the motor is a function of the motor torque $T_m$ , motor efficiency ${\eta }_m$ , and motor frequency speed $f_m$ , as given by Bossmanns and Tu ²³

$$H_m=2\pi f_m{T}_m\frac{1-{\eta }_m}{{\eta }_m}$$

(42)

where the motor efficiency ${\eta }_m$ is estimated based on the maximum efficiency, and the coefficients reflect motor load and speed. The distribution of motor heat generation between the stator and rotor is a function of rotational speed, motor torque, and motor slip, and can be written as

$$H_r=H_m\frac{f_{\mathrm{sl}}}{f_{\mathrm{sy}}},H_s=H_m-H_r$$

(43)

where $H_r$ and $H_s$ represent the heat distributed to the motor rotor and stator, respectively, and $f_{\mathrm{sy}}$ is the synchronous frequency of the motor field (i.e. the set speed). The slip frequency in equation (43) is

$$f_{\mathrm{sl}}=100\frac{f_{\mathrm{sy}}-f_m}{f_{\mathrm{sy}}}$$

(44)

Table 1.

AC motor specifications.

Collet	Output power	Rotary speed	Frequency speed	Phase	Pole	Voltage	Rated current	Cooling	Torque	Motor weight
ER-11	700 W	24,000 r/min	400 Hz	3ϕ	2P	180 V	4.5 A	Air	0.03 kg m	1.65 kg

Heat dissipation

Heat convection depends on the ambient fluid and solid body velocity. We consider the heat convection of a ball rolling in lubricant, the convection of a spindle rotating on its axis, and convection of the housing in ambient air.

Heat convection of a ball in lubricant

The convection coefficient of a ball at a certain orbital velocity is calculated using the Nusselt number. However, this number depends on the relationships among the Prandtl, Reynolds, and Grashof numbers. The Prandtl number is $\mathrm{Pr}=(\mu c_p/k_l)=(\eta \rho c_p/k_l)$ , where $c_p$ , $k_l$ , and $\rho$ are the specific heat, thermal conductivity, and density of the lubricant, respectively. The dynamic viscosity of the lubricant is $\mu =\eta \rho$ .

We used an oil-mist lubricant, for which the heat convection coefficient of a ball is ¹⁵

$$h_b=1.8283k_l\mathrm{Pr}{\left(\frac{{\omega }_c}{2\eta }\right)}^{1/2}$$

(45)

Heat convection of the inner ring

The Reynolds number for the inner ring in lubricant is

$${\mathrm{Re}}_i=\frac{\omega D_i}{\mu }$$

(46)

For Re _i  < 4 × 10⁵, as in Harris, ¹⁷ the convection coefficient of the inner ring is

$$h_i=\frac{0.19k_l({\mathrm{Re}}_i^2+{\mathrm{Gr}}_i^2)}{D_i}$$

(47)

where Gr _i is the Grashof number

$$\mathrm{G}{\mathrm{r}}_i=\frac{\varsigma g(T_s-T_{\infty })D_i^2}{{\mu }^2}$$

(48)

where $T_{\infty }$ (°C) is room temperature, $T_s$ (°C) is the temperature of the inner ring surface, g is the acceleration caused by gravity, and $\varsigma$ is the thermal expansion coefficient.

Heat convection of the spindle

The rotation of a spindle shaft in quiescent air was given in Özerdem; ²⁴ radiation is ignored so that the Nusselt number becomes

$${\overline{\mathrm{Nu}}}_s=0.318{\mathrm{Re}}_s^{0.571}$$

(49)

The spindle convection coefficient is then

$$h_s=\frac{k_{\mathrm{air}}{\overline{\mathrm{Nu}}}_s}{D_s}$$

(50)

Here, the spindle diameter $D_s$ varies by the position of the bearings. Therefore, we consider the values of the coefficient at different spindle positions.

Heat convection produced by the housing

The heat convection coefficient of housing in air at room conditions is derived by measuring temperature changes^15,17

$$h_h=23\mathrm{\Delta }{T}_h^{0.25}$$

(51)

where $\mathrm{\Delta }{T}_h$ is the change from the initial temperature. This is calculated using the experimental temperatures shown in Figure 5: the temperature is highest at the front of the housing, and the temperature is similar between the middle and rear parts of the housing. Four heat sources (four NSK 7904C bearings) are located in the front. Two heat sources (two NSK 7001C bearings) and one heat source (the inbuilt motor) are located at the rear and middle parts, respectively, of the spindle housing. In addition, at a specific spindle rotational speed (below 10,000 r/min), theory indicates that the heat generated by the bearings and motor is approximately equal. For example, at 5000 r/min (see Figure 6(a) and (b)), the heat generated by an NSK 7904C bearing, an NSK 7001C bearing, and the motor is 8.273, 7.032, and 10.76 W, respectively. However, the temperature of the middle housing is lower than that of the other two housing regions, because the middle part housing has the largest volume and surface area.

Figure 5.

Transient temperature of the experimental spindle housing at a rotational speed of 5000 r/min.

Figure 6.

(a) Heat generation by the motor and heat distribution to the stator and rotor, and (b) the heat generated by the bearings.

Figure 6(a) shows that heat is generated in the motor of the spindle system, principally in the stator.

Methodology and the ANSYS finite element model

Figure 7 shows the workflow for computing thermal effects on bearing stiffness in our spindle system. The entire procedure can be divided into two main stages: the first one is the computation of the thermal behavior of the spindle bearing system using an ANSYS finite element model; then, we use the thermally affected elements identified in the first stage to solve the bearing stiffness matrix by employing the Newton–Raphson method in the MATLAB environment.

Figure 7.

Block diagram of the model.

In the first stage, the boundary conditions for the ANSYS steady-state model include heat generation, thermal dissipation, and the ambient temperature. The total heat generated is calculated by reference to bearing geometry, lubricant viscosity torque, orbital rotational speed, and the preloads and spin moments (equation (34)). The heat generated in the motor, and the distribution thereof, is shown in Figure 6(a). Heat dissipation is determined by considering the convection coefficients of balls in lubricant $h_b$ (equation (45)), inner ring surfaces $h_i$ (equation (47)), the spindle shaft $h_s$ (equation (50)), and the housing $h_h$ (equation (51)). We measured temperature increases in the front, rear, and middle parts of the housing relative to the initial temperatures. The housing convection coefficients were then defined using equation (48).

After solving the thermal behavior, the thermal results were loaded into the ANSYS static structural model. For fitting a sphere through the coordinates of a set of points, at least four points on the spherical surface were collected. In this study, the coordinates of six nodal points on six regions of the ball surface were exported from the ANSYS static structural model to fit the radial changes (see Figure 8). For example, P₀₁(x₀₁, y₀₁, z₀₁) is the initial coordinate and P₁(x₁, y₁, z₁) is the final coordinate, after thermal factors are considered. As the coordinates are known, the ball radius can be computed using the sphere-fitting technique; the radius is approximated using six points. The thermal expansion ${\psi }_b$ between the initial and final radius is measured. As the ball radius increases, a more amount of heat in the bearing is generated and leads to a consequent thermal expansion. Figure 9 outlines the process of estimating the bearing radial expansion from the initial bearing parameters to the values at the steady-state condition; this is repeated until the radial change is so small that minimal heat generated. Thus, the steady-state thermal condition is attained, yielding temperature changes and amount of the inner ring-shaft, outer ring-housing thermal expansion. These parameters are then used to solve local equilibrium equations (equation (12)) and the global equilibrium (equation (22)).

Figure 8.

Nodal point coordinates before and after thermal displacement: (a) six-nodal point area, (b) initial coordinates, and (c) coordinates after deformation.

Figure 9.

The steady-state thermal diagram.

Previous studies^7,9,11,15,25 assumed that heat was transmitted equally from a ball to the bearing inner and outer rings. In fact, as balls rotate in a spindle system, heat is generated by friction torque. All nodal points on a ball surface contribute heat that affects thermal behavior via convection. Therefore, by fitting the ball radius using the coordinates of the set of points defined above (see Figure 8), thermal boundary conditions (the extent of internal heat generation) and physical conditions (changes in convection coefficients) are updated (see Figure 9); this overcomes a major limitation of previous thermal resistance networks.

Experimental system

Figure 10 shows the instruments used for data acquisition, the spindle simulator, and a cross-sectional view of the experimental spindle system. The front and rear regions have four NSK 7904C bearings and two NSK 7001C bearings, respectively. The stator, rotor, and air cooling systems of the motor are included. Temperature was measured using platinum PT100 resistance temperature detectors (RTDs) with three-wire interfaces (accuracy: ±0.5°C). RTDs were placed at the front, rear, and middle of the housing. Two more RTDs measured outer ring temperatures at the front (NSK 7904C) and rear (NSK 7001C) bearings. The RTDs were connected to an NI PXIE-4357 module of an embedded PXIE-8135 controller that sent data to a personal computer.

Figure 10.

The experimental system.

Results

Bearing geometry

The geometry of angular contact ball bearings is shown in Figure 11(a). Friction torque between the balls and raceways generates heat. The ball diameter $D_a$ , the bearing inner groove radius $r_{\mathrm{ig}}$ , and the bearing outer groove radius $r_{\mathrm{eg}}$ all affect thermal behavior significantly. We measured these parameters by cutting NSK 7904C and NSK 7001C bearings and measuring cross-sections using VMS-2010G machine, manufactured by Rational Precision Instruments Co., Ltd (resolution, 0.5 µm). The points used to fit the radii are shown in Figure 11(b) and (c). We also fit $r_{\mathrm{ig}}$ and $r_{\mathrm{eg}}$ values. Table 2 lists the specifications of the ball bearings. Apart from $D_a$ , $r_{\mathrm{ig}}$ , and $r_{\mathrm{eg}}$ , the other specifications were taken from datasheets.

Figure 11.

(a) Bearing geometry and measurement of ball diameters: (b) points measured and (c) fitting.

Table 2.

Bearing specifications.

Ball diameter, $D_a$ (mm)	Inner groove radius, $r_{\mathrm{ig}}$ (mm)	Outer groove radius, $r_{\mathrm{eg}}$ (mm)	Pitch diameter, $d_m$ (mm)	Width of bearing ring, $W$ (mm)	Bore diameter, $D_i$ (mm)	Outer ring diameter, $D_e$ (mm)	Number of ball, n	Initial contact angle, ${\alpha }_0$ (°)
NSK 7904C angular contact ball bearing parameters
4.695	2.4415	2.488	28.5	9	20	37	13	15
NSK 7001C angular contact ball bearing parameters
4.6915	2.4395	2.4865	20	8	12	28	10	15

Thermal behavior of the spindle

The ANSYS steady-state thermal model is shown in Figure 12. The inputs were the heat distribution in the motor and rotor, heat generation within bearing balls, heat convection coefficients, and the initial (room) temperature. The specific boundary conditions of the ANSYS steady-state model at a spindle speed of 5000 r/min are listed in Table 3.

Figure 12.

The ANSYS steady-state model.

Table 3.

Boundary conditions.

Boundary conditions of steady-state thermal model (spindle rotational speed 5000 r/min)
Internal heat generation-motor stator (W/m3)	67,406
Internal heat generation-motor stator (W/m3)	19,893
Internal heat generation-front ball (W/m3)	12.11 × 1e6
Internal heat generation-rear ball (W/m3)	23.58 × 1e6
Convection front ball (W/m2 °C)	1.2121
Convection rear ball (W/m2 °C)	1.1698
Convection front inner ring (W/m2 °C)	44.374
Convection rear inner ring (W/m2 °C)	44.374
Convection motor air cooling (W/m2 °C)	32.5
Convection front housing (W/m2 °C)	35
Convection middle housing (W/m2 °C)	32.5
Convection rear housing (W/m2 °C)	33.5
Ambient temperature (°C)	23

Figure 13 shows the simulated temperature changes at a spindle speed of 5000 r/min (see Figure 13(a)), and the transient temperatures measured at the bearing outer rings (see Figure 13(b)).

Figure 13.

(a) Temperatures at a rotational speed of 5000 r/min simulated using ANSYS and (b) transient temperatures measured experimentally at the bearing outer rings.

Figure 13(a) shows that the temperature was approximately 33°C and 32.5°C at the outer rings of the front and rear bearings, respectively. Compared to the initial temperature of 23°C, the values are in good agreement with the experimental results shown in Figure 13(b). We then varied the spindle rotational velocity from 1000 to 10,000 r/min. Using the boundary conditions listed in Table 3, the process of the thermal analysis for different spindle speeds was repeated and the results were shown in Figure 14. Figure 14 illustrates that the temperature changes in the simulated and experimental results were very similar.

Figure 14.

Simulated and experimental temperatures of the outer ring surfaces (front bearings) at different spindle rotational speeds.

After obtaining temperature data, ball radial thermal expansion was calculated by importing the body temperatures yielded by the ANSYS steady-state thermal model into the ANSYS static structural model, followed by fitting of the nodal point coordinates and computation of ball radii (see Figures 8 and 9). In the steady-state condition, the radial thermal expansion of balls in the front bearings is shown in Figure 15. Radial expansion of the inner ring-shaft and outer ring-housing was computed using equations (5) and (6), respectively.

Figure 15.

Thermal expansion of the front bearings.

Bearing stiffness analysis

Changes in radial, axial, angular, and coupling stiffness values by spindle rotational speed are shown in Figure 16. The off-diagonal terms in equation (23) are negligible if there is only axial-internal preloading. The nine dominant terms include the radial stiffness values ( $K_{\mathrm{xx}}$ (N/m) and $K_{\mathrm{yy}}$ (N/m)), axial stiffness value ( $K_{\mathrm{zz}}$ (N/m)), angular stiffness values ( $K_{{\theta }_x{\theta }_x}$ (Nm/rad) and $K_{{\theta }_y{\theta }_y}$ (Nm/rad)), and coupling stiffness values ( $K_{x{\theta }_y}$ (N/rad), $K_{y{\theta }_x}$ (N/rad), $K_{{\theta }_{x}y}$ (Nm/m), and ${K_{{\theta }_y}}_x$ (Nm/m)). However, given the symmetry of the stiffness matrix, the values of the five relatively large stiffness coefficients are presented in Figure 16. In general, as rotational speed increases, high-speed spindle bearings soften. Table 4 provides a more specific comparison by considering reductions in stiffness coefficients at a spindle speed of 10,000 r/min relative to the static condition. All dominant terms in the bearing stiffness matrices for both the front and rear bearings decreased. For the front bearings, the decreases in coupling stiffness coefficients (26.01%) were slightly greater than that in the axial (25.78%) and angular stiffness values (25.78%); the decrease in radial stiffness was the smallest (25.23%). Rear bearings exhibited similar tendencies. In particular, at 10,000 r/min, the reduction in rear bearing stiffness coefficient (∼4.88%) was much lower than that in the front bearings (∼26%), explained by the centrifugal force shown in Figure 17. The centrifugal forces acted on balls of the NSK 7904C bearing were much greater than those on balls of the NSK 7001C bearing; stiffness decreases more markedly.

Figure 16.

Elements of the bearing stiffness matrix according to the spindle rotational speed: (a) NSK 7904C and (b) NSK 7001C bearings.

Table 4.

Reductions in stiffness coefficients at 10,000 r/min compared to the static stiffness values.

Bearings	Decreases in stiffness (%)
	$K_{\mathrm{xx}}$	$K_{\mathrm{zz}}$	$K_{{\theta }_x{\theta }_x}$	$K_{x{\theta }_y}$	$K_{y{\theta }_x}$
NSK 7904C	25.23	25.78	25.78	26.01	26.01
NSK 7001C	4.53	4.94	4.94	4.88	4.88

Figure 17.

Increases in the centrifugal force on the front and rear bearings as the spindle spools up.

The relationship between centrifugal force and the decrease in angular contact ball bearing stiffness was explained earlier. ²² As the spindle spools up, the balls experience centrifugal force (see Figure 1); the contact between the ball and the bearing inner ring is pushed up and out from the center of the inner groove center. This causes the increase in the inner ring contact angle and the decrease in the outer ring contact angle, which are shown in Figure 18. The increase in the inner ring contact angle results in the decrease in the inner ring contact force; the increase in outer ring contact force is attributable to the decreased outer ring contact angle (see Figure 19). Finally, the balls tend to move away from the inner ring, and the stiffness coefficients decrease.

Figure 18.

The contact angles at the inner and outer rings of (a) the NSK 7904C bearing and (b) the NSK 7001C bearing.

Figure 19.

The contact forces at the inner and outer rings of (a) the NSK 7904C bearing and (b) the NSK 7001C bearing.

Behavior differed markedly according to whether thermally induced bearing stiffness was considered. Figure 20 shows the difference in stiffness ( $\mathrm{\Delta }K$ values) according to whether thermal effects were incorporated into the calculations. The difference is remarkable and increases as rotational speed rises, because thermal expansion increases the contact deformations ${\delta }_i$ and ${\delta }_e$ (equations (3) and (4)). Greater deformation increases the bearing contact forces $Q_i$ and $Q_e$ (equation (7)), ultimately increasing the bearing stiffness (equation (24)). Stiffness would vary less if bearings were fabricated from a material with a high Young’s modulus.

Figure 20.

(a) Difference in stiffness according to whether thermal effects are considered and (b) the rises of the stiffness differences as speed increases.

For validation, Figure 21 compares the thermally affected radial stiffness values derived in this study for the NSK 7014C angular contact ball bearing to those of Jorgensen and Shin. ¹⁵ With the preload input $(F_z=540\mathrm{N})$ as used in the work of Jorgensen and Shin, and using the similar thermal network (i.e. oil-mist lubrication), good agreement is evident. The tendencies are equivalent; the stiffness lines begin to decrease at 2000 r/min and the stiffness values vary as rotational speed increases because the thermal-related coefficients change, as do the distributed heat sources.

Figure 21.

Comparison of our radial stiffness measurements to those of Jorgensen and Shin (NSK 7014C bearing).

Conclusion

We built a thermal model to estimate heat generation by, and thermal distribution to, the ball bearings in a high-speed spindle system. We combined steady-state thermal and static structural models on an ANSYS workbench. The temperature and thermal expansion of bearing balls were calculated by reference to changes in nodal coordinates, and the extent of heat generation was then re-calculated. We simulated and experimentally verified outer ring temperatures. This study improves our understanding of thermal behavior in spindle systems with complex geometry.

We constructed a five-by-five bearing stiffness matrix that included radial, axial, angular, and coupling stiffness. The dominant matrix terms decreased dramatically as spindle rotational speed increased, and thermal effects increased stiffness. We validated the model by comparing the radial stiffness values generated herein to those of a previous study.