Next Article in Journal
Multi-Granulation Picture Hesitant Fuzzy Rough Sets
Previous Article in Journal
Designing a Supermarket Service Robot Based on Deep Convolutional Neural Networks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

School of Mechanical Science and Engineering, Huazhon University of Science and Technology, Wuhan 430074, China
*
Author to whom correspondence should be addressed.
Symmetry 2020, 12(3), 361; https://doi.org/10.3390/sym12030361
Submission received: 5 February 2020 / Revised: 24 February 2020 / Accepted: 27 February 2020 / Published: 2 March 2020

Abstract

:
With the application of finite element method on structure design and engineering analysis more and more widely, this paper presents a response surface model hybrid artificial bee colony method to optimize the thermal boundary conditions in finite element thermal analysis of a machine tool spindle to improve its finite element simulation precision. Initially, the thermal experiment and finite 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 artificial 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 finite element method of a machine tool system, and the simulation accuracy has been greatly improved.

1. 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 experimental method can only measure the temperature of limited points in the temperature field of a spindle unit, while the simulation method can comprehensively analyze the distribution of the temperature field of spindle without blind angle. In recent years, finite element method (FEM) has been frequently used by researchers to study the thermal problems of machine tools. Zhao et al. [8] studied the thermal errors of a turning center spindle at 2000 rpm by FEM, and recorded the thermal errors change at the spindle nose and temperature change at the selected points. Ultimately, the results show that the simulation results are satisfying to replace the experiment results. In Liu [9] and Ma’s [10] research, they considered the influence of thermal contact resistances (TCRs) on the temperature field & thermal deformation (TF&TD) of spindle system in the FE simulation, and verified the reliability of simulation results by experiments. Satisfactory results have been achieved in the above studies.
However, without a single exception, these studies are all carried out at relatively low spindle speeds and the TBCs in FE analysis are obtained directly from empirical formulas. According to Chen’s research, ref. [11] the calculated heat power of bearing by empirical formulas will exceed the actual value with the increase of rotational speed. The commonly used empirical formulas summarized 70 years ago for calculating bearing heat generation power are applicable to medium and low speed [12]; it is recommended to be reevaluated with the increase of bearing speed. In the past few years, the manufacturing technology of bearings has been constantly improved and the manufacturing accuracy has been improved too. For example, high-speed bearing is widely used in the machining of spindle nowadays, whose rolling element is made of ceramic material, which is lighter in weight, smaller in centrifugal force, and lower in heat than traditional steel balls, and its value of dn can reach 4 × 10 7 . It seems too ideal to estimate the heating power of bearings based on the past heat generating power of bearings. Moreover, the performance of bearings with different precision grades differs greatly, including heat generation power. Therefore, the empirical values of bearing generation power and convective heat transfer coefficients (CHTCs) will make spindle temperature’s simulation results different from the actual ones. In addition, the health status of spindles are different from each other, the CHTCs on the surfaces of the spindle are calculated on the basis of many assumptions, which will also have an impact on the results of analysis. In the study of the thermal characteristics of a lathe spindle, for the uncertainty of TBCs, Li et al. [13] introduced an inverse heat conduction theory to optimize the CHTCs of spindle system, and then get a satisfactory simulation results of thermal characteristics of spindle in ANSYS software. In addition, they also proposed a mean impact value method to select the thermal key points in the spindle system. Zhang et al. [14] introduced a biogeography optimization algorithm to optimize the heat transfer coefficient of a motorized spindle, and got an accurate predication of thermal deformation of motorized spindle using the experimental data of the spindle surface temperature. The average prediction error of thermal deformation of spindle is 0.72 μ m. Compared with the traditional prediction model of thermal deformation of spindle, the model has higher accuracy. Tan et al. [15] used a surrogate assisted differential evolution method to seek more accurate CHTCs in thermal analysis of spindle; then, the optimized CHTCs were brought into the finite element model of spindle, which improves the simulation accuracy of the finite element model, and its validity was verified by experiments. However, in the above optimization analysis, only CHTCs are considered, and the unreliability of thermal power of spindle bearing calculated by empirical formula was not considered. This being taken into account in this paper because of the thermal power of bearing is a very important boundary condition in the TD and TF analysis of spindle.
Since the simulation accuracy of spindle’s thermal characteristics can be improved by optimization algorithms. Based on the previous studies, a 9-variable, 1-output experiment was designed by a central composite design (CCD) method in this research, and an response surface model (RSM) is designed to establish the explicit expression between TBCs and the simulation errors; then, the response surface function is employed as the fitness function of ABC and the optimized solutions are obtained. Finally, the simulation accuracy of TF of spindle unit is improved, and the good effectiveness of the proposed RSM and ABC method is verified. The method proposed in this research can be used not only to assist in the design stage of spindle to reduce design cost and design time, but also to help engineers to find the thermal sensitive points of MT, which can be used to establish thermal error compensation model of MT.

2. 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.

2.1. 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
H = 1.047 × 10 4 × M × n
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 = M 1 + M v
M 1 can be expressed as
M 1 = f 1 × F β × d m
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
f 1 = z × ( F s C s ) y
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.
F β = F α 0.1 × F r F s = X 0 × F r + Y 0 × F α
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:
M v = 10 7 × f 0 × ( v 0 × n ) 2 3 × d m 3 , v 0 n 2000 M v = 160 × 10 7 × f 0 × d m 3 , v 0 n 2000
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).

2.2. 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 n 1 and A n 2 , 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 u a i r × λ a i r d i
In Equation (7),
N u a i r = 0.133 R e 2 3 P r 1 3 R e = u a i r ρ d i μ a i r u a i r = π d i n 60 P r = C p μ a i r λ a i r
where h f is the forced CHTCs between the rotating units ( A f 1 A f 5 ) and the ambient air, λ a i r 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:
d s = d 1 L 1 + d 2 L 2 + d 3 L 3 + d 4 L 4 L 1 + L 2 + L 3 + L 4
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.

2.3. 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.

2.4. 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.

3. 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, [ A f 1 A f 2 A f 3 A f 4 A f 5 A n 1 A n 2 H 1 H 2 ], 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.

3.1. 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:
Y = β 0 + i 1 k β i x i + i j β i j x i x j + ε ε N ( 0 , σ 2 ) , i = 1 , 2 , · · · , k , j = 1 , 2 , · · · , k
where Y denotes the dependent variable; ε denotes the normal random error; k is the number of design variables; β 0 , β i , and β i j denote the undetermined coefficients, respectively; β i denotes the linear effect of x i , β i j denotes the interaction effect between x i and x j , β i i 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:
f T e m p ( x ) = ( T 1 T 1 ) 2
Therefore, the corresponding values of min ( f T e m p ) are the optimal approximation values of f T e m p . The parameter modification problem of spindle FE model can be transformed into the optimization problem of formula (12):
F i t n e s s = m i n ( f T e m p )
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.
The expansion of the second order RSM can be expressed as:
f T e m p ( x 1 , x 2 , · · · , x 9 ) = β 0 + β 1 x 1 + β 2 x 2 + β 3 x 3 + β 4 x 4 + β 5 x 5 + β 6 x 6 + β 7 x 7 + β 8 x 8 + β 9 x 9 + β 11 x 1 2 + β 22 x 2 2 + β 33 x 3 2 + β 44 x 4 2 + β 55 x 5 2 + β 66 x 6 2 + β 77 x 7 2 + β 88 x 8 2 + β 99 x 9 2 + β 12 x 1 x 2 + β 13 x 1 x 3 + β 14 x 1 x 4 + β 15 x 1 x 5 + β 16 x 1 x 6 + β 17 x 1 x 7 + β 18 x 1 x 8 + β 19 x 1 x 9 + β 23 x 2 x 3 + β 24 x 2 x 4 + β 25 x 2 x 5 + β 26 x 2 x 6 + β 27 x 2 x 7 + β 28 x 2 x 8 + β 29 x 2 x 9 + β 34 x 3 x 4 + β 35 x 3 x 5 + β 36 x 3 x 6 + β 37 x 3 x 7 + β 38 x 3 x 8 + β 39 x 3 x 9 + β 45 x 4 x 5 + β 46 x 4 x 6 + β 47 x 4 x 7 + β 48 x 4 x 8 + β 49 x 4 x 9 + β 56 x 5 x 6 + β 57 x 5 x 7 + β 58 x 5 x 8 + β 59 x 5 x 9 + β 67 x 6 x 7 + β 68 x 6 x 8 + β 69 x 6 x 9 + β 78 x 7 x 8 + β 79 x 7 x 9 + β 89 x 8 x 9
Assuming that x 1 2 = x 10 , x 2 2 = x 11 , · · · , x 9 2 = x 18 , x 1 x 2 = x 19 , · · · , x 1 x 9 = x 26 , x 2 x 3 = x 27 , · · · , x 2 x 9 = x 33 , x 3 x 4 = x 34 , · · · , x 3 x 9 = x 39 , x 4 x 5 = x 40 , · · · , x 4 x 9 = x 44 , x 5 x 6 = x 45 , · · · , x 5 x 9 = x 48 , x 6 x 7 = x 49 , · · · , x 6 x 9 = x 51 , x 7 x 8 = x 52 , x 7 x 9 = x 53 , x 8 x 9 = x 54 . Then, Equation (13) can be rewrittrn as a multiple linear regression model:
f ( x 1 , · · · , x 9 ) = β 0 + β 1 x 1 + β 2 x 2 + · · · + β 9 x 9 + · · · + β 51 x 51 + β 52 x 52 + β 53 x 53 + β 54 x 54
The total number of the test is n = 156 , and then the RSM can be rewritten as matrix form
Y = X β + ε
where
Y = y 1 y 2 · · · y n ; X = 1 x 11 x i 1 1 x 12 x i 2 1 x 1 n x i n
β = β 1 β 2 · · · β n ; ; ε = ε 1 ε 2 · · · ε n ;
According to the standardized variables shown in Table 6, the coefficient vector β of the RSM can be obtained. By substituting β into Equation (13), the RSM can be obtained. The determination coefficient of f T e m p is R 2 = 0.89 , which reflects the outstanding fitting performance of the RSM model.
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.

3.2. 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:

4. 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.

4.1. 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.
  • Step1: Initialize system parameters and the population of solutions X i , and X i = [ x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 ] i
  • Step2: Evalute the population, which is to take the value of X i into f T e m p and compare the fitness value.
  • Step3: Cycle = 1
  • Step4: repeat
  • Step5: Produce new solutions ti for the employed bees by using v i j = x i j + ϕ i j ( x i j x k j ) 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 i t i / n = 1 N P f i t n .
  • Step8: Produce the new solutions v i j 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.
  • Step10: Determine the abandoned solution for the scout, if it exists, and replace it with a new randomly produced solution X i by x i j = x m i n j + r a n d [ 0 , 1 ] ( x m a x j x m i n j ) .
  • Step11: Memorize the best solution achieved so far.
  • Step12: cycle = cycle + 1.
  • Step13: until cycle = MCN.
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.

4.2. 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.

5. 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.

Author Contributions

Conceptualization, L.Z.; formal analysis, L.Z.; methodology, L.Z., J.X.; writing—original draft, L.Z., J.X.; wrinting—review and editing, J.X., T.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research is supported by the National Natural Science Foundation of China (Grant Nos. 51575202, and 51675204), Science Challenge Project (Grant No. TZ2018006-0102-01), and the National Science and Technology Major Project of China (Grant No. 2018ZX04035002-002).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Jedrzejewski, J.; Kwasny, W. Modelling of angular contact ball bearings and axial displacements for high-speed spindles. CIRP Ann. 2010, 59, 377–382. [Google Scholar] [CrossRef]
  2. Kim, J.D.; Zverv, I.; Lee, K.B. Thermal model of high-speed spindle units. KSME Int. J. 2003, 2, 10. [Google Scholar] [CrossRef] [Green Version]
  3. Yan, B.; Yan, K.; Luo, T.; Zhu, Y.; Li, B.Q.; Hong, A.J. Thermal coefficients modification of high speed ball bearing by multi-object optimization method. Int. J. Therm. Sci. 2019, 137, 313–324. [Google Scholar] [CrossRef]
  4. Li, Y.; Zhao, W.; Lan, S.; Ni, J.; Wu, W.; Lu, B. A review on spindle thermal error compensation in machine tools. Int. J. Mach. Tools Manuf. 2015, 95, 20–38. [Google Scholar] [CrossRef]
  5. Liu, H.; Miao, E.; Zhuang, X.; Wei, X. Thermal error robust modeling method for cnc machine tools based on a split unbiased estimation algorithm. Precis. Eng. 2018, 51, 169–175. [Google Scholar] [CrossRef]
  6. Li, Q.; Li, H.L. A general method for thermal error measurement and modeling in cnc machine tools’ spindle. Int. J. Adv. Manuf. Technol. 2019, 103, 2739–2749. [Google Scholar] [CrossRef]
  7. Liu, K.; Liu, Y.; Sun, M.; Li, X.; Wu, Y. Spindle axial thermal growth modeling and compensation on cnc turning machines. Int. J. Adv. Manuf. Technol. 2016, 87, 2285–2292. [Google Scholar] [CrossRef]
  8. Haitao, Z.; Jianguo, Y.; Jinhua, S. Simulation of thermal behavior of a cnc machine tool spindle. Int. J. Mach. Tools Manuf. 2007, 47, 1003–1010. [Google Scholar] [CrossRef]
  9. Liu, Z.; Pan, M.; Zhang, A.; Zhao, Y.; Yang, Y.; Ma, C. Thermal characteristic analysis of high-speed motorized spindle system based on thermal contact resistance and thermal conduction resistance. Int. J. Adv. Manuf. Technol. 2014, 76, 1913–1926. [Google Scholar] [CrossRef]
  10. Ma, C.; Yang, J.; Zhao, L.; Mei, X.; Shi, H. Simulation and experimental study on the thermally induced deformations of high-speed spindle system. Appl. Therm. Eng. 2015, 86, 251–268. [Google Scholar] [CrossRef]
  11. Chen, G.C.; Wang, L.Q.; Le, G.U.; Zhen, D.Z. Heating analysis of the high speed ball bearing. J. Aerosp. Power 2007, 22, 163–168. [Google Scholar]
  12. Harris, T.A.; Kotzalas, M.N. Rolling Bearing Analysis, 5th ed.; CRC/Taylor & Francis: Boca Raton, FL, USA, 2007. [Google Scholar]
  13. Li, Y.; Zhao, W.; Wu, W.; Lu, B. Boundary conditions optimization of spindle thermal error analysis and thermal key points selection based on inverse heat conduction. Int. J. Adv. Manuf. Technol. 2016, 90, 2803–2812. [Google Scholar] [CrossRef]
  14. Zhang, L.; Gong, W.; Zhang, K.; Wu, Y.; An, D.; Shi, H.; Shi, Q. Thermal deformation prediction of high-speed motorized spindle based on biogeography optimization algorithm. Int. J. Adv. Manuf. Technol. 2018, 97, 3141–3151. [Google Scholar] [CrossRef]
  15. Tan, F.; Wang, L.; Yin, M.; Yin, G. Obtaining more accurate convective heat transfer coecients in thermal analysis of spindle using surrogate assisted differential evolution method. Appl. Therm. Eng. 2019, 149, 1335–1344. [Google Scholar] [CrossRef]
  16. Fu, C.B.; Tian, A.H.; Yau, H.T.; Hoang, M.C. Thermal monitoring and thermal deformation prediction for spherical machine tool spindles. Therm. Sci. 2019, 23, 2271–2279. [Google Scholar] [CrossRef]
  17. Li, Y.; Wei, W.M.; Su, D.X.; Zhao, W.H.; Zhang, J.; Wu, W.W. Thermal error modeling of spindle based on the principal component analysis considering temperature-changing process. Int. J. Adv. Manuf. Technol. 2018, 99, 1341–1349. [Google Scholar] [CrossRef]
  18. Yin, Q.; Tan, F.; Chen, H.X.; Yin, G.F. Spindle thermal error modeling based on selective ensemble bp neural networks. Int. J. Adv. Manuf. Technol. 2019, 101, 1699–1713. [Google Scholar] [CrossRef]
  19. Liu, T.; Gao, W.G.; Zhang, D.W.; Zhang, Y.F.; Chang, W.F.; Liang, C.M.; Tian, Y.L. Analytical modeling for thermal errors of motorized spindle unit. Int. J. Mach. Tools Manuf. 2017, 112, 53–70. [Google Scholar] [CrossRef]
  20. Incropera, F.P. Fundamentals of Heat and Mass Transfer; John Wiley Sons: Hoboken, NJ, USA, 2011. [Google Scholar]
  21. Yan, K.; Hong, J.; Zhang, J.; Mi, W.; Wu, W. Thermal-deformation coupling in thermal network for transient analysis of spindle-bearing system. Int. J. Therm. Sci. 2016, 104, 1–12. [Google Scholar] [CrossRef]
  22. Jiang, S.Y.; Mao, H.B. Investigation of variable optimum preload for a machine tool spindle. Int. J. Mach. Tools Manuf. 2010, 50, 19–28. [Google Scholar] [CrossRef]
  23. Jedrzejewski, J.; Kowal, Z.; Kwasny, W.; Modrzycki, W. High-speed precise machine tools spindle units improving. J. Mater. Process. Technol. 2005, 162, 615–621. [Google Scholar] [CrossRef]
  24. Hinkelmann, K.; Kempthorne, O. Principles of Experimental Design; Wiley: Hoboken, NJ, USA, 2007. [Google Scholar]
  25. Galtchouk, L.I. On sequential estimation of parameters for linear regression with martingale noise. In MODA 5 Advances in Model-Oriented Data Analysis and Experimental Design; Atkinson, A.C., Pronzato, L., Wynn, H.P., Eds.; Physica-Verlag HD: Heidelberg, Germany, 1998; pp. 105–114. [Google Scholar]
  26. Yang, J.; Honavar, V.G. Feature Subset Selection Using a Genetic Algorithm; Springer: Boston, MA, USA, 1998. [Google Scholar]
  27. Gozali, A.A.; Fujimura, S. Localization Strategy for Island Model Genetic Algorithm to Preserve Population Diversity; Computer and Information Science; Springer: Cham, Switzerland, 2018. [Google Scholar]
  28. Pant, M.; Thangaraj, R.; Abraham, A. A New PSO Algorithm with Crossover Operator for Global Optimization Problems; Springer: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
  29. Karaboga, D.; Akay, B. A comparative study of artificial bee colony algorithm. Appl. Math. Comput. 2009, 214, 108–132. [Google Scholar] [CrossRef]
  30. Karaboga, D.; Basturk, B. A powerful and efficient algorithm for numericial function optimization: Artificial bee colony (abc) algorithm. J. Glob. Optim. 2007, 39, 459–471. [Google Scholar] [CrossRef]
Figure 1. The mechanical structure of spindle.
Figure 1. The mechanical structure of spindle.
Symmetry 12 00361 g001
Figure 2. Setting of convection heat transfer coefficients.
Figure 2. Setting of convection heat transfer coefficients.
Symmetry 12 00361 g002
Figure 3. Bearing rotating speed and its heat power.
Figure 3. Bearing rotating speed and its heat power.
Symmetry 12 00361 g003
Figure 4. Forces acting on a bearing ball and force equilibrium graph.
Figure 4. Forces acting on a bearing ball and force equilibrium graph.
Symmetry 12 00361 g004
Figure 5. Experiment set up.
Figure 5. Experiment set up.
Symmetry 12 00361 g005
Figure 6. Temperature rise process of spindle.
Figure 6. Temperature rise process of spindle.
Symmetry 12 00361 g006
Figure 7. Temperature field of spindle (a) heat load; (b) TCRs between bearing inner ring and shaft; (c) TCRs between bearing outer ring and housing.
Figure 7. Temperature field of spindle (a) heat load; (b) TCRs between bearing inner ring and shaft; (c) TCRs between bearing outer ring and housing.
Symmetry 12 00361 g007
Figure 8. Thermal characteristics of spindle (a) temperature field; (b) thermal deformation of shaft.
Figure 8. Thermal characteristics of spindle (a) temperature field; (b) thermal deformation of shaft.
Symmetry 12 00361 g008aSymmetry 12 00361 g008b
Figure 9. RSM between X 6 X 9 and F T e m p .
Figure 9. RSM between X 6 X 9 and F T e m p .
Symmetry 12 00361 g009
Figure 10. Flowchart of artificial bee colony’s application on TBC optimization.
Figure 10. Flowchart of artificial bee colony’s application on TBC optimization.
Symmetry 12 00361 g010
Figure 11. Optimization process.
Figure 11. Optimization process.
Symmetry 12 00361 g011
Figure 12. The optimized variables.
Figure 12. The optimized variables.
Symmetry 12 00361 g012
Figure 13. Thermal characteristics of spindle (a) temperature field; (b) thermal deformation of shaft.
Figure 13. Thermal characteristics of spindle (a) temperature field; (b) thermal deformation of shaft.
Symmetry 12 00361 g013
Table 1. Parameters of the spindle bearings.
Table 1. Parameters of the spindle bearings.
Bearing Type7008C7009C
C s (kN)15.919.3
d m (mm)5460
Preload force (N)290340
Table 2. TBCs of spindle.
Table 2. TBCs of spindle.
Spindle Speed (rpm)Heat Transfer Coefficients (W/m 2 C) H 08 (W) H 09 (W)
A f 1 A f 2 A f 3 A f 4 A f 5 A n 1 A n 2
200042.243.344.032.840.69.79.76.89.2
400066.468.369.752.164.59.79.721.329.2
600085.489.491.468.384.59.79.741.857.2
8000102.3108.3110.782.7102.49.79.767.592.4
10,000123.3126.6128.495.9118.89.79.797.7133.8
12,000139.2143145.2108.3134.29.79.7132.4181.3
14,000154.4158.5160.8120148.79.79.7171.1234.4
16,000168.7173.3175.7131.3162.69.79.7213.6292.7
18,000182.5187.4190142175.89.79.7259.9356.1
20,000195.8201203.9152.3188.69.79.7309.7424.4
Table 3. TCRs of spindle joints.
Table 3. TCRs of spindle joints.
Joint SurfaceThermal Resistance (m 2 C/W)
Bearing inner ring/shaft 1.46 × 10 4
Bearing outer ring/housing 5.73 × 10 4
Rolling elements/raceways 6.46 × 10 4
Other contact surfaces 8.46 × 10 4
Table 4. Material parameters of spindle unit.
Table 4. Material parameters of spindle unit.
Spindle Structure (45)Bearing (GCr15)Bearing (Coolants (air))
Density (kg/m 3 )780080001.1796
Thermal conductivity (W·m 1 · K 1 )60.540.11.0069
Young’s Modulus (GPa)210200
Poisson’s ratio0.280.3
Linear expansion coefficient (K 1 ) 1.2 × 10 5 1.3 × 10 5
Table 5. Comparison of simulated and experimental temperatures with initial empirica TBCs.
Table 5. Comparison of simulated and experimental temperatures with initial empirica TBCs.
Temperature PointsExperimental ( C)Simulated ( C)Error ( C)Error (%)
T127.462.535.1128.1
Table 6. Variable parameters and their levels.
Table 6. Variable parameters and their levels.
Number−1.732−1011.732
x 1 33.455.485.4115.4137.4
x 2 37.459.489.4119.4141.4
x 3 22.151.491.4131.4160.7
x 4 16.338.368.398.3120.3
x 5 15.244.584.5124.5153.8
x 6 14.79.714.718.4
x 7 14.79.714.718.4
x 8 7.221.841.861.876.4
x 9 5.227.257.287.2109.2
Table 7. CCD method of the spindle system.
Table 7. CCD method of the spindle system.
NumberDesign Variable (Normalized) f Temp ( x )
x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9
1111−1111−1136.2
2−11−1−1111−1126.9
···········
155−11−1−1−1111138.6
1561−111−1−1−11−122.2
Table 8. Comparison of simulated and experimental temperatures with optimized TBCs.
Table 8. Comparison of simulated and experimental temperatures with optimized TBCs.
Temperature PointsExperimental ( C)Simulated ( C)Error ( C)Error (%)
T127.428.71.34.7

Share and Cite

MDPI and ACS Style

Zhang, L.; Xuan, J.; Shi, T. Obtaining More Accurate Thermal Boundary Conditions of Machine Tool Spindle Using Response Surface Model Hybrid Artificial Bee Colony Algorithm. Symmetry 2020, 12, 361. https://doi.org/10.3390/sym12030361

AMA Style

Zhang L, Xuan J, Shi T. Obtaining More Accurate Thermal Boundary Conditions of Machine Tool Spindle Using Response Surface Model Hybrid Artificial Bee Colony Algorithm. Symmetry. 2020; 12(3):361. https://doi.org/10.3390/sym12030361

Chicago/Turabian Style

Zhang, Leilei, Jianping Xuan, and Tielin Shi. 2020. "Obtaining More Accurate Thermal Boundary Conditions of Machine Tool Spindle Using Response Surface Model Hybrid Artificial Bee Colony Algorithm" Symmetry 12, no. 3: 361. https://doi.org/10.3390/sym12030361

APA Style

Zhang, L., Xuan, J., & Shi, T. (2020). Obtaining More Accurate Thermal Boundary Conditions of Machine Tool Spindle Using Response Surface Model Hybrid Artificial Bee Colony Algorithm. Symmetry, 12(3), 361. https://doi.org/10.3390/sym12030361

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop