A dynamic thermal-mechanical model of the spindle-bearing system

This paper presents a dynamic thermal-mechanical model to investigate the thermal characteristics in a spindle-bearing system. In this model, transient thermal analysis, static structure analysis and calculation of the boundary conditions are conducted as a solution loop. The transient boundary conditions, such as bearing stiffness, bearing heat generation and thermal contact conductance, are calculated with the appropriate formulas and solution methodology. The thermal feedbacks, which are seldom considered in the previous studies, are calculated in details. In order to validate the prediction accuracy, thermal equilibrium experiment is conducted on a test rig of spindle-bearing system. The predictions of the proposed model, such as bearing preload, temperature and thermal displacements, are in close agreement with the experiment results. The comparisons of the proposed model with two traditional simulations show that, the thermal feedbacks on the boundary conditions and thermal contact conductance are of great importance to the real-time estimation of the thermal characteristics. The proposed model provides a practical method to improve the prediction accuracy. It could also be generalized in other mechanical systems to investigate the dynamic thermal characteristics.


Introduction
The development of the modern machine tools proceeds towards high precision and speed, and highlights the importance of thermal effects on mechanical systems.The thermal error accounts for almost 40-70 % of the error sources in a machine tool (Han et al., 2012) The spindle-bearing system is a core component of the machine tool.The performance of the spindle-bearing system is of great significance to the machine tool.Due to coupled interactions between the thermal and mechanical behaviors, it is extremely difficult to forecast the thermal characteristics of the spindle-bearing system accurately.As a result, only a few researches in thermalmechanical analysis are available for designing the spindlebearing system.
There are generally two schools in the thermal-mechanical analysis of the spindle-bearing system.One school analyzes the thermal effects with the statistical models, such as regres-sion (Yang, 2003;Wu and Kung, 2006) and neural network (Mize and Ziegert, 2000;Guo et al., 2010).The statistical model can be considered as a "black box".The inputs of the model are operating parameters or temperatures, and the output is thermal deformation.Based on fitting analysis of the experiment results, the model reflects nonlinear relationship of thermal effects.Owing to the high efficiency of calculation, the statistical model is usually used in the thermal error compensation.However, the models could not be used to identify the thermal characteristics without the experiment results.
The other school focuses on the theoretical methods to identify the thermal characteristics.Finite element model (FEM) (Ma et al., 2015;Xu et al., 2007) is the most common computing methods.With the help of softwares such as AN-SYS and ABAQUS, the temperature and the thermal deformation are calculated upon the basic laws of heat transfer and stress-strain.Bossmanns andTu (1999, 2001) presented one Stress-strain field under thermal effects Step 1 Step 2 Step 3 Figure 1.Thermal-mechanical analysis.
of the pioneering works on the thermal-mechanical model of spindle-bearing system.The friction efficiency of the built-in motor was utilized to predict the heat generation.The estimation of the heat generation in the bearing was based on the research of Harris (1991).Jin et al. (2012) improved the calculation model of heat generation in the bearing by considering various operating conditions.Due to different bearing configurations of the spindle-bearing system, the uneven axial expansions of the housing and shaft was calculated in the research of Li and Shin (2004).Jiang and Mao (2010) predicted the temperatures of a motorized spindle based on FEM.At high speed range, the spindle preload was determined by the constraint of temperature rise of bearings.At low speed range, the spindle preload was resolved by the fatigue life of bearings.Lee et al. (2015) analyzed the relationship between thermal deformation and vibration during rotation of the spindle.Takabi and Khonsari (2013) developed a lumped-mass model to analysis the transient temperatures in the spindle-bearing system.
For the theoretical methods based on FEM, there are two limitations in the previous studies.Firstly, the thermal contact conductance (TCC) of the joint surface is often neglected (Li and Shin, 2004;Lee et al., 2015;Jiang and Mao, 2010) or modeled as fractal method (Ma et al., 2015;Xu et al., 2007).The TCC will induce the decrease of temperature between two rough surfaces.It is of great importance to the thermal-mechanical analysis.Some researchers modeled the TCC with fractal method.However, the fractal parameters could not be estimated as the functions related to the traditional roughness parameters such as rootmean-square roughness.The fractal parameters of fractal dimension and fractal roughness are calculated by the structure function method after precise measurement (Majumdar and Tien, 1990) The difficulty in estimation of fractal parameters restricts its use in practical application.Secondly, insufficient attentions have been paid on the influence of thermal feedbacks on the boundary conditions.In general, the thermal-mechanical analysis is performed as the following steps presented in Fig. 1.The step 1 is the calculation of the temperature field with the boundary conditions.The step 2 is the procedure from the temperature field to the stressstrain field caused by thermal expansion.The step 3 is the calculation of thermal feedbacks on the boundary conditions caused by the variations of temperature field and stress-strain fields.Though there are a large number of researches in the area of thermal-mechanical analysis, most of them only involved the first two steps.For the thermal feedbacks, the previous studies (Li and Shin, 2004;Takabi and Khonsari, 2013) mainly considered the uneven axial expansions of the housing and shaft.However, the thermal-induced variations of the variables including bearing dimensions, bearing stiffness and thermal contact conductance are not considered.
This paper aims to propose a dynamic thermal-mechanical model of the spindle-bearing system to investigate the transient thermal characteristics, and is organized as follows: In Sect.2, the solution procedure of the proposed model is first introduced.The proposed model is established based on the loop calling between the softwares of MATLAB and Ansys Parametric Design Language (APDL).In addition, appropriate formulas and solution methodology are conducted to calculate the time-varying boundary conditions.The stiffness and heat generation of ball bearings are calculated based on the quasi-static mechanical model of ball bearing.The Newton iteration is adopted to solve the mechanical model of the bearing.The TCC is predicted by the elastic-plastic contact model.The TCC varies with the contact pressure of the joint surface.The convection coefficient is calculated with the Nusselt number.Section 3 presents a test rig of spindlebearing system with measurable and adjustable preload.The bearing axial force, temperature and thermal displacement are measured to validate the prediction accuracy of the simulation.Section 4 compares the simulation results with the experiment results.The simulation results of the proposed model are in close agreement with the experiment results.Compared with two traditional simulation models, the proposed model achieves remarkable improvements in prediction accuracy.And finally the conclusion is given in Sect. 5.

The solution procedure
In order to investigate the transient thermal characteristics in spindle-bearing system, the solution procedure shown in Fig. 2 is applied.The solution loop contains thermal analysis, structure analysis and the calculations of the boundary conditions.The thermal analysis is used to calculate the transient temperature variation at each time increment.As the time in- crement is set as 10 s in this paper, the temperature field could not achieve equilibrium in such a short time.As a result, the transient thermal analysis is used in the solution procedure.For the structure analysis, the variation of stress-strain field induced by the thermal expansion is calculated.From a micro perspective, thermal expansion is caused by the interatomic spacing change of the material due to temperature change, namely the stress-strain field changes in real time according to the variation of the temperature field.So the set time increment is long enough for the stress-strain field to reach the static state.Static structure analysis, instead of transient structure analysis, is conducted to improve the solving efficiency.The software of MATLAB is used to calculate the transient boundary conditions under thermal effects.AN-SYS(APDL) and MATLAB are called circularly during the simulation.Section 2.2-2.4 presents the calculations models in MATLAB.

Geometric relations analysis
The mechanical model of ball bearing under thermal effects is established in this paper.The geometry relations in the ball bearing are shown in Fig. 3.The rotating coordinate system, (XY Z), is fixed in the space where Z-axis is coincident with the axis of the inner race.The local coordinate system, (tnb), is used to describe the motion of the rolling ball.Figure 3b presents the displacement of the j th rolling ball.As the mechanical model of ball bearing is calculated under outer raceway control theory, the outer raceway groove curvature center (Point A) is unchanged under static state and rotating.j is the distance between the j th rolling ball center and raceway groove curvature center under rotating.A 1j and A 2j are the distances between outer and inner raceway groove curvature centers in n-axis and b-axis respectively.According to the geometric relations, the dimensions of ij , oj , A 1j and A 2j in Fig. 3b could be calculated as Eqs.( 1) and ( 2) Where D is diameter of the rolling ball and f is curvature radius coefficient of raceway.δ j is contact deformation between race way and rolling ball.rb is thermal expansion of the rolling ball radius.
Where X , Y and Z represent the displacements of the inner race.θ X and θ Y are the rotational angles around the X-axis and Y -axis.or_n and ir_n are the dimensional changes of outer and inner raceways in n-axis of coordinate system (tnb) and could be achieved from the structural analysis of APDL.In general, or_n and ir_n will be affected by three factors including thermal expansion, contact force and centrifugal force.
Based on the Pythagorean Theorem in Fig. 3b, Eq. ( 4) could be achieved:

Velocity analysis
Ball velocities are required in force analysis and calculation of heat generation.The velocities of the bearing could be decomposed as Fig. 4. In Fig. 4a, U is the spinning axis of the rolling ball.The angle β and β are the ball rotation angles in the n − b and n − t plane.Though β varies with the rotation speed, the assumption of β = 0 still provide reasonable predictions of β and angular velocities (Jorgensen, 1996).For the ball bearing in spindle-bearing system, orbital angular velocity ω m and angular velocity of the spinning axis ω R could be calculated as Eqs.( 5) and ( 6) respectively:  Where ω is the angular velocity of inner race.
In Fig. 4b, ω si and ω so are the spin angular velocities of the rolling ball on the inner and outer race respectively, Eq. ( 8) could be achieved according to the velocity resolution:

Force and moment analysis
The state of the rolling ball will be affected by the contact forces Q j , centrifugal force F cj and the friction forces caused by the gyroscopic moment M gj as shown in Fig. 5.
Considering the force equilibrium of the rolling ball, Eq. ( 9) could be achieved.
Where λ ij and λ oj are the influence parameters.If the bearing is under the outer raceway control theory (Harris, 1991), λ ij = 0 and λ oj = 2.The forces and the moments could be calculated as: The force and moment balance equations of the inner race are expressed as Eq. ( 13): Where F x , F y , F z , M x and M y are the external forces and external moments imposed on the inner race.F z is composed of bearing preload and cutting force.The preload will be affected by the uneven axial thermal deformations of shaft and housing as shown in Fig. 6.
Based on the bearing configurations, the preload under thermal effects is calculated as: Where iz and oz are the thermal expansions of the spindle and the housing respectively.PL init is the initial bearing preload.K z denotes the axial stiffness of the bearings in pairs and calculated as:

Heat generation
The heat power is the product of torque and angular velocity, and the total torque of the ball bearing is a sum of three torques: the torque due to the external load M F , the torque due to the viscous friction M v and the torque due to the spinning of the rolling ball on the bearing race M s .
The torque due to the external load M F is calculated as: Where F β is the dynamical equivalent load of the bearing.f 1 is defined as f 1 = η(F s /C s ) y .F s is the static equivalent load.C s is the basic static load rating.η denotes the coefficient related to the nominal contact angle of the bearing.
The torque due to the viscous friction M v is defined by the empirical equation: Where f 0 depends on the lubrication form.n is the revolution speed of the spindle.v 0 is the kinematic viscosity of the lubricating oil or grease and v 0 is calculated by the viscositytemperature curve.
The spinning torque M s is calculated as: Where µ is the friction coefficient and ranges from 0.06-0.07(Jones, 1959).c is the eccentricity of the ellipse c = 1 − b 2 /a 2 .a and b denote the semi-major and semi-minor axes of the elliptic contact area, and both of them are decided by the curvature radius and contact force between the raceway and the rolling ball.ψ 2 (c) is the complete elliptic integrals of the second kind and calculated as: The Newton iteration is adopted to solve the mechanical model of the bearing.And the flow chart is shown in Fig. 7.

Thermal contact conductance model
The spindle-bearing system contains a large number of joint surfaces conformed by two rough surfaces.The presence of the surface roughness results in the imperfect contact which severely affects the temperatures of the contact surfaces.As the bearing is the main heat source in the spindle-bearing system, the joint surfaces around the bearing are of great importance to the thermal analysis.
The thermal contact conductance (TCC) is used to evaluate the heat flux resistance caused by the joint surface.The influential factors of the TCC include surface roughness, contact pressure, thermal conductivity and material hardness.The elastic-plastic contact model of TCC proposed by Sridhar and Yovanovich (1996a) Where σ and m are the effective RMS roughness and effective absolute mean asperity slope respectively.k = 2k 1 k 2 /(k 1 + k 2 ) is the effective thermal conductivity of the material.P is the contact pressure.H ep is the elastic-plastic microhardness of the softer material.ε * c is the dimensionless contact strain and defined as ε * c = 1.67(mE/S f ).E and S f are the effective elastic modulus and flow stress of the material respectively.
However, m is unspecified in some experiments.Lambert and Fletcher (1997) correlated σ and m with experiments, and the estimation for m is given as Eq. ( 22).m = 0.076(σ × 10 6 ) 0.52  (22)   The relationship between the contact pressure and the elasticplastic hardness is given as (Sridhar and Yovanovich, 1996b): Where H * B = H B /3178, H B is the Brinell hardness.And Eq. ( 23) is applicable to the material with the hardness ranging from 1300 to 7600 Mpa.
The contact pressures of the joint surfaces will be affected by the interference value.Due to the uneven thermal expansion of the contact components, the contact pressure is not a constant during the simulation.
In addition, the TCC between the rolling ball and the raceway is calculated as (Nakajima, 1995): Where ψ 1 (c) is the complete elliptic integrals of the first kind.The TCC varies with the contact force Q ij and Q oj .

The convection model
The temperatures of the surfaces exposed to the air will be affected by the heat transfer of convection.The heat flux q of convection is decided by the convection coefficient and temperatures: Where T air is the measured environment temperature.The convection coefficient h could be classified into two types: the forced convection coefficient and the free-convection coefficient.The forced convection coefficient for the surfaces of the rotating components can be estimated as: Where k air is the thermal conductivity of the air.d is the equivalent diameter of the rotating component.Nu is the Nusselt number calculated as (Zhang et al., 2013): Where Pr air = 0.703, v air = 15.06×10−6 m 2 s −1 and n is the rotating speed.
For the stationary surfaces such as the housing and the flange, the free-convection with the natural air takes place.The convection coefficient is set as h = 9.7 W m −2 K −1 .

Experiment
In order to validate the proposed model, a test rig of spindlebearing system (Fig. 8) is designed and constructed.Four PT100 temperature sensors are used to measure the temperatures.Two sensors are fixed on the test rig as shown in Fig. 8c and the other two sensors are used to measure the temperatures of environment and hydraulic oil respectively.The axial thermal deformation is measured by the displacement sensors of laser CCD and eddy current.In order to measure the axial force of the bearing, the strain gauges (KYOWA KSN-2-120-E4) are stuck on the sleeve next to the front bearing.In order to eliminate the thermal-induced strains of the sleeve in the measurement, the strain gauges are configured as halfbridge.The relationship between the output voltage v 0 and the strain ε could be presented as Eq. ( 31).The variation of the axial force is proportional to the strain as Eq. ( 32).After the calibration of the strain gauge, the variation of the bearing axial force could be measured.
Where K s is the sensitivity coefficient of the strain gauge.E is the elasticity modulus.ξ is the poisson ratio.A is the sectional area of the sleeve.
The hydraulic cylinder is used to adjust the initial preload of the bearings.The spindle-bearing system in machine tools usually works under the fixed-position preload.Therefore, the experiment is conducted under the same condition.The process of preload adjustment is carried out as following: (1) Adjust the hydraulic pressure to the target value.(2) Use the circular locking nuts in Fig. 8c to fix the axial displacement of the bearings.(3) Unload the hydraulic pressure.

Comparison
The thermophysical parameters of materials in the simulation are obtained from the research of Zhang et al. (2013).
Bearing parameters are from the measurement results listed in Table 1.
The experiment under the initial preload of 5300 N is taken as an example to investigate the thermal characteristics.The test rig starts at the steady state.The spindle motor rotates at 500 RPM for 8400 s and then stays still.
Besides the proposed model, another two simulation cases are presented as comparisons.Case 1 is the simulation without considering the TCC.For the model simplification, many previous studies (Li and Shin, 2004;Lee et al., 2015;Jiang and Mao, 2010) neglected the TCC.In such a case, AN-SYS will set the TCC of the joint surface with interference fit as a large value by default.In Case 2, the dynamic updates of boundary conditions only involve the uneven axial expansions of housing and shaft.The bearing radial dimensions, bearing stiffness and TCCs are considered as constants.These variables are calculated at the initial state of the experiment.Such simulation (Hongqi and Shin, 2004;Tak-abi and Khonsari, 2013) was commonly used in the thermalmechanical analysis of the spindle-bearing system.
As there is no cutting force imposed on the shaft, the measured strain of the sleeve is totally caused by the variation of bearing preload.The measured and simulated results of bearing preload, temperatures and thermal displacements are presented in Fig. 9.The maximum absolute error (MAX) and the root mean square error (RMSE) are used to evaluate the prediction accuracies as presented in Table 2.
For the bearing preload, the simulation result of the proposed model is in close agreement with the experiment result.The overall variation of the bearing preload could be divided into four stages: 1. Stage1: 0-350 s The bearing preload in the experiment increases dramatically from 5300 to 5782 N in 350 s.There are two possible reasons for the increase of bearing preload.One reason is the increasing rotating speed.Figure 10a is the relationship between the axial displacement and the axial force for NSK7220C without considering thermal effects.The initial preload is loaded without rotation (n = 0 RPM).As the test rig works under fixed-position The axial force (N) The axial displacement (m) ∆Z, n=0 r min -1 ∆Z, n=1000 r min -1 ∆Z, n=2000 r min -1 ∆Z, n=3000 r min -1 Variation with time  preload, the preload will rise with the increasing of rotating speed in order to ensure the bearing axial displacement unchanged.However, the increase of bearing preload induced by rotating speed is limited under low rotation speed and heavy load, and the increase is only 14 N under the operating condition of 5300 N and 500 RPM.
The other reason is thermal effect.Most of the heat remains in the bearing at the starting stage.The rapid radial expansion of the bearing changes the forces and motion state of the bearing.Figure 10b presents the relationship between the bearing axial displacement and the axial force at different times.Under the same axial displacement, the axial force increases dramatically.

Stage2: 350-8400 s
In this stage, the bearing preload decreases slowly.The heat in the bearing has been gradually transferred to other components.The interference of the joint surface between inner race and shaft is larger than that between outer race and housing.The difference in TCCs results in the uneven thermal expansions of shaft and housing.Due to the O configuration of the bearings, the axial displacement is reduced.The bearing preload decreases with the variation of the displacement.

Stage3: 8400-8750 s
There is no more heat generation in the bearing to support the high temperature.The influence of free convection on temperatures is significant.With the sharp drop of the temperature in the bearing, the bearing preload decreases dramatically.

Stage4: 8750 s till the end
The uneven thermal expansions of shaft and housing dominate the variation of bearing preload.The bearing preload increases slowly.
For Case 1, the heat in bearing races will be conducted to the housing and shaft at the same level.The ranges of temperatures in housing and shaft are close to each other.The gap between axial expansions of the housing and shaft is reduced.As a result, the preload of Case 1 in Fig. 9a is smaller than that of the proposed model.For Case 2, the bearing radial dimensions, bearing stiffness and TCCs are considered as constants.As presented in Fig. 9a, Stage1 and Stage3 in the preload variation disappear.
As an important variable, the bearing preload directly affects the calculation of bearing stiffness and heat generation.Eventually, the deviation in preload will lead to the prediction errors of temperature and thermal displacement.The RMSE of the thermal displacement for the proposed model is 3.818 × 10 −6 m.Compared with Case 1 and Case 2, the improvements are 17.22 and 13.64 % respectively.Other thermal characteristics of the proposed model are shown in Fig. 11.The TCCs of the bearings in Fig. 11a and b vary with the contact pressure of the joint surfaces.Due to the large interferences between inner race and shaft, the TCCs in Fig. 11a are quite large.The overall variation ranges of them are within 10 %.The joint surfaces between the outer race and housing are under the transition fit.Therefore, the TCCs are much smaller.In addition, the TCC of the front bearing varies more dramatically than the TCC of the rear bearing in Fig. 11b.The multilayer structure of the rear bearing constrains the radial thermal expansion of the hydraulic cylinder.The contact pressure of the rear bearing is larger than that of the front bearing, which results in the larger value of TCC for the rear bearing.Appropriate adjustment in the interference of joint surface and feasible structural design are effective ways for the temperature control.
The stiffness in Fig. 11c is the axial stiffness of two bearings in pairs (Eq.15).The ratio of minimum stiffness to maximum stiffness is 0.8718, namely, the thermal effects bring in a thirteen percent decrease in stiffness.

Conclusions
This paper presents a dynamic thermal-mechanical model of the spindle-bearing system.Based on the loop calling of the softwares of ANSYS and MATLAB, the proposed model provides an effective way to simulate the interactions between the boundary conditions and the FEM analysis.
In order to investigate the thermal feedbacks, a spindlebearing system with measurable and adjustable bearing preload is designed and constructed.Though the bearing preload varies all through the operation process, the proposed model achieves accurate simulation results.In addition, two traditional simulation models are conducted as comparisons.The comparisons reveal that TCC and thermal feedbacks on boundary conditions are non-ignorable for the accurate simulation.Besides the thermal drift at the end of the shaft, the thermal effects bring in a thirteen percent decrease in axial stiffness, which has a negative effect on the performance of the system.
An accurate thermal-mechanical model is the basis for the thermal performance analysis and thermal optimization.The proposed model provide a practical and accurate way to predict the transient thermal characteristics.It is believed that the employed modeling technique is applicable to a wide range of mechanical systems.
Figure 3. Geometric relations in ball bearing.(a) Angular position of rolling ball.(b) Displacement of the j th rolling ball.

Figure 4 .Figure 5 .
Figure 4. Velocities of ball bearing.(a) View of inner race.(b) View of rolling ball.

Figure 10 .
Figure 10.Relationship between axial displacement and axial force.(a) Without considering thermal effects.(b) Considering thermal effects.

Figure 11 .
Figure 11.Thermal characteristics of the proposed model.(a) TCCs between inner race and shaft.(b) TCCs between outer race and housing/hydraulic cylinder.(c) Axial stiffness of bearings in pairs.

Table 2 .
Prediction accuracy of the models.