An improved thermal performance modeling for high-speed spindle of machine tool based on thermal contact resistance analysis

The temperature field distribution has a significant influence on structural performance, thermal deformation, and thermal error compensation. To improve the prediction accuracy of the temperature distribution of the spindle system, a comprehensive model considering the contact thermal resistance (TCR) of the interfaces was established to analyze the thermal performance of the high-speed spindle system in the present work. An elastoplastic contact model was used to calculate the contacting areas and loads of interfaces, which were employed to establish the contact thermal resistance model of the primary interfaces of the spindle, such as bearing rings and tool holders. Based on the TCR parameters, a finite element analysis (FEA) model was proposed to analyze the temperature distribution of the spindle system. And a temperature test experiment was set up to verify the accuracy of the FEA model. The results show that the relative errors of representative test points were all less than 5%, which means the established model can appropriately reflect the temperature field distribution of the spindle.


Introduction
As the core functional component of precision machine tools, the high-speed spindle system is crucial for the implementation of high-speed and high-accuracy machining [1,2]. However, with the improvement of running speed, the influence of vibration and thermal becomes more and more significant. Especially, the error is generated by thermal deformation. It is shown that the errors induced by heat account for nearly 40-70% of the total error of a precision machine tool [3]. Therefore, it is crucial to study how the heat influences the accuracy of the machine tools.
Usually, a thermal-mechanical coupling model is used to analyze the temperature distribution and deformation of the main parts of machine tools. Those models can be used in both the design and working stage. During the design stage, the thermal-mechanical model can predict the deformation of the error of spindle tools under the action of specific conditions [4][5][6]. Those models can also be used for the thermal error compensation [7,8] and thermal error control techniques [9,10] in the working stage. The distribution of temperature is the foundation of the thermal-mechanical coupling model. There are usually two foundation methods to solve the temperature field, the thermal grid method [11][12][13], and the finite element method [14][15][16].
For high-spindle system, the heat generated by the friction of bearing and built-in motor flowed through parts causes linear and nonlinear thermal expansions, which distorts associated structures and produces deformations. What is more, the spindle system is complicated system composed of many components with many kinds of interfaces. Their thermal and mechanical behaviors are not only determined by these components but also by the characteristics of various interfaces [17]. When the heat flowing contact areas, such as interfaces between the bearing and the shaft, the bearing and the bearing support, and the ball screw and nut, the temperature changes are always nonlinear and produce nonlinear thermal deformations, which adversely affect the stiffness of interfaces [18]. Therefore, the thermal error and contact thermal resistance (TCR) on the interface are among the most important objects to consider for designing machine tools.
To acquire the value of TCR, both theoretical and experimental methods are used. Compared with the experimental method [19,20], the theorical methods [21][22][23] got more scholars' attention because of their broader applicability. There are three key parameters to calculate the TCR: pressure distribution, really contact area, and thickness of contact region. All those parameters and their relationship can be deduced by the rough surfaces contact theories. For a given load condition, the real contact area and the displacement approach amount can be confirmed by the statisticsbased method [24,25] or fractal-based method [26,27], and then get the value of TCR. So, the contact models determine the validity and convenience of TCR calculation. No matter which method is used, the basic work deals with the contact problem between a single asperity and a flat rigid surface. Three-stage contact models, elastic, elastoplastic, and plastic stage, are now wildly used [28,29].
Based on an advanced contact method, which analyzes the relationship between actual area and deformation of contacted rough surfaces under the specific load, a comprehensive model for the thermal study was proposed in the present work. Together with thermal boundary conditions, such as heat generating, and thermal convection, the proposed model took three kinds of typical interfaces (bolt joints, ring joints, and tapered joint) into account. The TCR of those interfaces were acquired according to the distribution of pressure of joints under different loaded characteristics. And experiments were carried out on a vertical machine center for 5.5 h running and 1.5 h stopping to verify the proposed thermo-mechanical model by comparing the temperature and thermal error data.

Rough surface contact model
The contact of two rough surfaces can be simplified as a rough surface in contact with a rigid plane. Based on the Hertz contact theory and statistical method, Greenwood and Williamson proposed a contact analysis model, named as GW model, which assumed that the height of asperities obeyed a certain distribution (usually Gaussian distribution) and that asperities had the same radius of curvature [30]. The deformation for a single asperity experienced elastic or plastic stage in this model. However, researchers believed that there was a transitional phase between those two stages. Through the method of finite element analysis, Kougt and Etsion discovered that the asperity undergoes an elastic-plastic mixed deformation stage before it began to enter the plastic deformation stage [31].
In the present work, the contact behavior of asperities would be divided into three stages: elasticity, elastoplastic, and plasticity. When the deformation is smaller than critical elastic deformation ω e , the asperity is in the stage of elasticity. When the deformation ω is bigger than critical elastic deformation ω e , the asperity is in the stage of elasticity, elastoplastic deformation. And when the deformation ω is bigger than 110 * ω e , the asperity is in the stage of totally plasticity. The critical elastic deformation ω e is where k μ is the average contact pressure coefficient [32], related to Poisson's ratio of softer materials ν, and k μ = 0.46 45 + 0.3141ν + 0.1943ν 2 .
According to the Hertz contact theory, the relationship between the elastic contact area and the elastic contact load with the deformation ω can be established as following.
. E 1 and E 2 are the shear modulus of the material of contact surfaces. ν 1 , ν 2 are the Poisson's ratio of the material. R is the average radius of all asperities.
2. ω e < ω < 110*ω e where κ 1 (ω) and g(ω) are the polynomial expression which can be constructed as the function of the independent variable ω [33]. And that is where H is the hardness of the softer material in the two contact surfaces.

Average contact pressure of interfaces
The real contact area is often determined by the contact pressure. Different methods should be taken to deal with different types of contact surfaces. The contact pressure of the bolted joints can be obtained by dividing the given bolt preload by the nominal contact area, which can be easily acquired from the geometrical dimensions. Another two types of joint, spindle-handle and shaft-ring, with different load condition should be processed, respectively.

The interference fit contact pressure
As shown in Fig. 1, the contact pressure of the joint surface of the interference fit, such as shaft and inner rings of bearing, is determined by the interference.
The average contact pressure could be calculated according to the following formula where E 1 and E 2 are the modulus of the material of rings; Δ denotes the interference fit; μ 1 and μ 2 denote the Poisson's ratio of the bearing inner ring and shaft journal, respectively; D denotes the contact diameter; D 3 denotes the inner diameter of the shaft journal; and D 1 denotes the outer diameter of the bearing inner ring.
The joint surface of the spindle shank is a tapered shape, and the contact load is determined by the broach force. As shown in Fig. 2, under the action of axial broach force F p , the force balance equation of the shank with taper θ iswhere N is the normal pressure of contact surface and F f is the friction force between spindle and shank, respectively.
where N is the normal pressure of contact surface and F f is the friction force between spindle and shank, respectively.
According to the relationship of normal pressure and friction, the above equation can be rewritten as where μ is the coefficient of friction. For the most cases, this coefficient can be assumed to be in the range of 0.12 to 0.2.
So, the average surface pressure of the cone is

Thermal contact resistances
As shown in Fig. 3, actual contact surfaces are always pointcontact or an uneven small area contact. Compared with perfect contact surfaces, there is an additional thermal resistance during the heat conduction process, which is the contact thermal resistance (TCR) [34]. TCR is related to the actual contact state and can be calculated based on contact theory. TCR is fundamentally decided by heat transfer coefficient, height of gap and really contact area, which can be formulated as where L δ and A r are related to contact conditions (such as material, pressure, and roughness), and can be calculated by rough surfaces contact theories.
To simplify the analysis, the thermal contact conductance (TCC), which is the reciprocal the of TCR, is used to calculate the TCR of the interfaces.
For a single asperity, L δ is the distance from smooth rigid plane to average reference plane, and A r is the hertz contact area. So, the TCC of a single micro-contact asperity in different contact situations will be where z is the height of the asperity relative to the average line.
According to Sect. 2.2, the total TCC of the whole joint H e can be calculated by 3 FEA model

Geometric model
The thermal performance of the spindle system was analyzed by the finite element software ANSYS Workbench. As shown in Fig. 4, substructures such as spindle, box, test mandrel, and bearings were built and assembled in ANSYS. The spindle was supported by two pairs of bearings mounted on the spindle box, and the way arranged bearings that the front had three rows and the rear two rows. The test mandrel was fastened to the holder by pulling force provided by discspring. Hex dominant method was used to mesh geometries, and the number of the nodes and elements were 85, 313 and 25, 788, respectively. The connection between parts of the spindle system was defined by contact pair in ANSYS Workbench. The TCR of interfaces can be represented by the parameter "H" in this item. There were three kinds of interfaces: plane joints, bearing joints, and spindle-holder interface. The bearing joints included two types: the joint between inner ring and shaft, and the other was the joint between outer ring and box. The major materials are shown in Table 1.

Heat generation of bearings
The primary heat generation of the system was caused by the cutting process, and the friction of balls and races [35]. Assumed that most cutting heat was taken away by coolant and chips, and the heat generated by bearings was the dominant cause of temperature changing.
The heat generated by a bearing can be computed as where H f is the heat generated power (W), n is the rotating speed of the bearing (rpm), and M is the total frictional  torque of the bearing (N.mm). The total frictional torque M consists of two parts. One is the torque M 1 due to applied load and the other one is the torque M 2 due to viscosity of lubricant. And, Here, f 1 is a factor related to the bearing type and load, p 1 is the bearing preload (N), and d m is the mean diameter of the bearing (mm).
Here, f 0 is a factor related to bearing type and lubrication method, and v 0 is the kinematic viscosity of the lubricant.

Thermal convection coefficient
Rotating parts exposed to the air, such as the spindle, the taper shank, and the test rod, caused the air near those parts to be wound, thereby causing forced convection heat exchange. The calculation between the rotating body and the air using the Nuer Shelter equation [36]. The heat transfer coefficient for convection h is defined as Here, λ a is the thermal conductivity of the ambient air, N ua is the Nusselt number, and d is the equivalent diameter of the cylinder.
N ua is computed from the Reynolds number R e , and the Prandtl number P r is based on different convection conditions. For this research, the following equation is used and, where ω is the angular velocity of the rotating body, and v f is the kinematic viscosity of the air.
The thermophysical properties of dry air at 20 °C were as follows: λ air = 0.0259 W/(m K), v f = 15.06*10 −6 m 2 /s, Pr air = 0.703. For free convection around stationary surfaces, the heat transfer coefficient for convection h = 9.7 W/(m 2 K) was used in this paper.
From the above data combined with the size of the rotating shaft of the machine tool spindle, and the heat transfer coefficient of the inner and outer surfaces of the rotating body at each experimental speed can be obtained as shown in Table 2.

FEA result
According to the thermal boundary conditions given by Sect. 3, the spindle, bearing, and spindle box's transient temperature field distribution was obtained by the established FEA model. There were two stages in this procedure. In the first stage, which lasted for 5.5 h, the spindle was running, and the heat generated by bearings flowed to other parts. What is more, forced convection occurred on the surfaces of rotating parts, such as the node of the spindle and the test mandrel. The second stage lasted 1.5 h to simulate the natural heat dissipation process after the spindle stopped. At this stage, surfaces had a natural heat transfer, and there was no heat source anymore.
The temperature distribution after the spindle speed was operated at 6000 r/min for 5.5 h (at the end of stage I) was shown in Fig. 5b. It can be seen that the parts with apparent temperature rising were mainly concentrated in the region of the two sets of bearings, and the temperature in the region away from the bearing was nearly unchanged because of the heat source passing.
The bearing friction was the primary heat source of the spindle system, and the internal heat dissipation was difficult. So, the temperature rising was most obvious near bearings. Specifically, the spindle temperature near the front bearing group was up to 55 °C, and the temperature of the spindle near the rear bearing group was near to 48 °C. The total temperature rising of the front and rear bearing sets was 27.04 °C and 21.56 °C, respectively. Although the test mandrel was not far away from the front bearings, the temperature hardly increased. Because surfaces exposed to the outside had natural heat exchange or forced convection with the surrounding air, the heat from bearings was then dissipated into surroundings immediately. So, the temperature rise was relatively small. Figure 5b shows the temperature variation of a small region near the front bearing. From the variations of the temperature-time curves, it was apparent that there were three periods throughout the time frame, which were marked as I, II, and III, respectively. The temperature increased rapidly during period I, and the temperature of this point raised by nearly 24 °C. In period II, the temperature changed slowly, lasting almost 3 h until the spindle stopped, which indicated that the system tended to reach thermal equilibrium. In period III, the temperature began to decrease, and the rate of decline was faster than later. Figure 6 shows the experimental measuring scheme and equipment of the spindle system. The transient temperature field of the spindle was measured by 18 thermocouples and a set of the data acquisition system. During the experimental procedure, the thermocouples were attached to the stationary surface of the spindle head, and the spindle was run for 5.5 h, then stopped and cooled for 1.5 h. During the whole process, the temperature data of thermocouples were collected every 10 s. Figure 7 shows the variation of temperature over time, which presented a nonlinear relationship. In the first period, the temperature began to rise immediately, and then, its rise began to slow until tending to thermal balance in the second period. The temperature of the whole system was the same as ambience before the machine worked, and the heat generated by bearing conducted to other parts resulting the temperature rising once the machine tools working. As the increasing of the temperature, heat conduction got slower due to the change of temperature gradient. In the third period, the temperature began decrease until the ambient temperature. That was because the bearings no longer heated up, and air convection continued the outside surfaces. However, several abnormal phenomena were also observed. The temperature of some sensors had a slight rise and then started to fall. When the spindle stopped working, the temperature of the inner parts was higher than that of the surfaces exposed to the air, which means there was an uneven temperature field, and the heat continued flowing to the outside surfaces. What is more, the situation was not the same for all the 18 sensors. For example, as nearing bearings, t1 ~ t4 had higher temperatures than other sensors. The sensors away from the heat generation site, such as t17 and t18, had lower temperatures.

Discussion
To verify the simulation results, the experimental data of the spindle was compared with the corresponding position simulation temperature data, and the result was shown in Fig. 8. The simulation results were consistent with the trend of the experimental results.
The four observation points with the different temperatures rising at the front end of the spindle box were compared with the data from the simulation analysis at the corresponding position. Here, we compared two sets of data: the temperature of spindle stopping and the final temperature. As shown in Table 3, at the moment of spindle stopping, the relative error of simulation with TCR was no more than 3% for the four observation points. But the relative error of simulation without TCR varied between 5 and 10%, which were higher than that of simulation with TCR. Because the TCR prevented heat conduction outward, less heat was transferred to the surface. So, the surface temperature of the spindle of simulation with TCR will be lower. Due to uncertain factors such as processing or assembly of parts of the spindle, the actual TCR was slightly higher than that of theoretical analysis results, which resulted in the temperature of the experiment result were all lower than that of two kinds of the simulation method.
After the spindle stopped, the temperature of the simulation dropped faster than the experiment result, which can also be explained by the TCR. During the spindle stopping period, the temperature change of test point #1 was 10.7 °C, 14.14 °C, and 17.49 °C for the experiment, simulation with TCR, and simulation without TCR. Similar situations can be observed in test points #3, #10, and #11. So, the temperature raised slowly and declined slowly during the existence of the TCR. What is more, the relative error of simulation with TCR was higher than that of simulation without TCR in test points #1 and #3, which seems the improved model had lower accuracy. However, the temperature change of simulation with TCR was still smaller than that of simulation without TCR, which meets the above temperature changing conclusion.

Conclusion
An FEA method to analyze the thermal performance of the spindle system, which considered the effect of primary interfaces, was proposed based on the elastoplastic mixed contact model. The model can be used to accurately predict the temperature distribution of the spindle at the design stage without making prototypes, which was very useful for speeding up product development. Experiments verified the validity of the modeling method. Based on the results and analysis, the following conclusions were obtained.
1. Based on the elastic, elastoplastic, and plastic contact model, the TCR of a single point was analyzed by theoretical method. The mean TCR of primary interfaces was calculated according to the surface topography and contact load. The result showed that the value of TCR was mainly affected by the actual contact area of rough surfaces. So, the contact algorithm of rough surfaces was crucial for calculating TCR. 2. An improved simulating model was proposed to analyze the temperature distribution of the spindle box. In addition to considering the thermal convection conditions and heat generation of bearing, the proposed model focused on the influence of the primary interfaces' TCR.
The results showed that the temperature changed quickly initially, and then the temperature will gradually achieve balance. 3. To verify the accuracy of the proposed model, an experiment was carried out to measure the temperature variety when the spindle was running. The comparing results showed that the relative error of the proposed model was below 5% for the spindle stopping temperature. Mean-while, the relative error of the model without considering the influence of TCR was 9.6%. The comparing results confirmed the accuracy of the proposed model.
Author contribution Bing Fang conceptualized the idea, designed the methodology, and undertook data curation, investigation, formal analysis, project administration, manuscript writing, reviewing, and editing; Tianqi Gu undertook writing, original draft, review, and editing; Mengna Cheng and Dapeng Ye undertook experiment, data curation and investigation, and acquired resources.