Thermal error prediction of ball screw feed system based on inverse heat transfer analysis

The thermal error of the ball screw feed system will reduce the positioning accuracy of the machine tool, and the main reason for the thermal error is the friction heat generated by bearings and screw nut in the feed system. Therefore, a time-varying thermal error prediction method of the ball screw feed system based on the inverse method is proposed in this paper. Based on data-driven and dynamic thermal network model, the heat source estimation method of ball screw feed system under unsteady thermal effect is established, and the thermodynamic equilibrium equation is deduced into explicit form for heat source load identification. Aiming at the common matrix ill-conditioned problem of load identification, the regularization algorithm is used to identify the heat source load, and the optimal selection method of regularization parameters is proposed based on Monte Carlo algorithm. Using the collected temperature experiment data and the position data of the moving nut, the dynamic heat source load, temperature field, and thermal error of the feed system under the actual working condition are predicted and analyzed by using the inverse method. Finally, the accuracy and effectiveness of the prediction method are verified by experiment. The inverse method proposed in this paper has great application potential for the prediction and estimation of heat source and temperature field of machine tools and various mechanical structures.


Introduction
As the key part of precision transmission and positioning, the feed system of CNC machine tool plays a vital role in the positioning accuracy of machine tool. Ball screw feed drive system has the advantages of high efficiency, high stiffness, and long life. It is widely used in CNC machine tools. However, during the working process of the feed system, friction heat leads to temperature rise, resulting in thermal deformation of the structure and affecting the tool positioning accuracy [1,2]. With the increasing demand for feed speed, the linear speed of ball screw can reach up to 100 m/ min, but the influence of thermal effect is more prominent [2]. Previous studies have shown that the precision machining error caused by thermal deformation accounts for 40 to 70% of the total error [3,4]. Therefore, the research on the thermal error of machine tool feed system has become an important topic in the field of ultra-precision machining.
The thermal error of feed system is a time-dependent nonlinear process caused by uneven temperature variation in mechanical structure. In view of the current research of scholars, there are two main methods to reduce thermal error, i.e., avoid thermal error and error prediction compensation. Error avoidance and error prediction compensation are the two main technical measures to address the thermal issues of the ball screw. Most of the research on error avoidance is to reduce the thermal error by optimization design, strengthening the manufacturing process and improving machine tool structure. The main method is to eliminate the potential thermal error of the system in the design stage by means of structural symmetrical design, active cooling system, and the application of low temperature expansion materials [5][6][7][8]. Although this kind of method is reliable, it has great limitations. It has a high economic cost and cannot break through the limitation of existing processing capacity, and has no adaptability to complex working conditions. Therefore, error avoidance is not an economical method and has not been widely used.
Because the error prediction compensation method has the advantages of relatively low cost and wide application range, most scholars are committed to developing various error modeling strategies. The common method of error prediction is to establish a thermal error model through the temperature and thermal error data obtained from experiments to describe the relationship between temperature and thermal error. Generally, various data-driven model are widely used to establish thermal error models, including time series method [9,10], gray theory [11], and neural network [12][13][14]. Miao et al. [9] used the time series method to establish a thermal error model, predicted the change of time series by studying the variables and their inference mechanism, gave greater weight to the data near the prediction, and increased the impact of short-term parameters on the model, so as to improve the prediction accuracy. Abdulshahed et al. [11] used cuckoo search algorithm to correct gray model, and then verified this new intelligent modeling method of gray model as an alternative optimization algorithm for thermal error compensation model. However, the above datadriven modeling method needs multiple sensors to achieve, which is relatively expensive. Moreover, the research process needs lengthy and time-consuming experimental analysis, so it is difficult to maintain the accuracy and robustness of the compensation model in practical application. If the obtained data is not overall, the established model will not be suitable for random circumstance.
Other researchers attempt to model thermal error through numerical analysis and simulation calculation, and study the thermal deformation mechanism. This kind of modeling method is also known as mechanism-driven method [15][16][17][18][19][20]. Wang et al. [16] proposed a hybrid modeling method for z-axis thermal error of machine tool. The heat transfer equation was established based on the finite difference method, and the model coefficients were optimized by genetic algorithm. The results showed that the model has better prediction accuracy. Liu et al. [17] explored a new model for predicting the thermal contact resistance based on rough surface morphology to improve the simulation accuracy of finite element transient thermal analysis of feed system. The multi-objective and multi-parameter optimization model based on response surface analysis was established, and the temperature field and thermal deformation of the feed system were simulated by ANSYS, which improved the prediction accuracy of the finite element simulation. Liu et al. [19] established the thermal error prediction model of machine tool feed system based on the finite difference method and presented the identification method of model parameters. The modeling method does not need too many measuring points and the optimal layout analysis of measuring points, and the prediction results are still excellent under actual working conditions. To sum up, the researchers had made abundant achievements in the thermal error modeling based on the mechanism-driven method. At present, most of the friction heat sources are still determined by simple empirical formula, and the dynamic change of the heat sources under the influence of thermal is not considered. The mechanismdriven is separated from the experimental data, so the modeling method is not perfect and the robustness is not good. Therefore, it is necessary to explore a new idea combining the advantages of data-driven and mechanism-driven methods to solve the thermal error problem of the machine tool feed system.
A century ago, with the rapid development of computers, inverse analysis has attracted more attention in the past few decades [21][22][23]. When the system parameters are physical quantities that are difficult to measure or cannot be measured, the inverse method technique is very useful for obtaining the information of physical quantities. The commonly used early inverse method techniques include least square method, Kalman filter, linear quadratic estimation. Wang et al. [20] estimated the heat flux during grinding using mathematical modeling inverse method based on experiment data, which reflected the law of heat flux with time by establishing a time-varying system. Li et al. [15] achieved effective estimation of thermal boundary conditions for machine tool feed system by combining the finite element method with optimization method. Liu et al. [22] proposed a new method to solve the nonlinear convection diffusion heat problem of two-dimensional rectangle and three-dimensional rectangle. The numerical solution of temperature was substituted into the nonlinear equation to identify the unknown heat source. The numerical simulation results showed that the method was very accurate to solve the unknown heat source function.
In general, the inverse method had been successfully applied to the estimation and determination of heat source and heat transfer coefficient for heat transfer problems. There are two widely used inverse methods: one is the inverse method which combines the identification calculation method with commercial software, and the other is the inverse method of mathematical modeling based on heat transfer theory. The inverse method is still rare in the research of machine tool thermal analysis. Therefore, aiming at the problem that the heat source of machine tool feed system cannot be accurately modeled and calculated, this paper proposes a thermal error prediction model of feed system based on inverse method. Based on the experimental data and dynamic thermal network model, the inverse method of heat source estimation under unsteady thermal effect is established, and the thermodynamic equilibrium equation is deduced into explicit form for heat load identification. The dynamic heat source, temperature field, and thermal error of the feed system under actual working conditions are predicted and analyzed by inverse method based on temperature data of surface measuring points. The experimental results show that the proposed model considers the influence of the dynamic change of heat source and can more accurately describe the thermal error under working conditions than the traditional prediction model.

Model description
The two-dimensional structure diagram of the z-axis ball screw feed system of CNC machine tool is shown in Fig. 1. The heat transfer process of the feed system is heat conduction from the heat source to other parts, convective heat transfer between the heat transfer boundary and the surrounding environment, and finally reaches heat balance with working time. Due to the limitation of sensor placement, it is difficult to measure the temperature of friction contact surface, so it is difficult to accurately predict the internal temperature field of the feed system. In addition, there is a complex coupling relationship between friction heat and thermal deformation of the moving pair, so it is difficult to predict the temperature distribution.

Thermal network model
In this paper, the heat transfer model of the ball screw feed system is established by using the thermal network method. The modeling principle of the thermal network method is based on practical needs and convenience of measurement. The heat transfer model is divided into several units, and each unit is replaced by a temperature node. According to the heat conduction relationship between nodes, the heat transfer equation of the whole system is established, and the transient temperature field of the system is obtained by solving the temperature of each node. Due to the structural complexity of the ball screw feed system, it is necessary to reasonably arrange the thermal network temperature nodes of the feed system and divide the structural characteristics of the feed system. According to the principle of structural symmetry of the feed system, the typical temperature nodes are selected to establish the thermal network model. The ball screw feed system includes two fixed heat sources and a moving heat source of screw nut. The thermal network node arrangement of the heat source is shown in Fig. 2. The thermal loads of node 19 and node 20 are equal, accounting for one-half of the total heat flow of the bearing, and node 68 and node 69 are equal, accounting for one-half of the total heat flow of the bearing at the motor end. The node distribution of the screw shaft shall be consistent as far as possible, and the mass and volume of each node shall be the same, so as to facilitate modeling calculation in the later stage. The node description of the feed system is shown in Table 1. The thermal network model is similar to the circuit diagram, but the temperature difference of each node in the system is the power of heat flow transfer, thermal resistance is the resistance of heat transfer, heat capacity is the ability to absorb heat, and heat flow is the form of heat transfer. Based on the node layout in Fig. 2, the thermal network model of the feed system is established, as shown in Fig. 3. The following assumptions are made for the model before modeling and calculation: 1. Ignore the effect of the temperature rise of the feed system on the environment temperature, and the environment temperature remains unchanged. 2. The thermal conductivity and heat capacity of the material are constant. 3. Thermal radiation is not considered:

Heat transfer
The resistance encountered in the heat transfer process is called thermal resistance, which represents the heat transfer capacity. Before modeling, it is necessary to calculate the thermal resistance of all nodes. The calculation formula of thermal resistance for the feed system is expressed as follows.  where d i and d o are inside and outside radius of the installed object, A is the contact area, k I is the thermal conductivity, and L is the length of the object. 2. For the thermal conductivity of the feed shaft, it can be equivalent to several overlapping hollow cylinders in the radial direction, and the axial and radial thermal resistances are calculated as where L x is the axial distance between nodes, and k shaft is the thermal conductivity coefficient of shaft material.
(1) (ln( 3. Thermal contact resistance R b is related to ball parameters and preload. The thermal resistance of the ball is expressed as where k b and k r are the thermal conductivity, and a and b are the long/short axis of contact ellipse formed. 4. The thermal resistance of the inner ring in contact with the shaft is expressed as Here A is the apparent contact of two contact surfaces, A * R is the actual contact area, and L g is the void thickness. k d1 and k d2 are the thermal conductivity of two objects respectively. 5. The calculation method between the outer ring and the bearing housing is expressed according to the literature query.
where L S is the gap thickness, and k f is the thermal conductivity of the medium. k D1 and k D2 are the thermal conductivity of two objects respectively. 6. The contact thermal resistance of convection heat transfer in the process of heat conduction is calculated where N u is the Nusselt number, K L is the thermal conductivity of the fluid, and A is the contact area. h is the heat transfer coefficient, and the calculation method is not discussed in detail here.

Heat balance equation
The principle of the energy balance method is shown in Fig. 4. The thermal network model generates heat transfer between nodes, and each node is connected to multiple nodes. Different nodes exist convective heat transfer or heat conduction, and in order to establish a general mathematical model of thermal resistance network, any node is identified as The total number of its associated nodes is n, the thermal resistances of node i is defined asR in , and the convective heat transfer resistance of node is defined asR air i , if the thermal network node is an internal node of the structure, and the thermal resistance Internal node of the feed shaft of the preload end bearings 9~16 Feed shaft surface nodes for preloaded end bearings 17,18 Inner ring node of the bearings 19, 20 Ball node of the bearings 21, 22 Outer ring node of the bearings 23~32 Bearing housing nodes 33, 34 Bearing seal nodes 35~38 Bearing cover nodes 39, 40 Bearing sleeve nodes 41 Rolling body nodes of the moving nut 42~44 Moving nut raceway nodes 47~49 The outer surface node of the moving nut 50~57 Internal node of the feed shaft of the motor end bearings 58~65 Feed shaft surface nodes for the motor end bearings 66, 67 Inner ring node of the bearings 68, 69 Ball node of the bearings 70, 71 Outer ring node of the bearings 72~83 Bearing housing nodes 84, 85 Bearing seal nodes 86~91 Bearing cover nodes 92~93 Bearing sleeve nodes 94~N Screw shaft nodes Air Ambient temperature  where Ξ i is the set of nodes related to node i.
Compared with the traditional thermal resistance network model, this paper proposes to model the heat transfer analysis by adding the heat capacity term, and the formula is expressed as follows: For the first-order differential equation in Eq. (8), the implicit scheme is used for discrete processing (Crank-Nicolson Method), and the transient first-order difference processing is carried out for the temperature rise rate of the node i.
where T k+1 i is the temperature of the node i at the time k, T k+1 i is the temperature of the node i at time k + 1, and Δt is the time step.
The temperature change is taken as the average value of the time interval, which is expressed as Substituting Eqs. (10) and (11) into Eq. (8), Eq. (12) is obtained by rearrangement.
In the resistance capacitance thermal network of the feed system, determine the total number of nodes as N, and set the vector x to represent the temperature of each node.
The matrix form of numerical calculation can be obtained from Eq. (12).
In Eq. (14), the coefficient matrices A and B are square matrices, and their diagonal elements can be written as follows: The non-diagonal elements are related to the thermal resistance between the nodes. If the thermal resistance between the ith node and the ilth node is R j , it can be written as follows: If there is no thermal resistance connection between the nodes, the corresponding elements of the A and B matrices are 0. The matrix P is an N-dimensional column vector, and when node i is a structural surface node with convective heat exchange with the air, thenp i = 1∕R air i ; otherwise, p i equals zero.
The matrix Q is the heat flux vector, and the heat source input nodes are arranged according to node serial number, as shown in Fig. 2. The matrix L is the heat source mapping matrix, and if there is no heat flux input to the network node, the corresponding row of the L matrix is 0. If the corresponding node in the thermal network is a heat source node and corresponds to the jth element of matrix Q, the value of the jth element of matrix L corresponding to the row is 1 and the values of the other elements are 0. The coefficient matrices A and B are symmetrical, and the matrix A is strictly diagonally dominant, which can be obtained from Eq. (14).

Heat source identification
For the thermal network model of feed system, it is assumed that the system is tested in a constant temperature environment, the ambient temperature under working conditions remains unchanged, and D is the fixed vector. If the heat transfer system starts from the ambient temperature, the initial value of each node is the ambient temperature.
If the sensor detects m temperature points and its detection value is y l (l = 1, 2, ⋯ , m) , then the detection vector is expressed as For a given time step, the calculation relationship of the detection temperature is expressed as where R is a matrix of m × N. If the detection value y i corresponds to the lth temperature node in the network, then the other elements are 0. Multiply matrix R by Eq. (22) to obtain Create the z vector: where H is the block matrix.
For the calculation of thermal network system, the system temperature reference boundary is determined first. It is assumed that the heat transfer system works in a constant temperature environment and the ambient temperature remains unchanged during the working process, and the ambient temperature is selected as the reference temperature. The system numerical calculation Eq. (28) is a zero initial problem, which is expressed as: For a constant vector C, when a node exists in the heat source, its value is the product of the value of the Q i allocated by the node and the time step; when the node has no thermal source, its value is 0.
Through the above modeling method, the Z can be deduced for a given system temperature field. The ill-condition degree of matrix H is estimated according to the condition number of matrix H. When the condition number of the matrix is too large, the matrix of the system is an ill-conditioned matrix. The judgment formula is expressed as follows: If the condition number is too large, it is not feasible to use the method of direct inversion identification. The heat source load is strongly dependent on the collected data, and weak interference will also affect the results, so it is difficult to identify the accurate heat load. There are measurement errors, truncation errors, and other errors in the detection data of each test point of the system. The system temperature field equation containing the errors can be expressed as: where e is the error vector.
In order to overcome the interference effect of the illconditioned matrix on the inverse problem, the regularization method is selected to solve the problem. The regularization method is a common method for solving ill-conditioned problems in linear algebra theory, including Tikhonov regularization, various iterative regularization, and wavelet regularization. When solving the thermal inverse problem of the temperature field, it is assumed that the following equation is satisfied between the accurate data and the collected data containing error disturbances. This paper adopts Tikhonov regularization method to solve it.
where α represents the regularization parameter. The solution of the inverse problem becomes to find the minimum value of the objective function, and the regularization solution Q , is the solution of the inverse problem. In order to find the minimum value of the objective function, the derivative of J(Q, ) is expressed as follows: If the above formula is equal to 0, the regularization solution can be derived to satisfy the following equation: The explicit unique solution of Q , can be obtained:

Regularization method
The above derivation process gives the explicit solution for identifying the thermal load by Tikhonov regularization, and the solution process is very important to the regularization parameters. Although l-curve method can select regularization parameters, the calculation in the process of drawing l-curve is cumbersome and the efficiency is low, and it is very inconvenient to use. Therefore, this paper proposes a parameter selection method based on Tikhonov regularization, which can select regularization parameters more effectively through identification-contrast-optimize (ICO). The specific implementation process of ICO method is shown as follows: 1. Identification Through estimation, the regularization parameters are appropriately selected, and the regularization algorithm is used to bring them into the matrix equation group to identify the thermal load. The identified thermal load is substituted into the heat transfer model; the matrices A, B, D, and E are solved; and the temperature field is solved iteratively by Eq. (32) to obtain the temperature data.

Contrast.
The temperature response obtained by iterative solution and the experiment temperature response signal measured in real time are calculated according to Eq. (40). The calculation is completed until the regularization parameter satisfies the output termination condition

Optimize.
When the regularization parameters do not meet the termination conditions, it indicates that the parameters are not accurate enough and need to be reasonably optimized and adjusted to meet the accuracy requirements of the solution. In this paper, the regularization parameters are simulated and determined by the Monte Carlo method. The Monte Carlo method is proposed based on probability and statistics. Assuming that the lower and upper bounds of the regularization parameters are LM and HM respectively, and then a random function is used to generate k sets of regularized sample values: Each group of regularization parameters is substituted into the first step for identification calculation, and the estimated results are brought into the second step for comparison. Perform K sets of simulation iterative calculations, and search for the minimum value of W in K sets of simulation results; if the calculation accuracy is met, The selection method of the regularization parameters is based on the established thermal network model, using direct method and inverse method to calculate the error of the temperature response of the experiment data and the identification results, and this paper gradually optimizes the regularization parameters. When the computer simulation results meet the accuracy requirements, the search iteration process is ended. This method is adapted to practical needs. Figure 5 shows the numerical calculation process of this study.

Dynamic contact node processing method
Due to the reciprocating movement of the moving heat source, the contact nodes of the screw shaft change with time and position, which increases the difficulty of thermal  Guide  T1  T2   Screw   T4   T3   P2  P3  P4  P5  P6  P7  P8  P9 P10 P11   T5  T6  T7  T8  T9  T10  T11  T12  T13 T14 T15 Worktable Worktable (a) Location of temperature measurement point (b) Location of positioning error measurement point Fig. 9 The temperature rise curves of different nodes at different feed rates analysis of the ball screw feed system. Therefore, the thermal resistance of the nodes on the shaft surface changes dynamically with the movement of the screw nut. Therefore, this paper proposes a processing method for dynamic contact nodes. Figure 6 shows the relationship between stroke position and the heat transfer of screw nut heat source, and the length of nut is L N . It can be seen from Fig. 6 that the residence time of the screw nut on each network node is different, and the residence time of the middle part is longer, so the heat transfer is higher than that on both sides. Here, the time of moving nut from left to right is Δt , and the residence time of any node of the feed screw shaft in contact with the heat source is Δt i , then the proportion of the residence time of the node in contact with the heat source can be expressed The heat conduction time between the node and the air is calculated Thus, the heat transfer equation of any node of the stroke can be expressed as where T Q is the temperature of the heat source node, and R Q is the thermal resistance between the node and the heat source.
By rearranging Eq. (46), it can be rewritten as Therefore, it can be understood that the thermal resistance has changed, and the Eq. (47) is redefined as According to the above method, the H matrix can still be derived according to Eq. (31). The change of thermal resistance of dynamic contact node is shown in Eq. (48), which solves the problem of dynamic contact node for the thermal network model.

Experiment
In order to verify the effectiveness of the thermal source identification method of the feed system, the experiment system is shown in Fig. 7. The experiment was measured at a spindle of 0r/min, and the experiment object was the z-axis feed system of the CNC machine tool. The experiment was carried out in a constant temperature laboratory, which can ensure that the ambient temperature did not change. The parameters of z-axis ball screw feed system are shown in Table 2.
In this experiment process, no-cutting load, the ambient temperature is 25 °C, and the distribution of temperature measuring points and positioning error measuring points of the experiment system is shown in Fig. 7. In order to explore the temperature rise law of the heating parts of the feed system, thermocouple was used to detect the temperature rise data. According to the measuring criteria of ISO 230-2, the key point temperatures of the designated points T1, T2, T3, and T4 were measured by thermocouples every 0.5 s. During the test, the infrared thermal imager was used to measure the surface temperature of the lead screw, and the temperature of the measuring points T5-T15 were recorded 10 min. The worktable reciprocated at a steady rate between T5 and T15 measuring points, and its stroke was 0-420 mm. In these experimental measurements, the laser interferometer was used to measure the positioning error of designated points P1-P11, and the point was recorded every 42 mm. When the lathe had been heated for 10, 20, 30, 40, 50, and 60 min, the temperature and positioning error of each point were recorded by infrared thermal imager and laser interferometer. During the test, FOCAS function was used to monitor the position of screw nut and the torque current information of drive motor through the whole Fig. 10 Temperature field distribution of screw shaft based on simulation calculation 1 3 Fig. 11 Distribution of test points for temperature and positioning error on Z-axis process using the developed computer acquisition system (Fig. 8). The computer acquisition system of sampling interval is 0.5 s, and experiment parameters are shown in Table 2.

Results and discussion
In order to verify the effectiveness of the load identification method, the load identification process shown in Fig. 5 is used for modeling and solving. The experiment point T1 is equivalent to nodes 29 and 30 in the thermal network model, T2 is equivalent to nodes 82 and 83 in the thermal network model, and T3 is equivalent to node 49. The temperature data collected by the experiment are imported into the load identification program as input signal. The temperature acquisition data need to be pre-processed before solution.
The thermocouple acquisition period of the test system is 0.5 s, and the temperature sampling period in the load identification program can be adjusted appropriately. After identifying the simulation calculation, the node temperature data as shown in Fig. 2 can be obtained. The temperature rise curves of different nodes at different feed rates are shown in Fig. 9, and the thermal network model can simulate the internal temperature field of the whole system. When the feed rate is 5 m/min, the temperature field distribution of the whole lead screw shaft is shown in Fig. 10. Figure 11 shows the comparison between the experiment temperature rise and the identification results under the working conditions. From the diagram, it can be seen that the identification prediction results are in good agreement with the experiment data, and the temperature field gradually reaches a stable state after the rapid increase of time. It is proved that when the heat generation and heat dissipation reach equilibrium, the temperature field is close to the thermal steady state. At the same time, it can be seen from the identification results that heat source estimation method has a very good effect on the temperature field prediction. Therefore, using the experiment data to identify the thermal load to predict the temperature field of the feed system is closer to the real working conditions than directly giving the fixed thermal source load for thermal analysis. The results show that the accuracy of the thermal source load plays a key role in the thermal analysis of the feed system, and the accuracy of the thermal load is related to the accuracy of temperature field simulation. From the above results, the feasibility of the proposed inverse problem solving method for thermal analysis is verified, and the method can well predict the temperature field of the mechanical structure. Figure 12 shows the calculation results of the heat generation rate identified by the inverse method. When the feed system is started, the heat generation rate reaches a peak value, and then the heat generation rate decreases rapidly in a short period of time and finally tends to be stable. The reason for the decrease may be the change of lubricant of kinematic viscosity and preload. According to the above results, it is proved again that the heat generation rate is not a fixed value, and the heat generation rate changes dynamically with the Fig. 12 Heat source identification results of the different rates working time of the feed system. The feed system is tested without external force, and it can be guessed that the friction heat is related to the mechanical structure, lubricating oil, preload, and its dynamic temperature rise. Therefore, it is proved that the mechanical structure has complex thermal mechanical coupling relationship under working conditions.
In the whole process of the experiment, PC104 computer was used to collect the torque current of the servo motor and the position data. According to the corresponding relationship between torque current and time, the corresponding relationship between torque current and worktable position was obtained through interpolation calculation. As shown in Fig. 13, comparing the average torque current at the same position for different running times, it can be seen that the servo motor torque current decreases as the continuous running time of the system increases. The experimental results of torque current are in agreement with the predicted results of heat generation rate in this paper, which proves the correctness of the proposed method. At the same time, it can be seen that when the servo motor moves at a uniform feed rate, the relationship between the average torque current and the position of the worktable is approximately parabolic; that is, the middle is low and the two ends are large. With the increase of the running time of the feed system, the position of the parabola gradually decreases, but tends to stabilize with the extension of time.
The thermal error of the screw shaft can be predicted by calculating the temperature field distribution of the screw shaft. The screw shaft of the feed system is elongated due to the nonlinear temperature rise of the screw shaft caused by the frictional heat. The axial thermal error of the feed system is expressed as [18].
The temperature rise of each node is obtained by the identification calculation method, and the thermal elongation of any node during the working period of the system is where l j is the length of the j node; ΔT j is the temperature rise, and the thermal expansion coefficient is  Figure 14 shows the thermal error of each measuring point of the screw shaft under different feed rates. The thermal error of the feed system at different running time approximately shows a linear variation relationship with the position of the worktable, and increases gradually with the increase of the system running time, but the increase amplitude gradually decreases. When the running time exceeds to 30 min, the thermal error of the worktable with the position changes gradually tends to be stable. Comparing the thermal error experiment data with the simulation results, the maximum error occurs when the feed speed is 15 m/min. It can be seen that the thermal error results are in good agreement with the experimental results. The above results show that the heat source loads determine the accuracy of the thermal error and temperature field, and the accurate boundary conditions of the heat transfer process can make the predicted temperature field and thermal error closer to the real working conditions. At the same time, the boundary conditions of heat transfer system are determined by identifying heat source heat generation rate based on experimental data, and the simulation results of this method are more consistent with the actual heat transfer process. The heat generation rate identification method of the feed system based on the thermal network is more practical, and the method combining data drive and mechanism drive to predict the temperature field and thermal error results can improve the accuracy of the results.

Conclusions
Based on the thermal error mechanism of the feed system, this paper proposes a time-varying heat source prediction method of the feed system based on the inverse method. This method can be used to calculate the temperature distribution and thermal error of the feed shaft with time without using finite element method or any iterative simulation: 1. According to the geometric and physical characteristics of the ball screw feed system, the thermal resistance parameters of the ball screw feed system are determined, and the resistance capacitance heat conduction model is established. Based on the heat flow equilibrium relationship and the central difference principle, the thermodynamic equilibrium equation is deduced into an explicit expression, and a numerical calculation algorithm for thermal load identification under unsteady thermal effect is proposed. 2. Aiming at the common ill-conditioned problems of load identification, the regularization algorithm is used to solve it, and the selection method of the regularization parameters is proposed. The dynamic heat source, temperature field, and thermal error of the feed system under actual working conditions are predicted and analyzed by using the temperature data collected from the experiment and the inverse heat transfer model. It is found that the heat source of the feed system is a dynamic heat source varying with time. 3. By comparing the numerical simulation results with the experiment data, the practicability and effectiveness of inverse heat transfer analysis method is verified. The proposed method can predict the heat source and temperature field of the feed system. The results show that the method can be used for dynamic heat source identification of complex mechanical structures and has good prediction accuracy.

Availability of data and material
The data cannot be shared openly, to protect study participant privacy.

Code availability
The code cannot be shared openly, to protect study participant privacy.

Declarations
Ethical approval and consent to participate Not applicable.