Development of a Digital Model for Predicting the Variation in Bearing Preload and Dynamic Characteristics of a Milling Spindle under Thermal Effects

Abstract

The spindle tool is an important module of the machine tool. Its dynamic characteristics directly affect the machining performance, but it could also be affected by thermal deformation and bearing preload. However, it is difficult to detect the change in the bearing preload through sensory instruments. Therefore, this study aimed to establish a digital thermal–mechanical model to investigate the thermal-induced effects on the spindle tool system. The technologies involved include the following: Run-in experiments of the milling spindle at different speeds, the establishment of the thermal–mechanical model, identification of the thermal parameters, and prediction of the thermal-induced preload of bearings in the spindle. The speed-dependent thermal parameters were identified from thermal analysis through comparisons with transient temperature history, which were further used to model the thermal effects on the bearing preload and dynamic compliance of the milling spindle under different operating speeds. Current results of thermal–mechanical analysis also indicate that the internal temperature of the bearing can reach 40 °C, and the thermal elongation of the spindle tool is about 27 µm. At the steady state temperature of 15,000 rpm, the bearing preload is reduced by 40%, which yields a decrease in the bearing rigidity by approximately 16%. This, in turn, increases the dynamic compliance of the spindle tool by 22%. Comparisons of the experimental measurements and modeling data show that the variation in bearing preload substantially affects the modal frequency and stiffness of the spindle. These findings demonstrated that the proposed digital spindle model accurately mirrors real spindle characteristics, offering a foundation for monitoring performance changes and refining design, especially in bearing configuration and cooling systems.

1. Introduction

The spindle tooling system is a critical component of the machine tool, and its structural performance directly impacts processing capability and machining precision. Particularly noteworthy is the impact of thermal errors resulting from heat generation during high-speed operation, which significantly affects the machining performance of a machine tool [1]. Consequently, thermal deformation compensation technology has been extensively developed and implemented in machine tool spindles to enhance machining accuracy [2,3,4,5,6]. This technology includes a thermal error prediction model and real-time compensation through online error detection. Ouafi et al. [2] proposed a thermal error compensation model based on multi-point temperature measurements. Through this compensation method, errors induced by spindle temperature rise were reduced from 19 μm to less than one μm. Yang et al. [3] utilized five sets of displacement gauges to measure radial and axial deviations of the milling tool, accounting for thermal deformation. With these measurements, they established a mathematical relationship between thermal error and temperature rise, subsequently integrating the model into the Siemens 840D controller for error compensation. The experimental results demonstrated an impressive 86% improvement in accuracy. Pahk et al. [4] developed a thermal error model by incorporating the spindle thermal error measuring system with multiple linear regression and neural network techniques. Their findings indicate a significant improvement in machine tool accuracy by approximately 4–5 times, following real-time compensation. In a study conducted by Chen et al. [5], four points exhibiting high sensitivity to temperature rise were identified as pivotal locations for constructing a thermal error prediction model. With the model, the compensation system was devised to control displacement variation and successfully applied in milling operation. Peng et al. [6] proposed a physical-based model that accounts for the influence of thermal bending deformation on spindle axial thermal errors. Experimental verification revealed superior performance compared to typical data-driven models. Further, various thermal error modeling algorithms, such as Bi-Directional Long-Term and Short-Term Memory (Bi-LSTM) deep learning [7,8], Gaussian Process Regression (GPR) [9], Back Propagation (BP) and Radial Basis Function (RBF) neural network [10,11] were proposed to enhance the prediction accuracy. It is worth noting that thermal error models developed for specific working conditions may lack adaptability to changes in environmental and working conditions, potentially leading to a decline in prediction accuracy [12].

In other respects, the dynamic frequency response of the milling spindle plays a pivotal role in determining the machining stability of the milling machine [13,14,15]. Essentially, both the static and dynamic characteristics of the spindle system are heavily influenced by factors such as bearing specifications, bearing configuration and arrangement, and preload amount [16,17,18,19,20]. For instance, Harris [16] noted that increasing the initial preload of the bearings can enhance the spindle’s rigidity and increase damping in the spindle system to mitigate machining chatter. Similarly, Kim et al. [17] demonstrated that higher bearing preload can ensure spindle operation accuracy under temperature rise. Ciou and Lee [20] highlighted that while higher bearing preload can guarantee rotational accuracy and dynamic spindle performance, it can also lead to temperature rise and higher vibration amplitudes. The preload status of bearings and contact stiffness at rolling interfaces are subject to change due to thermal deformation, consequently affecting the vibration characteristics of the spindle system. Bian et al. [21] further discovered that the thermal-induced preload of bearings during high-speed rotation correlates with the initial preload amount. These effects were found to influence the spindle dynamics to a different extent, depending on the specific bearing arrangements [19].

Considering the effects of thermal-induced bearing preload with the temperature rise, Jiang and Mao [22] and Hwang et al. [23], respectively, proposed the concept of adjustable bearing preload and implemented it in the spindle as a replacement for the traditional fixed preload method. They conducted a Finite Element modeling approach to further explore the impact of bearing preload on the rigidity of the spindle shaft as the temperature rises. In general, the thermal effects on the stiffness of the spindle bearing system have been extensively investigated using thermal-mechanical analysis based on analytical approaches and the Finite Element method [24,25]. For instance, Zahedi et al. [24] introduced a simplified thermal–mechanical coupling model for the spindle-bearing system, simplifying the spindle housing and shaft as six-degree-of-freedom Timoshenko beam elements. This model was used to analyze the deformation and temperature increase in the spindle at various speeds, revealing a decrease in bearing stiffness and modal frequency with higher speeds. Truong et al. [25] studied the impact of thermal effects on bearing stiffness using the Newton–Raphson technique, incorporating thermal radial expansion and centrifugal expansion of rolling balls into their analytical model through Finite Element modeling. The results indicated a significant reduction in all dominant bearing stiffness coefficients at high speeds compared to static stiffness values. Additionally, integrating bearing stiffness into the Finite Element model showed that the dominant frequency and vibration amplitude of the spindle tool decreased with rising temperatures [26].

Consequently, the thermal-induced variations in the structural characteristics of the spindle system were expected to have an impact on the dynamic behavior of the milling machine and its machining accuracy [27,28,29,30,31]. Gao et al. [27] developed a Finite Element model of a spindle using the Timoshenko beam theory to analyze the thermal-induced bearing preload and stiffness. They discovered that the temperature rise significantly affected the critical cutting depth of the milling system, depending on the initial preload and spindle speed. Liu et al. [28] proposed a thermal–mechanical coupling approach based on the Finite Element method to evaluate the tool point frequency responses under thermal effects during high-speed operation. The findings indicated that the thermal–mechanical coupling effect had a significant influence on the stability lobe of the milling tool. Additionally, Li et al. [29] examined the thermal-induced preloads on combined bearings and their relationship with bearing stiffness, which were incorporated into spindle dynamic models to demonstrate the thermal effects on operational stability. Recent experimental studies by Zhang et al. [32] also revealed noticeable variations in the modal parameters of the hobbing spindle system with changing thermal conditions.

The digital twin model has become a crucial component in Industry 4.0 in recent years and is recognized as a technological strategy for integrating the virtual and real worlds [33,34,35]. According to a study by Fuller et al. [33], the majority of digital twin research has primarily focused on the manufacturing field rather than other fields such as smart cities, maritime and shipping, healthcare, and aerospace. Segovia and Alfaro [36] define the digital twin as a means of mapping the physical world into the digital realm to enable optimization, monitoring, simulation, prediction, diagnosis, and control of physical processes within a virtual space. In 2017, Cai et al. [37] integrated manufacturing information and real machine sensing signals to predict cutting accuracy and enable remote machining monitoring. Liu et al. [38] also developed a thermal–mechanical coupling model to analyze the thermal deformation of a motorized spindle unit and improve the spindle cooling system. Xiao and Fan [39] developed a digital twin application for analyzing the thermal characteristics of a motorized spindle and found that the prediction results were accurate up to 95%.

The aforementioned studies clearly showed that the temperature increase in the spindle bearing system significantly affects the bearing stiffness and the dynamics of the spindle tool during machining, leading to a potential decrease in machining accuracy. It was observed that most studies utilized an analytical approach combined with the Finite Element method. Thermal expansion of rolling bearings obtained from Finite Element analysis was integrated into the analytical solution to calculate the bearing stiffness matrix using numerical difference techniques. The analytical solution allows for the estimation of the thermal impact on bearing stiffness and spindle dynamics. However, the relationship between the speed-dependent transient thermal behavior and thermal properties of the spindle tool was not established, which is crucial for analyzing the thermal–mechanical interaction of spindle dynamics under different operating conditions. Additionally, the simplified geometrical model used in the analytical solution may not accurately depict the temperature distribution in the milling spindle and the thermal expansions of all spindle components, such as the bearing raceways and spacers. The clearance or interferences among these components substantially affect the contact force and contact stiffness at rolling balls with raceways. Therefore, inaccurate calculation of thermal expansions may consequently affect the static and dynamic characteristics of the spindle tool. The changes in frequency responses of the milling spindle due to time-varying thermal effects have not been adequately measured.

On the other hand, for a spindle tool system, the variations in dynamic characteristics due to changes in bearing stiffness and preload caused by temperature rises are notably crucial for high-speed machining in machine tools. Error compensation through temperature sensing and predictive mathematical modeling for spindle deformation is a well-established technology in practical use. Nevertheless, there remains a lack of effective methods for measuring the variations in bearing stiffness and dynamic properties of milling spindles resulting from temperature rises. In order to overcome the challenge of implementing a preload sensing system in a spindle tool, this research aimed to develop a thermal–mechanical model of a milling spindle system to estimate changes in bearing preload and dynamic characteristics of the milling spindle at varying operating speeds. Initially, a thermal–mechanical analysis model of the milling spindle tool system was created to analyze transient thermal behaviors. This model incorporated thermal parameters such as bearing heat generation and convection coefficients, which were validated through comparisons with experimental measurements from a run-in test of a milling machine operating at different speeds. By utilizing this model, the thermal deformation, bearing preload, and vibration frequency of the milling spindle under various speeds were calculated, providing valuable insights into the machining stability during operation. The milling spindle model can be utilized to simulate the interaction between thermal and mechanical characteristics under different operating conditions, serving as a valuable tool for enhancing machining performance in the presence of thermal effects.

2. Theoretical Background

2.1. Thermal Modeling

The approach of the thermal model allows a steady and transient simulation that considers the spindle parts, the bearings, and their thermal–mechanical interaction. The Finite Element method (FEM) was utilized to analyze the thermal characteristics of the spindle-bearing system, especially to simulate the temperature field and thermal deformation under thermal loads from heat dissipation in the spindle under operation. Thermal–mechanical analysis, heat generation, and heat dissipation were considered the crucial boundary conditions for thermal analysis in the abovementioned models. Additionally, the contact status in the bearing was also included in the digital model because of the difference in heat transfer capability between the continuum and joint surface heat transfer [40]. Also, the axial and radial stiffness of the bearing has an apparent influence on thermal deformation [41]. Therefore, it is important to take into consideration the change in contact stiffness at the rolling interface between rolling balls and raceways when creating the spindle bearing model.

2.1.1. Heat Generation of Ball Bearing

The aim of this study is to develop a thermal–mechanical coupling model to analyze variations in bearing preload and dynamic characteristics of spindle tools. The key is to identify an appropriate friction model capable of accurately depicting heat generation within bearings. To date, numerous friction models have been proposed with varying levels of detail in defining friction moments for calculating power losses. However, they may yield disparate results in different scenarios. Among the most commonly used models are the Palmgren, Harris, and SKF models [42,43,44]. In the Palmgren and Harris models, friction torque is divided into two components: viscous friction moment and load-dependent friction moment. The former depends on hydrodynamic losses due to internal friction in the bearing lubricant, while the latter is determined by the contact characteristics of the rolling element and raceway under external loads. However, variations in contact status due to changes in loading conditions and surface irregularities may occur, which further affect the friction moment and heat generation in the bearing [45,46,47,48]. Research by Liu et al. [45] found that friction moment fluctuates due to surface waviness during bearing operation. Rodionov et al. [47] investigated the relationship between machining errors of ball-bearing raceways and friction torque, concluding that raceway surface shape errors primarily influence fluctuations in bearing friction torque. Additionally, fluctuation of friction torque and heat generation due to contact surface waviness can be effectively reduced by improving the surface finish of bearings [48]. In other respects, the Palmgren model and Harris model offer an efficient way of modeling and calculating friction torque. In addition, minimal effort is required to implement parameters in the Finite Element thermal model. The friction model has achieved popular acceptance as an accurate method to construct the thermal model of a high-speed spindle [38,39,49,50]. Therefore, the Palmgren model and Harris model could be an appropriate selection for dealing with heat generation. Heat loss from bearing caused by bearing friction moment can be calculated through the empirical formula as follows [44]:

$$H=1.047\times {10}^{-4}\times M\times n$$

(1)

In this case, H is the heat value with the unit of energy in the form of Watts (W), n is the rotating speed of the bearing (rpm), and M represents the sum of the moments in (N mm), which can be contributed by mechanical friction torque (M_l) and viscous friction torque (M_v).

M_l is determined by bearing type load factor f₁, bearing load F_β, and the diameter of ball bearing d_m. The factor f₁ is determined by the static equivalent load f_s as a function of axial load (F_a) and radial load (F_r) components, basic static load rating (C_or), and constant of z and y. While F_β is the bearing load with axial and radial components. F_a and F_r are load components derived from the bearing preload; C_or, X₀, and Y₀ are provided by the manufacturer’s data [49]; and constants z and y are available in [43]. The following equation expresses the formula for M_l.

$$M_l=f_1·F_{\beta }·d_m= z {\left(\frac{X_0 F_r+Y_0 F_a}{C_{or}}\right)}^y·\left(F_a-0.1 F_r\right) ·{ d}_m$$

(2)

M_v is correlated with fluid friction. Although the complexity of calculating the forces in viscous friction is related to the elastohydrodynamic lubrication, the viscous friction torque can be empirically determined according to Harris’s theory [43], which can be expressed as Equation (7).

$$\begin{gathered} Mv={{10}^{-7}·f}_0 {\left(v_{0}n\right)}^{\frac{2}{3}} ·d{m}^3, for v_{0}n \ge 2000 \\ Mv={{10}^{-7}·f}_0 ·d{m}^3, for v_{0}n \le 2000 \end{gathered}$$

(3)

where n is the rotational speed in rpm, and v₀ is the kinematic viscosity of the lubricant in mm²/s. Standard values of the kinematic viscosity of lubricating grease commonly used in spindle bearing are reported as 10–32 mm²/s at 40 °C [49,51], which also changes with operating temperature. f₀ is a factor related to the type of bearing with the condition of bearings arrangements and lubrication type, in which f₀ = 2 for single angular contact bearing with grease lubrication and f₀ = 4 for duplex arrangement and with lubrication using grease [43].

2.1.2. Convective Heat Transfer Coefficient

Heat transfer coefficients can be appropriately defined, but due to the complexity of the heat problem, which is directly related to the system characteristic, this process can be challenging. In the milling spindle, the predominant methods of heat transfer are conduction and convection. The convective heat transfer of the entire spindle unit includes two types of convections: free (natural) and forced convection. The natural convection coefficient for all stationary components exposed to air under standard room circumstances is determined based on the variations in air temperature [52]. The equation is expressed as in Equation (4), where h_h is the convection coefficient in W/m²·°C and ΔT_h is air temperature changes in °C.

$$h_h=23 ∆{T_h}^{0.25}$$

(4)

Basically, the spindle head stock structure and the outer surface of the spindle housing mainly interact with ambient air by natural convection with a heat transfer coefficient given as 9.7 W/m²·°C [53]. A rough approximation of the forced convection coefficients of the rotating components can be estimated as below.

$$h_f=\frac{N_u·\mathsf{\lambda } }{D_s}$$

(5)

In which $N_u$ and $\mathsf{\lambda }$ is the Nusselt number and the thermal conductivity of the air in unit W/m·°C. $D_s$ is the equivalent diameter of the rotating components. Details for calculating this value are available in reference [54]. The Nusselt number can be obtained using the formula [55].

$$N_u=0.133·R_e^{\frac{2}{3}}· P_r^{\frac{1}{3}} ; { R}_e=\frac{u·D_s }{\mathsf{\nu } }$$

(6)

where $P_r$ is the Prandtl number of the fluid, which is given as 0.703 for air at room temperature, and R_e represents the Reynolds number of the air flowing over the rotating surface. $\mathsf{\nu }$ is the kinematic viscosity factor of the fluid, which is given as 15.06 × 10⁻⁶ m²s⁻¹ for air. u is the velocity of the rotating surface (m/s). In addition, when considering the heat transfer between the rolling balls and the lubricant, the following equation can be used to estimate the forced convection coefficient $h_b$[56].

$$h_b=0.0322·{\lambda }_b·P_r^{\frac{1}{3}}·{\left(\frac{u_s}{D_m{·\nu }_b}\right)}^{\frac{1}{2}}$$

(7)

in which λ_b and P_r are the heat conductivity in unit W/m·°C and the Prandtl number of the lubrication grease. u_s is one-third of the shaft/housing surface velocity, and D_m is the diameter of the bearing’s outer surface.

2.2. Bearing Configuration

As mentioned previously, the tool spindle must have excellent dynamic and static rigidity in the practical CNC tool machining to ensure the machining quality and accuracy. The bearing system in the spindle system is used to support the rotating shaft. It is rigid enough to withstand the cutting load through appropriate preloading, and the rotating shaft can perform precise rotational motion. Bearings for spindles are usually arranged in paired combinations of angular ball bearings, which can produce rigidity in both axial and radial directions. A preloading mechanism was installed inside the bearing to meet the high rigidity requirements to generate a preload load. The commonly used preloading methods are axial position preload and constant pressure preload, as shown in Figure 1 [49].

For the position preloading method, the axial preload was determined by the axial clearance of the inner and outer spacers, which were mounted between the front and rear bearings. Similarly, the interference fit between the bearing inner ring and the spindle mandrel causes the radial preload. Taking the DB arrangement of the bearing combination shown in the figure as an example, reducing the length of the inner spacer generates an axial clearance between the inner and outer rings, increasing the preload. Conversely, reducing the outer spacer length reduces the preload amount. The function of preload is not only to eliminate the clearance of bearing components but also to reduce the noise and vibration of the bearing. Moderate preloading can also improve the bearing rigidity. Still, it also changes the bearing contact angle, increases the deformation of the balls and the bearing ring seat, causes more significant wear, and affects the bearing life. Therefore, the preload is usually adjusted according to the recommended value in the technical manual of the bearing manufacturer.

2.3. Preload and Axial Clearance of Spacers

For the spindle module used in this study, axial positioning preload is usually used to generate axial and radial rigidity, as shown in Figure 2. The preload amount depends on the axial clearance set for the spacers, which can usually be selected according to the technical manual provided by the bearing manufacturer [49].

When the inner ring is compressed and deformed by an axial displacement “δ_a”, the center point of the arc surface of the ring shifts to the right by an equal amount. At this time, the contact angle shift of the ball is α, and the elastic deformation along the normal contact line is δ. As shown in Figure 2b,c, the relationship between the contact force Q of each ball and the axial preload F_a is:

$$Q=\frac{F_{a }}{Z Sin \alpha }$$

(8)

where Z is the number of balls, α is the contact angle at contact interface geometry. In the unloaded condition, it can be calculated based on geometry parameters using Equations (9) and (10).

$${\alpha }^o={Cos}^{-1}\left(1-\frac{P_d}{2A} \right)$$

(9)

$${A=r_o+r_i-D}_b$$

(10)

where r_o, r_i, D_b, P_d in Equations (9) and (10) are the radius of inner raceway curvature, outer raceway curvature, ball diameter, and diametric clearance, respectively. A ball-raceway normal force Q depends on this contact angle α. The normal force can be projected to radial and axial force denoted as Q_r and Q_a, respectively, as depicted in Figure 2 and written in Equations (11) and (12).

$$Q_r=Q Cos \alpha$$

(11)

$$Q_a=Q\mathrm{Sin}\alpha$$

(12)

Since the Q value depends on the contact angle, any geometrical change at contact geometry could alter normal force. Figure 2c shows the shifted inner raceway position due to rail misalignment δ_a or by force F_a. This will change the “A” value in Equation (9) to become B × D_b + δ_n. Here, B = f₀ + f_i − 1, in which f₀ and f_i are form factors, the ratio of curvature radii of the outer and inner raceways to ball diameter, respectively. Since $Q=K_h{{\delta }_n}^{1.5}$ and ${\delta }_n=B·D_b\left(\frac{\cos {\alpha }^o}{\cos \alpha }-1\right)$, then substitute to Equation (12) becomes the following equation.

$$\frac{F_a}{Z{K}_h{(B·D_b)}^{1.5} }=Sin \alpha {\left(\frac{Cos {\alpha }^o}{Cos \alpha }- 1\right)}^{1.5}$$

(13)

The above equation is derived from the geometrical change caused by the shift in the position of the inner raceway. For a given axial preload F_a, the actual contact angle α can be solved iteratively. Geometrical change may happen because of factors such as rail misalignment, external force, and preload loss by surface wearing. These factors commonly appear in the application of linear bearing. The normal contact stiffness K_n at the contact surface was estimated by using Hertz contact theory [57,58]. In this context, the symbol Q represents the contact force, δ denotes the elastic deformation occurring at the contact site, and K_h represents the Hertz contact stiffness coefficient, which is defined by the geometric dimensions of the ball groove or raceway and the material qualities of both contacting components.

$$Q= {K_h \delta }^{3/2}$$

(14)

$$K_n = \frac{dQ}{d\delta }= \frac{3}{2} {K_h \delta }^{1/2}$$

(15)

The normal contact stiffness K_n at the contact surface can be obtained as written in Equation (15). This statement clarifies that the contact stiffness in the normal direction is influenced by the contact force, which is related to the axial preload of the bearing pairs in the spindle. Consequently, any external perturbation applied to the spindle components would alter the bearing stiffness, thereby affecting the natural frequency of the system.

2.4. Modal Characteristics of Spindle Tool

Generally, the modal parameters of the mechanical system can be obtained from the frequency response function (FRF) by conducting the tapping test on the system. Figure 3 shows the typical frequency response function, which indicates that the system has several vibration modes with significant amplitude. The frequency with maximum dynamic deflection can be regarded as the primary mode that affects the system dynamics. Essentially, the dynamic characteristics or modal parameters of a milling spindle change with the preload amount. Therefore, by analyzing the change in this modal parameter, we can easily observe whether the bearing preload changes. A single-degree-of-freedom system can be used to describe the vibration behavior of a specific vibration mode, which can be represented by the following function.

$$m\frac{d^{2}x}{d^{2}t}+c\frac{dx}{dt}+kx=F\left(t\right) or (m{s}^2+cs+k)x(s)=F\left(s\right)$$

(16)

where m, c, k is the mass, damping coefficient, and modal stiffness of the system, respectively.

The relationship between the displacement and excitation force can be expressed as the transfer function,

$$FRF=\mathsf{\Phi }\left(s\right)=\frac{x\left(s\right)}{F\left(s\right)}=\frac{1}{({ms}^2+cs+k)}$$

(17)

The frequency response under cyclic excitation can be mathematically represented as a function of the excitation frequency, as demonstrated below.

$$\mathsf{\Phi }\left(i\omega \right)=\frac{x\left(i\omega \right)}{F\left(i\omega \right)}=\frac{{\omega }_n^2}{k}\frac{1}{({\omega }_n^2-{\omega }^2+i2\xi {\omega }_n\omega )}=G\left(\omega \right)+iH(\omega )$$

(18)

In the above, ω_n is the modal frequency, and ξ is the damping ratio of the system, which can be related to the system parameters by the formula $\xi =c/2\sqrt{km}$.

G(ω) and H(ω) are the real and imaginary parts of the frequency response, which can be expressed by the following equations, respectively.

$$G(\omega )= \frac{{\omega }_n^2}{k}\frac{({\omega }_n^2-{\omega }^2)}{{({\omega }_n^2-{\omega }^2)}^2+{\left(2\xi {\omega }_n\omega \right)}^2}$$

(19)

$$H\left(\omega \right)= -\frac{{\omega }_n^2}{k}\frac{2\xi {\omega }_n\omega }{{({\omega }_n^2-{\omega }^2)}^2+{\left(2\xi {\omega }_n\omega \right)}^2}$$

(20)

$$\left|\mathsf{\Phi }\left(i\omega \right)\right|=\sqrt{G^2\left(\omega \right)+H^2\left(\omega \right)}$$

(21)

$$\mathrm{\varnothing }={\tan }^{-1}\frac{H\left(i\omega \right)}{G\left(\omega \right)}={\tan }^{-1}\frac{-2\xi {\omega }_n\omega }{({\omega }_n^2-{\omega }^2)}$$

(22)

These modal parameters can be assessed by utilizing the real and imaginary components of the frequency response function of the spindle system. When the spindle system has a temperature growth, the bearing preload is affected by thermal deformation [50]. Therefore, changes in the dynamic characteristics of the milling spindle system can be detected by examining the frequency response function through vibration experiments.

3. Architecture of Digital Model of Spindle Tool

As mentioned above, for high-speed milling spindle, the change in bearing preload generated by thermal deformation from temperature rise cannot be directly measured in real-time. Due to the rise in temperature, the axial and radial deviations of the main shaft and its components are easily affected by the assembly accuracy of the machine structure. Using the sensing element for online monitoring and measurement is also difficult. Therefore, the purpose of this research is shown in Figure 4. By establishing a digital model of the spindle tool system, the thermal deformation and variation in the preload state of the spindle tool can be predicted based on the operation conditions and temperature measurements as input parameters, which are imported into the analysis model for subsequent analysis. Creation of the digital thermal–mechanical coupling model was performed following the procedures below:

• Establishment of the thermal–mechanical model with characterized thermal parameters. Initially, a series of spindle run-in experiments were carried out. Based on the thermal–mechanical model of the spindle tool with bearing modules, the temperature rise history was predicted and compared with the measured results. From the comparisons, the thermal parameters, such as the heat loss of bearings and heat convection coefficients, were appropriately calibrated from initial defaults by an analytical approach.
• Establishment of thermal–mechanical model with mechanical characteristics. The dynamic characteristics, such as the modal parameters of the spindle tool, were assessed by conducting the tapping tests on the milling machine, which was used to calibrate the preload state of the spindle bearing model at the initial amount.
• Verification of the model. Additional experiments were conducted to assess the temperature rise history and dynamic characteristics at a steady state after run-in operation. Thermal-induced deformation and preload variation in the spindle tool were predicted based on an analysis model with imported temperature loadings form thermal analysis. The predicted bearing preload and dynamic characteristics under thermal effects were compared with experiments.

Figure 4. Architecture for establishing the digital model of the spindle tool system.

Figure 4. Architecture for establishing the digital model of the spindle tool system.

4. Assessments of Thermal and Mechanical Behaviors

4.1. Run-In Temperature Rises Experiments

In establishing the digital model of the milling spindle system for thermal–mechanical coupled analysis, thermal parameters such as bearing heat flow, heat conduction properties of components, and convection coefficients are essential factors affecting the thermal behavior of the spindle in operation. The thermal parameters can be appropriately defined by an analytical approach based on the mechanical specification of the spindle system and the environmental condition. Still, they should be appropriately identified according to the thermal behavior assessed from practical conditions in operation. To this end, we conducted a series of run-in temperature rise experiments to assess the transient thermal behavior of the spindle system. The time histories of temperature growth under different operation conditions were assessed and employed to validate the digital model of the spindle tool system.

As shown in Figure 5, three temperature sensors were, respectively, mounted on the outer surface of the spindle housing near the rear bearing (T1), front bearing (T2), and machine base (T3). Another sensor was mounted on the machine frame (T4). During the experiment, the environment was controlled at a room temperature of 20 °C (±1 °C). The spindle was operated at speeds of 3000, 6000, 9000, 12,000, and 15,000 rpm, respectively. For each run-in operation, the temperature at the spindle housing gradually increased with time and reached a stable thermal state after around 3.5 h. The time histories of temperature rise at each sensor were recorded during experiments. After reaching a steady state in temperature, the spindle was stopped completely to conduct vibration experiments and assess the frequency response function of the spindle tool.

Figure 6 illustrates the history of temperature rise in the spindle tool at different points at various speeds. It can be seen from Figure 6 that the temperature rises at spindle housing near bearings are higher than at other points. The temperature rise at spindle housing near the front and rear bearings are 25.1 and 27.2 °C, 26.9 and 27.3 °C, 28.30 and 25.45 °C, and 28.34 and 28.58 °C for spindle speeds of 3000, 9000, 12,000, and 15,000 rpm, respectively. The temperature rises at the spindle housing at around 5.5 to 9.3 °C depending on the spindle speed at 3000, 9000, 12,000, and 15,000 rpm.

4.2. Impact Vibration Test

In this section, the dynamic characteristics of a milling tool are assessed by conducting the impact vibration test on the milling machine. The experiment configuration with sensor points on the milling machine with the arrow points to the applied force of the impact test is shown in Figure 7a. A solid carbide steel bar with a length of 100 mm was fixed in the tool holder. The spindle headstock was positioned at medium height, a distance of 295 mm between the tool end and the working table. Accelerometers were, respectively, mounted at the tool holder and steel bar in the X and Y orientations. The hammer was employed to strike the tool end on the opposing side sequentially to induce vibration within the spindle tool system. The frequency response function of the tool holder was derived from the recorded Fast Fourier Transform (FFT) spectrum. The modal parameters, including the damping ratio and dynamic stiffness related to the main vibration modes, were evaluated based on the measured frequency response functions.

The measured vibration response of the milling tool is illustrated in Figure 7b, which is expressed in terms of the dynamic compliance at the end of the tool holder as a function of frequency. It is found that there are some peaks occurring at the frequency of 176, 390, 792, 1070, and 1420 Hz, in which the modes below 800 Hz are associated with the structure vibration motion, and the higher modes more than 1000 Hz are mainly dominated by the vibration of spindle shaft with tool holder. The vibration mode at 1070 Hz shows a significant compliance of 0.0405 µm/N. This mode is associated with the bending vibration of the spindle shaft and is considered the primary mode affecting the machining stability of the milling tool. The modal parameters of these modes are identified by the curve fitting approach and summarized in

Table 1. Modal parameters of the spindle tool holder.

Mode	Modal Frequency (Hz)	Dynamic Compliance (μm/N)	Modal Stiffness (N/μm)	Damping Ratio (%)
1	176	0.0447	403.91	5.16
2	380	0.0167	680.60	6.48
3	675	0.0192	629.22	6.44
4	1070	0.0405	519.32	3.38
5	1420	0.0151	548.66	5.29

, which indicates that the predominant mode has a modal stiffness of 519.32 N/µm and a damping ratio of 3.54%.

4.3. Variation in Modal Parameters due to Temperature Rise

The frequency response functions of the spindle tool under the influence of thermal deformation are shown in Figure 8, which compares the difference in the frequency response of the spindle after the temperature rises at different speeds. It is seen that the temperature rises at different speeds, eventually affecting the dynamic response characteristics of the spindle tool to different extents. The temperature rise after the operation affects the main modal frequency and dynamic compliance associated with the spindle mode, which is around 1000 Hz. The modal parameters corresponding to this critical spindle mode are listed in

Table 2. Modal parameters measured after run-in operation.

Spindle Speed (rpm)	Frequency (Hz)	Damping (%)	Compliance (μm/N)	Modal Stiffness (N/μm)	Dynamic Stiffness (N/μm)
Before run-in	1070	3.380	0.0405	519.32	24.69
3000	1030	2.447	0.0407	483.43	24.57
6000	1015	2.483	0.0404	401.99	24.75
9000	1012	2.490	0.0447	361.03	22.36
12,000	1010	2.495	0.0476	338.25	21.01
15,000	1010	2.505	0.0520	311.59	19.23

and illustrated in Figure 9 for comparison. At higher speed operation, the modal frequency of the milling tool decreases, but the deflection increases. For example, after 15,000 rpm, the frequency decreases from 1070 Hz to 1010 Hz by 6%, and the dynamic compliance increases from 0.0405 μm/N to 0.052 μm/N, about an increment of 28.4%, corresponding to a decrement of dynamic stiffness of about 22%. Such changes after the low-speed operation are relatively small, but this shows that the thermal temperature rise under high-speed operation affects the dynamic characteristics of the spindle tool. As found, when the spindle temperature increases, the natural vibration frequency decreases, and the dynamic deflection increases, which implies a softening of the spindle system with increasing speed [24,25]. In addition, as shown in Figure 10, the bearing preload positively affects the modal frequency and modal stiffness with a correlation coefficient (R > 0.9), in which the bearing preload was predicted by thermal–mechanical model, reported in Section 6. This also indicates that the preload of the bearing and the rigidity of the spindle shaft are affected by deviating from the original values due to the increase in temperature. Essentially, the decrease in the dynamic stiffness due to the temperature rise results in a decrease in the machining stability of the spindle tool [25], which is expected to affect the machining performance of the milling machine.

5. Creation of a Digital Model with Thermal–Mechanical Characteristics

5.1. Finite Element Model of Milling Spindle

Figure 11 shows the milling machine (TC510, L. K. Machinery corp., Taichung City, Taiwan) and spindle system (DDS BT-30, POSA Spindle, Taichung City, Taiwan). The milling machine is composed of a vertical column structure, spindle headstock, feeding mechanism, and spindle tool module. The spindle headstock was driven to move along the column in the vertical direction. The high-speed spindle installed in the headstock was designed with a spindle shaft supported by the front and rear bearing pairs coded 7008C, arranged in a DB configuration, respectively. Referring to the bearing manufacturer [49], the geometry information of the angular ball bearings 7008C is given as follows: contact angle (15°), inside diameter (40 mm), outside diameter (68 mm), width (15 mm), ball number (18), and ball diameter (7.89 mm). Also, the basic dynamic and static load rating is 20.6 kN and 15.9 kN, respectively. The preload is an essential factor dominating the rigidity and the dynamic characteristics of the spindle tool, which basically can be preloaded at different amounts using spacers with different clearances, as shown in Figure 12. The clearance refers to the discrepancy in the axial dimensions of the outer and inner ring spacers without considering the radial clearance. In this study, a standard clearance of around 30 µm between the inner and outer ring spacers was made to generate a medium–high preload on bearings, approximately 540–590 N, according to the recommendation in reference [49].

The Finite Element model of the machine structure is shown in Figure 13. In order to consider the heat transfer across bearing components in thermal analysis and the subsequent mechanical analysis for calculating the contact forces of rolling elements, bearings and inner/outer spacers were included in the spindle unit. The bearing was modeled using 3D elements, including inner and outer races and balls. Basically, the bearing preload was determined by the interferences or clearances among these components. Therefore, these mating components were assumed to be in rough contact modes with appropriately defined interferences, which allowed for modeling the preload conditions generated by the change in interference or clearance among these mating components due to thermal expansion. In Finite Element computation, the contact force can be obtained by integrating the nodal force at each nodal point over the entire contact surface between rolling balls and raceways. The modeling of bearing in this way ensures the spindle model is more realistic as a physical unit and adaptable to the change in thermal boundary conditions.

The structural components of the system were meshed using an eight-node hexahedron and a ten-node tetrahedral element. The system consisted of a total of 879,752 elements and 3,064,098 nodes. This study mainly focused on investigating the thermal effect of the spindle tool. Therefore, the model considered here was composed of the spindle model and spindle headstocks, neglecting the vertical column structure and the feeding mechanism. With this model, the contact force at the raceways of bearings was analytically calculated by static analysis. The axial preload of bearing was then obtained as around 542–559 N when the interference between the inner spacer and inner raceway of bearings was set at 30 μm, which is comparable to the expected value of the assembled spindle. The result ensures the correct setting of the contact statute between the rolling balls and raceways in bearings for subsequent thermal–mechanical modeling. As shown in Figure 13c, the localized contact stress of the raceways was distributed in an elliptical area with a Hertzian contact pattern. On the other hand, in modal analysis, the contact stiffness at the rolling interface of the ball bearing was calculated as 347 N/μm under the axial preload of 550 N according to the Hertz contact theorem, which has axial and radial rigidity of 180 and 672 N/μm, respectively.

The milling spindle tool model was further validated by comparing the modal characteristics obtained from impact vibration tests with the computation of the Finite Element model. Figure 14 shows the fundamental vibration modes of the spindle tool headstock, including (1) pitching or nodding motion of the spindle headstock (146 Hz), (2) swinging vibration motion of the headstock about Z axis of the vertical column (160 Hz), (3) yawing vibration of the spindle headstock (370 Hz), (4) the bending vibration of the spindle shaft with tool holder at 1065 Hz, (5) higher bending mode of spindle shaft with tool holder and tail coupling at 1425 Hz. The linear guide modules on the feeding mechanism between the headstock and vertical column structure mainly dominate the first two modes. The fourth and fifth modes at higher frequencies are mainly dominated by the rigidity of the spindle bearings. It is noted that the predicted modal frequencies associated with the spindle modes agree well with those obtained from impact vibration tests. This obviously indicates that the stiffness of the spindle bearing was well implemented in the model based on the proposed modeling approach, which enabled the spindle model to exhibit the vibration characteristics of the milling machine.

5.2. Thermal–Mechanical Analysis and Calibration of Thermal Parameters

The thermal–mechanical coupling analysis was carried out to predict thermal behaviors under various operation conditions. The thermal parameters based on the previous studies were introduced as the boundary of the model under a thermal environment to validate the FE model of the spindle tool. With this, the temperature rise in the spindle tool can be predicted, and hence the mechanical characteristics with the thermal effects can be modeled. Prior to the thermal–mechanical analysis, the determination of the correct thermal parameters and boundary conditions of the spindle tool model is a prerequisite for creating the digital model that can mimic the thermal characteristics of the physical spindle tool system. Essentially, the parameters were initially obtained based on the analytical and empirical formulas in Section 2.1.1 and then imported into the digital model in a parametric mode from the thermal analysis.

In this study, the spindle bearing was lubricated with grease (SKF LGLT 2) with thicker (Lithium soap) and synthetic PAO base oil, which has a kinematic viscosity of 18 mm²/s at 40 °C [49,51]. The heat loss of the bearing was calculated as 17.6, 97.5, and 222.0 W for spindle rotational speeds of 3000, 9000, and 15,000 rpm. The forced convection coefficient generated by the lubricating greases with rolling balls and inner rings in bearing was calculated as 83.3 to 186.3 W/m²·°C for a speed of 3000 to 15,000 rpm, and the forced convection coefficient due to the rotation of the shaft was estimated as 101.6 to 301.9 W/m²·°C, changing with the increase in rotation speed from 3000 to 15,000 rpm. The temperature rises in the spindle under various speeds were obtained through thermal analysis. The thermal parameters were identified by comparing the time variation in temperature rise obtained from the experimental measurements with the results obtained from the thermal analysis. In this process, the results of the thermal analysis are affected by the heat loss and thermal convection coefficients introduced in the thermal model, which essentially were obtained through empirical calculations, but they are calibrated to fit the prediction with the measurements by considering the variation in the viscous characteristics of the lubricant, which is the main factor contributing to the heat loss and heat convection transfer in bearing [40].

The lubricant grease (SKF LGLT 2, NLGI grade2) used in bearing has a kinematic viscosity rated as 18 mm²/s at operating temperature (40 °C) and as 4.5 mm²/s at temperature (100 °C). The temperature-dependent viscosity can be expressed as an exponential form based on studies [40,59], which can be stated as the equation ν = 45.357 × e^−0.023T. Therefore, in transient thermal modeling, thermal parameters were imported into the digital model in the parametric model. The heat loss and heat convection coefficient adapt to change with the temperature of the bearing when the spindle is operated at different speeds. With consideration of the variation in lubricant viscosity, the predicted transient temperature was shown to be comparable to the measurements.

5.3. Validation of Digital Model

Figure 15 shows the transient temperature rise in the housing at the front bearing of the spindle tool under speeds of 3000, 9000, 12,000, and 15,000 rpm, which were compared with the measured temperature rise. As observed, the temperature rise histories predicted under different spindle speeds agree with the measurements. The observation verifies the accuracy of the digital model of the spindle tool system in predicting thermal behaviors.

With the modeling of thermal behaviors, the temperature rises in bearings at a steady thermal state were predicted, as shown in Figure 16 and

Table 3. Comparisons of the temperatures at bearing sites between measurements and predictions.

RPM	Measured Temp. at Bearing Housing (°C)		Predicted Temp. at Bearing Housing (°C)		Predicted Temp. at Ball Bearing (°C)	Temperature Rise (°C)
	Front	Rear	Front	Rear	Rear
3000	27.17	25.10	26.63	25.64	29.5	9.5
6000	27.96	28.07	27.34	27.25	31.5	11.5
9000	27.03	27.00	27.91	27.57	34.5	14.5
12,000	28.34	28.45	28.45	28.75	37.9	17.9
15,000	28.35	28.58	28.45	28.62	40.7	20.7

. The bearing temperature is around 29.4, 35.3, 37.9, and 40.8 °C for spindle speed at 3000, 9000, 12,000, and 15,000 rpm, respectively, corresponding to the temperature rise of about 9.4, 15.3, 17.9, and 20.7 °C. The temperature distributions within the spindle at a steady thermal state were predicted and illustrated in Figure 17.

Results of run-in experiments and thermal modeling have shown that the temperature increases with the spindle speed change, and the value of bearing heat generation ranges from 23 W at 3000 rpm to 198 W at 15,000 rpm. From the validation of the thermal model of the milling spindle, the thermal characteristics, such as bearing heat loss and forced convection coefficients with change in spindle speed, can be appropriately identified, as listed in

Table 4. Thermal parameters.

Spindle Speed (rpm)	Bearing Heat Loss (W)	Forced Convection around Spindle Shaft (W/(m2.°C)	Forced Convection at Bearing Rotating Surface (W/(m2.°C)
3000	23	101.6	83.3
6000	52	162.4	117.8
9000	87	213.7	144.3
12,000	138	259.6	166.6
15,000	198	301.9	186.3

. The thermal parameters can be accurately related to the spindle speed by curve fitting of the data, as shown in Figure 18. Equations for bearing heat loss Q_b (w) and forced convection coefficient H_f (W/(m²·°C) can be expressed as the function of speed S (rpm) shown in Equations (23) and (24), respectively.

$$Q_b=6.0\times {10}^{-7}·S^2 + 0.0034·S + 7.8$$

(23)

$$H_f=0.4508·S^{0.6767}$$

(24)

6. Application of Thermal–Mechanical Model

6.1. Prediction of Thermal Deformation

A typical deformation pattern of the spindle tool with the headstock under a stable thermal state is illustrated in Figure 19, which shows the maximum deformation occurring at the tool tip. The axial deformation of the spindle tool due to temperature rise under speeds of 3000, 6000, 9000, 12,000, and 15,000 rpm is predicted at around 11.4, 14.7, 18.0, 23.1, and 27.2 μm, respectively.

The thermal deformation of the spindle tool is more significantly affected by the increasing spindle speed. This phenomenon was clearly illustrated in Figure 20, in which the thermal deformation increases nonlinearly with spindle speed. The higher the speed of the spindle in operation, the higher the temperature rise inside the bearing, hence inducing more axial elongation of the spindle tool, which was believed to affect the dimensional accuracy of the workpiece. The axial deformation δ_t (µm) of the spindle tool can be well related to the spindle speed S (rpm) by the following equation.

$${\delta }_t=9.3376\times e^{7\times {10}^{-5}\times S}$$

(25)

6.2. Prediction of Thermal-Induced Bearing Preload

As revealed in the previous section, there is an axial deformation due to the rise in temperature under high-speed operation. The thermal elongations of the spindle shaft, spindle housing, and the inner and outer spacers are different, which in turn change the clearance of the inner and outer spacers between the front and rear bearings, as shown in Figure 21. The axial deformation will change the contact state, such as the contact zone and contact force between rolling balls and raceways, altering the bearing preload amount and rigidity.

In this case, the bearings were arranged in a double-paired DB arrangement, in which the negative clearance between the outer ring spacer and the inner ring was made to generate the preload in bearings. According to the thermal–mechanical coupled analysis, the axial preload of bearing under thermal deformation can be calculated based on the contact force of the rolling balls on the raceways. It is noted that the spindle shaft with the inner spacer has more elongation than the outer spacers with spindle housing. Hence, it will lessen the interference between rolling balls and raceways and reduce the preload on rolling balls. The variation in bearing preload with changing of the spindle speed is shown in Figure 22a, which shows the decrease in preload with increasing speed in a nonlinear way. The fitted expression of preload P_r (N) with respect to spindle speed S (rpm) is given as

$$P_r=1920· S^{-0.191}$$

(26)

The initial preload of the spindle bearing is about 550 N, while it decreases to 365 N when operating at a speed of 6000 rpm, corresponding to a reduction of 35.6% in preload. If the spindle speed increases to 15,000 rpm, the bearing preload will decrease to 325 N, about a reduction of 40%. The reduction in bearing preload substantially reduces the contact stiffness and bearing rigidity by 16%. This implies that the dominant vibration frequency of the spindle with significant deformation will be affected by a decrease in temperature rise under high-speed operation, as was observed in the experiments conducted. The frequency responses of the spindle tool under different speeds were further predicted by harmonic analysis based on the proposed model, as illustrated in Figure 22b. It is noted that the predicted dynamic compliances at the dominant frequency around 1000 Hz decrease with the increase in spindle speed. The frequency responses affected by thermal effects on bearing preload are well comparable to those observed in experiments (Figure 8).

The current results clearly demonstrate that the digital model can effectively mimic the thermal–mechanical behavior of the milling spindle tool under different operation speeds. The change in bearing preload and the modal parameters dominating the machining stability under thermal effects can also be accurately predicted by using the digital model instead of the experiment measurements. As a whole, the model developed with thermal–mechanical coupling effects and thermal parameters identified from experiments can provide the basis for the improvement of the mechanical performances of the milling spindle system. Some implications in practical application can be drawn as follows: (1) Design of bearing configuration. Front and rear bearings in the spindle can be configured in different configurations, such as front-to-front or back-back arrangements [49]. As known, bearing preload was susceptible to decrease or increase due to thermal expansion of the spindle components. For different bearing configurations, the influence extent of such thermal effects on the dynamic characteristics and machining stability of the milling spindle can be evaluated following the proposed model, which can help engineers make appropriate evaluations at the design phase of the spindle. (2) Initial preload setting. Preload was normally set at a specific amount according to the machining purpose, light preload for high-speed finish machining or heavy preload for rough machining, but it causes a rise in bearing temperature and variation in bearing preload to different extents. With this model, an appropriate evaluation of the variation in dynamic characteristics of the spindle under temperature rise can help to determine the optimum preload.

7. Conclusions

This study aimed to develop a comprehensive thermal–mechanical model for the spindle tool system to investigate transient thermal behavior and its impact on the dynamic characteristics of the milling machine. The temperature rise histories and the dynamic characteristics of the spindle were first measured by conducting run-in tests on a milling machine and by tapping tests, which were used as model verifications. Some conclusions are drawn below:

• Experimental findings revealed that temperature rises in the spindle unit affect the modal frequency and dynamic compliance to change with the increase in spindle speed, indicating a decrease in bearing stiffness due to thermal effects at higher speeds.
• Thermal parameters played a crucial role in ensuring the effectiveness and accuracy of the model in mimicking the thermal–mechanical behavior of the milling spindle. The study successfully identified temperature-dependent thermal parameters by comparing results from transient thermal analysis with experimentally measured temperature rise histories. Notably, thermal parameters varied with temperature rise. In addition, the thermal-induced bearing preload can also be predicted by the digital model without the need for sensor implementation in the spindle.
• Furthermore, results from thermal–mechanical coupling analysis highlighted that thermal expansion of spindle components and headstock led to a decrease in bearing preload with increasing speed. At a speed of 15,000 rpm, bearing preload exhibited a reduction of approximately 45% from its initial state, consequently increasing the dynamic compliance of the spindle tool, as confirmed in experiments.
• These findings demonstrate that the spindle model developed in this study accurately reflects the characteristics of real milling spindles. Moreover, this model, incorporating a milling spindle with a machine frame, provides a valuable foundation for monitoring changes in spindle performance and improving the design of both spindle and machine frame structures. This is particularly relevant for optimizing bearing configuration, feeding mechanisms, and cooling systems to mitigate thermal effects.