# Effects of Machine Tool Spindle Decay on the Stability Lobe Diagram: An Analytical-Experimental Study

Authors: Omar Gaber; Seyed M. Hashemi
Journal: Shock and Vibration
 (2016)

---

## Abstract

An analytical-experimental investigation of machine tool spindle decay and its effects of the system’s stability lobe diagram (SLD) is presented.  A dynamic stiffness matrix (DSM) model for the vibration analysis of the OKADA VM500 machine spindle is developed and is validated against Finite Element Analysis (FEA). The model is then refined to incorporate flexibility of the system’s bearings, originally modeled as simply supported boundary conditions, where the bearings are modeled as linear spring elements. The system fundamental frequency obtained from the modal analysis carried on an experimental setup is then used to calibrate the DSM model by tuning the springs’ constants. The resulting natural frequency is also used to determine the 2D stability lobes diagram (SLD) for said spindle.  Exploiting the presented approach and calibrated DSM model it is shown that a hypothetical 10% change in the natural frequency would result in a significant shift in the SLD of the spindle system, which should be taken into consideration to ensure chatter-free machining over the spindle’s life cycle.


---

## Body

## 1. Introduction

A great number of airframe structural components are manufactured by high speed milling, where problems can arise related to the instability in the process, dimensional errors in the work pieces, and even breakage of the tools. The instability of the process, a vibration phenomenon known as chatter, appears in the high removal rate roughing, as well as in the finishing of low rigidity airframe sections. Machine chatter has been researched since the early 20th century [1]. It has been well established that chatter is dependent on the dynamic behavior of machine tool structure and is directly linked to the machine vibrational characteristics [2], which is often expressed in terms of the frequency response function (FRF) of the system at the tool point. As the machine vibrational characteristics change, the machine stability behavior changes.In the mid-1990s, Altintas and Budak [3] presented an analytical form of the so-called two-dimensional (2D) stability lobe theory for milling. The two-dimensional stability lobe theory deals with the stability of solutions for dynamical cutting systems, which usually stands for the spindle speed and axial depth of cut. As a function of these two cutting parameters, the border between a stable cut (i.e., chatter-free) and an unstable one (i.e., with chatter) can be visualized in a chart known as stability lobes diagram (SLD). A stable cut is defined as cut where the tool tooth displacement decreases from one pass to the next dampening the effect of the initial deflection. If the displacement increases or stays the same causing a wave pattern on the part, the cut is unstable. Numerous effective experimental and analytical techniques have been developed to establish SLD and to predict stable processes in recent decades (see, e.g., Altintas and Weck [4]). A more recent review of methods of obtaining stability lobe diagram in high speed milling operation is presented by Palpandian et al. [5]. 3D stability lobes were also later established considering radial depth of cut as another parameter [6]. Both of these stability lobe theories can help to select the appropriate cutting parameters of the spindle speed and axial depth to avoid chatter in machining processes. Attempts have also been made to integrate the dynamical behavior variation of the part with respect to the tool position in 3D lobes construction, with application to thin-walled structure milling, in order to determine optimal cutting conditions during the machining process [7].The stability lobes calculation requires the dynamic parameters of the system, namely, stiffness, natural frequency, and damping ratio of the workpiece for each natural mode. Thevenot et al. [7] used FEM-based numerical models to determine the system’s stiffness terms, but the damping ratio cannot be easily calculated numerically, and therefore it is usually evaluated through experimental measurements. Ertürk et al. [8] modeled the machine tool spindle with a set of springs and dampers to simulate bearing behavior and to predict tool tip FRF. Cao and Altintas [9] investigated the effect of preload applied to a bearing on the overall system natural frequency and showed that the greater the preload, the higher the system’s natural frequency.Many researchers have studied the stability through machine behavior, assuming a rigid work piece. The tool tip transfer function is then elaborated through models or experimental approaches. In addition, most of the previous models reported in the open literature have been developed assuming that spindle-tool set dynamics do not change over the full spindle speed range. However, this assumption needs to be reconsidered in high speed machining, where gyroscopic moments and centrifugal forces on both bearings and spindle shaft induce spindle speed dependent dynamics changes. Furthermore, the change of spindle system dynamics has not been accounted for in most existing stability studies. Few studies were also reported considering the flexibility of work piece [10].In summary, the methods of obtaining stability lobe diagram (SLD) can be divided in three main categories, namely, experimental, semianalytical, and analytical approaches.The aim of experimental methods is to obtain SLD by conducting a series of experiments on work piece by machining it using a milling machine tool; while machining at a certain depth of cut along the tool path, forced vibrations turn into self-excited vibrations, causing the milling process to become unstable, that is, chatter onset. This procedure is used in various experiments and is repeated for various depth of cut and spindle speed combinations.In semianalytical methods, most of the parameters required to obtain stability lobes are calculated analytically. The modal parameters of spindle/tool-holder/tool systems, however, are obtained experimentally (e.g., using tap test), from the resulting FRF of the system. The system’s parameters can then be used to calculate its SLD.The analytical approaches (see, e.g., [3]) aim to obtain the transfer functions of structure, which is required in the stability model, and to eliminate the need for series of costly experimental impact testing at various points on a work piece. The dynamic characteristics of the entire (spindle/tool-holder/tool) system, contributing in the transfer functions of structure, have been shown to depend on a large number of factors, including holder characteristics, spindle shaft geometry and drawbar force, and the stiffness and damping provided by the bearings. Thus the analytical approaches, complemented with numerical methods such as FEM, eliminate the need for experimentation and save the time and cost involved in determining the stability lobes.The aim of this paper is to present a semianalytical stability technique, developed to incorporate the spindle’s dynamic behavior variations in the stability lobes diagram (SLD). The change in the spindle’s dynamic behavior, also referred to as aging, is generally caused by system’s bearings wear, translated through a reduction in the system’s natural frequencies. Exploiting and adapting the dynamic stiffness matrix (DSM) method [11], the model of the spindle system is generated, where the boundary conditions (BC) at the bearings are originally enforced using simple supports. The spindle DSM model is also validated against conventional FEM-based simulations generated in the educational version of commercial software ANSYS V13. Once the correctness of the DSM method has been established, the simple supports BC are replaced with linear spring elements to incorporate the inherent flexibility of the bearings. The experimental (working) fundamental frequency of the spindle system, determined from FRF data obtained from tap testing, is subsequently used to adjust the spring stiffness constants, leading to a calibrated dynamic stiffness matrix (CDSM) [12] model of the system. Using the system’s CDSM model and the experimental FRF data, the SLD can then be determined. Exploiting the proposed semianalytical method, it will then be possible to evaluate updated SLD and optimized machining parameters, should one know the variation of system’s stiffness (i.e., changes in the fundamental frequency) over the spindle life cycle.The application of the proposed model is demonstrated through an OKADA VM500 machine spindle, where the shift in the SLD resulting from a simulated 10% reduction of the system’s fundamental frequency and its effects on the machining parameters are investigated. The spindle, initially examined while mounted on the original machine tool, was then installed on a bench top fixture to carry out further experimentations.


## 2. Theory and Governing Equations

In what follows, the differential equations governing the bending-bending (BB) vibrations of a spinning beam segment are briefly discussed. Following mainly the theory presented by Banerjee and Su [11], the effect of torsional stiffness is assumed to be large enough so that the torsional vibrations can be ignored (see, e.g., [12]). Figure 1 shows the beam in a right-handed rectangular Cartesian coordinates system. The beam length is L, the mass per unit length is m=ρA, the polar mass moment of inertia per unit length is Iα, the principal axes bending rigidities are EI for both planes, and the torsional rigidity is GJ.Figure 1
Degrees of freedom of a spinning beam segment.The differential equations governing the bending-bending (BB) vibrations of a spinning beam segment can be written as [11, 12](1)EIyyu′′′′-mΩ2u+mu¨-2mΩv˙=0,(2)-2mΩu˙-EIxxv′′′′-mv¨+mΩ2v=0,where u and v are displacement of a point, P, at an arbitrary cross section at z, in the X and Y directions. The cross section is allowed to rotate or twist about the OZ axis. The resulting loads are then found to be in the following forms, written for shear forces as(3)Sx=EIxxu′′′,Sy=EIyyv′′′and for Bending Moments as(4)Mx=EIxxv′′,My=-EIyyu′′.Assuming simple harmonic motion, (5)uz,t=Uzcos⁡ωt,vz,t=Vzsin⁡ωt,where ω is the circular frequency of oscillation and U and V, respectively, are the amplitudes of u and v. Substituting (5) into (1) and (2), they can be rewritten as(6)EIyyU′′′′-mΩ2+ω2U+2mΩωV=0,(7)EIxxV′′′′-mΩ2+ω2V+2mΩωU=0.Now we introduce ξ=z/L and D=d/dξ, which are nondimensional length and the differential operator. We substitute them back into (6) and (7): (8)D4-mΩ2+ω2L4EIyyU+2mΩωL4EIyyV=0,(9)D4-mΩ2+ω2L4EIxxV+2mΩωL4EIxxU=0.Combining the above equations into one 8th order, the resulting differential equation in terms of W, satisfied by both U and V, is then written as(10)D8-λx2+λy21+η2D4+λx2λy21-η22W=0,where(11)λx2=mω2L4EIxx,λy2=mω2L4EIyy,η2=Ω2ω2.The solution of the differential equation is sought in the form (12)W=erξ.Substituting (12) into (10) leads to(13)r8-λx2+λy21+η2r4+λx2λy21-η22,where(14)r1,3=±α,r2,4=±β,r5,7=±iα,r6,8=±iβ,(15)α2=12λx2+λy21+η2+λx2-λy221+η22+16λx2λy2η2(16)β2=12λx2+λy21+η2-λx2-λy221+η22+16λx2λy2η2. Based on the roots presented in (14), solutions for U and V can then be written as (17)U=A1sin⁡αξ+A2cos⁡αξ+A3sinh⁡αξ+A4cosh⁡αξ+A5sin⁡βξ+A6cos⁡βξ+A7sinh⁡βξ+A8cosh⁡βξV=B1sin⁡αξ+B2cos⁡αξ+B3sinh⁡αξ+B4cosh⁡αξ+B5sin⁡βξ+B6cos⁡βξ+B7sinh⁡βξ+B8cosh⁡βξ. Based on Euler-Bernoulli beam bending theory, the corresponding bending rotation, that is, slope, about X and Y axes, Θx and Θy, respectively, is given by (see [12] for further details):(18)Θx=dVdz=-1LdVdξ,Θy=dUdz=1LdUdξ.By doing similar substitutions in load equations (3) and (4), one finds(19)Sx=EIyyL3d3Udξ3,Sy=EIxxL3d3Vdξ3,Mx=EIyyL2d2Udξ2,My=EIxxL2d2Vdξ2.Substituting the boundary conditions into the governing equations, enforced at z=0 (ξ=0) and z=L (ξ=1), leads to (20)δ=BR,where (21)δ=U1V1Θx1Θy1U2V2Θx2Θy2T,R=R1R2R3R4R5R6R7R8T.Substituting similarly for the force equation one finds(22)F=AR,where (23)F=Sx1Sy1Mx1My1Sx2Sy2Mx2My2T.Finally, the frequency-dependent, dynamic stiffness matrix (DSM) of the spinning beam segment, K(ω), can be derived by eliminating R, where the force amplitude is related to the displacement vector by(24)F=Kδand F=KBR=AR, where assuming KB=A finally leads to(25)Kω=AB-1.The element matrices are then assembled to find the global dynamic stiffness matrix, Kω, and the system boundary conditions are applied. Setting force vector F=0 (for free vibrations) then leads to the following nonlinear eigenproblem:(26)Kωδ=0,where the solution then consists of searching for system natural frequencies, ω, satisfying (27)Kω=0.


## 3. Numerical DSM Results

The DSM formulation was used to model an OKADA VM500 vertical milling machine spindle [13] (see Figure 2 for the schematic structure of the spindle). The spindle has a total of 10 bearings, with four sets of angular contact bearings and two other ball bearings, and two belt assemblies, which provide the necessary driving force from the spindle motor. The spindle dimensions and geometric/sectional data, modeled using CATIA V5, were taken from [13]. For the sake of modeling convenience, the tapered and curved areas of the spindle and tool were modeled as uniform sections. The spindle along with the tool-holder and tool was modeled as a single piecewise uniform, stepped system, as shown in Figure 3, where the connection between the tool and heat shrink tool-holder is assumed to be rigid, as instructed by the tool-holder manufacturer. In the case of system studied here, the grip force is 5000 to 12000 lbf, depending on the extent of the interference fit, that is, 0.001–0.002 in. The rigid connection assumption is justified, as the cutting forces required for aluminum are normally much less than the nominal grip forces, that is, 2000 N or 450 lbf of maximum reference force, as reported in [2]. Also the tool-holder is made of hardened steel and therefore the deflection of this connection is minimal.Figure 2
Schematic diagram of an OKADA VM500 spindle [13].Figure 3
Spindle model.The spindle was modeled using the DSM method presented in the previous section, where only one element per uniform segment (i.e., a total of 26 elements) was used and simply supported boundary conditions were applied at the bearings locations. The material properties of tool steel, elastic modulus ofE=210 GPa and mass density of ρ=7850 Kg/m3, were used. Modal analysis was carried out and the fundamental flexural natural frequency of the spindle system was found to be 2303 Hz. To validate the spindle DSM model, frequency data were compared with those obtained from a 154-element FEM model created in ANSYS V13 software [14]. The BEAM188 element used is one-dimensional (1D) line element and has 2 nodes and 6 degrees of freedom (3 translational and 3 rotational). Comparing the DSM and FEM frequency results shows minor differences, with less than 3% error, confirming the correctness and accuracy of the DSM model, making it attractive for practical modeling applications.


## 4. Experimental Modal Analysis

The experimental modal analysis results, obtained using tap testing, were used to determine the system’s nominal fundamental frequency. Experimental Apparatus (Figure4) includes the following components:(1)
352A21 (Light, 0.8 g) accelerometer,(2)
086C04 (5000 N Hammer) Hammer,(3)
SIM3 Module Photon+ Data Acquisition,(4)
MetalMax Software (TXF).The accelerometer was attached to the edge of the tool-holder using wax. The spindle was struck 10 times and the average FRF graphs were generated and the experimentally evaluated frequency response function (FRF) data were measured for the spindle system while mounted on the machine (see Figures 5 and 6). As can be observed from FRF graphs, the most predominant natural frequency measured experimentally is found to be 1722 Hz, corresponding to the lowest peak of the imaginary part of the FRF graph.Figure 4
Experimental equipment.Figure 5
Spindle experimental frequency response function (FRF) (real part).Figure 6
Spindle experimental frequency response function (FRF) (imaginary part).


## 5. Experimentally Calibrated CDSM Results

The system’s fundamental frequencies were found to be much higher than that of the experimentally measured fundamental natural frequency, that is, 2303 Hz (from DSM) and 2367 Hz (from FEM) versus 1722 evaluated experimentally. The large difference (34%) between numerical and experimental results can be associated with the bearings flexibility. As mentioned earlier in this paper, the spindle DSM and FEM models were both initially equipped with simply supported (frictionless pin) boundary conditions at the bearing locations, that is, infinite lateral stiffness. A calibrated dynamic stiffness matrix (CDSM) model of the system was created by implementing spring boundary conditions (i.e., finite lateral stiffness) in the spindle’s DSM model (Figure7). The spring constant, KS, was then varied to calibrate the model to the experimentally evaluated fundamental frequency of the spindle system (see Figure 8). As it is clear from Figure 8, as the spring stiffness, KS, is increased, the natural frequency asymptotically approaches the natural frequency value for the spindle with simply supported boundary conditions, that is, bearings with zero flexibility, or infinite stiffness.Figure 7
Spindle with bearings replaced with springs.Figure 8
System fundamental frequency versus bearing equivalent spring constant for a low speed spindle (in log scale).To achieve the system’s experimental natural frequency of 1722 Hz, the required bearing stiffness was found to beks=1.28×109 N/m (see Figure 8). Furthermore, the spindle’s critical speed was also evaluated using both uncalibrated DSM and calibrated CDSM models and is found to be 138,464 RPM and 103,142 RPM, respectively (See Figure 9). Both values, however, are well above the operational RPM of the spindle; the suggested manufactures operational speed is 6000 RPM.Figure 9
Spindle natural frequency versus spindle RPM; high torque spindle.


## 6. Stability Lobes

The stability lobes were generated using the experimental FRF data obtained. The peak picking method was used to determine the equivalent stiffness at the tool tip [2]. Consider again the initial FRF data, shown in Figures 5 and 6. In the imaginary graph (Figure 6), the frequencies labeled at C correspond to the minimum peaks of the FRF. The corresponding frequency at this peak represents the dominant natural frequency of the system, ωn. The equivalent damping ratio at each mode can be calculated by using the two frequencies in the real part (ωA and ωB):(28)ϛ=ωB-ωA2ωn.The peak value of FRF in the imaginary part (Figure 6) is used to calculate the equivalent stiffness values at each mode:(29)Kexp⁡=-12ϛD.Using above Kexp⁡, one can then generate stability lobes diagram (SLD) using methods described in [2] (see Figure 10, the graphs in red color, where damping ratio of ϛ=0.029136 and the equivalent stiffness value of Kexp⁡=3.2289E+07 N/m have been used to generate the SLD).Figure 10
Stability chart an OKADA VM500 spindle using 2 different frequencies.Furthermore, exploiting the calibrated CDSM method, it is possible to evaluate updated SLD and optimized machining parameters, should one know the variation of system’s stiffness (i.e., changes in the fundamental frequency) over the spindle life cycle. This reduction can correspond to a combination of bearing degradation and loss of bearing preload due to age and/or damage. For the spindle system at hand, however, the authors were not able to find any data publicly available on the time history and reduction of the system’s fundamental frequency versus service time. However, to demonstrate the application of the presented CDSM and the process, let us consider a hypothetical 10% reduction in the system’s fundamental frequency.As presented in [2], the stability lobes are generated based on a spindle structure modeled as a single-DOF system, subjected to end force. However, the spindle models investigated here consist of continuous DSM elements (or large number of beam FEM elements) and the supporting bearings represented by multiple linear spring elements. Therefore, one needs to first convert the rotating system to a spring-damper system, in order to evaluate the equivalent stiffness of the entire updated system, Ktheo. In what follows, the process of evaluating Ktheo is briefly presented.Exploiting the graph given in Figure8 (the system’s natural frequency versus spring constant), the required bearing stiffness, Ks, in the calibrated DSM model to achieve a fundamental frequency of 1550 Hz (i.e., 90% of the nominal experimental frequency) is first evaluated (ks=7.72×108 N/m). The updated calibrated dynamic stiffness matrix of the system is then converted to a static stiffness matrix by linearizing the nonlinear eigenvalue problem of the system, as described by Hashemi (1998) [15]. Based on this method, the static stiffness of the system can be evaluated at zero frequency (i.e., no vibration: ω=0), or in the neighborhood of any given frequency [16].Let us consider a very small natural frequency,ωL, in the vicinity of zero (ω≈0). Adding a very small frequency increment, Δω, to ωL, leads to (30)ωU=ωL+Δωand an average frequency value, ω-, defined as(31)ω-=ωU+ωL2.Linearizing the updated system’s nonlinear eigenproblem (26), in the vicinity of ω-, leads to(32)Kω-=K-ω-2M,where Kω- is the system DSM and K and M are the static stiffness and mass matrices, respectively. Subsequently, the mass matrix of the system is written as(33)M=-dKω-dω-2and substituting (33) back into (32), the system static stiffness matrix is obtained as (34)K=Kω-+ω-2Mω-.Once the static stiffness matrix of the system, [K], has been evaluated, one can then apply any force to the node representing the tool tip in the spindle model and calculate the deflection, x, in the same direction as the applied force, Fx, using the following expression:(35)U=K-1F.This deflection, x, is then used to calculate the theoretical tool stiffness, as follows:(36)Ktheo=Fxx.The updated system’s parameters are then used to develop the revised stability lobes diagrams (SLD). The updated stability lobes, together with original SLD, are shown in Figure 10 graphs in blue, where damping ratio of ϛ=0.029136 and the equivalent stiffness value of Ktheo=2.90601E+07 N/m have been used to generate the SLD. As it is clear from Figure 10, a 10% reduction in the system’s fundamental frequency leads to a significant shift in the stability lobes. If this shift is not accounted for and the machining parameters are not updated, a chatter-free cutting program can start yielding chatter as the spindle ages and/or gets damaged. A typical 10% reduction in the system’s natural frequency over the spindle’s lifetime has been previously observed and reported in an earlier work by the authors [17]. Furthermore, it is worth noting that an FEM model (see, e.g., [14]) could be developed and calibrated for the spindle modeling purposes and subsequently used to produce the system’s FRF. However, the aim of the present study was to introduce the dynamic stiffness matrix (DSM) approach to assess the system’s dynamic behavior. The DSM method, formulated based on the continuous beam model and exact within the limits of the theory, is validated against FEM and used to ultimately evaluate Ktheo, representing the system’s overall stiffness, treated as an equivalent single-DOF system. It is to be noted that one could equally use conventional FEM to achieve the same goal.Finally, the reason for not using the FRF directly and entirely is the extremely large amount of data points in the FRF file. The goal of the present study was to use a quick approach to generate stability lobes that would benefit someone working as a machine operator in a manufacturing facility, with limited computation capabilities, that is, the CNC machine programmers. The stability lobes generated in this research, however, were confirmed against ones directly generated from FRF graphs. Additionally, as also suggested in [18], including the higher/additional modes obtained from the FRF does not increase the accuracy of the stability lobe chart significantly.


## 7. Conclusion

An experimentally calibrated dynamic stiffness matrix (DSM) formulation was presented and used to model an OKADA VM500 spindle. Unlike conventional FEM, the modeling method commonly used in this field, the presented element DSM is exact within the limits of the theory. As a result, in comparison with FEM, the DSM model requires much fewer elements to achieve the desired precision, that is, one element per uniform segment. The system’s natural frequency found experimentally was used to calibrate the spindle model, equipped with bearings represented as linear spring elements, and to generate spindle’s stability lobes. This natural frequency was then reduced by 10% to simulate a bearing degradation, or loss of preload due to spindle aging and/or damage. Incorporating the reduced frequency into the calibrated DSM, it was utilized to calculate the spindle’s equivalent stiffness value, used to generate updated stability lobe diagram. It was found that a 10% reduction in frequency leads to a significant shift in the stability lobes, which should be taken into consideration to ensure that a machine spindle runs chatter-free over its entire life cycle.



---

*Source: 4390185-2015-12-29.xml*