Transient thermal error modeling of a ball screw feed system

The thermal behaviors of the ball screw feed system significantly affect the machining accuracy of the machine tool, which is widely valued by users and manufacturers. In the current ball screws thermal error models, the moving heat source at nuts is simplified to a fixed heat source to simplify the calculation, which can decrease the model’s prediction accuracy. A novel thermal error modeling method for the ball screw feed system is proposed by this paper based on the FEM. In the proposed method, the moving heat source at the nut is loaded into the model by changing the grid distribution at different time steps. Meanwhile, a series of thermal experiments are carried out to collect the thermodynamic data with the machine tool’s feed drive system operating at a different speed. And then the multi-objective optimization method is introduced to obtain the heat generated by the nuts and the convection coefficient between the air and the screw. The experimental results show that the thermal error of the ball screw feed system is different when the feed system runs in different intervals, even if the temperature rise at measuring points is the same. In the final, experiments under different conditions are carried out, and the proposed model is accurate enough to predict the thermal error of the ball screw feed system.


Introduction
With the increasing demand for machining accuracy and efficiency, the control of thermal error has been more and more valued by machine tool manufacturers and users [1]. Though the feed system with full closed-loop control has better thermal stability, the manufacturing cost of the machine tool will be significantly increased due to the application of grating rulers. Therefore, the ball screw feed system with semiclosed loop control is still irreplaceable. When semi-closed loop control is adopted, the thermal deformation of the ball screw will directly decrease the machine tool's positioning accuracy, which can also cause the machining accuracy to decrease. Therefore, various thermal models have been proposed by many scholars to predict the thermal behavior of the ball screw feed system and to improve positioning accuracy by some optimization design.
There are two main methods of modeling: databasebased methods and model-based methods [2,3]. Recently, various intelligent learning algorithms represented by neural networks trend to be applied in the machine tools' thermal error modeling, with the successful application of them in various fields. The single-directional multi-layered neural networks with error back-propagation (MLP), radial basis function neural networks (RBF) and Kohonen networks were adopted by Rojek et al. to establish the ball screw feed system thermal error model respectively. And then experiments were performed to validate the accuracy of the neural network modeling method for the ball screw feed system, and the prediction accuracy of the three neural network models was also compared [4]. A new inverse random model was proposed by Li et al. to predict thermal error. The randomness of influencing factors was taken into consideration by this model through the combination of stochastic theory, genetic algorithm and radial basis function neural network (RBFNN) [5]. The model took the temperature at the nut and the bearings as the input and was validated by the experiment under a single working condition. Elman neural networks (ENs) were employed by Yang and Xing [6] to carry out the thermal error modeling of the high precision feed system. The differential evolution (DE) algorithm was used to optimize the initial weights and thresholds of the ENs. The feeding system's thermal error is reduced from 1.73 to 0.88 μm with the online compensation. In addition to the neural networks, the multiple linear regression [7], support vector machine (SVM) [8], and grey rough set theory [9] are also widely adopted in the feed system thermal error modeling. The temperature at the nut and bearing where the data is easier to collect is usually taken as the input to these data-based models, due to the difficulty of measuring the screw temperature distribution directly. However, when the nut is running at different intervals on the feed system, it may turn out the different thermal error, although the temperature rise at the measuring point is the same, which is caused by the different temperature distribution of the lead screw. In this case, the data-based model cannot predict the feed system thermal error accurately.
When the model-based methods are applied for modeling, analytical or numerical method can be chosen based on the model's complexity [10]. An analytical model was established by Shi et al. to predict the thermal error of the ball screw feed system based on thermodynamic parameters calibrated by the experiments [11]. Feng and Li et al. [12] proposed an analytical transient temperature model of the ball screw feed system in the warm-up and cooling period respectively. And then the thermal error of the system can be predicted based on the thermal model, which was validated by the cutting compensation experiments. Finite difference methods (FDM) and finite element methods (FEM) are currently widely used numerical modeling methods for machine tool's thermal errors [13]. The feed system thermal error was divided into two parts by Liu et al. including internal heat source induced error and environment temperature-induced error. And the linear regression method and finite difference methods were adopted respectively to establish thermal error models [14]. Monte Carlo (MC) simulation-integrated FEM method was used by Li and Zhao et al. to calculate the heat generated by the heat source of the feed system. And then a numerical model for thermal error of the lead screw was put forward based on FDM [15]. Mian et al. [16] established a thermal deformation model of machine tools based on FEM, which was used to explore the influence of ambient temperature on its accuracy. Min and Jiang [17] developed an integrated thermal model with the aid of the FEM to analyze the temperature distribution of a ball screw feed drive system. The thermal contact resistance between the bearing and its housing was taken into consideration in this model. A modified lumped capacitance method (MLCM) was developed by Zhu et al. [18,19] based on thermal behavior model of the ball screw system to estimate its thermal error and the effectiveness of the air cooling system, which was validated by three different work conditions [20].
The heat generated by the nut is a moving heat source in the feed system, which is difficult to be loaded in the simulation model accurately and is also one of the key factors to affect the accuracy of the thermal error model. In the above models, the moving heat source at the nut is transferred to the uniform heat flux loaded on the whole stroke, the position of which does not change with time. Thus, predicted errors of these model can be induced, because the boundary conditions are inconsistent with the actual situation.
A transient thermal error model of the ball screw feed system is developed by this paper based on the FEM, to overcome the existing models' shortcomings. A novel method to realize the moving heat source at the nut loading is proposed in this paper by changing the grid at each calculation time step in the proposed model. Meanwhile, the multi-objective optimization method is applied to calibrate the heat generated by the nuts and the convection coefficient between air and the screw. Finally, the model is validated by series of experiments.
The structure of this paper is arranged as follows: The ball screw feed system thermal error model based on FEM is established in Section 2. Meanwhile, the method for loading the moving heat source at nut into the model is also introduced in Section 2. Section 3 briefs the calculation methods for the boundary conditions of the established model. The heat generated by the nuts and convection coefficient between air and the screw is calculated based on the multi-objective optimization method in Section 4, and then a series of experiments under different conditions are conducted to validate the proposed model. In the final, the conclusions and prospects of this study are offered in Section 5.

Thermal error modeling of ball screw feed system
The screw can be simplified to a one-dimensional rod model, because of its large aspect ratio and simple symmetrical structure, as shown in Fig. 1. In the proposed model, the q 1 , q 2 and q 3 are the heat flux generated by the front bearings, nut and rear bearings respectively, which will be input the screw along its radial direction as the second boundary. Meanwhile, the heat q 2 is a moving heat source and reciprocates along the axis of the screw at the speed of v. The convective heat transfer, with convective coefficient h, is conducted between the rest of the screw cylindrical surface and the air. The simplified one-dimensional heat conduction equation of the ball screw is as follows: The boundary conditions and initial conditions are written as where, k, ρ and c p are the thermal conductivity, density and specific heat capacity respectively, T(x,t) is the temperature at the x position of the screw at time t and T evn is the environmental temperature. Γ qi (i = 1,2,3) represents the boundary of the heat flux input, while Γ h is the boundary of convective heat transfer with air.
The analytical solution of the equation is difficult to be obtained, due to the loading position of heat flux q 2 in this model changing with time. Therefore, the FEM is applied in this paper for the screw temperature modeling, as it has higher solution accuracy and is more suitable for irregular collection areas than the FDM.

Screw temperature modeling based on FEM
The screw temperature field T(x,t) can be approximated by where N i is the shape function, n is the number of the element nodes and T i is the time-dependent node temperature.
The following equation can be obtained by integrating Eq. (1) in a finite element based on the Galerkin Method.
The matrix form can be written as follows: Two-node one-dimensional linear element is selected in the proposed model, as shown in Fig. 2; the shape function of which can be expressed by The differential of the temperature function T(x,t) to the space coordinate x be expressed as follows: Then, the capacity matrix [C] e , stiffness matrix [k] e and load vector{f} e of each element can be calculated as follows: where l is the length of the element, A and P are the area and perimeter of the element's cross section.
The equation systems can be obtained by assembling the element matrix, as follows.
The one-dimensional model of the ball screw feed system The above equation can be rearranged as follows: The backward difference method (θ = 1) is applied in this paper. Then, the transient temperature distribution of the screw can be predicted based on Eq. (13), after the initial conditions and thermal boundary conditions given.

Moving heat sources loading and thermal error modeling
The position of the heat flux generated by the nut is constantly changing when the feed system is running. A novel method proposed by this paper to realize the loading of the moving heat source at the nut by changing the grid distribution on each time step. The screw is divided into five parts along its axial direction as shown in Fig. 3. where x(t) is the position of the nut at different times, L z is the width of the bearing, L n is the width of the nut and N i (i = 1,2,…,5) is the number of elements in each part. The different positions of the nut x(t) can be calculated at each time step. And then the mesh of the model will be updated at the corresponding moment. The complex work condition of feed system can be taken as the combination of its single reciprocating motion in a certain interval, which is defined by a work condition vector a = [p 1 , p 2 , v, t] T , where p 1 , p 2 , v, t are the starting point, ending point, feed rate and running time respectively under the single work condition. And then the x(t) in the current time step will be calculated through the vector a.
The calculation process of the x(t) is described by Fig. 4. The function mod() is the modulo operation that returns the remainder or signed remainder of a division. And the floor() is the function that can give as output the greatest integer less than or equal to the input. As shown in Fig. 4a, if the expression mod floor v • t∕ p 2 − p 1 , 2 equals to 1, the nut will run in the negative direction of the x axis at time t. Otherwise, the nut will go along with the positive direction of the x axis, which is revealed in Fig. 4b.Then, the position x(t) of the nut on the screw can be written as follows: The loading location of boundary condition in the model is The node temperature at the previous time step is required by the backward difference method. However, the grid of the model is changing, because location of q 2 , x(t), is a function of simulated time. Therefore, the linear interpolation method is introduced to calculate the temperature of the current nodes (t + Δt) at the previous time t, as shown in Fig. 5. 1) The initial temperature of the screw T(x,0) will be given.
And the time t is set to 0.
2) The position of the nut on the screw x(t) at the time t can be calculated by Eq. (14). And then the mesh of the model can be updated.
3) The screw temperature at time t i can be solved after the boundary conditions are calculated by Eq. (15) and add them into the model. 4) If t is less than the total simulated time t all , t is accumulated by one-time step Δt, and then the process is returned to step 2). The position of the nut x(t) will be recalculated. And the grid and the node temperature of the model will be also updated. 5) The above loop will repeat until t is equal to t all .
The complex working conditions of the feed system can be regarded as a combination of single work condition, which can be expressed in the form of a matrix.
Each column in matrix A is a single work condition vector. The screw temperature of the previous work condition can be used as the initial condition for the calculation of the next working condition. And then the ball screw temperature field distribution under complex working conditions can be obtained.
Having the screw temperature calculated, the thermal error will be predicted. There will be different equations to calculate the thermal extent of the screw, due to the different bearing installation form. When the screw was installed with both sides of the screw fixed, the thermal error of the screw can be calculated by Eq. (18). If one side of the screw is fixed, and the other side of it is supported, Eq. (19) will be used to predict the screw feed system's thermal error.
where α is the thermal expansion coefficient.

Heat generated by bearings
The bearing heating is caused by the friction loss and rolling resistance between the ball and the raceway in the contact area. An empirical formula for calculating bearing friction torque was proposed by Palmgren [21], which is composed of two parts, including the friction torque M 1 generated by the viscosity of the lubricant when the bearing is idling and the friction torque Mv generated by the load action independent of the speed. The bearing heating is related to the friction torque of the contact area and rotation speed,

Heat generated by nuts
Similar to the bearing heating, the heat generated by the screw nut pair is mainly caused by the friction between the ball and the nut groove: where n is screw rotation speed, M is the total friction torque, M a is the torque that overcomes axial force and cutting force, M b is the sum of resistance torque, F a1 is the axial force loaded on the nut, F a0 is the nut preload, P h is the lead of screw and η is transmission efficiency.

Heat dissipation calculation
The forced convection will be conducted between the screw surface and air, when the feed system is running. The convective coefficient can be expressed as follows, according to the Nusselt criterion.
where N u is the Nusselt number, λ f is the thermal conductivity of the air and L is the hydraulic diameter.
The Nusselt number can be calculated as follows [22]:

Experimental setup and testing results
A series of experiments were carried out on the X-axis ball screw semi-closed loop feed system of a gantry milling machine in this paper, as shown in Fig. 7. The 7006C/DF angular contact ball bearings with back-to-back installed The dimensional parameters and material properties of the screw are given in Tables 1 and 2 respectively. As shown in Fig. 8, four thermal sensors are arranged around the feed system to measure the front/rear bearing, nuts and environment temperature. Meanwhile, the laser interferometer is adopted to measure the positioning accuracy of the feed system at different times, and then the thermal error can be calculated.
The temperature rise at each test point is shown in Fig. 9, when the feed system runs with the work conditions in Table 3.
During the test, the environmental temperature changed less than 1 ℃, the influence of which on the positioning accuracy can be ignored. The feed system can reach the thermal steady state after running for nearly 3 h. The temperature rise at the front bearing house is the highest, which is closed to 3 ℃. While, the temperature rise at the nut is the lowest, around 2 ℃ revealed in Fig. 9. Figure 10 shows the thermal error of the feed system under different work conditions. When the feed system runs under different work conditions, although the temperature rise at the measuring point is roughly the same, the positioning error is completely different. As shown in Fig. 9a, b, when the nut runs on the front half of the  screw (0-200 mm and 0-300 mm), thermal error at end of the screw can reach 10 ~ 15 μm. The thermal extension of the screw is larger on the front half of the screw. So the thermal error tends to be concave-down along the axial direction, while the changing trend of thermal error is the opposite when the nut is running on the rear half of the screw (400-600 mm and 400-800 mm) as revealed in Fig. 9c, d. The maximum thermal error can reach around 20 μm under work condition 3 and 4.

Model parameter identification
The accuracy of the boundary condition value is one of the key factors to guarantee the accuracy of the thermal model. However, the convective coefficient between the screw surface and its surrounding air is difficult to be accurately calculated with simple empirical formulas, due to the thread grooves on the screw surface. In addition, the heat generated by nut is also difficult to be accurately expressed based on the simple empirical formula, because it can be affected by various factors such as preload, abrasion. Therefore, the parameter identification method was adopted by this paper to calibrate the convective coefficient and the heat generated by the nut based on the experimental data of the above work conditions. The mean square error (MSE) of the simulated value and the experimental value is applied to construct the objective function: where mse i is the MSE under the i th working condition, m is the total number of samples, ε e (k) is experimental thermal error at the k th sampling, ε sim (k) is the simulated thermal error at the k th sampling. The convective coefficient h and the heat generated by nut q 2 were taken as the optimized variables. The estimated value of h and q 2 will be calculated first based on these empirical formulas, and then a more reasonable solution domain for them will be set, in order to improve efficiency accuracy of the optimization. And then, 50 Design Points of Experiments (DOE), within the range of 50 ≤ h ≤ 60, 690 ≤ q 2 ≤ 720 , are generated through the Optimal Latin Hypercubes (OLHCs) method [23]. And the moving least squares method [24] (MLS), expressed by Eq. (28), is introduced to map optimized variables and objective functions which is shown in Fig. 11. Then, the original FE model will be replaced so that there will be a more efficient model in the calculation of optimization.
where, a(x) is the function of curve correction coefficient; p(x) is the basis function. Compared with least squares method, more smooth fitting curves will be obtained by the MLS. Meanwhile, the segmentation fitting can be avoided because of the introduction of weighting function W(s) as follows: In the final, the Pareto front of h and q 2 was obtained by the multi-objective optimization method, as shown in Fig. 12.
The boundary condition values in the model are listed in Table 4.

Validations under different operating conditions
The feed system runs with the different work conditions that are listed in the Tables 5 and 6 to validate the proposed model.
The proposed thermal error model is coded by Python3 and runs on the common PC (CPU: Intel Core i5-9400F, RAM:16G). The calculation time of the model with different element numbers and time step Δt is listed in Table 7, when the work condition 5 was used as input.
Although the grid needs be updated at every time step in the proposed method, high computational efficiency still can be guaranteed due to the small amounts of elements. And the computational time is much less than the simulated time (3 h). The number of elements is 236, and the time step Δt of 1 s is applied for the following simulation to ensure the calculation accuracy and efficiency. Figure 13 illustrates the experimental results and the simulated results of the feed system, which runs with the   working condition listed in Table 5. The maximum thermal error of the feed system can reach up to 14 μm, if the feed system runs according to the work condition 5.
The maximum prediction error of the proposed model is within − 4 μm. The maximum thermal error of the feed system under work condition 6 is nearly 14 μm. The model's prediction error is within 5 μm, as shown in Fig. 13b. While the prediction error of the model also is within 4 μm, when the feed system runs with work condition 7, which is revealed in Fig. 13c. The ratios of prediction error to total thermal error of the feed system from work condition 5 to work condition 7 are listed in Table 8.
In the transient simulation of combined work conditions, the feed system runs at the feed speed of 0.09 m/s for 0.5 h in each interval, according to the working conditions listed in Table 6. The laser interferometer is used to collect thermal error every 0.5 h.
The matrix A 1 and A 2 are respectively created as the input of the proposed model, and the simulated value of the thermal error would be obtained in the final. The simulated results under combined work condition are illustrated in Fig. 14. And the simulated data at 0.5 h, 1.0 h and 1.5 h are extracted respectively to be compared with the experimental results to validate the thermal error model proposed in this paper.
As shown in Fig. 15, the maximum thermal error of the feed system with the work condition 8 varies between 11 and 18 μm. And the prediction errors of the proposed model are 0 ~ − 3 μm, − 2 ~ + 2 μm and 0 ~ 2 μm respectively.      model are − 2 ~ + 2 μm, − 2 ~ + 1 μm and 0 ~ 5 μm respectively. The ratios of prediction error to total thermal error at different time in the transient simulation are listed in Table 9.
To sum up, the ball screw feed system thermal error model proposed in this paper has high accuracy under various working conditions, the prediction accuracy of which is more than 80%. However, the prediction accuracy of the model decreases, when the feed system is running on a larger stroke, what is revealed in Figs. 13b and  16c. And the predicted value is lower than the actual value, which may be caused by the actual heat generation increasing during runtime of the feed system, because the lubrication status gets worse in the case of large stroke running.

Conclusion
A transient thermal error model of the ball screw feed system is developed in this paper based on FEM. The loading of the moving heat flux at the nut is realized by changing the grid of the model at each time step to predict the thermal error of the feed system under complex work conditions. In addition, the multi-objective optimization method is applied in this paper for calibrating the boundary conditions of the model to improve the accuracy of the model, because of the difficulty to calculate the convection coefficient and the heat generated by the nut accurately. The core conclusions of the study are as follows: (1) Under different work conditions, even if the temperature rises at the measuring point are approximately close, the thermal errors of the feed system can also be completely different due to the difference in the screw temperature distribution. (2) The proposed model is reliable and accurate to predict the transient thermal error of the ball screw feed system under different work conditions, the prediction accuracy of which is more than 80%.
In further investigations, the thermal behavior of the feed system that continuously starts and stops will be explored. Meanwhile, the change of the heat generated by the nut during the feed system running needs to be taken into  consideration to improve the accuracy of the model. In addition, the proposed model can be further embedded in the NC system and compensate the thermal error online to improve the stability of the machine tool's positioning accuracy.
Author contribution DS: methodology, data curation, formal analysis, writing-original draft, review and editing; YL: conceptualization, supervision, funding acquisition, writing-review and editing; WZ: conceptualization, supervision, writing-review and editing; HZ: data curation, validation, writing-review and editing.
Funding This work was financially supported by the Key-Area Research and Development Program of Guangdong Province (Grant No. 2020B090927002). The first author would like to thank the China Scholarship Council (CSC) for financial support and the scholarship award.
Data availability All data generated or analyzed during this study are included.
Code availability Not applicable.

Declarations
Ethical approval Not applicable.

Consent for publication
The authors consent that the work entitled as "Transient thermal error modeling of a ball screw feed system" for possible publication in International Journal of Advanced Manufacturing Technology. The authors certify that this manuscript is original and has not been published in whole or in part nor is it being considered for publication elsewhere.