Thermal Characteristics of Spindle System Based on the Comprehensive E ﬀ ect of Multiple Nonlinear Time-Varying Factors

: The thermal characteristics of the spindle system for CNC machine tools are in ﬂ uenced by multiple factors which are nonlinear and time-varying. In this paper, a nonlinear time-varying thermal characteristics solving model for the spindle system was established based on the numerical solution method. Through theoretical deduction and data ﬁ tt ing, mathematical models of nonlinear time-varying factors including the friction torque generated by lubricants, convective heat transfer coe ﬃ cient, and coolant and ambient temperature are constructed. The temperature and displacement of the spindle system at each time step are solved by considering the comprehensive e ﬀ ect of multiple nonlinear time-varying factors. And the actual temperature and axial deformation data of the spindle system are obtained through thermal characteristics detection experiments. By comparing solution results a ﬀ ected by multiple nonlinear time-varying factors and by non time-varying factors with experimental data, it can be concluded that the nonlinear time-varying thermal characteristics model has advantages in re ﬂ ecting the trend of numerical changes and the accuracy of result solving over a method considering non time-varying factors and the solution values of temperature a ﬀ ected by multiple nonlinear time-varying factors are almost consistent with detection values and the relative errors are all within ±3%. The relative error of axial deformation between the value solved by the model and the detection value is close to − 1%. This conclusion demonstrates the rationality and accuracy of the thermal characteristics solving model and the construction of non-linear time-varying factors. This study is of great signi ﬁ cance for exploring the thermal characteristics of the spindle system and improving CNC machine tool performance in depth.


Introduction
The technical level of CNC machine tools directly affects the development of the manufacturing industry.The accuracy and stability of CNC machine tools are crucial for product quality [1].The spindle system, whose accuracy is one of the core indexes to measure the performance of a CNC machine tool, is the end effector during the machining process of the workpiece [2].As a complex assembly composed of spindle, bearing supports, lubrication, cooling, and other component structures, its machining error is the result of the comprehensive effect from multiple factors related to the components and their interrelationships [3].Heat-related factors account for 40% to 70% of all factors which Citation: Lin, X.; Deng, X.; Zheng, J.; Academic Editor: Iztok Golobič cause accuracy loss of the spindle system [4,5].Starting from the factors affecting the thermal characteristics to study their comprehensive effect on the temperature and deformation of the spindle system is an effective method which can improve the performance of CNC machine tools.
In recent years, there have been numerous studies on the thermal characteristics and their influencing factors on the spindle system.Lee et al. [6] researched the changes in temperature and deformation of the spindle system due to the influence of spindle speed by the finite element method and thermal analysis method.The relationship between thermal deformation and vibration was analyzed and modeled, which was verified experimentally.Raja et al. [7] developed a coupled fluid-thermal model for spindles.The model can simulate the fluid-structural conjugate heat transfer of the spindle system.They measured the minimum deviation of the model through experiments as 7.6%.Marez et al. [8] developed an approach based on TFs which is used for effective thermal error modeling for machine tools.This approach provides insight into the share of each source in the total machine thermal error through a combination of linear parametric models.Kaftan and Wegener et al. [9] studied the adverse effects of internal and external heat sources on the complex non-symmetric structure of machine tools.They developed a novel method based on artistic intelligence that compensates for thermal errors associated with hidden boundary condition changes to address the drawbacks of traditional methods.Tan et al. [10] studied the thermal characteristics of spindle bearings under preloading and achieved nonlinear prediction of thermal characteristics.The nonlinear changes in contact angle, contact force, preload, and stiffness were analyzed by considering the effects of contact thermal resistance, bearing parameters, lubricant viscosity, and time-varying temperature on the heat source.The proposed method significantly reduced the computational workload.Liu and Ma et al. [11] proposed a closed-loop iterative modeling method.The heat generation of bearings and built-in motors, convection coefficient of bearing joints, and thermal contact resistance were calculated by modifying the heat source and thermal boundary conditions in each calculation step.By considering the comprehensive effects of lubricant viscosity changes and bearing thermal preload, the bearing heat generation was corrected.Xiang et al. [12] proposed a data-driven prediction approach which can establish a dynamic linear model of spindle thermal error.The serious contradictions in traditional prediction methods have been resolved by predicting current thermal errors through historical temperature data without physical mechanism information.Wu et al. [13] established a thermo-mechanical coupling analysis model for the spindle bearing system of machine tools.Through this model, they obtained the relationship between bearing preload and frictional heat generation, as well as between coolant and system thermal balance.Brecher, Wenkler, and Ihlenfeldt et al. [14,15] also explored the thermal related influencing factors of CNC machine tools in the research of thermal error modeling.
After the CNC machine tool is started, its state and surrounding environment exhibit nonlinear time-varying characteristics.The comprehensive effect of multiple changing factors is an important cause for the poor thermal characteristics and accuracy loss of CNC machine tools.At present, according to the research of many scholars, the research on factors such as heat sources, bearing loads, and coolant is relatively in-depth.The nonlinear and time-varying effects of these factors on the spindle have also been studied separately.However, there are few reports on solving the thermal characteristics of the spindle system under the comprehensive effects of multiple nonlinear time-varying factors such as heat sources, cooling conditions, surrounding environment, and heat transfer capability.Therefore, this paper establishes a numerical solving model for the nonlinear time-varying thermal characteristics of the spindle system.Based on the time-varying variables of the model, mathematical models of nonlinear time-varying factors such as friction torque generated by lubricant, convective heat transfer coefficient, and coolant and ambient temperature are constructed through theoretical derivation and data fitting.The time-varying temperature of the spindle system is solved taking into account multiple factors.Based on the time-varying value of the temperature, the deformation of the spindle at each moment is calculated.This paper provides a theoretical basis for the thermal characteristics study of spindle systems under complex working conditions.

Establishment of the Nonlinear Time-Varying Thermal Characteristics Model
During the operation of CNC machine tools, the heat transfer process of mechanical components in the spindle system follows the control equation shown in Equation (1) [16].
where ρ is the material density of the component, CP is specific heat, T is temperature, t is time, k is thermal conductivity, and S is source terms.
Based on the finite volume method, this study will use hybrid grids to partition the control volume of the physical model for the spindle system.For each control volume, the central difference scheme and time implicit format are introduced.The pressure-based solver of Fluent 2022 R2 software will be used for solving.The least square cell-based and second-order upwind are adopted for spatial discretization.Therefore, within the allowable range of discretization error and numerically controllable, Equation ( 1) can be expanded as where i t P T and 1 i t P T  are the average temperature of the local control volume P at the ith and i−1th time step, VP is the volume of P, M is the number of adjacent control volumes to P,  is the average temperature of control volume N adjacent to P at the i-th time step, kPN is the harmonic mean of thermal conductivity between control volume P and N, is the distance between the center of P and N, APN is the heat transfer surface area between P and N, s is the intensity of the heat source, i t  is the i-th time step size, and The heat flux on both sides of the interface between the spindle system and the external fluid zone in contact with the system is equal under the state of fluid-solid coupling heat transfer [17].According to Fourier's law and Newton's flow equation, Equation (3) can be obtained.
where L is the distance along the vector normal to the interface, h is the convective heat transfer coefficient, and T  is the temperature difference between the fluid and wall surface.
Assuming the temperature of fluid in contact with the spindle system is TF, Equation (4) can be established within the allowable range of discretization error and numerically controllable.
where kp is the thermal conductivity of P, i t W T is the temperature of the wall surface at the i-th time step, is the distance between the center of P and the wall surface, and .F h is the convective heat transfer coefficient of the fluid.
In components, there is obviously Equation (5), as follows.
It is assumed that the thermal conductivity of the spindle system material is constant.The nonlinear time-varying thermal characteristics of spindle systems should be considered from heat generation and heat dissipation.The heat generation is reflect in the intensity of the heat sources whose values at the i-th time step can be expressed as i t s in the control equation.And the heat dissipation is usually related to the external ambient conditions and the capabilities to transfer heat outward, which is mainly reflected in the convective heat transfer coefficient and the temperature of the fluid in contact with the system.Their values at the i-th time step can be expressed as i t F h and i t F T in the control equation.By combining Equations ( 2) and ( 6), the nonlinear time-varying thermal characteristics model shown in Equation (7) where J is the number of convective heat transfer surfaces of P, and APF is the convective heat transfer surface area of control volume P.
The heat transfer equation shown in Equation ( 7) is the model for the nonlinear timevarying thermal characteristics of the spindle system.According to the model, i t s , i t F h , and i t F T are the variables which can affect the value of the temperature of control volume P and its adjacent control volume N.They are usually nonlinear and time-varying.As the definite solution conditions for the temperature of the spindle system, the models of nonlinear and time-varying factors related to heat source intensity, convective heat transfer coefficient, and the temperature of fluid in contact with the system should be constructed and substituted into the solution.

Time-Varying Factors of System Heat Source
For non-heating component in the spindle system, ( ) 0 s t  and for components that generate heat, s(t) is obtained according to Equation (8).
where Q(t) is the time-varying function of heat flow from the heating component, and V is the volume of the heat source component.The heating of the motor and bearings is the main cause of temperature rise in the spindle system [18,19].Due to the installation of the motor outside the spindle box, less heat is transferred to the spindle components through the mounting surface of the box.The bearings heating is due to the viscosity of lubricants and the torque generated by bearing loads.The heat flow Q can be expressed by Equation ( 9) [11].
where M0 is the friction torque generated by lubricant viscosity, M1 is the friction torque generated by the load on the bearing, and n is the spindle speed.
During the rotation of bearing, the viscosity of the lubricant changes with the increase in temperature.In other words, if oil lubrication is used, the viscosity voil is a time-varying factor related to temperature Toil, as shown in Equation (10).
The spindle rotates at a speed of n.M0 generated by lubricant viscosity exhibits timevarying and nonlinear characteristics.When , Equation ( 11) can be obtained according to references [20].
where 0 ( ) M t is the time-varying function of the friction torque generated by lubricant viscosity, f0 is the coefficient determined by bearing structure and lubrication, and dm is the middle diameter of bearing.
From Equations ( 9) and ( 11), it can be seen that the heat generated by the lubricant will change the viscosity of the lubricant, whose change will also affect the heat generation in turn.The two mutually constrain each other through the intermediate variable Toil (t).The value of Toil (t) is obtained by replacing i t F T and iteratively calculating based on the nonlinear time-varying thermal characteristics model shown in Equation (7).The coupling relationship of voil, Toil, and M0 makes them exhibit time-varying and nonlinear characteristics.

Time-Varying Factors of Heat Dissipation Outside the System
The heat dissipation outside the spindle system mainly involves convection heat transfer with air forcefully and naturally.The air characteristics affected by ambient temperature are nonlinear and time-varying due to the long time from startup to thermal equilibrium of the spindle system, which causes both the convective heat transfer coefficient and the temperature difference inside and outside the system to change simultaneously.In other words, the construction of time-varying factor models for the external heat dissipation of the system should include ambient temperature and convective heat transfer coefficient.
Forced convection heat transfer mainly occurs on the rotating surface of the spindle system.When the Reynolds number Reair and the Prandtl number Prair satisfy Reair < 4.3 × 10 5 and 0.7 < Prair < 670, according to reference [21], the forced convection heat transfer coefficient f ( ) h t can be calculated using Equation (12).
where Nuair(t) is the Nusselt number that varies over time, λair(t) is the time-varying function of air thermal conductivity, and ds is the equivalent diameter of rotating surfaces.According to references [22,23], the forced convection heat transfer coefficient f ( ) can be derived by expressions of the Nusselt number as follows: where ρair(t), µair(t), and Prair(t) are the time-varying functions of density, kinematic viscosity, and the Prandtl number for air, respectively.For spindle segments with different diameters, the equivalent diameter ds can be calculated by Equation ( 14) [24]: where l is the total length of spindle, and di and li are the diameter and the length of the ith spindle segment.
For natural convection heat transfer surfaces, the time-varying function of the natural convection heat transfer coefficient hn(t) can be derived according to Equation ( 15) [20].
where g is the standard gravitational acceleration, βair(t) is the time-varying function of air thermal expansion coefficient, D is the feature size, and ΔT (t) is the temperature difference between ambient air and the heat transfer wall, which is obtained using the following equation: where is the time-varying function of ambient temperature, and   W T t is the time-varying function of the heat transfer wall surface and its value determined by Equation (5).

Time-Varying Factors of Internal Cooling Inside the System
The convective heat transfer of the spindle system not only occurs externally with the air but also internally with the coolant such as cooling oil.Similar to the time-varying factors of external heat dissipation, the time-varying characteristics of internal factors such as convective heat transfer coefficient hoil and cooling oil temperature Toil need to be considered.
Due to the spiral-shaped slender channel in the cooling jacket of spindle system, the flow of cooling oil is mostly laminar.This paper assumes that the oil density remains constant.The nonlinear time-varying function of the convective heat transfer coefficient hoil of cooling oil can be constructed as Equation ( 17) based on reference [25].
where oil ( ) P C t is the time-varying function of oil-specific heat, ρoil is the density of oil, dc is the hydraulic diameter of cooling channel, and λoil(t) is the time-varying function of oil thermal conductivity.lc is the channel length.
The nonlinear time-varying function oil ( ) T t of the cooling oil temperature is obtained by fitting the data detected by the experiment.

Experimental Platform Construction and Data Detection
This paper conducts experimental research on the G1160 CNC machine tool.This CNC machine tool adopts No. 32 cooling lubricating oil, whose cooling channel is spiral shaped.The motor drives the spindle through a synchronous belt.The spindle is equipped with 7014 C angular contact ball bearings with DT installation method.The internal structure of the G1160 CNC machine tool spindle system is shown in Figure 1.The spindle runs at a speed of 6000 r/min.The intelligent thermal characteristic detection and compensation instrument for the CNC machine tool spindles is adopted to detect data of the spindle system.The working temperature of this instrument is 4-50 °C, and the sampling frequency is 5 s per time.The system temperature and ambient temperature are detected through a PT100 temperature sensor manufactured by Heraeus in Germany.It has a 16 bit A/D converter and a measurement accuracy of 0.4%.The cooling oil temperature Toil of the CNC machine tool is detected and read using the RCO-15PTS oil cooler digital temperature display instrument which is equipped with the CNC machine tool and produced by Ruike Refrigeration Plant Co., Ltd, Dongguan, China.To detect the spindle deformation, a BT40 inspection rod with a diameter of 40 mm and a length of 300 mm is clamped onto the spindle.The axial deformation of the inspection rod is detected by the intelligent thermal characteristic detection and compensation instrument for CNC machine tool spindles, using the CPL230 sensor produced by LION PRECISION in America (Minneapolis, MN, USA).It includes a high-precision capacitive displacement probe and a compact multi-channel driver.As a non-contact sensor, its detection range is 125-375µm and the working temperature is 4-50 °C.The RMS resolution of the sensor is 18 nm.The construction of an experimental platform for the data detection of the G1160 CNC machine tool spindle system and the instruments are shown in Figure 2. Taking into account the structure of the spindle system, as well as the heat source and heat dissipation, the temperature sensors numbered T1-T8 were arranged according to Figure 3.The temperature values detected by the T1-T8 sensor are Ti, i = 1,2… The temperature data of the spindle end, spindle bearings, flange, spindle body, both sides of the spindle box, and spindle box support plate are detected separately.A temperature sensor numbered T9 is placed in the air to detect the ambient temperature Tair simultaneously.The sensor used to detect the axial displacement Zs of the inspection rod is numbered Z.The settings for sensor positions are shown in Figure 3 and the description of the sensor positions is provided in Table 1.

The Solution of Thermal Characteristics Based on Multiple Time-Varying Factors
It can be seen from Equations ( 10) and ( 11) that the friction torque generated by the lubricant is affected by the characteristics of kinematic viscosity.According to reference [26], the nonlinear time-varying characteristics of the No. 32 oil can be described using Equation (18).
The torque M1 generated by the load on the bearing can be calculated according to Equation ( 19) [27].
where f1 is a coefficient determined by the structure and load of the bearing, and P1 is the calculated load for the bearing torque.For angular contact ball bearings, f1 and P1 can be calculated according to Equation (20) [28].
where F0 and C0 are the equivalent static load and rated static load of the bearing, and Fr and Fa are the radial and axial load of the bearing.F0 is calculated according to Equation (21).
where X0 and Y0 are equivalent static load coefficients.If P1 obtained from Equation ( 20) is less than Fr, then P1 = Fr.The load on the bearing of the spindle installing the inspection rod mainly includes the bearing preload Fp, the synchronous belt compression spindle force FQ, and the gravity G of the rotating component.Fp provided by the enterprise is 700 N. G can be obtained by summing up the masses of each component.FQ is calculated according to Equation ( 22) [29].
where KA is the operating condition coefficient, Pm is the nominal power, d1 is the diameter of the small pulley, and n1 is the small pulley speed.According to Figure 4, the ambient temperature rise range is only 2.9 °C during the detection period.The air property parameters in the convective heat transfer coefficient are all monotonic functions with air temperature as the independent variable.Equation ( 23) is constructed to determine the sensitivity of multiple air property parameters to the convective heat transfer coefficient.Low-sensitivity property parameters can be set as constants to simplify calculations.( ) where Sg is the sensitivity, h  is the extreme difference of the convective heat transfer coefficient affected by a certain property parameter within the detection temperature range, h is the time averaged value of the convective heat transfer coefficient during the experimental period, g(Tmax) and g(Tmin) are the peak values of a certain property parameter at the maximum and minimum detection values of experimental temperature, respectively, and te is the end time of the experimental period.According to reference [30], the natural convection heat transfer coefficient is much smaller than the forced convection heat transfer coefficient.Considering the actual radiation heat dissipation, the equivalent natural convection heat transfer coefficient is set to 9.7 W/(m 2 •°C).
Tair(t) and Toil(t) are the important nonlinear time-varying factors affecting the thermal characteristics of the spindle system.The fitting degree of their mathematical models for experimental data can affect the accuracy of the calculation results.The fitting degree is evaluated using root mean square error (RMSE) and R-square.For this study, RMSE was calculated according to Equation ( 24) [31].
where n is the number of data samples, and T(ti) and ( ) are the detected value of the experiment and the fitting value of the mathematical model at the i-th time point.
The calculation of R-square is shown as Equation ( 25) [32].
where ( ) i T t is the average value of the detected temperature.The rational function fitting method is used to fit the detected data of Tair and Toil.Both the Numerator degree and Denominator degree are set to 1.The fitted mathematical models are shown in Equation (26) and Equation ( 27 [33], the values of each parameter at different temperatures can be obtained by linear interpolation.According to Equation ( 23), the sensitivity solution results of air property parameters are shown in Table 3.The small variation in air property parameters, especially Prair, due to the small range of ambient temperature rise detected in this study can be seen from the table.The SST k-Omega turbulence model was used for solving in this study.The solution process adopts the SIMPLE strategy.In transient calculations, the values of time-varying factors and their parameters at each time step will be solved and updated by the solver in two ways: the time-varying factors are directly solved by substituting time steps into their time functions or the time function of property parameters; the factor values for the next time step are solved and updated based on the temperature values obtained from the current time step.It is worth noting that the different materials of different components and the different thermal characteristics of the same component at different positions can affect the waste heat dissipation after the last operation of the CNC machine tool.Therefore, there is a certain gradient in the initial temperature field of the spindle system, for which the initial detected temperature varies at different detection positions.Due to the focus of this study on the process of multiple factors affecting thermal characteristics, before solving, the spindle system is initialized with 7.4 °C, which is the lowest initial temperature value of each sensor detection position.

Temperature Analysis of Spindle System
The temperature distribution on the surface and inside of the spindle system under the comprehensive effect of multiple nonlinear time-varying factors at 7000 s is shown in Figure 5.It can be seen that the high-temperature area is mainly distributed at the positions of the upper and lower bearings.Due to their influence, the temperatures of the spindle body, pulley, cooling jacket, and the root of the inspection rod are also relatively high.The temperature peak of the entire system is 30.8 °C, which is mainly near the installation position of the bearings.The high temperature value extends from this position and gradually decreases.The spindle body maintains around 20 °C because of its position in the middle of the upper and lower bearings.Due to the installation surface near the lower bearing group, the temperature at the root of the inspection rod is about 22 °C or above.The temperature of partial areas on the front, left, and right sides of the box is increased through the thermal conductivity of the connecting rib plates.The temperature of the points numbered T1  ~T8  , which correspond one-to-one with the positions of sensors numbered T1~T8 on the spindle system, is solved at each time point.The temperature values are i T  , I = 1,2…The curves of the experimental detection value and the solution value at each point are plotted.The relative error i  is calculated according to Equation ( 28), which is drawn to evaluate the accuracy of the model for solving the nonlinear time-varying thermal characteristics of the spindle system considering the comprehensive effects of internal and external factors.
From the comparison shown in Figure 6, it can be seen that the experimental detection values and the solution values gradually increase over time.Their overall trend remains consistent.The relative error gradually decreases at each moment.Except for the several initial moments with relatively large relative errors, as time progresses, the relative error values of most points are within ±5% after 3000 s.At the final moment, the relative error of each point is below ±3%, and more than 60% of the points even reach within ±1%.Relative error/% By analyzing the temperature curves of each point, there may be four reasons for the existence of relative error.Firstly, there are significant differences in the initial temperature values of each detection point which, to some extent, affect the temperature changes in the initial stage.Also, due to the small initial temperature values, the relative error is relatively large.However, it can be seen that the difference in initial temperature has a relatively short impact time.The overall fit between the solution value and the detection value of the temperature is relatively high.Secondly, Tair detected by the sensor T9 for detecting ambient temperature has difficultly describing the actual complex temperature field around the spindle system.Due to the interference and gradient at different positions, the temperature inside and outside the system is different, which results in a difference in convective heat transfer coefficients.Thirdly, the oil temperature detected by the oil cooler digital temperature display instrument experiences a temperature drop during the entire pipeline transportation process due to unstable heat dissipation at different moments, which to some extent affects the solution results.Finally, there are complex variations in the characteristics of sealed components and contact thermal resistance between multiple components, as well as other minor time-varying factors that have not been considered.These changes affect the system temperature to varying degrees as the system temperature rises.Especially during the stage of rapid temperature rise, the actual temperature values are higher than the solution values and the relative error curve exhibits certain fluctuations.However, overall, the model for solving the nonlinear time-varying thermal characteristics of the spindle system has high accuracy.

Axial Displacement of the Spindle
The temperature field solution data at different time points are loaded onto the spindle system structure as thermal load conditions.In addition, the bearing preload and BT40 spindle tension, with the values provided by the enterprise, are loaded as the mechanical conditions.The transient solution for the deformation of the spindle system has been completed.And the axial displacement s Z  of the inspection rod is simultaneously calculated.The total deformation of the spindle system at the final moment is shown in Figure 7.The relative error δs is calculated according to Equation (29).The comparison between the solution value s Z  and the experimental detection value s Z of axial displacement at the end of the inspection rod, as well as the relative error s  curve between the two, are shown in Figure 8.The installation surfaces of the spindle system are located on the back of the box.At the final moment, the total deformation of the spindle system shown in Figure 7 continuously accumulates as the distance from the installation surface increases.The maximum deformation occurs at the end of the inspection rod, which is 77.6 µm.
It can be seen from Figure 8 that both the detection value and the solution value of axial displacement at the end of the inspection rod gradually increase, and the trend of change is almost consistent throughout the entire solution period.Except for the initial few time points where the relative errors are large because of the small displacement values, the relative error of displacement after about 10 min is all within 10% and gradually approaching −1%.Due to the relative error in solving the temperature of the spindle system, the solution value of axial displacement at the end of the inspection rod is also smaller than the actual displacement.But over time, the values of the two gradually approach.This also fully reflects that the establishment of the thermal characteristics solving model and the construction of multiple nonlinear time-varying factors for the spindle system has high rationality and accuracy.

Comparison of Thermal Characteristics Solution with Non Time-Varying Factors
The traditional methods often consider changes in the heat source conditions of the spindle system but rarely comprehensively consider the external environment, internal cooling conditions, and the heat transfer ability caused by their temperature during the actual operation of the CNC machine tool.Therefore, the ambient temperature and cooling oil temperature of the spindle system are set as non time-varying values in this study.In order to compare the thermal characteristics value at the last moment, separately, the ambient temperature and cooling oil temperature are set to the detected 12.3 °C and 16.4 °C, respectively.The thermal characteristics of each point in the spindle system are solved.The temperature comparison of the detection values and solution values affected by timevarying factors and by non time-varying factors at the final moment, as well as their relative errors, are shown in Figure 9.It can be seen from Figure 9 that except for T5, the relative errors of the solution values affected by time-varying factors are closer to 0 than the solution values affected by non time-varying factors.This indicates that the proposed solution method considering the comprehensive effect of time-varying factors has more advantages in solving the accuracy of temperature in the spindle system.
The axial displacement curves of the detection values and solution values affected by time-varying factors and by non time-varying factors, as well as their relative error curves are shown in Figure 10.It can be seen that due to the constant ambient and cooling conditions, the axial displacement curve affected by non time-varying factors takes a relatively short time to approach a steady state, about 2500 s.However, in reality, the steady state is not achieved even at 7000 s.The change trend of axial displacement is not as close to the actual axial displacement rise, as solved by considering the effect of time-varying factors.From the perspective of calculation accuracy, at 7000 s, the relative error of the solution value calculated by non time-varying factors is −5.13%The value obtained by time-varying factors is −1.24%.The above analysis indicates that the nonlinear time-varying thermal characteristics model has advantages in reflecting the trend of numerical changes and the accuracy of result solving.

Conclusions
In this paper, a model for solving the nonlinear time-varying thermal characteristics of the spindle system is established based on numerical solution methods.And the nonlinear time-varying models of internal and external factors are constructed through theoretical derivation and data fitting.The temperature and deformation of the spindle system are solved taking into account the comprehensive effect of multiple factors.By comparing solution results affected by multiple nonlinear time-varying factors and by non time-varying factors with experimental data, the following conclusion can be drawn: the nonlinear time-varying thermal characteristics model has advantages in the solution of temperature and deformation for the spindle system over the method considering non time-varying factors; the established model for solving the nonlinear time-varying thermal characteristics of the spindle system has high rationality in solving the temperature field values of the spindle system; the construction method of mathematical models for nonlinear timevarying factors including friction torque generated by the lubricant, convective heat transfer coefficient, and coolant and ambient temperature is correct; and the solution results for the temperature and axial displacement of the spindle system by substituting multiple nonlinear time-varying factors into the nonlinear time-varying thermal characteristics model are almost consistent with the actual values.The relative errors of temperature and axial displacement are within ±3% and close to -1%, respectively.The thermal characteristics study of the spindle system based on the comprehensive effect of multiple nonlinear time-varying factors provides a theoretical basis for solving the thermal characteristics of spindle systems under complex working conditions.

Figure 1 .
Figure 1.Internal structure of the spindle system.

Table 1 .
Description of sensor positions.of Tair and Toil used for fitting at each time point are shown in Figure 4.

Figure 4 .
Figure 4. Detected ambient temperature and cooling oil temperature curves.

Figure 5 .
Figure 5. Temperature distribution on the surface and inside of the spindle system.

Figure 6 .
Figure 6.Temperature and relative error curves of each point.

Figure 7 .
Figure 7.The total deformation of the spindle system.

Figure 8 .
Figure 8.The axial displacement and the relative error curve.

Figure 9 .
Figure 9. Temperature comparison of 8 points at the final moment.

Figure 10 .
Figure 10.Comparison of axial displacement and relative error curves.

Table 2 ,
), respectively.arethe fitted mathematical models of ambient temperature and cooling oil temperature.The RMSE and R-square of a ir ( ) which reflects a high degree of fit in a ir( )

Table 2 .
RMSE and R-square of a ir ( )

Table 3 .
The sensitivity of property parameters.