A general thermal model of machine tool spindle

Abstract

As the core component of machine tool, the thermal characteristics of the spindle have a significant influence on machine tool running status. Lack of an accurate model of the spindle system, particularly the model of load–deformation coefficient between the bearing rolling elements and rings, severely limits the thermal error analytic precision of the spindle. In this article, bearing internal loads, especially the function relationships between the principal curvature difference F(ρ) and auxiliary parameter n_δ, semi-major axis a, and semi-minor axis b, have been determined; furthermore, high-precision heat generation combining the heat sinks in the spindle system is calculated; finally, an accurate thermal model of the spindle was established. Moreover, a conventional spindle with embedded fiber Bragg grating temperature sensors has been developed. By comparing the experiment results with simulation, it indicates that the model has good accuracy, which verifies the reliability of the modeling process.

Introduction

A study on the thermal model and improvement of the computer numeric control (CNC) machine tool spindle accuracy is the most important pursuit for researchers in precision manufacturing field. Moreover, obtaining an accurate model of the spindle is beneficial to analyzing the bearing running status and the thermal error. However, the fundamental prerequisite is to analyze the spindle bearing internal action mechanism.

There will be a large amount of heat when the spindle is running at a high speed. And the heat generated by the bearing is the main heat source of the spindle. To calculate the heat source of the bearing, the bearing internal loads’ analysis has to be done. According to Harris’ ¹ research, the internal loads of bearing have complex relations with the contact angles between the balls and raceways, as it exits many nonlinear equations with coupled parameters; many researchers, such as C Jin et al., ² T Xu et al., ³ and C Ma et al., ⁴ have done the bearing internal loads’ analysis in the recent years. But they did not show the detailed calculation steps, especially the calculation of load–deformation coefficient, which is coupled with the contact angles between the ball and raceway and is different at each ball position causing ball centrifugal force and gyroscopic moment; and semi-major axis a and semi-minor axis b of the spinning contact area between the balls and raceways are not given clearly either.

Only when the bearing internal action mechanism is analyzed, the bearing heat generation can be received using the data of bearing internal loads, and the spindle temperature field can be solved. Based on the bearing load analysis, the variation in the spindle temperature field machine tools can be analyzed based on finite element model (FEM), such as C Ma et al., ⁴ JK Choi and DG Lee, ⁵ M Mori et al., ⁶ H Zhao et al., ⁷ E Creighton et al., ⁸ D Chen et al., ⁹ WS Wang et al., ¹⁰ Q Yang et al. ¹¹ and S Jiang and H Mao. ¹² Based on FEM method, the accurate result can be obtained, but the result is based on the accurate calculation of heat source, not the method itself, and the heat sources and boundary conditions of the spindle system have been precisely calculated in this article.

In this article, a conventional spindle was studied; the calculations of bearing internal loads, heat generation, and heat sinks of the spindle are described in detail. Moreover, the specific relationships between the principal curvature difference F(ρ) and auxiliary parameter n_δ, contact ellipse semi-major axis a, and semi-minor axis b to accurately calculate the ball spinning friction torque on the contact areas were given out precisely, and also FEM method is used to solve the thermal model. In addition, the fiber Bragg grating (FBG) temperature sensors were selected to accomplish the measurement, which was the new sensing material in monitoring the mechanical system dynamic.^13,14 The nomenclatures in this article are shown in Table 1.

Table 1.

Nomenclature.

Db  diameter of the ball rb  radius of the ball M  total frictional torque of the bearing i   refer to inner ring d  bearing inside diameter Z  number of balls wi rotation speed of bearing element fo outer ring groove curvature radius coefficient f1 inner ring groove curvature radius coefficient n  bearing rotating speed, r/min $w_{\mathrm{si}(o)}$ rotation speed of ball on the contact areas

$F_a,F_r$ external loads r radius of the raceway groove curvature ${\alpha }^0$  free contact angle o refer to outer ring D bearing outside diameter dm pitch circle diameter Σρi sum of balls and inner ring principal curvature Σρo sum of balls and outer ring principal curvature $F(\rho )$  principal curvature difference

Bearing internal loads’ analysis

The angular contact ball bearings are the most popular bearing type used in the high-speed spindle for machine tools, which are also utilized in this article. For a thermal model of spindle system, the accurate heat sources have to be received. For the generated heat of the angular contact ball bearing, the internal loads of the bearings have to be obtained first, which are related to the bearing generated heat. Modeling the bearing internal loads is based on the quasi-static mechanical analysis of the bearing, which contains the ball centrifugal force and bearing gyroscopic moment, which is the classical and most popular method, and also applied in this article. For angular contact ball bearing, simplified schematic diagram is shown in Figure 1.

Figure 1.

Simplified schematic diagram of angular contact ball bearing.

When loads and rotation speed were applied to the bearing, the ball positions will be changed. First, the force balance equations can be established. The force applied to every rolling element is depicted in Figure 2.

Figure 2.

Ball loads at angle ψ_j.

According to Figure 2, the ball force balance equations can be obtained

$$\begin{array}{l}{Q}_{ij}\sin {\alpha }_{ij}-Q_{oj}\sin {\alpha }_{oj}\\ -\frac{M_{gj}}{D}({\lambda }_{ij}\cos {\alpha }_{ij}-{\lambda }_{oj}\cos {\alpha }_{oj})=0\end{array}$$

(1)

$$\begin{array}{l}{Q}_{ij}\cos {\alpha }_{ij}-Q_{oj}\cos {\alpha }_{oj}\\ +\frac{M_{gj}}{D}({\lambda }_{ij}\sin {\alpha }_{ij}-{\lambda }_{oj}\sin {\alpha }_{oj})+F_{cj}=0\end{array}$$

(2)

where $Q_{\mathrm{ij}}={\kappa }_{\mathrm{ij}}{{\delta }_{\mathrm{ij}}}^{1.5}$ , $Q_{\mathrm{oj}}={\kappa }_{\mathrm{oj}}{{\delta }_{\mathrm{oj}}}^{1.5}$ , and ${\kappa }_{\mathrm{ij}/\mathrm{oj}}=2.15\times {10}^5(\sum {\rho }_{i/j}{)}^{-1/2}(n_{\delta }{)}^{-3/2}$ . n_δ is the auxiliary parameter which is calculated by the function with principal curvature difference F(ρ) obtained through polynomial curve fitting; the contrast of different curve fitting methods is shown in Figure 3; the dots show the original data. ¹⁵ And the specific equation F1(x) is as follows

$$\begin{array}{cc}\multicolumn{1}{c}{F1(x)=}& -2837.329x^9+12437.693x^8\hfill \\ \multicolumn{1}{c}{}& -22725.233x^7+22396.826x^6\hfill \\ \multicolumn{1}{c}{}& -12868.75x^5+4345.347x^4-824.367x^3\hfill \\ \multicolumn{1}{c}{}& +77.79x^2-2.726x+1.004\hfill \end{array}$$

(3)

Figure 3.

n_δ curve fitting methods’ contrast.

For a specific bearing, when the contact angles are known, the principal curvature difference F(ρ) is decided. So using equation F1(x), n_δ can be obtained, which can be used to calculate the bearing load–deformation coefficient. Particularly, when the principal curvature difference F(ρ) of the contact area is in the range of 0.85–0.999, which are the appropriate values for angular contact ball bearing. This calculation method can avoid the complicated process to build the relationship with the contact angles and show their relationships clearly.

For steel ball bearings in this article, the centrifugal force and gyroscopic moment acting on a ball are calculated as follows^2–4

$$F_{\mathrm{cj}}=2.26\times {10}^{-11}\mathrm{dm}{D}^3{n}_{\mathrm{mj}}$$

(4)

$$M_{\mathrm{gj}}=4.5\times {10}^{-12}{n}_{\mathrm{bj}}{n}_{\mathrm{mj}}\sin {\beta }_j$$

(5)

$$\begin{array}{c}{n}_{\mathrm{bj}}=\mathrm{dm}/D*n*\hfill \\ \frac{(1-D\cos {\alpha }_{\mathrm{ij}}/\mathrm{dm})(1+D\cos {\alpha }_{\mathrm{oj}}/\mathrm{dm})}{(1-D\cos {\alpha }_{\mathrm{ij}}/\mathrm{dm})\cos ({\alpha }_o-{\beta }_j)+(1+D\cos {\alpha }_{\mathrm{oj}}/\mathrm{dm})\cos ({\alpha }_i-{\beta }_j)}\hfill \end{array}$$

(6)

$$n_m=\frac{n*(1-D\cos {\alpha }_{\mathrm{ij}}/\mathrm{dm})\cos ({\alpha }_{\mathrm{oj}}-{\beta }_j)}{(1+D\cos {\alpha }_{\mathrm{oj}}/\mathrm{dm})\cos ({\alpha }_{\mathrm{ij}}-{\beta }_j)+(1-D\cos {\alpha }_{\mathrm{ij}}/\mathrm{dm})\cos ({\alpha }_{\mathrm{oj}}-{\beta }_j)}$$

(7)

If outer raceway control is approximately applied at a given ball location

$${\beta }_j=\frac{\sin {\alpha }_{\mathrm{oj}}}{\cos {\alpha }_{\mathrm{oj}}+D/\mathrm{dm}}$$

(8)

Second, the relative positions of the ball and inner and outer raceway groove curvature centers can be found, as shown in Figure 4. It is assumed that the outer raceway groove curvature center is fixed in space and the inner raceway groove curvature center moves relative to the fixed center. As observed from Figure 4, the contact angles of the ball and inner and outer raceway are dissimilar. The relative axial and radial deformations between the inner ring and the outer ring are δ_a and δ_r, respectively. A_1j, A_2j and X_1j, X_2j are the auxiliary parameters. δ_ij and δ_oj are the deformations of the inner and outer raceway at angle ψ_j, respectively.

Figure 4.

Relative positions of ball center and raceway groove curvature centers.

The auxiliary parameters A_1j, A_2j, X_1j and X_2j and the equilibrium equations of the entire bearing are defined as follows

$$A_{1j}=\mathrm{BD}\sin {\alpha }_0+{\delta }_a$$

(9)

$$A_{2j}=\mathrm{BD}\cos {\alpha }_0+{\delta }_r\cos {\psi }_j$$

(10)

$$\begin{array}{c}\cos {\alpha }_{\mathrm{oj}}=\frac{X_{2j}}{(f_o-0.5)D+{\delta }_{\mathrm{oj}}},\sin {\alpha }_{\mathrm{oj}}=\frac{X_{1j}}{(f_o-0.5)D+{\delta }_{\mathrm{oj}}}\hfill \\ \cos {\alpha }_{\mathrm{ij}}=\frac{A_{2j}-X{}_{2j}}{(f_i-0.5)D+{\delta }_{\mathrm{ij}}},\sin {\alpha }_{\mathrm{ij}}=\frac{A_{1j}-X{}_{1j}}{(f_i-0.5)D+{\delta }_{\mathrm{ij}}}\hfill \end{array}$$

(11)

$$F_a-\sum _{j=1}^Z(Q_{\mathrm{ij}}\sin {\alpha }_{\mathrm{ij}}-\frac{{\lambda }_{\mathrm{ij}}{M}_{\mathrm{gj}}}{D}\cos {\alpha }_{\mathrm{ij}})=0$$

(12)

$$F_r-\sum _{j=1}^Z(Q_{\mathrm{ij}}\cos {\alpha }_{\mathrm{ij}}+\frac{{\lambda }_{\mathrm{ij}}{M}_{\mathrm{gj}}}{D}\sin {\alpha }_{\mathrm{ij}})\cos {\psi }_j=0$$

(13)

where F_a and F_r are the applied loads on the bearing.

For calculating bearing internal loads, α_i, α_o, δ_i, δ_o, X₁, X₂, A₁, A₂ at each ball position, inner ring, and outer ring relative deformations δ_a and δ_r should be determined, combining with nonlinear equations (1), (2), and (9)–(13), which can also be solved by Newton–Raphson method. After obtaining all the parameters, the bearing working status can be received.

The specific parameters of 7012AC bearing are shown in Table 2. The bearing without radial load is calculated in Figures 5–7.

Table 2.

The specific parameters’ values of the bearing.

Parameter name	Value	Parameter name	Value	Parameter name	Value
d (mm)	60	D (mm)	95	Bearing width (mm)	18
Z	19	$D_b$ (mm)	10.7	$d_m$ (mm)	77.5
α 0	25	fo	0.525	fi	0.515

Figure 5.

${\alpha }_i$ , ${\alpha }_o$ versus different axial loads and spindle speed.

Figure 6.

$Q_i$ , $Q_o$ versus different axial loads and spindle speed.

Figure 7.

Axial deformation versus different axial loads and spindle speed.

From Figure 5, it can be observed that the contact angles α_i between the balls and the inner raceway are comparatively bigger than the contact angle α_o between the balls and the outer raceway in the same bearing axial load and speed conditions; the smaller the axial load, the bigger the difference between the contact angles at the same bearing rotation speed. And in a ball bearing, depending on the contact angles, ball gyroscopic moments and ball centrifugal forces can be of significant magnitude such that the inner raceway contact angles tend to increase and the outer raceway contact angles tend to decrease. From Figure 6, the contact load Q_i between the balls and the inner raceway is relatively smaller than the contact load Q_o with the same bearing axial load and speed; the contact loads Q_i and Q_o are symmetrically distributed on the sides of static situation with different axial loads. And the axial deformations tend to become similar to the static station with large axial loads at big bearing speed, which can be observed in Figure 7.

Spindle thermal analysis

Bearing heat generation

In the machine tool operations, the major heat sources include the heat generated by the cutting process and the heat from the bearings. It is assumed that the majority of cutting heat is taken away by coolant, and therefore, the heat generated by bearings is the dominant heat causing the thermal deformations. After receiving bearing working data, the heat generated by a bearing can be computed by the following equation ²

$$H_f=1.047\times {10}^{-4}\mathrm{nM}$$

(14)

where n is the rotating speed of the bearing, M is the total frictional torque of the bearing, and H_f is the generated heat. The frictional torque M is a sum of three torques: (1) the torque due to applied load, M₁; (2) the torque due to viscous friction, M₀; and (3) the torque due to spinning motion at contact area, M_si(o).

The torque due to the applied load

The torque due to the applied load can be empirically approximated by the following equation^2–4

$$M_1=f_1{P}_1{d}_m$$

(15)

where f₁ is a factor depending on the bearing design and relative load, and for single-row angular contact ball bearing, f₁ = 0.0013(p₀/c₀)^0.33; P₁ is the dynamical equivalent load, and for angular contact ball bearings, P₁ = F_a − 0.1F_r; $d_m$ is the pitch circle diameter, and d_m = 0.5(D + d).

The torque due to viscous friction

The torque due to viscous friction is related to bearing type, rotating speed, and lubricant type. For bearings that operate at moderate speeds and under non-excessive load, the viscous friction torque can be empirically expressed as follows^2–4

$$\begin{array}{c}\multicolumn{1}{c}{\nu n\ge 2000,M_0={10}^{-7}{f}_0(\nu n{)}^{2/3}{d}_m^3}\\ \multicolumn{1}{c}{\nu n\ge 2000,M_0=160\times {10}^{-7}{f}_0{d}_m^3}\end{array}$$

(16)

where f₀ is connected to bearing type and lubrication mode, and for single-row angular contact and grease lubrication ball bearings, f₀ = 2; v is the movement viscosity of the lubricant under operating temperature, and the movement viscosity of grease lubrication is 30.5 cst at 20°C, 18 cst at 40°C, 10 cst at 60°C, and 1 cst = 10⁻⁶ m²/s.

Spinning torque at contact area

The spin friction torque M_si(o) generated by ball’s spinning motion can be expressed as^2,4

$$M_{s(i/o)}=\frac{3\mu ·Q_{i/o}·a_{i/o}·E_{i/o}}{8}$$

(17)

So the total friction torque due to ball’s spinning motion can be obtained as

$$M_s=Z{M}_{s(i)}+Z{M}_{s(o)}$$

(18)

For calculating M_s, the elliptical semi-major axis and eccentricity ratio should be computed. Also, the relationships between the contact ellipse semi-major axis a, semi-minor axis b, and principal curvature difference F(ρ) were obtained by curve fitting; the specific results are shown as F2(x), F3(x) (Figure 8); the original data were also cited in Okamoto ¹⁵

$$\begin{array}{cc}\multicolumn{1}{c}{F2(x)=}& 9.218\times {10}^5{e}^{-{\left(\frac{x-1.024}{0.0073}\right)}^2}+15.88e^{-{\left(\frac{x-1.02}{0.02696}\right)}^2}\hfill \\ \multicolumn{1}{c}{}& +0.9293e^{-{\left(\frac{x-0.977}{0.008799}\right)}^2}\hfill \\ \multicolumn{1}{c}{}& +1.287e^{-{\left(\frac{x-0.9626}{0.0143}\right)}^2}+0.632e^{-{\left(\frac{x-0.9416}{0.0171}\right)}^2}\hfill \\ \multicolumn{1}{c}{}& +0.4765e^{-{\left(\frac{x-0.915}{0.03154}\right)}^2}\hfill \\ \multicolumn{1}{c}{}& +4.673\times {10}^{13}{e}^{-{\left(\frac{x-44.13}{7.86}\right)}^2}+0.1264e^{-{\left(\frac{x-0.8622}{0.6216}\right)}^2}\hfill \end{array}$$

(19)

$$F3(x)=\frac{(-203.8x^2-179.8x+388.4)}{(x^5-791.8x^4+1241x^3-1024x^2+216.1x+383)}$$

(20)

Figure 8.

The specific curve fitting methods’ contrast results.

Through equations (19) and (20), the contact ellipse semi-major axis and semi-minor axis can be decided by the principal curvature difference F(ρ); taking full advantage of this method can save large amount of computing time. According to the bearing 7012AC, the results of three different friction torque and heat generation are shown in Figures 9–11.

Figure 9.

Spinning friction torque of each ball.

Figure 10.

Heat generation caused by viscous friction.

Figure 11.

Heat generation caused by applied load friction.

Heat sinks in spindle

The spindle heat sinks mainly include the convection between the shaft and air; the ambient air; the thermal resistance between bearing outer ring and housing, bearing inner ring, and shaft; and the thermal conduction between the components.

First, the convection η between the shaft and air can be defined as ¹²

$$\eta =\frac{\mathrm{Nu}·k_{\mathrm{air}}}{d}$$

(21)

where Nu = 0.133Re^2/3Pr^1/3, d is the shaft diameter, Re = V × d/v_air, P_r is the Prandtl number of air, V is the speed of air on the shaft surface, v_air is the viscosity coefficient of air, and k_air is the air conduction coefficient. The specific result is shown in Figure 12.

Figure 12.

The convection between shaft and air under different spindle speed.

Second, a free convection coefficient can be assumed for ambient air around stationary surfaces such as the spindle housing. The coefficient can be obtained as a = 9.7 W/(m² K) according to Jiang and Mao. ¹²

Third, the thermal resistance between the bearing outer ring and housing is the reciprocal number of the thermal conduction coefficient number Π between the clearances. ¹² And the thermal conduction coefficient number $\Pi$ is

$$\Pi =\frac{1}{\frac{h_{\mathrm{ring}}}{{\lambda }_{\mathrm{ring}}}+\frac{h_{\mathrm{gap}}}{{\lambda }_{\mathrm{air}}}}A$$

(22)

where h_gap is the average of the clearance, h_ring is the thickness of the ring, λ_ring is the thermal conduction coefficient of the bearing ring, λ_air is the thermal conduction coefficient of air, and A is the contact area. According to the specific values, the thermal in this article is 1439.5 W/m² K of the shaft and bearing outer ring and 908.75 W/m² K of the housing and the bearing inner ring. All the needed boundary conditions for bearing temperature analysis are conformed so far.

Spindle temperature field simulation and measurement

Spindle temperature field simulation

The general and popular way for spindle thermal analysis is using the commercial FEM software, such as ANSYS. ANSYS 14.0 was selected for solving the spindle thermal field and bearing temperature rise in this article. The spindle system is simplified and reconstructed by computer-aided design (CAD) software. Specially, the mesh size is refined at the bearings to guarantee the results’ accuracy. After loading all the boundary conditions, the spindle field is obtained at a spindle speed of 1000 r/min and axial load of 0 N; the specific results are shown in Figure 13.

Figure 13.

Spindle temperature field simulation results under 2000 r/min and axial 0 N.

The simulation results show that the highest temperature area is at the position of bearings, especially at the spindle front end. And because of the mounting position of the experiment setup, there is no convection on the bottom face of the spindle housing; the bearing temperature has uneven temperature field, but the difference can be ignored.

The experimental apparatus is shown in Figure 14. The experimental apparatus is mainly a spindle system. The coupling connects the servo driver and the shaft to avoid the motor’s generated heat flow into the spindle housing; the milling force is uploaded to apply the axial loads; the housing is fixed on a vibration isolation platform. Figure 15 shows the specific experiment setup.

Figure 14.

Experimental apparatus.

Figure 15.

The specific experiment setup.

The bearing outer ring temperature measurement is based on FBG temperature sensor, which has unique advantages in small installation space and harsh environment.^13,14 The specific layout positions of FBG are shown in Figure 16. There is a ring groove at the contact position of the housing and the bearing outer ring; FBG temperature sensors are arranged in the substrates, which are printed by three-dimensional (3D) printer. Thermal conductive grease fills the gap between FBG and the bearing outer ring.

Figure 16.

FBG sensors’ installation locations.

After the installation of FBG sensors, they should be calibrated, since the sensitivity of FBG will be affected by thermal expansion of the thermal grease and its own bending. The whole bearing slider was placed in a constant temperature environment to calibrate the FBG sensors. The specific results are shown in Figure 17. From Figure 17, every FBG temperature sensor has a very good linear character to the temperature signal. The average slopes of the fitting curves are 10.72, 10.75, 10.27, 10.9, 10.44, and 10.96 separately.

Figure 17.

The calibration results of FBG sensor.

Spindle bearing temperature measurement

Bearing temperature rise is the most concerned parameter in the machine tool spindle system, which is related to the running state of bearing, spindle working life, and spindle thermal error; the measurement of bearing temperature is helpful to ascertain the bearing running characteristic and analyze the spindle thermal error. So monitoring the bearing temperature plays a critical role in the work process of spindle. So the bearing temperature was selected to verify the model built in this article. Figure 18 shows the measurement results of FBG temperature sensors at spindle speed 2000 r/min and axial load 0 N. From Figure 18, it can be observed that the thermal balance is achieved at 150 min; when the spindle stopped running at 230 min, the bearing temperature has small increase, caused by the decrease in shaft and air convection.

Figure 18.

The specific measurement temperature of FBG sensors.

After thorough analysis of the spindle heat generation and heat sinks in the spindle system, the measurement and simulation results of the bearing temperature under 1000 r/min, 2000 r/min, 3000 r/min, 0 N, 1000 N axial loads, and different ambient temperatures have no much difference with error smaller than 1°C. The specific contrast results are shown in Figure 19.

Figure 19.

The contrast results between the simulation and measurement results.

Conclusion

Nowadays, the spindles can be divided into two groups: conventional spindle and motorized spindle; the latter is developed from the former, but the core units—bearings—have not been changed. So a study on the traditional spindle can push for a resolution for the problems on the machine tool spindle. In this article, the thermal model of a conventional spindle is developed; the bearing internal loads between the balls and raceways, the heat generation, and heat sinks of a conventional machine tool have been put in detail further. Especially, the load–deformation coefficient between the raceways and the balls is explained clearly. Also, the functions between the principal curvature difference F(ρ) and the auxiliary parameter n_δ, semi-major axis a, and semi-minor axis b have been clearly given out. These functions can directly be utilized in the future to calculate the bearing internal load distribution and the spinning torque at the contact areas, reducing the whole calculating time. After receiving the bearing simulated temperature, the embedded FBG temperature sensors have been used to measure and verify the simulation bearing temperature. The contrast results show that the thermal model caught up in this article can be used for spindle state analysis.