On the Free Vibration Modeling of Spindle Systems : A Calibrated Dynamic Stiffness Matrix

The effect of bearings on the vibrational behavior of machine tool spindles is investigated. This is done through the development of a calibrated dynamic stiffness matrix (CDSM) method, where the bearings flexibility is represented by massless linear spring elements with tuneable stiffness. A dedicatedMATLAB code is written to develop and to assemble the element stiffnessmatrices for the system’smultiple components and to apply the boundary conditions.The developedmethod is applied to an illustrative example of spindle system. When the spindle bearings are modeled as simply supported boundary conditions, the DSM model results in a fundamental frequencymuch higher than the system’s nominal value.The simply supported boundary conditions are then replaced by linear spring elements, and the spring constants are adjusted such that the resulting calibrated CDSMmodel leads to the nominal fundamental frequency of the spindle system. The spindle frequency results are also validated against the experimental data. The proposed method can be effectively applied to predict the vibration characteristics of spindle systems supported by bearings.


Introduction
The booming aerospace industry and high levels of competition have forced companies to constantly look for ways to optimize their machining processes.Cycle time has been a major concern at various industries dealing with manufacturing of airframe parts and subassemblies.When trying to machine a part as quickly as possible, spindle speed or metal removal rates are no longer the limiting factors but rather the chatter that occurs during the machining process.It is well established that chatter is directly linked to the natural frequency of the cutting system, which includes the spindle, shaft, tool, and hold combination.
The first mention of chatter can be credited to Taylor in 1907 [1]; however, the first comprehensive investigation of the problem was published by Arnold, in 1946 [2].Arnold's experiments were conducted on the turning process, where he theorized that the machine could be modeled as a simple oscillator and that the force on the tool decreased as the speed of the tool in relation to the work piece increased.Gurney and Tobias later theorized the now widely accepted belief that chatter is caused by wave patterns traced onto the surface of the work piece by preceding tool passes [3].The phase shift of the preceding wave to the wave currently being traced determines whether there is any amplification in the tool head movement.If there exists a phase shift between the two tool paths, then the uncut chip cross-sectional area is varied over the pass.Thomas and Shabana [4] showed that the cutting force is dependent on the chips cross-sectional area, and, as a result, a varying cutting force is produced.This was demonstrated through an illustrative example of a grinding machine, modeled as a one-degree-of-freedom mass-spring system, as opposed to an oscillating system.The spring-mass system is also the widely used modeling theory for how a vibrating tool should be characterized today.
Prior to 1961, the machining procedures were regarded as steady state and discrete processes [5].This erroneous idea led to the creation of overly heavy and thick-walled machines, believed to provide high damping to the tooltip forces (assumed to be static).However, one must note that machining is a continuous, dynamic process, with constantly fluctuating tooltip forces.Moreover, while being identical in construction, the dynamics and response of similar machines may be slightly different due to structural imperfections, imbalances, and so forth.Therefore, the specific characteristics and dynamics of each machine must also be 2 Shock and Vibration taken into consideration when performing chatter-related calculations [6].The structural modes of the machine determine the frequency and the direction that the tool is going to vibrate at [7].Eisele [5] and Peng et al. [8] stressed that designers must investigate the mode-shapes, weak points, bearing clearances, and self-inducing vibratory components of the machine to reduce chatter.Butlin and Woodhouse also investigated the number of structural modes required to produce accurate results [9].They showed that low order models, incorporating only two modes, are sufficiently accurate to model the machines.
In 1981, Tlusty and Ismail published one of the first papers that documented the nonlinearity of the vibratory system occurring during chatter [10].Kondo et al. [11] added to Tlusty and Ismail's investigation the behaviour of the tool after the onset of chatter.An important study, carried out by Lee et al. in 1989, helped bring a more complete theory of chatter to light [12], where they looked into the response of the work piece to chatter.They modeled the tool as a typical two-degree-of-freedom spring-mass system and the work piece as an Euler beam.When the responses of both the tool and the work piece were incorporated into the typical spring-mass equations, the results obtained far surpassed the accuracy of previously used equations.They observed a large increase in surface finish quality while machining, even with the use of long end mills.
While a comprehensive set of mathematical models have already been derived to explain chatter using spring-mass and Euler beam systems, which have led to dramatic increases in surface finish quality, chatter is still an issue in the machining industry as the surface tolerances keep increasing.From the 1990s on, this drove researchers to use more refined mathematical models to characterize the machines [13].With the rise of computer aided calculations and simulations, it has become increasingly easier to perform more complex analyses of the machining process [14].Tool wear is an oftenoverlooked factor that contributes to chatter.With the aid of more powerful computers, this variable can now be included in simulations.Clancy and Shin [15] investigated the wear of turning tools and how that affects the stability of the system.It was found that as the tool becomes worn, its limits of stability increase.As a result, stability lobes have now become a function of tool wear.Li et al. [16] determined that the coherence function of two crossed accelerations in the bending vibration of the tool shank approaches unity at the onset of chatter.A threshold needs to be set [17] and then detected using simple mechanism to alert the operator to change the machining conditions and avoid increased tool wear.
In most of the existing stability prediction methods, a frequency response function (FRF)-acquired using an impact test-is required to extract the system fundamental frequency, used to perform the calculations.In order to perform accurate impact tests, the machine would need to be shut off and allowed to cool down and the coolant drip off.This process could lead to several hours of down time.An offline method of obtaining the FRF could greatly benefit manufacturing companies by eliminating the downtime needed for the impact tests.While stability lobe graphs have become more accurate, the industry is trending towards computer simulations and numerical analyses.With computing power becoming relatively cheap and readily available, this method has become feasible.Providing the most accurate results due to its constantly updated stability calculations, this is the area of research that engineers will most likely progress towards.Lee et al. [18] suggested that the whole system comprising machine, tool and work piece could be modeled using finite element method (FEM) analysis.A computer simulation would be able to predict the FRF during all phases of the machining process.As the part is machined and becomes thinner, its response to vibration changes dramatically.Most existing research papers assume that the FRF remains constant throughout the whole process.However, a constantly updated FRF would allow for accurate, real time stability calculations.
In an attempt to numerically model the entire spindle systems for the vibration modeling and analysis and to incorporate the effects of bearings, Cao and Altintas [19] modeled spindle bearing and machine tool systems for virtual simulation of milling operations, where they presented a general, integrated model of the spindle bearing and machine tool system, consisting of a rotating shaft, tool holder, angular contact ball bearings, housing, and the machine tool mounting.The proposed model, also verified experimentally, allows virtual cutting of a work material with the numerical model of the spindle during the design stage and predicts bearing stiffness, mode shapes, frequency response function (FRF), static and dynamic deflections along the cutter and spindle shaft, and contact forces on the bearings with simulated cutting forces before physically building and testing the spindles.The study showed that the accuracy of predicting the performance of the spindles requires integrated modeling of all spindle elements and mounting on the machine tool.Jones' bearing model, which considers the bearing balls and rings as elastic parts, and the Hertzian contact theory were used to calculate the contact force and displacement.Patwari et al. [20] presented a detailed systematic procedure for experimental and analytical modal analysis techniques, to evaluate structural dynamics of a vertical machining centre.Abuthakeer et al. [21] carried out the dynamic characteristics analysis of high speed motorized spindles.More recently, Delgado et al. [22] used FEM to study thenonlinear dynamic behavior of a high speed machine tool spindle.
It is known that FEM predictions and the numerical results are often different from test and experimental results.However, based on the model updating, corrections can be made on the FEM models by processing records of test results.Motthershead and Friswell [23] presented a survey of model updating in structural dynamics.Later, Mares et al. [24] investigated the model updating using robust estimation, and more recently Mottershead et al. [25] presented a basic introduction to the most important procedures of computational model updating, including tutorial examples and a large scale model updating example of a helicopter airframe and discussed the sensitivity method in FEM model updating.
To the authors' best knowledge, except authors' previous works [26,27], no studies on the calibrated/updated dynamic stiffness matrix (DSM) models of spindle systems have appeared in the open literature.The DSM method [28] provides an analytical solution to free vibration problems, achieved by combining the coupled governing differential equations of motion of the system into a single higher-order ordinary differential equation.Seeking the most general solution to the resulting higher-order ODE and enforcing the end displacements and loads through extensive matrix manipulations lead to the DSM of the structural element.Unlike limited nodal DOFs in conventional FEM [29], the DSM-based element is modeled as a continuous system with infinite number of DOFs within the element.The DSM formulation has been proven to provide exact frequency (within the limits of the theory) for various beam configurations and beam-structures, with the use of a single continuous element per segment (see, e.g., [28,[30][31][32][33]).The objective of the present study is to determine the natural frequencies/vibration characteristics of machine tool spindle systems by developing the corresponding dynamic stiffness matrix (DSM) [31,32] and the proper boundary conditions.These results would then be compared to the existing results for a common cutting system to validate/tune the model developed.The coupled bending-bending (B-B) vibration of a spinning beam is investigated.A MATLAB code is developed to assemble the DSM element matrices for multiple components and to apply the boundary conditions (BC).The bearings are first modeled as simply supported (S-S) frictionless pins, which are then replaced by linear spring elements to incorporate the flexibility of bearings into the model.In comparison with the existing data on spindle's fundamental frequency, the bearing stiffness coefficient,   , is then varied to achieve a calibrated (or updated) dynamics stiffness matrix (CDSM) vibrational model.The proposed formulation can also be extended to include torsional degreeof-freedom (DOF) for further modeling purposes.

Problem Description and Governing Equations
Computer numeric control (CNC) machines are becoming the norm in manufacturing plants worldwide.These machines are generally 3-, 4-, or 5-axis, depending on the number of degree-of-freedoms the device has.Having the tool translating in the , , and  directions accounts for the first three degrees of freedom.Rotation about the spindle axes account for any further DOF.A typical spindle contains the motors that rotate the shaft connected to the tool holder, the bearings, the tools, and all the mechanisms that hold the tool in place (see Figure 1 for a sample spindle configuration).
As mentioned earlier in this paper, the fundamental natural frequency of this system dictates the limitation of the cutting parameters.If these limits are surpassed, the system tends to vibrate which causes chatter.Also, as the bearings wear, the system characteristics and spindle's natural frequency change.
In what follows, the derivation of the governing differential equations for the coupled bending-bending vibration of a spinning beam for a rotating shaft is first briefly discussed.Bearings restrain the spinning motion of the spindle system in five degrees of freedom, that is, displacement in , , and  directions and the rotation in  and  directions.The beam torsional rigidity, , is large enough so that the torsional vibrations can be ignored [31,32].Figure 2 shows a cylindrical beam in a right-handed rectangular Cartesian coordinates system.The beam length is denoted by , mass per unit length is  = , polar mass moment of inertia per unit length is   , and the principal axes bending rigidities are EI for both planes.At an arbitrary cross section along -axis,  and V are displacement of a point  in the  and  directions.The cross section is allowed to rotate or twist about the  axis.The position vector r of the point  after deformation is given by (refer to Figures 2 and 3) where i and j are unit vectors in the  and  directions.The velocity of point  is given by The kinetic and potential energies of the beam ( and ) are given by ( Using the Hamilton principle in the usual notation state, where  1 and  2 are the time intervals in the dynamic trajectory and  is the variational operator.Substituting the kinetic and potential energies of the beam in the Hamilton principle, collecting like terms and integrating by parts, the following differential equations are obtained: which govern the bending-bending (B-B) vibration of a beam, coupled by the spinning speed.The resulting loads, also force boundary conditions at free ends, are then found to be in the following forms, written for shear forces as and for bending moments as Assuming simple harmonic motion, and substituting (8) into (5), they can be rewritten as where  is frequency of oscillation and , , and Φ are the amplitudes of , V, and .Introducing  = / and  = /, which are nondimensional length and the differential operator, and back substitution into (9) lead to

Dynamic Stiffness Matrix (DSM) Solution.
As suggested by Banerjee and Su [31,32], the above equation (10) can be combined into one 8th-order differential equation written as where "" represents both lateral displacements  and  and The general solution of the differential equation is sought in the form Substituting ( 13) into (11) leads to where  1,3 = ±√,  2,4 = ±√,  5,7 = ±√,  6,8 = ±√, and Based on Euler-Bernoulli beam bending theory, the corresponding bending rotation, that is, slope, about -axis and -axis, Θ  and Θ  , respectively, is given by (see the appendix for further details) By doing similar substitutions in load ( 3) and ( 4), one finds Substituting the boundary conditions into the governing equations, enforced at  = 0 ( = 0) and  =  ( = 1), leads to where Substituting similarly for the force equation one finds where The frequency-dependent, dynamic stiffness matrix (DSM) of the spinning beam, K(), can then be derived by eliminating R. We also know that the force amplitude is related to the displacement vector by Explicit matrix forms of A, B, and R are given in the appendix.

Classical Finite Element Method (FEM).
For validation and comparison purposes, a conventional finite element formulation of the problem and a MATLAB-based FEM code were developed to carry out the free vibration analysis of the spindle system [33].The Galerkin method of weighted residuals [29,33] was employed to develop the integral form of the governing equation (10).Performing integration by parts on the resulting integral equations leads to the weak integral form of equations, which also satisfy the principle of virtual work.The system is then discretized using beam elements with four DOF per node (i.e., two lateral displacements,  and , and two slopes,   and   , eight degrees of freedom in total).Introducing cubic Hermite-type polynomial approximations [29,33] for bending displacements in the element weak integral form, the classical finite element formulation is developed.This formulation results in constant (or static) element stiffness, mass, and coupling matrices, which are then assembled to find the global matrices.Finally, enforcing the boundary conditions leads to a linear eigenvalue problem, which is solved to extract the system natural frequencies and modes.There has also been an attempt to develop a dynamic finite element (DFE) [34] formulation for the problem in hand, using frequency-dependent interpolation functions, and the preliminary results were found to be promising [35].

Application of the Theory and Numerical Results
The DSM formulation was first validated against various analytical and numerical frequency data available in the literature [31,32] for various uniform simple beam configurations.The DSM frequency data were also compared with those obtained from commercial FEM-based software ABAQUS and the effect of various element types on the convergence was also examined.In addition, the conventional FEM formulation [29,33] was used for further comparisons; however, as the in-house FEM-based code and commercial software both led to similar results, the FEM results have been omitted here for the brevity.A typical spindle system model (see Figure 4) was then constructed using multiple segments.As the nominal fundamental frequency was available for nonspinning spindle, therefore the effect of spinning speed was taken out by simply setting Ω to zero in FEM models.However, one should note that setting the spinning speed exactly to zero would cause the DSM formulation to collapse.This is due to the fact that the DSM formulation pivots around combining the two expressions (10) into a single on (11), through the coupling terms involving Ω [31,32].However, an extremely small number (Ω ≅ 0) can be used to model the stationary spindle.
The entire system was assumed to be made from the same material, that is, tooling steel.Initially it was also assumed that the system is simply supported at the bearing locations, that is, the bearings modeled as frictionless pins, rigidly attached to the spindle casing with zero flexibility, allowing only for rotations in all directions.This system was modeled using the DSM method described above, as well as finite element analysis (FEA) in ABAQUS commercial software to evaluate the first six natural frequencies of the spindle system.Several different element types and mesh sizes were examined to ensure the convergence (Table 1); beam element B33 uses a 2-node cubic beam (similar to the inhouse FEM code), C3D20R is a general-purpose quadratic brick element, and C3D10M is a general-purpose tetrahedral element (see Table 1 for results).Unlike the conventional FEM formulation, the DSM method results in exact solution for all natural frequencies of a spinning uniform beam [31,32].Therefore, as the spindle system in hand is made of 12 uniform sections, only a 12-element piecewise-uniform (stepped) DSM model was required.As expected, the results using beam element type B33 yielded the best agreement with the DSM results, as the 3D and shear effects are neglected in both.However, the fundamental frequency values were found to be different from the nominal value, that is, 1393 Hz compared to 1000 ± 50 Hz, as provided by the manufacturer.
In order to address this issue, the DSM and FEM models were calibrated to take into account the flexibility of the bearing systems.To this end, the simply supported boundary conditions were modified and replaced by simple linear spring elements, attached to the spindle shaft in  and  directions, representing general flexibility of the bearing systems (see Figure 4).The bearings were all assumed to be identical, with the same spring coefficient,   .The spring stiffness value,   , was then varied and a calibrated model equipped with spring elements of constant   = 2.1×10 8 N/m was found to result in a fundamental bending frequency of 1004 Hz, closely equivalent to the system's nominal value (Figure 5).The effect of spinning speed, Ω, on the system's fundamental frequency was also investigated.In the absence of spinning speed, Ω = 0 the system exhibits uncoupled behavior (see (7) and ( 8)), undergoing two separate bending displacements in  and  directions.As a result, for a circular shaft with   =   , the modal analysis results in repeated frequency values, appearing in pairs (Table 1).By introducing the spinning speed, Ω ̸ = 0, the bending natural frequencies in   and  directions become coupled, changing the repeated natural frequencies to two distinct ones (see Table 2), with one decreasing and approaching zero (Figure 6) and the other increasing from their uncoupled values.Further increase in the spinning speed would make one of the bending frequencies go to zero; that is, whirling instability, also reported and discussed in [31,32].The first critical spindle speed, in the case of calibrated model, was found to be Ω cr = 6 × 10 5 RPM, which is well above the operating speed of the spindle, that is, 3.5 × 10 4 RPM (see Figure 7).In order to further validate the proposed CDSM model, the tap testing method was used to experimentally collect the frequency response function (FRF) of the spindle system in hand, while being installed on its corresponding machine tool.A 1-inch diameter blank tool, with a 2-inch protrusion, in a typical shrink fit tool holder was used.The tool and holder were inserted into the spindle.The spindle was set to be horizontal (i.e., neutral) position and an accelerometer (352A21, 0.8 g) was then attached to the edge of the tool holder.The spindle was struck 10 times in the  direction using a hammer (086C04, 5000 N) and   the average FRF graph was generated (SIM3 Module Photon + Data Acquisition).Similarly, the spindle was also tested along the  direction.The nonspinning spindle's natural frequencies extracted from the FRF graphs are presented in Table 3, alongside the CDSM results.As can be seen, the maximum difference between the first three CDSM flexural  natural frequencies and the experimentally evaluated data are, respectively, 8%, 0.07%, and 14.3%, in  direction, and 8.4%, 0.06%, and 13.3%, in  direction.The average difference is found to be around 7.4%, which is within the acceptable range.

Discussion and Conclusion
It has been established that the service life changes the vibrational behavior of spindle systems, that is, reduced natural frequency over the time, associated with the bearings wear.Therefore, in order to generate an updated spindle model, it is essential to include the effect of the bearings in the system.An analytical model of a multisegment spinning spindle, based on the dynamic stiffness matrix (DSM) formulation and exact within the limits of the Euler-Bernoulli beam bending theory, was developed.The beam exhibits coupled bending-bending (B-B) vibration and, as expected, its fundamental coupled frequency was found to decrease with increasing spinning speed.The bearings were included in the model using two different models: rigid, simply supported, frictionless pins, and flexible linear spring elements.The simply supported bearing model was found to yield an excessively rigid system with a much higher natural frequency than the nominal value of the system.Replacing the simply supported boundary condition with linear spring elements, the spring stiffness,   , was then calibrated so that the fundamental frequency of the system matched the nominal value, that is, a more realistic and accurate representation of the spindle system.The accuracy of the proposed calibrated dynamic stiffness matrix (CDSM) method was confirmed through comparison with both numerical and experimental results.When compared to the conventional finite element method (FEM), the CDSM frequencies were found to best match with those obtained from a cubic beam finite element (B33 in ABAQUS).Moreover, in comparison with experimental data, the CDSM results were found to be within acceptable accuracy; the average error for the first three bending natural frequencies is less than 8%.
It is worth noting that the proposed CDSM can also be fine-tuned to match spindle's experimental fundamental frequency, instead of the nominal value.Further research is underway to improve the CDSM by using more accurate material properties for different segments of the spindle.The developed calibrated model could then be exploited in an attempt to take into account the effects where

Figure 3 :
Figure 3: Degrees of freedom of the system.
frequency (Hz) Bearing equivalent spring constant K (N/m) DSM Simply supported ANSYS Spindle manufacturer's designed natural frequency (1000 Hz)

Figure 5 :
Figure 5: Varying spring constant versus system natural frequency using DSM and FEA.