Thermal error modeling of spindle and dynamic machining accuracy reliability analysis of CNC machine tools based on IA and LHSMC

E-mail addresses: Machining accuracy reliability as a key index of CNC machine tools is seriously influenced by the geometric and thermal errors. In the paper, a spindle unit thermal error modeling and machining accuracy reliability analysis method is proposed. By analyzing the heat generation mechanism, a thermal error model was developed to describe the thermal deformation of the electric spindle. Based on the immune algorithm (IA), the heat generation power and the heat transfer coefficient were optimized, and the thermal error was obtained by finite element thermal-mechanical coupling. By adopting the multi-body system theory (MBS), a dynamic machining accuracy model was put forward including the geometric and thermal errors. Based on the Latin hypercube sampling Monte Carlo method (LHSMC), a machining accuracy reliability analysis method was proposed to characterize the machining accuracy reliability considering the geometric and thermal errors. The method was employed to a machine tool, and the experimental results indicate the verification and superiority of the method. Highlights Abstract


Introduction
CNC machine tools are favored by various enterprises as core equipment owing to their high-speed and high-accuracy characteristics. Machining accuracy reliability is an important criterion, which reflects the capability of machine tools to achieve the desired requirements. The machining accuracy is directly expressed in the dimensional error of the workpiece in the work condition. Geometric and thermal errors are the major contributors to variations in machining accuracy [21,23]. Among the total error sources, the occupancy rate of thermal error is 40%-75% [14]. The higher machining accuracy reliability of machine tools, the greater their competitiveness, which is a factor that manufacturers need to focus on [4]. Thus, it is crucial to consider the impact of thermal error on the machining accuracy reliability. With the increased requirements for high precision and efficiency in manufacturing technology, the effective prediction of thermal error is becoming increasingly important in studying machining accuracy [30]. Hence, how to characterize the reliability of machining accuracy under the effect of thermal error is critical to evaluate the machining capability of machine tools [18].
The temperature variation of the motor and bearings are the major cause of thermal error [15]. Since it is difficult to be calculated effectively, scientific modeling methods to predict the thermal error are necessary [13]. The thermal error modeling process has the following steps: analyzing the spindle temperature rise field; establishing the connection between the critical region of temperature and thermal deformation field; collecting the critical region temperature to predict the thermal error [24]. Recently, many scholars have devoted themselves to the research of thermal error modeling methods. Zhao et al. calculated the convection heat transfer coefficient on the spindle surface, the temperature and deformation fields of the spindle were derived based on the finite element method and validated experimentally [31]. Hou et al. analyzed the coupled thermal deformation of machine tool components, and a multi-objective genetic algorithm was applied to derive a high robustness thermal error model [6]. By using the improved particle swarm optimization (IPSO) method, Li et al. developed a thermal error model, and the higher modeling efficiency and accuracy were validated by experiments [11]. Shi et al. introduced a thermal error model by Bayesian neural network, and used fuzzy cmeans (FCM) clustering analysis to select the sensitive points of temperature. The model has higher prediction accuracy compared to BP neural network [20]. Uhlmann and Hu established a three-dimensional finite element model to predict the thermal characteristics of electric spindles, quantified the heat transfer process internal to the spindle, and validated it under consideration of complex boundary conditions [22]. Zheng and Chen studied the thermal performance of angular contact bearings, analyzed the heat transfer process of the bearing sub source, and developed a comprehensive thermal grid model to predict the temperature rise of the bearings considering the effect of various factors such as constraints and assembly relationships [33]. To effectively forecast the transient temperature and thermal deformation fields of the spindle unit, Yang et al. developed a coupled thermalstructural model, and compared with the TCP thermal errors obtained by the regression model to validate the feasibility of this model [25].
Most of the above-mentioned thermal error modeling methods were proposed based on empirical formulas, which deviate from the actual situation and decrease predictive performance. In addition, the traditional temperature-thermal error and multivariate-thermal error models did not study the formation mechanism of thermal error of the spindle unit in depth, which leaves room for improvement in the application of these two methods [16]. Hence, how to construct a spindle thermal error model with higher prediction accuracy is an urgent problem.
Meanwhile, geometric and thermal errors are the direct cause of the machining accuracy deviation from the desired requirements [5]. Thus, it is imperative to perform machining accuracy modeling. To this end, there are a lot of efforts made by many researchers to develop modeling methods, such as HTM [10,7,26,5], screw theory [32,17], D-H method [9], Lagrange method [19], Spearman rank correlation method [3] and so on. Geometric errors can be measured using instruments, the translation geometric errors are determined using a laser interferometer, while the angle errors are determined using a ball-bar [34]. The machining accuracy requires consideration of multiple error sources, which is the key point of this paper's research.
Machining accuracy reliability can effectively indicate the capability of the machine tool to achieve the required machining accuracy [27]. So far, much work has been done by many researchers to characterize the reliability of machines. Cheng et al. developed the reliability and sensitivity model to analyze the machining accuracy of the machine tool based on Monte Carlo method, and verified with a three-axis machine tool [2]. Jiang et al. calculated the reliability of the electric spindle system based on the quasi-Monte Carlo method (QMC), used the Kriging method instead of the QMC method to make the calculation more efficient, and experimentally verified the applicability of the model [8]. Cai et al. developed the machining accuracy reliability and sensitivity models considering multiple failure modes based on the first-order and second-moment (AFOSM) and performed experimental validation [1]. The aforementioned scholars have conducted many works for the reliability of machine tools, which shows that the development of machine tool reliability analysis methods is critical. High accuracy reliability prediction methods can improve the competitiveness of machine tools [12]. Detailed descriptions of this aspect are also provided in this paper.
In the paper, a thermal error modeling of spindle unit and machining accuracy reliability analysis method for CNC machine tools is presented, which has the following advantages: The first advantage is the establishment the thermal error mod-1.
el of the electric spindle based on heat generation mechanism and IA. Besides, the heat generation power and heat transfer coefficient were optimized by IA. This model can reveal the inner mechanism of the thermal error and screen the key factors compared to the current empirical methods and temperature-thermal error neural network methods, which can provide the methodological support for predicting the spindle thermal error. The second advantage is the development of the machining 2.
accuracy model by MBS, which considers the comprehensive influence of geometric and thermal errors. This model extends the prediction of machining accuracy from considering a single error source of geometric errors to a multiple error sources of geometric and thermal errors; The third advantage is the proposal of the machining accuracy 3. reliability analysis method by LHSMC, which considers the multiple error sources and the randomness of errors. What's more, the model expands the machining accuracy reliability assessment from the traditional static category to the dynamic category, which is closer to the actual machining process and can provide a valid theoretical guide to analyze the machining accuracy reliability of machine tools.
The other contents are arranged as follows. Section 2 gives a description of the thermal error modeling and machining accuracy modeling process. In Section 3, the machining accuracy reliability model based on LHSMC is established. Then, a four-axis machine tool is utilized to validate this model through experiments, by comparing with AFOSM in Section 4. The conclusion of the paper is given in section 5.

Dynamic machining accuracy modeling considering
the thermal error of spindle unit

Heat generation mechanism analysis
As the core of the machine tool, the spindle unit contains many heat generating components. Rotating parts such as bearings and motor inevitably generate heat due to the friction under the high speed rotation. When the heat accumulates, the spindle unit generates the deformation. To decrease the large deformation induced by high temperature rise, the spiral cooling water jackets are placed outside the heat generating components to accelerate the heat dissipation process. When the electric spindle unit is working, the existence of temperature gradients internal to the electric spindle unit drives the heat transfer from the high temperature area to the low temperature area. The components in contact with the air transfer heat to the air in the form of thermal radiation. The heat transfer process between the coolant and the heat generating components carries away most of the heat, as shown in Fig. 1.
It is assumed that the effect of the change in ambient temperature is not considered. According to the law of energy conservation, the rela- where pro Φ is the component heat generation power, thr Φ is the heat conduction through the continuous cross section of the spindle unit, dis Φ is the coolant heat exchange power, sol c is the specific heat capacity coefficient, t is the heat transfer time, θ is the thermal conductivity of the spindle unit, and A is the area of the heat exchange surface.

Thermal error modeling based on IA
Since the electric spindle unit consists of many heat generating components, the process of thermal error modeling can be categorized as following: 1) Front and rear bearing heat generation power modeling The heat generation power of the bearings operating under the load is expressed as: where 1 Q is the heat generation of rotor, 2 Q is the heat generation of stator, 1 V is the volume of rotor, and 2 V is the volume of stator.

3) Coolant heat transfer coefficient modeling
Coolant is an important medium in the cooling process. The heat transfer process internal to the spindle unit contains three parts: the direct contact heat exchange between the rotating components and the ambient air; the heat exchange between the bearings and the coolant; and the heat exchange between the motor and the coolant. The heat transfer coefficient is expressed as follows: where Re is the Reynolds value, Pr is the Prandtl value, and a l is the cross-sectional circumference of the spindle unit.

4) Optimization of heat generation power and heat transfer coefficient
Generally, there is a deviation between the theoretical values and practice values acquired from Eq. (2) to Eq. (5). To narrow the gap between them, the optimal heat generation power and heat transfer coefficient are found based on IA. The iterative optimization process of IA is implemented by the operators, including the operator for evaluating affinity, the operator for evaluating antibody concentration, and the operator for calculating the degree of excitation and so on.
The equation for parameter optimization can be written as: are the heat generation power of the front bearing, rear bearing and motor, respectively; are the forced and natural convection heat transfer coefficient, respectively; are the proportionality factor of heat generation power and heat transfer coefficient, and are the deviation modification factor of heat generation power and heat transfer coefficient.
The vector form of optimized variables can be written as: The affinity indicates the binding strength of the immune cells to the antigens, and the antibody affinity function based on the Euclidean distance is expressed as follows: i j aff q q is the affinity between antibodies, , i k q and .
j k q are the k-th dimension of antibody i and the k-th dimension of antibody j, and r is the number of antibody dimension.
To solve the optimization problem in this paper, Eq. (8) can be transformed into: where F T , M T and B T are the front bearing, rear bearing and motor simulation temperature, respectively; and The antibody concentration reflects the quality of antibody population, and the high concentration indicates the presence of many similar individuals which inhibits the global optimization procedure. To guarantee individuals have the characteristics of diversity, the antibody concentration is defined as follows: Den q S q q L ff q q and S q q ff q q where L is the population size, (, ) i j Sq q is the similarity between antibodies, and s δ is the similarity threshold.
The incentive degree is an evaluation indicator for antibodies considering both affinity and concentration, and then the next generation of antibodies is screened. It can be described as: where () i sim q is the incentive degree of the antibodies i g , and a, b are calculation parameters depending on the actual situation.
The cloning process was performed by selecting the antibodies with the high incentives in the antibody population, which is represented as follows: where clone(q i ) is the set consisting of m cloned antibodies identical to q i , and m is the number of clones.
By updating the antibodies with low incentive degree in the population and replacing them with new antibodies generated randomly, the global search is achieved. The optimization process is shown in Fig. 2, which can be summarized as follows: the basic parameters are set such as the immune algebra G=200, the variation probability P m =0.7, and the incentive degree factor α =2; the initial range of variation of the variables k and b are set, and the random values of k and b are dynamically generated for the calculation of Eq. (6)；the values calculated by Eq. (6) are used in the finite element simulation, the simulated temperature values and the experimental temperature values are used for the calculation of Eq. (7); and the optimal values of k and b are obtained by circular iterations.
5) The finite element simulation of spindle unit On fluid-solid contact surface, the heat generation and the heat transfer processes of the motor and bearings can be described as: where u ϑ is the energy carried per unit mass,

Machining accuracy modeling based on MBS
MBS serves as a complete abstraction of complicated mechanical systems, which can abstract the four-axis machine tool to be a system with multiple independent bodies. The machining accuracy model is then established by determining the position relationship between each body. The kinematic model of the four-axis machine tool is shown in Fig. 3. Fig. 3(a) shows the overall structural model of the machine tool. Fig. 3 (b) shows the topological framework of the machine tool. Fig. 3(a)  The position change relationship between each two adjacent bodies is described by the homogeneous coordinate transformation matrices. The motion of any individual in the machine tool can be decomposed into two sub-motions, namely rotation around the axis and translation along the axis, so six errors are generated. Taking the Z-axis as an explanation, the matrix of error homogeneous transformation is written as follows: 16 1 The four-axis machine tool includes X, Y, Z translation axes and B rotation axis. So a total of 30 geometric errors are generated. They are given in Table 1.
The forming point of the tool is in its own coordinate system as: Similarly, the forming point of the workpiece is in its own coordinate system as:

Fig. 2. Optimization process of immune algorithm
Ideally, their forming point coordinates in their own coordinate systems will overlap, so the machining accuracy can be written as: 16 16 where p, s are the static state and motion state, and , , , However, in practice, the geometric and the thermal errors can both affect the machining accuracy. The homogeneous transformation matrices between adjacent bodies are given in Table 2.
In summary, the machining accuracy model in the state of actual machining can be written as: 16 16 16 are the homogeneous transformation matrices of the static error and motion error, respectively.

Preliminary
The machining accuracy reliability of the machine tool reflects its ability to achieve a specific function under specified conditions in a predetermined time period [26]. Here, it is regarded as a criterion to evaluate the merits of the four-axis machine tool. Besides, a machining accuracy reliability prediction method is presented based on LHSMC.
The process of LHSMC can be organized as follows: The sample size is determined by Latin hypercube sampling (LHS); the cumulative distribution curve is divided into equal intervals on the cumula- (2) One sample is extracted from each of the K subintervals as to variable V j . Then, only one random number is generated as to each interval, and it is taken as the representative value of the interval. As to the i-th interval, when the representative value u i is chosen randomly and the random number u is generated in U

3-4 Spindle
when it is selected at the center of the interval.
(3) The K sample values of variable V j are randomly sorted according to the ordinal number of the interval to which they belong, and they are placed together in the order of the variables. It is equivalent to constructing a sampling ordinal matrix [ ] ij K n r × = R with K rows and n columns. Besides, the variables are in the column order, and each column is a random and unique arrangement of the ordinal numbers 1, 2,..., K , and the K samples of each variable are arranged by numbers in the columns.
The sequential matrix R is randomly generated, and the columns which introduce the statistical correlations affect the simulation results. The Spearman coefficient of the order matrix R is employed to reduce the statistical correlation, described by [ ] S ijs n n q × = ρ . The Spearman rank correlation coefficients of the i-th and j-th columns can be obtained by: where ijs q is a symmetric matrix, and it is equal to the unit matrix I n in the case of uncorrelated columns. There are no columns with the same sort in R, so the matrix S ρ is positive definite.
The LHSMC reflects the variable distribution characteristics by using the few samples, which can effectively simulate the failure possibility of machine tools.

Machining accuracy reliability analysis
When a part of the machine structure or the whole exceeds the specified state and cannot operate in accordance with the desired functional requirements, the specified state is a limit state which is critical to judge the structure reliability [28]. To analyze the structure reliability, it is necessary to determine whether it has reached its limit state. This section studied the machining accuracy reliability under the stochastic uncertainty of errors, and the performance function is represented as: contains the n basic random variables affecting the machining accuracy reliability, 0 < η means the structure is in a failure state, and 0 = η means the structure is in a limit state.
Function η is a continuous random variable, and the failure possibility can be expressed as: where ( ) f z η is the possibility density function of η , and ( ) F z η (25) Then, the failure possibility can be described by the following equation: Up to now, a machining accuracy reliability analysis method for the machine tool based on LHSMC is formed. The proposal procedure of the method is illustrated in Fig. 4.

Validation of thermal error
The four-axis machine tool is taken as a research example, and its basic parameters are given in Table 3. Through analyzing the heat generation mechanism, the heat generation power and heat transfer coefficient were determined as the characteristic parameters related to the thermal error. Then, they were optimized by IA, and they were employed as the thermal load and boundary conditions of the steadystate thermal analysis of the electric spindle unit in finite element analysis software. Table 3. Basic parameters of the machine tool

Motor power 22kW
Movable distance in X-axis 1400mm Movable distance in Y-axis 1000mm Movable distance in Z-axis 900mm B-axis rotation 0°-360°T he electric spindle unit is assembled from many parts. Inevitably, there are small features inside the spindle unit. To decrease the complexity of the simulation, the spindle system was simplified by ignoring small chamfers, mounting holes and screw holes. The interference fit between parts is not considered. The bearing balls are regarded as rings for finite element simulation purpose. The simplified spindle unit is shown in Fig. 5(a), and Fig. 5(b) shows the half-section structure of the spindle unit, which shows the internal components and position relationship of the spindle unit. The three-dimensional model of the spindle unit is imported into the analysis software, and its finite element mesh model is obtained by the automatic mesh division method, as shown in Fig. 5(c) and Fig. 5(d). In addition, the material of the spindle unit structure is defined in the material library of analysis software.
The following conditions are set before performing the steady-state thermal temperature field simulation: the material parameters of each component of the electric spindle are set; the environmental and coolant inlet temperature are set at 20°C; the coolant flow rate is set at 5L/ min as a constant; the optimized heat generation power is applied to the motor and bearings as the heat source, and the optimized heat transfer coefficient is applied to the heat transfer surface as the boundary condition. The simulation process is in the steadystate thermal analysis module of the analysis software.
The temperature field distribution of spindle unit in the thermal steady state is shown in Fig. 6. Fig. 6(a) shows the temperature field distribution of main components of the spindle unit. The highest temperature and the lowest temperature of the body are 33.09°C and 20.52°C, respectively. The temperature in the middle of the cooling jacket is higher than the two ends, which is related to the direct contact between the cooling jacket and the heat source. Fig. 6(b) shows the temperature distribution of the heat generating parts such as rotor, stator and bearings of the spindle unit. Fig. 6(c) shows the internal temperature distribution of the spindle unit, which implies the internal heat accumulation in the thermal steady state.  Table 4. The two methods are both used to optimize the characteristic parameters. Table 4 indicates that the temperature rise of different components of the spindle unit was predicted based on IA with less error than PSO, which explains the thermal error modeling based on IA has better prediction effect. The temperature field model was imported into the finite element mechanics analysis module, and the fixed constraints were configured. The deformation of the electric spindle unit was obtained as shown in Fig. 7 and Fig. 8. Fig. 7(a) shows the total deformation distribution of the spindle unit, and the maximum deformation is 1.7446e-5m. The max deformation along the X-, Y-and Z-directions are 1.6799e-5m, 6.5465e-6m and 6.6051e-6m, as shown in Fig. 7(b)-(d). Fig. 8 shows the axial deformation of the core shaft is 1.3078e-5m.

Validation of machining accuracy
The machining accuracy model comprising the geometric and thermal errors was developed in the paper. To validate the model, the thermal deformation of the core shaft was extracted directly from the finite element simulation results as the thermal error. 30 geometric errors of the machine tool were measured by the relevant measuring instruments. Fig. 9 shows the platform and workpiece of the four-axis CNC machine tool. 30 measuring points are selected from the workpiece coordinate system, and their coordinates are: The machining accuracy of each point is calculated based on MBS, and the results are shown in Fig. 10. Fig. 10 describes the predicted and measured values of the machining accuracy in the X-, Y-and Z-directions, respectively, which illustrates that the machining accuracy with thermal error exhibits less deviations from the measured machining accuracy. Hence, the thermal error model established in the paper is verified.

Validation of machining accuracy reliability
The distribution characteristics of errors need to be determined before discussing the machine tools reliability. Normally, the geometric errors are regarded as Gaussian distribution [26]. These 30 measuring points are employed to predict the machining accuracy reliability. To predict the machining accuracy reliability at different coordinates, the Monte Carlo method was used as the simulation standard. The LHSMC and the AFOSM methods were used as the prediction methods, and the deviation of their reliability prediction values from Monte Carlo was shown in Fig. 11. Fig. 11 demonstrates that the prediction value of machining accuracy reliability obtained by LHSMC, AFOSM and Monte Carlo, respectively. It can be noted that the reliability prediction value of LHSMC is closer to Monte Carlo than AFOSM. The results illustrate that the LHSMC method in predicting the machining accuracy reliability has a higher accuracy compared with AFOSM.
To verify the machining accuracy reliability analysis method proposed in this paper, the coordinates covering the entire workpieces are selected according to the dimensions of the workpieces. A total of 300 workpieces are machined, as shown in Fig. 12. The failure possibility

. Thermal deformation nephogram of the spindle unit, a) the total deformation, b) the deformation in X-direction, c) the deformation in Y-direction, d) the deformation in Z-direction
is calculated as the experimental data of machining accuracy reliability. The results are shown in Table 5. Table 5 shows the predicted values of machining accuracy reliability based on different methods, and the comparison with the experimental values verifies that the machining accuracy reliability analysis method proposed in the paper has a higher accuracy than AFOSM.

Conclusion
CNC machine tools are high continuity equipment, and their machining accuracy reliability mainly is affected by the geometric and thermal errors. Therefore, a thermal error modeling and machining accuracy reliability analysis method is presented, and the detailed conclusions comprise as following: Based on the heat generation mechanism and IA, a thermal er-1.
ror model was established, which was used to estimate the thermal error in the actual machining process of machine tools; Applying the MBS theory, a machining accuracy model in-2.
cluding both geometric and thermal errors was established to derive the machining accuracy of machine tools; A machine accuracy reliability prediction method was put for-3.
ward considering the spindle thermal error. The effectiveness and superiority of the method were experimentally demonstrated. It can provide the methodological support for spindle unit thermal error modeling and dynamic machining accuracy reliability prediction of machine tools.