Method for Vibration Response Simulation and Sensor Placement Optimization of a Machine Tool Spindle System with a Bearing Defect

Bearing defects are one of the most important mechanical sources for vibration and noise generation in machine tool spindles. In this study, an integrated finite element (FE) model is proposed to predict the vibration responses of a spindle bearing system with localized bearing defects and then the sensor placement for better detection of bearing faults is optimized. A nonlinear bearing model is developed based on Jones' bearing theory, while the drawbar, shaft and housing are modeled as Timoshenko's beam. The bearing model is then integrated into the FE model of drawbar/shaft/housing by assembling equations of motion. The Newmark time integration method is used to solve the vibration responses numerically. The FE model of the spindle-bearing system was verified by conducting dynamic tests. Then, the localized bearing defects were modeled and vibration responses generated by the outer ring defect were simulated as an illustration. The optimization scheme of the sensor placement was carried out on the test spindle. The results proved that, the optimal sensor placement depends on the vibration modes under different boundary conditions and the transfer path between the excitation and the response.


Introduction
Due to the increasing demands from aerospace, automotive, die/mold and other industries, machine tools with higher speed, precision and reliability are required urgently. In a machine tool, the spindle directly affects the cutting ability of the whole machine tool, since it either carries cutting tools as in milling operations, or work-pieces as in turning. In general, there are four types of spindle depending on the type of drives used: belt drive, gear drive, direct drive and integrated (built-in) drive [1]. The belt-driven spindles and gear-driven spindles are widely used in conventional machining, while the direct drive spindles and integrated spindles can run at higher rotation speeds, from 15,000 rpm. Recently, Abele et al. [2] have provided a detailed review on the historical development, recent challenges and future trends of machine tool spindles. Mechanical design of spindles, modeling the thermal and dynamical behavior of spindle units, sensor and actuator integration, were addressed in detail. In the application, the condition monitoring and fault diagnosis are very important to maintain the performance of spindles. Unexpected failure of spindles can lead to severe part damage and costly machine downtime, affecting the overall production logistics and productivity [3].
Bearings play an important role in machine tool spindle systems. Compared with hydrostatic, aerostatic or magnetic bearings [4], rolling element bearings are still most commonly used today in the spindles, which can provide the required precision, load carrying capacity, and spindle speeds. The bearings are usually extremely reliable, but failure can still occur by improper operation, overloading and thermal seizure. The early stage of bearing damage is often characterized by a sizeable local defect (pits or spalls) on inner raceway, outer raceway or rolling balls. When this occurs, subsequent rolling over of the damage will produce repetitive shocks or short-duration impulses. Extensive work has been carried out on the detection and diagnosis of bearing faults over the past several decades, as reviewed in several papers [5][6][7]. There is considerable interest in diagnostics and prognostics of rolling element bearings based on vibration analysis and signal processing [8]. Many powerful diagnostic methods have been developed to detect bearing failures based on advanced signal processing techniques such as envelope analysis [9,10], wavelets [11][12][13], empirical mode decomposition (EMD) [14][15][16], spectral kurtosis [17][18][19][20], and information entropy [21,22].
Most studies have focused on feature extraction of the faulty bearing from vibration signals, however, just a few of works attempted to model bearing faults mathematically so that the vibration signals generated by the defects can be simulated and explained theoretically. McFadden and Smith [23] gave a pioneering model to simulate the vibration signals produced by localized defects of rolling ball bearings. The successive impacts generated by a defect during rotation were modeled as a periodic impulse series. Based on this work, Wang and Kootsookos [24] proposed a general model for faulty bearing signals, which contained effects caused by non-uniform load distribution, machinery induced vibration and measurement noise. In another aspect, Ho and Randall [25] improved McFadden and Smith's model by incorporating slight random variations in the spaces of the excitation pulses. They pointed out that the random fluctuation causes the harmonic frequency components to smear laterally and this viewpoint gave a good explanation of the vibration spectra observed in reality. The bearing fault signals were generated as a series of impulse responses of a single degree of freedom (DOF) system. Similarly, Brie [26] also concerned the non-strict periodicity of the impacts induced by the defect and a concept called quasi-periodic impulse train was proposed for the impact signals. The bearing is modeled as a mass-spring-damper system with time-varying parameters. Later, Antoni and Randall [27] found out that vibration signals generated by localized faults were not exactly quasi-periodic since the random slippages are non-stationary and consequently developed a stochastic model to describe the vibration signals generated by localized defects. They treated the repetitive impacts as a non-stationary process for the first time. The vibration signals were viewed as the response of a linear system which is time-varying and accommodate some degrees of stochasticity. Choudhury and Tandon [28] modeled the rotor-bearing system as a 3-DOF system, and the vibration response due to a localized defect under radial load conditions were predicted with the model.
In the above research, rolling element bearings are modeled linearly and usually simplified as single DOF systems. However, due to the nonlinear Hertzian force/deformation relationship, the varying compliance vibration effect, and lubricant film effect, rolling ball bearings show high non-linearity and time-varying characteristics during operation, so a linear model is not adequate to express the dynamic behavior of bearings. Recently, Sopanen and Mikkola [29,30] developed a dynamic nonlinear model of deep groove ball bearings, which included the effects of non-linear Hertzian contact deformation and elastohydrodynamic fluid film. Distributed defects and localized defects were both taken into consideration. Sawalhi and Randall [31,32] presented an integrated simulation model for a gearbox test rig through the incorporation of a time-varying, non-linear stiffness bearing model into a developed gear model. The incorporated bearing model has the capacity to simulate both localized spalls (inner race, outer race and rolling balls) and extended inner and outer race faults (rough surfaces). Similarly, Rafsanjani et al. [33] proposed an analytical model to study the nonlinear dynamic behavior of rolling ball bearing systems including surface defects. Mathematical expressions were derived for local defects of inner race, outer race and rolling balls. Sawalhi and Randall [8] simulated the acceleration time signal responses resulting from a rolling element entry into and exit from a typical spall, which has the potential to enable quantification of the fault size.
Bearings in rotating machines are mechanically coupled to supporting structures, and thus defect-induced transient signals are often masked by interfering signals or background noise [34]. When bearing defects occur, the generated impacts can cause system vibration at many frequencies from different structures. The measured vibration responses to the localized defects are affected by the transmitting media between defects and sensors. However, most existing fault models don't contain the effects induced by other components other than bearings.
In the machine tool spindle system, bearings work as a subsystem. Accurately predicting the vibration signals of the spindle with bearing defects remains a challenging task because of its complicated nonlinear behavior. In this paper, the machine tool spindle system is focused on and an integrated finite element (FE) model for a spindle-bearing system is proposed. The nonlinear bearing model is developed based on Jones' bearing theory [35], while the drawbar, shaft and housing are modeled as Timoshenko's beam. The bearing model is integrated into the FE model of drawbar/shaft/housing by assembling equations of motion. In this way, bearings are coupled with other spindle components. The dynamic model is validated on a test spindle and then used to simulate the vibration signals generated by localized defects of the bearing outer ring. Considering the influences of the transmitting media between defects and sensors, the dynamic responses at different positions are compared to obtain an optimal sensor placement strategy.
The rest of the paper is organized as follows: in Section 2, a dynamic model of a spindle-bearing system is given with experimental validation. Localized bearing faults modeling is presented in Section 3, followed by the optimization of sensor positions in Section 4. The conclusions are given in Section 5.

Dynamic Model of a Spindle-Bearing System
The machine tool spindle system usually consists of drawbar, shaft, housing, bearings and other accessories. A typical test spindle is shown in Figure 1, which has five angular contact ball bearings in an O-type configuration. This spindle was used in [36][37][38] to demonstrate the FE modeling process as well. The spindle drawbar/shaft/housing and the angular contact ball bearing are modeled separately. The drawbar, shaft, and housing are modeled as Timoshenko's beam, while the nonlinear bearing model is developed based on Jones' bearing theory which considers the centrifugal force and gyroscopic effects from rolling balls. The bearing model is integrated into the FE model of shaft/housing by assembling equations of motion. The pulley, clamping unit and other accessories are modeled using rigid disk elements. The spacers between bearings are modeled using bar elements. This FE modeling process can be extended to high-speed spindles of integrated structure, in which the motor is modeled using rigid disk elements. With the FE model of the spindle, the Newmark integration method is used to obtain the vibration responses numerically.

Drawbar/Shaft/Housing Modeling
The Timoshenko beam element (as shown in Figure 2) is used to model the drawbar, shaft and housing. The motion of each node q is described by three translational (δ x , δ y , δ z ) and two rotational (γ y , γ z ) degrees of freedom: where the torsional motion is not included. The equation of motion for the drawbar/shaft in matrix form is expressed by including the centrifugal force and gyroscopic effects as: (2) where K b and M b are the stiffness and mass matrices of beam elements, G b is the skew-symmetric gyroscopic matrix of the rotating shaft, M C is the mass matrix for the centrifugal force effect on the shaft, Ω is the rotating speed, and F b is the external force vector. The subscript b represents the beam. For the housing, Equation (2) is still applicable. The rotating speed Ω is set to zero as the housing is stationary. In the spindle system, disk and other accessories (e.g., clamping units, nuts) are commonly used, and they are modeled using rigid disk elements. The equation of motion is given as: where M d is the mass matrix of disk elements, G d is the gyroscopic matrix and F d is the external force vector.

Nonlinear Bearing Model
The geometric drawing of an angular contact ball bearing is shown in Figure 3. The radii of the inner and outer ring grooves are: , where f i and f o are curvature ratios of the inner ring and the outer ring, D is the diameter of the bearing ball. The angular position of the k th ball is: .
Similar to the shaft/housing, the bearing is modeled with 'bearing element'. Each 'bearing element' consists of two nodes-the inner ring node and the outer ring node. The motion vectors of bearing nodes are expressed by q i (the inner ring node) and q o (the inner ring node), respectively: Under the working condition, the bearing rotates at high speed with applied axial loads and/or radial loads and the relative positions between the inner ring, ball and outer ring change consequently. Based on Jones' bearing model, the geometry relations of the bearings in the X-Y plane are shown in Figure 4.    Figure 4, the following geometric equations can be obtained using the Pythagorean Theorem: where sin ∆ , cos ∆ , δ ik is the ball-inner raceway contact deformation, and δ ok is the ball-outer raceway contact deformation.
Considering the plane passing through the bearing axis and the center of a ball located at azimuth φ k (see Figure 3), the load diagram without consideration of noncoplanar friction forces is obtained, as shown in Figure 5. θ ik -inner ring contact angle; θ ok -outer ring contact angle; Q ik -inner ring contact force; Q ok -outer ring contact force; F ck -centrifugal force; M gk -gyroscopic moment.
From Figure 5, the equilibrium of forces in the horizontal and vertical direction are given as: , Combining Equation (7) and Equation (8), the following equations are obtained: where ε 1 , ε 2 , ε 3 and ε 4 are error variables. The Newton-Raphson iteration method is used to solve the Equation (9) and the values of X bk , Y bk , δ ik , δ ok can be found. Then the contact loads are calculated by Hertzian theory for spherical contact as: where K i and K o are the load-deflection coefficients. The forces , , , , T applied on the inner ring of the bearing are given as: .
Similarly, the forces , , , , T applied on the outer ring of the bearing can be obtained. The bearing stiffness matrix is obtained by calculating the derivative of the force , , , , T acting on the bearing rings with respect to the displacement , , , , T of bearing rings, namely: (12) The bearing stiffness matrix K B has the form: (13) where k xx is the axial stiffness, k yy and k zz are the radial stiffness, k θxθx and k θyθy are the angular stiffness, and other terms are coupled stiffness.

Finite Element Model of the Spindle-Bearing System
The finite element model of the spindle-bearing system is shown in Figure 6. The black dots represent nodes, where each node has three translational and two rotational degrees of freedom. The connections between the drawbar and shaft are modeled by spring elements. The displacement relations among the spindle shaft, bearings, and the housing correspond to the configuration and preload mechanism of bearings [39]. For the test spindle (as seen in Figure 1), the five bearings are preloaded by a hydraulic system. Therefore, The inner ring of the first bearing (B1), the outer ring of 3 x y x y x y the second bearing (A2), and the inner ring of the fifth bearing (B5) are fixed to the corresponding nodes of the shaft or housing. Other rings of the five bearings are sliding in axial direction. Figure 6. The finite element model of the spindle-bearing system.
By assembling the equations of the drawbar/shaft/housing and bearings, the following general nonlinear dynamic equation for the spindle-bearing system is obtained: (14) where are the mass, damping and stiffness matrices. C s is the structural damping matrix constructed from modal damping ratios identified experimentally. The bearing stiffness matrix K B depends on the displacement , among other factors. In turn, the displacement is affected by the system stiffness and the external force F. Therefore, the dependency of bearing stiffness matrix on displacement is the root cause of the nonlinearity in the spindle system. The details of the matrices (M, C, K) can be referred in [40].
If the external load is time-varying, then the system responses are time-varying as well and the equation of motion at time t + Δt is given as: . (15) Assuming that the displacement, velocity, and acceleration vectors at time t (t 0) are known, the time history responses of the spindle system (i.e., displacement t+Δt , velocity t+Δt and acceleration t+Δt ) at time t + Δt can be evaluated numerically using the Newmark time integration method.

Model Validation
The test spindle was hung using elastic strings as a free-free system. The impact forces were applied by a hammer on the spindle nose in the radial direction and the vibration responses at the opposite part of the nose were recorded by an accelerometer, as shown in Figure 7. The frequency response function (FRF) were measured by CutPro-MalTF® Fourier Analyzer.
In the simulation, some parameters must be identified. These parameters include joint stiffness between the drawbar and the shaft, and the system damping ratios. It is possible to obtain a good match between the simulation and the measurement by tuning these parameters. The joint stiffness was found to be 2 × 10 8 N/m in the radial direction. The modal damping ratios were borrowed from experimental data and used in the FE simulations. From the simulated FRFs, the first mode ( = 965 Hz) and fourth mode ( = 2,722 Hz) have been found to be the most dominant. Corresponding to these two

Mx + Cx + Kx = F
main modes, experimental model damping ratios of 1.6% ( = 955 Hz) and 2.9% ( = 2,745 Hz) were used in the simulation. Damping ratios of other modes were set as 3%.  The time-history response of the acceleration due to the impact force is investigated as well. In the simulation, the measured impact force shown in Figure 9 was input to the spindle-bearing model. This is a typical impulse force whose duration is about 0.24 ms. The acceleration time responses and their spectra at the front side and the rear side of the spindle housing are predicted in Figure 10(a,b). The phases of the two acceleration time responses are converse ( Figure 10(a)). From Figure 10(b), it can be seen that there are two distinct peaks in the frequency range of 0-5,000 Hz, which are corresponding to the two dominant modes of the spindle, i.e., the first mode ( = 965 Hz) and the fourth mode ( = 2,722 Hz). Therefore, the predicted acceleration time responses at both the front side and the rear side of the spindle housing are acceptable. The FE model of the test spindle is capable of predicting the acceleration time responses due to the excitation. If bearing defects occur in the spindle system, the defect in one surface of a bearing strikes another surface and then an impact force is generated, which may excite resonances in the spindle-bearing system. With the input excitation from an appropriate bearing fault model, the simulation of vibration response caused by the bearing defects is feasible.

Localized Bearing Defects Modeling and Vibration Simulation
Vibration signals measured on a spindle contain rich physical information about the operating conditions. When local defects (e.g., cracks, pits, spalls) exist in one of the bearing components, transient impact forces occur whenever such a defect on one surface strikes its mating surface. The impact forces will excite the vibration responses at the natural frequencies of bearing parts and housing structure. During the rotating process, a series of approximately equally spaced force impulses are produced. The repetition rate of the force impulses is equal to the characteristic frequencies of the bearing, i.e., ball passing frequency for the outer raceway, ball passing frequency for the inner raceway, and twice the ball spin frequency [41].
Similar to the impact force generated by a hammer test, the shape of the force impulse is modeled as a triangular form, as shown in Figure 11. In practice, the pulse form may not be of such a regular shape. The duration time ΔT can be determined by dividing the defect width (L) by the relative velocity (v r ) between the mating elements, i.e., ΔT = L/v r . The magnitude A of the impulse is affected by both the contact load and the angular velocity at the defective surface, which can also be modulated due to relative motion of the load zone. The accurate calculation of the magnitude A is beyond the scope of this paper, which will be addressed in the future work. To explain the localized bearing faults modeling process, the outer ring defect of the first bearing (HYKH61914) in the spindle is modeled as an example. The bearing parameters are listed in Table 1.
With the expressions given in [41], the outer race defect frequency is calculated as 298 Hz when the rotating speed of the spindle is 1,200 rev/min.

Type Pitch diameter (mm) Ball diameter (mm) Ball numbers Contact angle (°)
HYKH61914 85 6.35 32 25 The load is applied in the radial direction and there is a pit on the surface of the outer ring, as shown in Figure 12. When the rolling balls pass through the pit, the force impulses are generated. The size of the pit is assumed very small and just one rolling ball passes through the defective region at a certain time. The successive force impulses generated by the localized defect on the outer ring is shown in Figure 13. The time interval between two adjacent impulses is determined by the ball passing frequency for the outer raceway. In reality, there is a slight random fluctuation in the spacing between each force impulse because the load angle on each rolling element changes as the rolling elements enter and leave the load zone [25]. The impulse responses due to the outer ring defect are displayed in Figure 14. For each force impulse, a decay response is generated. The characteristic frequencies of the outer ring defect are shown in Figure 15. The high harmonics of the defect frequency in the range of 4,000-8,000 Hz are quite clear (Figure 15(a)). The envelope spectrum ( Figure 15(b)) is obtained by using the Hilbert transform technique. The defect frequency of the outer ring (298 Hz) and its harmonics are extracted clearly. For the spindle-bearing system, out-of-balance contributions at the shaft rotational frequency f r , is always evident as various sources (e.g., shaft mis-alignment) contribute to unbalanced rotation of the shaft. To include the out-of-balance effect, the vibration responses induced by force impulses are modulated in amplitude, as shown in Figure 16. The characteristic frequencies of the outer ring defect with the out-of-balance effect are shown in Figure 17. Besides the defect frequency of the outer ring at 298 Hz and its harmonics, the rotating frequency f r (20 Hz) is extracted as well in Figure 17  For measurements in practice, the noise problem can never be avoided. For enhancing the simulation towards more realistic data, the impulse responses are contaminated with a normal distribution N(0,1) noise such that: (16) where w is the noise level. As an illustration, the amplitude modulated impulse responses of the outer ring defect with noise contamination is shown in Figure 18. The characteristic frequencies of the outer ring defect with the out-of-balance and the noise contamination effects are shown in Figure 19. Due to the noise, the spectra are not so clear as in Figure 17, and the amplitudes at the characteristic frequencies decrease. It can be predicted that when the signal-to-noise ratio is very low, the characteristic frequencies might be merged. Consequently, advanced signal processing techniques become indispensable to extract the bearing fault features from the raw signals contaminated by noise.

The Optimization of Sensor Placement
For the condition monitoring and diagnosis of the spindle bearings, vibrations produced due to defects are usually measured using the velocity transducer or accelerometers. However, the measured vibration responses are affected by the transmitting media between defects and sensors. To obtain a better monitoring effect, the optimization of the sensor position is necessary.
In this study, the optimization scheme of the sensor position is carried out on the test spindle, as shown in Figure 20. The spindle is fixed on the machine tool at the flange by bolts. It is assumed that the defect occurs at the outer ring of the first bearing (HYKH61914). Since the spindle housing is stationary, seventeen measured points are chosen evenly on the surface of the spindle housing to detect the bearing fault. The measured point 7 is chosen on the fixed flange. The position of the measured point 4 is coincident with the defect bearing in axial direction. Next, the best measured point for detecting the bearing defect will be determined. ...
It is well known that the displacement of a vibratory system due to the excitation is related with mode shapes. In Figure 21, the first three mode shapes of the test spindle are displayed. It can be seen that, the flexural mode of the shaft is the dominant vibration mode. For the housing, the displacement at the rear is much larger than that at the front. From the aspect of the mode shapes, the point 17 at the rear of the spindle housing may be the best choice to detect bearing defects. On the other hand, it is better to monitor the vibration response signals at the nearest location to the defective bearing from experience [42]. From this point of view, point 4 may be the best position to mount the sensor, since it is the closest measured point to the defective bearing.
The acceleration responses due to the localized defect are simulated by using the dynamic spindlebearing model. The dynamic responses at different positions are simulated and then compared with each other to obtain an optimal sensor position. The peak-to-peak value of the vibration responses at each measured point is displayed in Figure 22. Since the measured point 7 is chosen on the fixed flange, the simulated vibration response is zero from the dynamic model. It is necessary to clarify that the vibration response at point 7 exists in reality and can be measured definitely. In our dynamic model, the boundary conditions are simplified and the flange is fixed rigidly on the ground. Therefore, the simulated vibration response at point 7 is zero. The simulation gives a different aspect: the amplitude of the vibration response at point 1 is the largest, which means that the point 1 is the best position to measure the vibration responses induced by the localized defect bearing. Neither point 4 nor point 17 will be chosen in practice. Although the point 4 is the nearest point to the defective bearing, it is also close to the fixed flange which lowers the vibration response. For point 17, although the displacement of the mode shapes is the largest, the vibration energy decreases largely due to the long transfer path between the excitation and the response. It can be concluded that, the optimal sensor placement depends on the vibration modes under different boundary conditions, and the transfer path between the excitation and the response.
The impulse responses of the outer ring defect at three typical measured points (1, 4 and 17) are shown in Figure 23(a-c). The impulse responses are equally spaced. It is clear that the amplitude of the impulse responses is different. Here, the noise is not included.  Next, the noise contamination is included to enhance the simulation towards more realistic data. The impulse responses with noise contamination at the three measured points (1, 4 and 17) are displayed in Figure 24(a-c). It can be seen that the equally spaced impulse responses are submerged in the strong white noise. The signal-to-noise ratios (SNRs) of the impulse responses with noise are calculated. For the three measured points (1, 4 and 17), the SNRs are −1.52 dB, −4.61 dB and −25.87 dB, respectively. The envelope spectra of the impulse responses with noise are shown in Figure 25. For the measured point 1, the defect frequency of the outer ring (298 Hz) and its harmonics are extracted clearly. The rotating frequency f r (20 Hz) is extracted as well (Figure 25(a)). For point 4 ( Figure 25(b)), the defect frequency of the outer ring and its harmonics can still be extracted, but the information of the rotating frequency is lost. The worst case is the point 17 shown in Figure 25(c), where none characteristic frequency is extracted because of the low signal-to-noise ratio. Therefore, for the test spindle, if bearing defects occur at the first bearing, the point 1 is the best position to mount the sensor. This case study illustrates the significance of sensor placement optimization. If the sensor is mounted at improper position, the bearing faults may not be detected.

Conclusions
In the current study, the vibration response simulation and sensor placement optimization of machine tool spindles were investigated using an integrated FE model. The localized bearing defects in spindles were modeled visually, and vibration responses generated due to the outer ring defect were simulated as an illustration. The results have shown that: (1) The FE model of the spindle is capable of predicting the acceleration time responses due to the excitation. With the input excitation from an appropriate bearing fault model, the simulation of vibration response caused by bearing defects in machine tool spindles is feasible. (2) The noise decreases the amplitudes at the bearing characteristic frequencies in the envelop spectrum. If the signal-to-noise ratio is very low, the characteristic frequencies are submerged and cannot be detected from sensor signals.