Lifetime analysis of motorized spindle bearings based on dynamic model

The traditional probabilistic-based lifetime evaluation methods for motorized spindles ignore the effects of load dynamic and structural difference. Therefore, we propose a dynamic model-based lifetime estimation method that combines these effects to improve the estimation results for motorized spindles, especially at the design stage. Given that bearing lifetime dramatically influences the reliability of motorized spindles, this paper also designs a shaft-bearing-toolholder based on a dynamic model to estimate the lifetime of bearing group. The proposed dynamic model closely resembles the actual structure of spindles and indicates the stiffness of bearings and contact surface conditional on the nonlinearity of inputting radial and axial forces. The stiffness model is verified by performing an experiment and a finite element analysis. The load applied to bearings is accurately calculated using the dynamic model. Afterward, the load is introduced to a well-known bearing lifetime model, and the lifetime of each bearing and bearing group is calculated. The bearing lifetime results obtained under preload, clamping force and cutting force conditions are then discussed.


Introduction
As a key part of machine tools, the performance of a motorized spindle greatly influences the high-quality and efficient production of CNC machine tools [1][2][3]. However, given infrequent failures and slow performance degradation processes, obtaining enough data to evaluate the lifetime of motorized spindles consumes much time. Moreover, different structures and operation loads lead to the significant heterogeneity of motorized spindles and necessitate an individual estimation of their lifetime. Bearing has a crucial role in motorized spindle systems, whose failures account for a large proportion of system failures [4][5][6]. Therefore, to improve the reliability of motorized spindles at the design stage [7], an accurate and efficient method for evaluating the lifetime of bearings in consideration of differences in the applied load should be utilized.
The lifetime of machine tools and their components are commonly evaluated using probabilistic models [8]. However, the estimation accuracy of these models depends on the amount of fault data whose collection is timeconsuming for long-life systems or components. To improve estimation accuracy, Yang [9] and Peng [10] proposed the small sample size method to relieve the demand for a large amount of data. To develop a highly accurate model that considers other factors, Yang [11] and Mu [12] considered the influence of load conditions on the reliability of spindle and machine tools. Some scholars also proposed using prior information and observation data to establish a degradation model for spindle failure. Peng [13] analysed the degradation of a spindle using the Bayesian method under time-varying degradation rates. Guo [14] proposed a Bayesian information fusion method to model the degradation process of a heavy-duty machine tool's spindle from multiple sources. Both statistical and degenerate models use a large amount of prior data to eliminate the influence of external factors on homogeneous structures. However, for vastly different structures, these models may generate inexact or misleading results. Through its transmission across different structures, the load applied to the component may differ even if the input force and component are the same. Such load difference can change the lifetime of components and thereby impact the reliability of the entire system. Moreover, the input force is dynamic during the operation of spindles, indicating that the load applied to the components is also dynamic. This dynamic load has a complex relationship with component lifetime.
Dynamics-based lifetime estimation methods have been proposed to better estimate the lifetime of products having different structures, dynamic forces and long lifetimes. These methods consider the structure of the system and can accurately calculate the dynamic load applied to components. An accurate load result is very useful in lifetime estimation. Another advantage is that these methods do not depend on lifetime data, thus enabling them to quickly evaluate the reliability of new products at the design stage. Fu [15] analysed the fatigue damage of flange and bolt on a wind turbine tower with vibration response. Nejad [16] analysed the long-term fatigue damage of the gear tooth root bending in the drivetrains of a wind turbine with dynamics. Cheng [17] introduced a wear model for the ball screw mechanism that considers dynamic fluctuations in loading. For a motorized spindle system, Li [18] used dynamic load to calculate the lifetime of the shaft on the motorized spindle yet ignored the influence of the joint surface structure on load. Zhang [19] analysed the influence of the preload and bearing angular deviation on spindle system lifetime based on a dynamics model yet merely focused on the effect of static load on a simple structure without considering the toolholder. To the best of our knowledge, a lifetime analysis of motorized spindles that considers both the shaft-bearing-toolholder structure and dynamic load is yet to be conducted.
To address this gap, on the basis of the structural characteristics of a motorized spindle, a multiple-degrees dynamic model of the shaft-bearing-toolholder structure is established in this paper. This model integrates the influence of axial and radial forces on bearing and contact surface stiffness and analyses the dynamic response of each bearing in the motorized spindle under different conditions of preloads, cutting forces and clamping forces.
The rest of this article is organized as follows. Section 2 establishes the dynamic model for the shaft-bearingtoolholder structure in a motorized spindle. Section 3 presents a novel method for evaluating the lifetime of the bearing group in a motorized spindle based on the dynamic model. Section 4 experimentally verifies the stiffness calculated by the dynamic model. Section 5 compares the lifetimes of each bearing and bearing group under different conditions of bearing stiffness, clamping force and rotating speed. Section 6 concludes the paper.

Dynamic model of motorized spindle
The internal structure model of a motorized spindle is shown in Fig. 1. Due to the shaft deflection, gyroscopic moment and nonlinearity of bearing stiffness, the bearing load motivated by cutting force is dynamic and nonlinear across different positions. Therefore, the response of bearings in various positions differs. In this case, a simplified dynamic model of a motorized spindle should be established to determine the relationship between cutting forces and bearing loads.
On the basis of its internal structure, a motorized spindle can be divided into three parts, namely, shaft, bearing and toolholder. This section establishes a shaft-bearingtoolholder dynamic model whilst considering the interaction amongst these parts.

Establishment of the dynamic model
The finite element and transfer matrix methods are most commonly employed in establishing a rotor dynamics model [20,21]. This paper uses the Timoshenko beam element in the finite element method to model the dynamics.
The proposed dynamic model of the shaft-bearingtoolholder system is shown in Fig. 2. This model has two front bearings and two rear bearings. Preload is added to the bearings to ensure adequate support stiffness. The contact surface matches the toolholder and shaft through the clamping force. The vibration characteristics of the whole shaft-bearingtoolholder system can be obtained when applying cutting force to the end of the toolholder.
According to Timoshenko beam theory, the finite element model of spindle bearing is established by [22]: where [M] is the mass matrix and [C] is the damping matrix; in this paper, proportional damping is used to define the total damping matrix. [G] is the Gyro moment matrix, Ω is the velocity of spindle rotation angular, [K] is the stiffness matrix for the shaft, {q} is the generalized displacement, [K b ] is the stiffness matrix of bearings, [K C ] is the contact stiffness matrix of the spindle contact surface and the toolholder contact surface, and {F} is the vector of cutting force, including distributed force and concentrated force.
The stiffness matrix of bearing [K b ] and the contact stiffness matrix [K C ] are treated as spindle support and components connection parameters that need to be established. The tooltip frequency response test is usually performed to measure bearing and contact stiffness, and this test relies on the performance of the detection equipment and requires repeated measurements of spindles with different structures [23][24][25][26]. To avoid this limitation and facilitate the prediction of spindle dynamic performance at the design stage, this paper uses the physical model to establish the bearing and contact stiffness. The physical model is known to have better accuracy in estimating bearing stiffness [K b ] compared with other models. On the basis of the Hertz contact theory quasi-static model, the five freedom analysis method is adopted to derive the load distribution and bearing deflection under any working condition [27].
Contact stiffness [K C ] refers to the elastic connection between the toolholder and shaft, which is a weak point that can influence the stiffness of the entire spindle system. Yoshimura [28] proposed a physical model that describes the relationship between contact pressure and contact stiffness. This model was recently developed in consideration of the influence of many factors, such as stress distribution, centrifugal force [29], roughness [30] or contact stiffness. However, the dynamic load, radial force and axial force on the joint surface are rarely considered in the dynamic model.
In practice, during the spindle rotation, the cutting force dynamically alternates, whilst the amplitude fluctuates. In addition, the dynamic axial and radial forces on the spindle constantly change the stress distribution on the contact surface, thereby changing the contact stiffness itself. Given that the damping in the system limits the amplitude of instantaneous load, focusing only on static load yields inaccurate prediction results that cannot reflect the actual load condition. Therefore, the influence of dynamic load on contact stiffness should be analysed.

Contact stiffness considering dynamic force
This subsection formulates the dynamic contact stiffness on the toolholder under the condition of dynamic axial and radial force.
Assume that the spindle-toolholder does not have any manufacturing error. As shown in Fig. 3, the toolholder is divided into N elements, one of which is subjected to radial force F cx and F cy and axial force F cz . Under the action of axial force F cz , the contact pressure of the contact surface is given by where F cla is the clamping force, S c is the area of the element contact surface, θ t is the taper angle of the toolholder and μ t is the friction coefficient.
The cross-section of the contact surface element is shown in Fig. 4. The resulting radial force F cxy on one element of the contact surface is produced by radial force F cx and F cy in directions X and Y. The angle between F cxy and the X is denoted by γ. Under radial force F cxy , the contact pressure on the one side increases, whereas that on other side decreases. Consider the contact pressure produced by the radial force on the contact surface as a cosine function distribution [31]. Afterwards, the increasing contact pressure P a1 (φ) and decreasing contact pressure P a2 (φ) dependent on angle φ are computed as: where P is the maximum contact pressure of the change pressure distribution.
Therefore, the rule of pressure change generated by radial force F cx and F cy along angle φ follows: if P a2 (φ)≥P cz , the contact surface in this angle loses contact. As a result, no pressure is generated, and the contact stiffness is equal to 0. When the contact pressure on one side completely disappears, the contact stiffness only exists on the other side.
After changing the contact pressure, the force equation in the γ direction is where S φ is the cone area per unit angle and μ is the friction coefficient.
The equilibrium equation under the resultant force F cxy is According to Eq. (6), the maximum contact pressure P can be obtained. Therefore, the change of unit contact stiffness on the contact surface [28] is where α and β are constant parameters dependent on the surface material properties.
Therefore, the contact stiffness of the contact surface is computed as the sum of normal contact stiffness along the X and Y, given by The contact stiffness between the two nodes in the dynamic model can be expressed as   6  6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  7  5 ð9Þ The complete toolholder model should include the contact stiffness of the toolholder-shaft and toolholder-tool. Given that the method for modelling the toolholder-tool is similar to that for modelling the toolholder-spindle, the toolholdertool modelling process is not described in this paper.

Motorized spindle bearings lifetime evaluation
Cutting force is the key load that affects the lifetime of a motorized spindle. By testing the cutting force under the common processing parameters for the selected spindle, the cutting force is measured, and the cutting force data under the stable amplitude state are selected as typical cutting force {F} inputs into the dynamic model Eq. (1) of the spindle to realize bearing lifetime evaluation. By calculating the displacement time history {q} of each node in the whole system, the node time history displacement of the bearing {q b } can be extracted from {q}. Afterwards, the time history load F b of each bearing in the spindle is equal to the product of the bearing stiffness matrix [K b ] and time history displacement {q b } computed as The load F b is a time series with variable amplitude load that can be changed into constant amplitude cyclic load via rain flow calculation. The 2D data of bearing forces are then obtained with different amplitudes and average values. According to the principle of the bearing force, the maximum load of each cutting force cycle is taken as the applied radial load F r of the bearing. Given that the preload of the bearing also has a significant impact on the bearing lifetime, the resultant force F 0 is computed as the sum of the preload and axial force in the cutting force. The equivalent dynamic load F m of the bearing is calculated as when F 0 /F r ≤e, X b =1. When F 0 /F r ≤e, X b is the rotation factor, and Y b is the thrust factor of the bearing.
According to the Palmgren-Miner rule, the equivalent dynamic load F m in the cyclic load in consideration of preload is transformed into the average payload P m , that is where N is a total number of cycles. For the ball bearing p=3, the lifetime of the single bearing under the specific cutting force according to the calculation method of bearing lifetime [32] is where a 1 is the lifetime parameter, a 2 is the lifetime coefficient of the bearing characteristics, a 3 is the lifetime coefficient of the working condition, n is the rotating speed, and C r is the rating radial load of the bearing.
Assume that the lifetime of other parts of a motorized spindle is infinite. The lifetime of the bearing group is close to that of the motorized spindle. Therefore, where e=1.1 and h is the number of bearings in the motorized spindle.
The overall flow chart of the paper is shown in Fig. 5 to clarify the calculation performed in this paper. Because it is difficult to obtain a large number of samples for the life test of motorized spindle, the accuracy of the results for long-life motorized spindle equipment cannot be easily verified using limited samples. However, the same method has demonstrated high accuracy for bearings. In this case, a motorized spindle dynamic model should be established to accurately estimate the lifetime of bearing groups. In this way, the lifetime of a motorized spindle can be evaluated at the design stage.

Dynamic model verification
The contact stiffness model is verified by performing a finite element analysis on the stress distribution of the cone element, and the static stiffness model of the motorized spindle is verified by performing the impact test and load experimental.

Stress distribution on the contact element
The contact surface of toolholder is divided into several elements, and the stress distribution of one element is analysed. The small end diameter of the element is 20mm, the taper of the toolholder is 7/24, the clamping force on the element is 80N, and the friction coefficient for the contact surface is 0.1. In the calculation, the radial forces range from 0 to 80N.
The von Mises' contact pressure under radial force is illustrated in Fig. 6, where the contact pressure is biassed to one side. The pressure increases on one side of the toolholder unit by increasing the radial force, whereas the pressure on the other side decreases. The centre of the tool holder shifts to the loading direction, which results in an uneven distribution of contact pressure on both sides. The theoretical value on the pressure-decreasing side differs from that obtained using the finite element method. In the theoretical calculation, the centre of the toolholder is assumed to be fixed, but in the finite element method, due to the effect of radial force, a relatively small displacement is observed between the centre of the toolholder and the spindle, which leads to a relative increase in displacement for the two contact surfaces on the pressuredecreasing side. As a result, the pressure calculated by the finite element method is less than the theoretical value. Meanwhile, on the pressure-increasing side, the error between the theoretical and simulation values is small. Whilst the corresponding error on the pressure-decreasing side is large, the overall trend is the same.
The change in contact stiffness under radial force is shown in Fig. 7. The contact stiffness increases linearly with the increasing radial force of the contact element. The radial force changes the contact pressure distribution on the contact surface, and the contact stiffness sharply increases on the contactpressure-increasing side. Given that the effect of axial load in cutting force on the toolholder is only to increase or decrease the clamping force in proportion, the influence of axial force is not discussed in this paper.

Verification of stiffness model of motorized spindle
To verify the dynamic model of a motorized spindle, an impact test is performed on the motorized spindle for assembling the tool and test rod. The model is JSZD105C produced by Shangdong Beat Precision. The shaft is supported by four angular contact ball bearings under 280N preload. The model of the toolholder is BT30. A PCB force hammer is used to excite the motorized spindle at the tip of the tool and the test rod, and an accelerometer is installed on the tip of the tool and rod. The experiment setup is shown in Fig. 8, and results of the amplitude-frequency response comparison between the experiment and simulation are shown in Fig. 9. The prediction   Given that the machining error, thread, sleeve and other details of the spindle are ignored in the modelling, the other natural vibration frequencies are not entirely consistent. Nevertheless, the prediction trend of the entire model is relatively consistent with the test results.
The stiffness of the motorized spindle is nonlinear with the load. To further verify the correctness of the model, the motorized spindle is loaded statically as shown in Fig. 10. The loading position is 175 mm extended from the spindle interface, and the measuring position is 90 mm and 230 mm away from the interface. To avoid errors resulting from the displacement of the fixed device, a displacement sensor is used to measure the displacement of the motorized spindle housing. The motorized spindle is then loaded from 0N to 2000N, and the results are shown in Fig. 11. The displacement result for the motorized spindle under constant stiffness is linear. In practice, the displacement in either the 90 mm or 230 mm positions is nonlinear, and the corresponding change rate decreases along with increasing load. Therefore, the variable stiffness model that considers nonlinearity closely reflects the experimental results.

Model analysis and discussion
To illustrate the feasibility and effectiveness of the proposed method, a motorized spindle in the machining centre is analysed as a case. In the spindle structure, the bearings from front to back are denoted by bearings 1 to 4. The rotating speeds under the two working conditions are set to 4000 r/ min and 2500 r/min. The acquisition equipment with a Kistler dynamometer brand collected cutting forces, as shown in Fig.  12 and Fig. 13.

Influence of preload and cutting force on bearing lifetime
The preload is the initial axial load on the bearings required by the motorized spindle to achieve high stiffness and increase rotational. The preload has a great impact on spindle stiffness, rotating precision and service lifetime.
To analyse the preload for the lifetime, the filtering of the cutting force data and the obtained cutting forces I and II are integrated into Eq.  accurately, the influence of preload on bearing load is eliminated. The bearing 1 at the front end of the motorized spindle bears the most load, and the load increases gradually along with increasing preload due to the fact that as the stiffness of the front bearing increases, the displacement of the shaft behind the bearing support position is restrained. As a result, the rear bearing bears less load. However, when the preload reaches a threshold value, the load on bearing 2 gradually increases probably due to the fact that the front-end bearing is located close to the rigid support, thereby increasing the gyroscopic moment effect of the spindle and the deformation of bearing 2. The rear bearings are not sensitive to increasing stiffness and almost unchanged under the influence of overall cutting force level. By converting the two measured cutting forces, their average values are 684N and 449N, respectively, and their ratio is 1.5 times. However, the load on bearing 1 is not within this proportion, and the ratio of the load to the cutting force of the front axle increases along with rigidity. Therefore, the effect of cutting force on the motorized spindle bearing is nonlinear. The effect of preload on the bearing load is shown in Fig.  15. The preload is converted into bearing load, which increases along with preload. When the preload is low, the load on the front bearing 2 decreases because bearing 1 bears more cutting force with an increased stiffness. Given that the other bearing loads are small, the preload gradually dominates the bearing loads, and the load increases linearly along with preload. Fig. 16 shows the lifetime of bearings and bearing groups in a motorized spindle under the influence of cutting force. Given that the position at the front end is equivalent to a cantilever beam, the front bearing 1 is subjected to the most load and has the shortest lifetime. With the support of bearing   1, the load of bearing 2 sharply decreases, and the lifetime of bearing 2 is extended as a result. As the preload gradually becomes dominant in the bearing load, the bearing lifetime initially increases and then decreases along with preload. Nevertheless, the lifetime of the bearing group continues to decrease. The lifetime curve shapes under different cutting forces are proportional. Fig. 17 shows the influence of different clamping forces on spindle vibration displacement amplitude under cutting force I. The preload of each bearing is 280N. With an increasing clamping force, the contact stiffness between the cone surface of the toolholder and the shaft increases, and the distance between them decreases, thereby improving the stability of the spindle under cutting force and reducing vibration displacement amplitude. Fig. 18 shows the lifetime of bearings under clamping force. With an increasing clamping force, the supporting stiffness of the toolholder increases, thereby reducing the inertial force caused by the bending of the toolholder under cutting force when the motorized spindle rotates. Therefore, the lifetime of front bearing 1 increases along with clamping force, and the lifetime of the other bearings slightly increases. In sum, the lifetime of the entire motorized spindle significantly increases.

Conclusion
This paper proposes a nonlinear dynamic spindle model for evaluating the lifetime of motorized spindle bearings in consideration of dynamic cutting force. The main conclusions are as follows: (1) A nonlinear dynamic model of a motorized spindle is established. The measured cutting force is used as the input force of the dynamic model to obtain the time history load on each bearing in the motorized spindle. Afterwards, the time history load is converted into average load via cycle counting and inputted into the bearing life model. The lifetime of each bearing and bearing group is then estimated. This method can reflect the load on motorized spindle bearings in actual working  established. The contact pressure distribution on the contact surface under radial and axial forces is analysed and verified using the finite element method. The change in contact stiffness under two types of load is then deduced. Results show that the contact stiffness is linearly proportional to radial and axial forces. In addition, the validity of the whole dynamic model is verified by performing impact and loading tests on the spindle. Test results show that the proposed dynamic model accurately calculates the bearing load. (3) The influence of preload and clamping force on bearing lifetime is analysed. The lifetime of the bearing group is mainly determined by bearing 1, which bears the heaviest load. The load on bearing 2 can be affected by front bearing stiffness and initially increases before decreasing. The lifetime of the rear bearing is generally determined by the overall level of cutting force. When the load on the bearing is small, the bearing lifetime is mainly affected by the preload. The contact stiffness increases along with clamping force, and the lifetime of the bearing group can be increased. (4) Given its nonlinear stiffness, the load on each bearing does not increase linearly with cutting force. According to the bearing lifetime model, the relationship between bearing load and bearing life is cubic. Therefore, considering the nonlinear stiffness of the motorized spindle dynamic model can provide a better understanding of the relationship between cutting force and bearing load and offer a highly realistic, and accurate result for bearing group lifetime can be offered.
Availability of data and material All the data presented and/or analysed in this study are available upon request from the corresponding author.
Code availability Not applicable.
Author contribution Jun Ying: background research, methodology, data curation, software, validation, writing (original draft) and editing. Zhaojun Yang: supervision, project administration and funding acquisition.