Obtaining More Accurate Thermal Boundary Conditions of Machine Tool Spindle Using Response Surface Model Hybrid Artificial Bee Colony Algorithm

: With the application of ﬁnite element method on structure design and engineering analysis more and more widely, this paper presents a response surface model hybrid artiﬁcial bee colony method to optimize the thermal boundary conditions in ﬁnite element thermal analysis of a machine tool spindle to improve its ﬁnite element simulation precision. Initially, the thermal experiment and ﬁnite element thermal analysis of the machine tool spindle with initials that were calculated by empirical formulas were conducted, respectively. Additionally, focusing on thermal boundary conditions, a response surface model is designed to establish the explicit expression between thermal boundary conditions and the simulation errors; then, an artiﬁcial bee colony algorithm is used to solve the mixed-variable optimization problems of a response surface model. Finally, the optimized thermal boundary conditions are brought into the ﬁnite element method of a machine tool system, and the simulation accuracy has been greatly improved .


Introduction
The spindle system, as the core component of machine tools, is crucial for machining. In high-speed milling, such as computer, communication, and consumer electronics (3C) product processing, mechanical spindles are widely used. However, the thermal conditions of spindle, which are caused by the heat flow from internal and external sources, seriously affect the position accuracy of machine tools and the dimensional accuracy of manufactured parts [1][2][3][4][5]. According to statistics, the responsibility of thermal error on the total geometrical error of workpiece is as high as 70%. In Li and Liu's research, they found that the thermal error of numerical control lathe can be as high as 80 µm when n = 2000 rpm and NC machining center can be as high as 50 µm when n = 8000 rpm, which is obviously unacceptable for NC machining [6,7]. Therefore, it is very important to study the thermal characteristics of spindle and reduce its influence on the machining precision of machine tools.
At present, there are mainly two methods to study the thermal characteristics of the spindle:experimental method and simulation method; analytical methods are rarely used. Because the mechanical structure of spindle is three-dimensional and its heat transfer and thermoelastic deformation are very complex, it is difficult to obtain the analytical solution of its thermal characteristics. Compared with the experimental method, the simulation method is more economical and has become a better choice for studying spindles' deformation as long as the thermal boundary conditions (TBCs) are correctly defined and the spindle structure is properly meshed. Especially in design stages, simulation methods can avoid costly design modifications based on experimental studies. In addition, the thermal sensitive points of MT, which can be used to establish thermal error compensation model of MT.

Initial Thermal Boundary Conditions (TBCs) of Spindle
The objective of this study research around is a BT30 mechanical spindle of a FANUC a-D14MiA machining center in our lab. The shaft of the spindle is a multi-diameter and hollow part. As shown in Figure 1, the front bearing of the spindle is a pair of angular contact ball bearings, 7009C, whereas the rear bearing is a pair of 7008C, with a double back to back (DBB) assembly way and axial positioning preload. It can bear radial load and axial load in two directions. Under the same internal clearance, the pre-tightening force and stiffness are twice as large as those of double back (DB) combination, and the ultimate axial load is also larger. The heat transfer between the spindle and the environment is mainly completed by forced convection on the surface of the rotating units and natural convection on the surface of the fixed units. Unlike motorized spindles, the motor of the mechanical spindle is set externally, and the angular contact ball bearings are the only heat source. When the spindle unit is in rotational state, (usually more than 20,000 rpm), the generating power of a single bearing can exceed 400 W. The direction of heat flow is shown by the red arrow in the figure. The first step in solving the temperature field of the spindle is to accurately set the TBCs of the spindle.

Bearing Heating Power
By studying the influence of bearing radial load, axial load, moment, lubricant viscosity, and bearing speed on the load distribution of a rolling element, Palmgren [12] summarized the empirical formula for the heating power of rolling bearing. The generated power of bearing friction can be calculated as where H is in W, M is in N.mm and using bearing speed n in units of rpm. To simplify the analysis method, Palmgren considered that the friction torque of bearings mainly consists of two parts, namely, the mechanical friction torque M 1 and the viscous friction torque M v ; M 1 mainly lies on the structural parameters of the bearing itself M 1 can be expressed as where f 1 is a factor related to the bearing type and load, F β is the bearing load, and d m is the average diameter of the bearing where F s is the static equivalent load, C s is the basic static load rating, y is 0.33, and z is 0.0013 for angular contact ball bearings.
where X 0 and Y 0 are values of the spindle bearing and equal to 0.5 and 0.46, respectively. In this research, the bearing group consists of two 7008C and two 7009C bearings installed back-to-back, and the spindle of the vertical machining is set vertically. For the spindle bearings, the inner ring and the shaft are under an interference fit. Hence, an assembly force acts on the bearings in the radial direction, whereas the radial resultant force is zero because the structure of the spindle is symmetrical. The bearings are under moderate preload, and the preload forces of the bearings are 340 N. The forces of gravity that can be withstood by each 7008C and 7009C bearing is 42 N, respectively. The axial loads that can be withstood by each 7008C and 7009C bearing are 382 N, respectively. According to the previous formulas and data, the values of M 1 are shown in Table 1. M v can be obtained by the following formula: where f 0 is a factor related to the type of bearing and lubrication method, v 0 is the kinematic viscosity of the lubricant at corresponding temperature, (mm 2 ·s), and n is the spindle speed (rpm). The spindle bearing is lubricated by Kyodo Yushi Multemp PS No.2 grease, which has a kinematic viscosity of 16 mm 2 /s at 40 • C. It is worth noting that, due to the different dynamic viscosity of the selected grease, bearing type, and axial size of spindle, the thermal elongation of different spindles vary greatly [16][17][18]. Equation (6) is applicable when the spindle speed exceeds 63 rpm and v 0 n ≥ 2000 . For the single-row angular contact bearing with grease lubrication, f 0 = 2, while for a pair of bearing, f 0 = 4. The bearing viscous friction torque can be obtained by the first formula of Equation (6).

CHTCs of Spindle
There are three forms of heat transfer in the spindle system: heat conduction, heat convection, and thermal radiation. Among the above heat transfer forms, heat convection between the spindle and the surrounding environment is the core factor affecting the spindle TF and TD. In this part, the dimensional analysis method is employed to calculate the spindle CHTCs. For the rotating units surfaces of the spindle, (A f 1 − A f 5 in Figure 2), moving at high speed in the air and producing relative movement between the surfaces and nearby air. For the convenience of calculation, the CHTCs of the rollers and the inner ring surfaces of front and rear bearing group are treated as unified values, and the convection coefficients are calculated as forced convection using Nusselt's Equation (7). The stationary units of the spindle, A n1 and A n2 , mainly interact with the surrounding air by natural convection, and the CHTCs is h n = 9.7 W/m 2 ·K, which is provided by reference [19]. This value is the same order of magnitude as the theoretical free convection coefficient on flat plane [20] h f = N uair × λ air d i In Equation (7), r R e = u air ρd i µ air u air = πd i n 60 P r = C p µ air λ air (8) where h f is the forced CHTCs between the rotating units (A f 1 − A f 5 ) and the ambient air, λ air is the thermal conduction of the ambient air; d i is the equivalent feature size of the forced convection surfaces. P r is the Prandtl number of air and R e is Reynolds number. For simplicity, the equivalent diameter of the outer face of shaft can be calculated using formula: The empirical calculation methods of the seven CHTCs occurring on the spindle and the bearing thermal load were shown in Table 2. As shown in Figure 3, the thermal power of bearing increases rapidly with the increase of the spindle speed and the relationship between them is nonlinear. This phenomenon can be easily explained by the principle of bearing kinematics. As shown in Figure 4, ignoring the impact of gyroscopic moment (D m N < 1 × 10 6 ), there are four forces acting on the bearing in total. F a is axial preload force, F i and F o are the acting loads between the rolling elements and the inner and outer rings of bearing, and F c is centrifugal force. These forces will form a closed force polygon when the system is in an equilibrium condition. When the bearing rotation speed increases, the centrifugal force of rolling elements increases too. To remain balanced, F i and F o must also be increased. As a result, the contact load between the roller elements and the raceway increases, the power to overcome friction force increases, and then the heat generating power of bearing increases. In addition, the forced CHTCs on the surface of the spindle rotating unit also increase rapidly with the increase of bearing rotation speed.

Thermal Experiment and Finite Element Thermal Analysis (FETA)
In order to research the thermal characteristics of spindle unit, firstly, the thermal experiment of spindle at n = 20,000 rpm is conducted, as shown in Figure 5. Two PT100 temperature sensors (T 1 and T 2 ) are fixed on the spindle front bearing surface and MT structure to measure the temperature change of the spindle and environment, respectively. The change of spindle and ambient temperature is shown in Figure 6. The temperature recording experiment lasted 18,000 s. It is worth noting that, before the experiment, the dynamic balance of the spindle should be normal; otherwise, the abnormal temperature rise of the spindle will occur, and the ideal experimental results will not be obtained. The environmental temperature in the finite element (FE) software was set to be 17.5 • C, the same as average temperature of environment during measurement. The experiment is set up in the climate chamber, and the temperature change is not significant, thus the impact of environment change is ignored in this research. However, when the ambient temperature changes greatly, its influence on the thermal characteristics of the spindle cannot be ignored. In this experiment, the bearing groups are the only heat source of the spindles, and its temperature rise is different from that of motorized spindles [21,22].  Secondly, we simulated the temperature field and thermal deformation of the spindle unit according to the TBCs obtained by empirical formulas before. It is worth noting that spindle and MT systems are composed of many units. It is not recommended to mesh it directly because this will lead to over dense meshing in some unimportant places, (e.g., bolts) and results in long calculation that can affect the accuracy of the results. To ensure the correctness of the analysis results, the 3D structure of the spindle should be simplified. In this research, the spindle bolt connection, chamfers, and small steps are ignored.To get satisfying FEA results, the spindle region with a larger temperature gradient is meshed to be more refined, such as in the regions near the bearings. There are 676,665 elements and 589,914 nodes of the FE model. The thermal loads of bearings are evenly distributed on the rolling elements and raceways, as shown in Figure 7a. In addition, the TCRs between bearing inner rings and shaft, between bearing outer ring and housing, are both set, as shown in Figure 7b,c, and the values of them can refer to Table 3. (By defining the material properties, surface morphology, and pressure on both sides of the joint, the TCRs of for the joints can be obtained.) In FE software ANSYS Workbench, the TCRs of corresponding contact surfaces can be assigned by its reciprocal value, thermal contact conductance (TCC). In previous TF analysis of spindles, the bearing groups are usually simplified as a thick-walled ring, and the thermal load of the bearing can only be added to the ring wall surfaces. However, according to the heating generation principle of bearings, the heat of bearings is mainly caused by the friction moment of the rolling elements to overcome the grease and the friction between the rolling elements and the inner and outer raceways and cages, so it is too ideal to add the thermal load of bearings directly to the ring wall surfaces. The material properties of the spindle system are set according to Table 4. When n = 3000 rpm, Tan [15] takes a certain type of electric spindle as the research object, and does not consider the instability of empirical formula calculation, the temperature simulation error of the front bearing surface (T1) of the spindle is as high as 54.37%. In fact, we also find that the simulation error is about 50% when the spindle is at n = 3000 rpm when the instability of empirical formula is ignored. However, with the increase of the rotating speed, the simulation error of the temperature field of the spindle will be larger and larger after the boundary condition without optimization is brought into the finite element model. Thus, it reaches 124.5% at n = 20,000 rpm. Therefore, it can be concluded that, with the increase of the spindle speed, the thermal boundary conditions of the spindle obtained by the empirical formula and the simplified equation will deviate from the actual boundary conditions seriously.

FETA Results
From Figure 8a, the results show that temperature of the outer surface of spindle near the front bearings is 62.5 • C. However, the experimental result is 28 • C when the spindle reaches thermal balance, as shown in Figure 6. Because the bearing is heat source, the highest temperature area is at the front bearing groups. In addition, due to the existence of TCRs between the contact surfaces, there is a significant temperature drop between the contact surfaces of the spindle temperature field. Figure 8b shows the axial TD of spindle shaft; the maximum deformation happens to the nose of it, nearly 75.5 µm. The comparison of the simulated and experimental results are shown in Table 5. Obviously, the FEA results deviate from the experimental results seriously, which indicates that the TBCs calculated by the empirical formulas are inaccurate. This can be attributed to the heat generation of bearing being a very complex problem. The calculation value of empirical formula is not always reliable, let alone the CHTCs' empirical formulas which are based on a series of assumptions. According to Chen's research, the calculation result of the formula will be far greater than the actual situation with the increasing of bearing rotational speed [11]. In addition, in the past few years, the manufacturing technology of bearings has improved significantly. Moreover, the performance of bearings with different precision grades and material differs greatly, including heat generation power [23]. For example, the calculated empirical heat power value of ceramic ball bearing is the same as that of a metal ball bearing at the same rotational speed. However, due to the lightweight quality of ceramic bearing rolling elements, small centrifugal force, its heat power is lower than a metal ball bearing at the same speed. Therefore, it will make the spindle temperature's simulation results different from the actual ones by directly bringing the values calculated by empirical formulas to the finite element model. In addition, the health statuses of spindles are different from each other, which will also influence the thermal analysis results of FEA. It is clearly seen that there exist significant differences between them and the maximum relative simulation error was up to 128.1%. Accordingly, it is essential to increase the simulation accuracy, and, in this research, we concentrate on optimizing the TBCs to achieve that goal.

Optimization of Thermal Boundary Parameters
According to the previous analysis, the accuracy of the FE model by directly substituting the parameters calculated by empirical formulas into the FE model is not satisfactory. This can be attributed to the structure of the spindle being simplified and the thermal boundary conditions (TBCs) set in FEM software were inevitably different from the actual situation. This problem is difficult to solve with traditional methods. However, it can be solved by establishing the mathematical model between simulation error and TBCs. Abstract the problem into an optimization problem, and then use intelligent algorithms to solve it. This section will describe how to use an RSM hybrid ABC algorithm to correct the TBCs of the finite element thermal analysis model and improve its simulation accuracy. In this research, the TBCs, , are abstracted as the optimization variables, [x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 ], of the ABC algorithm, and the simulation errors are set as the fitness value. The main steps are described as below: (1) Establishing the simplified digital model of spindle unit based on the design parameters and meshing it. (2) Setting the TBCs according to simplified empirical formula, such as bearing heat power, thermal contact resistances between the contact surfaces, and CHTCs of the spindle units' surfaces. The ambient temperature is 17.5 • C, and spindle rotating speed is 20,000 rpm. (3) According to the simulation error, correct the TBCs based on an RSM hybrid ABC algorithm. (4) Then, the optimized TBCs are substituted to the FEM for analysis and compared with the experimental results.

DE and RSM
To establish the explicit expression between TBCs and the simulation errors with RSM, we need to design an experiment (DE) first, and the experimental design method should be determined according to the characteristics of the objects. Ref. [24] If the experimental design is not reasonable, it will be difficult to get the ideal results. For example, for a 9-factor and 3-level experiment, if one experiment is conducted at each level of each factor, the total number of experiments is at least 3 9 = 19,683 times. It is a very time-consuming work to collect the data from experiments. To obtain the relationship between TBCs parameters, usually, the commonly used experimental methods in RSM are central composite design (CCD) and Box-Behnken design (BBD) method, which will greatly reduce the number of experiments and simplify the analysis of experimental data. For the BBD method, the number of experimental factors is usually 3 to 7. When the experimental factor increases, the CCD method is usually used. In this work, there are nine experimental factors in total, so the CCD method, which can reflect the intrinsic characteristics of the entire design space with fewer sample points, is employed to design the experiment. RSM is a method that replaces the implicit constraints in the original design problem by constructing explicit approximate expressions. Usually, the response surface approximation function used in engineering is a second-order model, although the higher the order of the polynomial, the higher the fitting degree between the response surface equation and the actual data. However, the higher the order, the greater the accumulation of rounding errors in the calculation process. Ref. [25] Therefore, when the order is too high, the accuracy of the equation will decrease, or even the reasonable results can not be obtained. The second-order response surface model can be expressed as below: where Y denotes the dependent variable; ε denotes the normal random error; k is the number of design variables; β 0 , β i , and β ij denote the undetermined coefficients, respectively; β i denotes the linear effect of x i , β ij denotes the interaction effect between x i and x j , β ii denotes the secondary effect of x i . N shows that ε obeys standard orthogonal distribution. In this research, the objective function is the simulation error of the simulated and experimental values of three temperature measuring points at a certain time (thermal steady state). Assuming that the simulated and measured temperatures can be expressed as T and T respectively, the objective function can be expressed as below: Therefore, the corresponding values of min ( f Temp ) are the optimal approximation values of f Temp . The parameter modification problem of spindle FE model can be transformed into the optimization problem of formula (12): The variable parameters to be optimized have been shown in Table 2. There are nine design variables in the RSM in total, as shown in Table 6, and the design variables [x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 ] have been converted to standardized levels. For such a five level, nine variable parameters, they will have 540 groups of experiments for full implementation in CCD, which is extremely time-consuming. To save time and ensure the reliability of the experiment, 1/4 of the experiment was carried out, and its response values are listed in Table 7.  Table 7. CCD method of the spindle system.
In this case, x 1 , x 2 , x 3 , x 5 , x 6 , x 7 , x 8 , x 9 , x 6 x 9 are significant model terms. According to the different influence degree of each design variable on the output variable, the two parameters with the greatest influence degree are selected. Based on the experimental design data points, the response surface contour of x 6 x 9 and the output variable are fitted, as shown in Figure 9. . RSM between X 6 X 9 and F Temp .

ABC
According to the obtained RSM, the TBCs and the simulation error are treated as the variables to be optimized and the objective function value respectively, thus the searching for optimal TBCs can be equivalent to a constraint optimization problem and various optimization methods can be used to solve it, such as Particle Swarm Optimization (PSO) and genetic algorithm (GA). PSO is the abbreviation of particle swarm optimization. It is a kind of stochastic optimization technology based on population. It was proposed by Eberhart and Kennedy in 1995. PSO mimics the swarm behavior of insects, herds, birds, and fish. These groups search for food in a cooperative way. Each member of the group changes its search mode by learning its own experience and the experience of other members. GA is a computational model simulating the natural selection and genetic mechanism of Darwinian biological evolution. It is a method to search the optimal solution by simulating the natural evolution. However, the performance of the conventional PSO and GA method will become sensitive to the mutation strategy and the associated parameters; in addition, it has premature convergence, dimension disaster, and can be easily trapped in local extreme value. Refs. [26][27][28] Moreover, these optimization methods usually need generous assessments of finite element thermal analysis, which is not economical because a single assessment of finite element thermal analysis is always time-consuming. ABC is an optimization algorithm based on the intelligent behavior of honey bee swarm, compared with PSO and GA, and the ABC algorithm has the advantages of simple operation, less control parameters, high search accuracy, and strong robustness [29,30].
It has been proved by Griewank function, Rastrigin function, Rosenbrock function, Ackley function and Schwefel function, whose value is 0 at its global minimum [30]. By comparing the optimization results of ABC, PSO, and GA algorithm for these functions, the global minimum values of the ABC algorithm are 0.000329, 0, 0.012522, 4.6×10 −11 , and 27×10 −9 , respectively, in 500 search cycles. The global minimum of PSO algorithm is 0.079393, 2.6559, 4.3713, 9.8499×10 −13 , and 161.87. The global minimum values of GA algorithm are 0.050228, 1.3928, 46.3184, 0.59267, and 1.9519, respectively. The comparison shows that ABC algorithm is better than GA and PSO algorithms. In this paper, the RSM hybrid ABC method is proposed to optimize the TBCs of the FE model, and its main steps are described below:

Comparison of FEA Results of Spindle before and after Optimization
In this section, the TBCs optimized by the proposed RSM assisted ABC method were given and brought into the finite element model of spindle unit. Then, compared with the previous analysis results which directly brought into empirical formulas, the effectiveness of the optimization method is proved.

The Optimized TBCs
As shown in Figure 10, the optimization process of ABC can be seen. After 42 iterations, the optimum solutions are obtained and given in Figure 11, reflecting the high efficiency of the algorithm. All the optimized TBCs deviate from the calculated value of empirical formula to different extents. x 1 , x 2 , and x 6 are much smaller than the empirical values, x 3 and x 4 are much larger than the empirical values, and x 5 , x 7 , x 8 , and x 9 are close but less than the empirical values. It is worth noting that, because of the randomness of the ABC intelligent algorithm, this group of solutions is not unique. When the algorithm is run repeatedly, other solutions can be obtained.
Initialize system parameters，NP=20,limit=100,maxCycle=500 Compare food positions and retain the best one • Step1: Initialize system parameters and the population of solutions X i , and Step2: Evalute the population, which is to take the value of X i into f Temp and compare the fitness value. • Step3: Cycle = 1 • Step4: repeat • Step5: Produce new solutions ti for the employed bees by using v ij = x ij + φ ij (x ij − x kj ) and evaluate them • Step6: Apply the greedy selection process for the employed bees. • Step7: Calculate the probability values P i for the solutions X i by p i = f it i / ∑ NP n=1 f it n . • Step8: Produce the new solutions v ij for the onlookers from the solutions X i selected depending on P i and evaluate them. • Step9: Apply the greedy selection process for the onlookers. where NP = 20, limit = 100, maxCycle = 150, and D = 9. D is the number of optimization parameters, where SN denotes the size of employed bees or onlooker bees, and maxCycle denotes the maximum cycle number. The flowchart is shown in Figure 10.

Comparision of Finite Element Thermal Analysis Results
Since the optimal TBCs have been searched, as shown in Figure 12, then we brought it into the FE model of spindle and calculated its TF&TD. Meanwhile, the comparison between the simulation results and the experimental results of the average temperature of T 1 point is shown in Figure 13, compared with the TD and TF in Figure 8, which ignored the calculation error of TBCs. Obviously, the optimized TBCs greatly improved the simulation accuracy of the FEM of spindle unit to 4.7%. As shown in Table 8, the correctness and effectiveness of the RSM assisted ABC method's application on TBCs optimization are verified.  Although the RSM assisted ABC algorithm improved the simulation accuracy of the FEM of spindle unit to some extent, there are still errors between the simulation value and the experimental value due to the environmental temperature change, thermal radiation, unreasonable structural simplifications, and other factors, such as nonlinearity of thermal deformation of materials and thermal radiation, which were not considered in this research. In addition, the FEM software itself has calculation errors, so it is very difficult to completely match the simulation value with the experimental value. Actually, it is impossible to completely eliminate the thermal error of the machine tool, but it is very meaningful to predict and compensate the thermal error by some methods and control the thermal error within the acceptable range.

Conclusions
This paper presented an RSM hybrid ABC method to optimize the TBCs in finite element thermal analysis of an MT system to improve its simulation accuracy. Firstly, the thermal experiment, FETA, of the spindle with initial TBCs which was calculated by empirical formulas was conducted and compared, respectively. It is found that the simulation error is notable. Secondly, focusing on TBCs, the RSM hybrid ABC method was introduced. RSM is used to establish the explicit expression between TBCs and the simulation error, and ABC algorithm is employed to look for the global optimal solution of TBCs. Finally, the optimized TBCs are brought into the spindle FEM of spindle, and the simulation accuracy of spindle unit is improved to 4.7%. This method can be used not only to assist in the design stage of the spindle to reduce design cost and design time, but also to help engineers find the thermal sensitive points and establish the thermal error model. The thermal error can be compensated in advance by the numerical control system. However, the thermal dynamic characteristics of the spindle are also worth studying. In the second section, we simply analyzed the relationship between the bearing contact load and rotation speed. With the increase of the rotation speed, the heat generating power of the bearing increases rapidly due to the work to overcome friction. In addition, due to the uneven heating of the components of spindle unit, different degrees of thermal deformation will occur, and the support stiffness of bearing will also change, thus changing the natural frequency of the spindle unit. When the working frequency of the spindle coincides with the natural frequency, the spindle will produce a large resonance, which will directly affect processing quality of product, and even endanger the personal safety.