Identification of a Machine Tool Spindle Critical Frequency through Modal and Imbalance Response Analysis

Corresponding Author: Khairul Jauhari, Department of Mechanical Engineering, University of Diponegoro, Semarang, Indonesia Email: khairul.j4uhari@gmail.com Abstract: The most important part for machine tool system is spindle system. The structural properties of spindle will affect to the productivity of machining process and the work pieces quality. Thus it is necessary and very important to know its spindle dynamic behaviors for avoiding forced vibration due to resonance. The Finite Element Method (FEM) has been adopted for obtaining spindle dynamic behaviors. FEM model and analysis of dynamic behaviors for the spindle system are similar to those developed in rotor-dynamic. The main purpose of this study, we implemented an ANSYS Parametric Design Language (APDL) program based on finite element method for obtaining full analysis of rotor dynamic in order to investigate the spindle dynamic behaviors. The programs efficiently performed the full analysis by determining the Campbell diagrams, critical speeds and response of imbalance due to mass unbalance at the cutting tool. Results show that the critical speeds were calculated previously are far enough from the spindle speeds operating range, then the spindle would not experience resonance and the maximum response of imbalance occurs on operating speed is also in acceptable limit. ANSYS Parametric Design Language (APDL) can be used by spindle designer as tools in order to increase the product quality, reducing cost and time consuming in the design and development stages.


Introduction
The most important components for machine tool system is the system of spindle. The dynamic properties of spindle directly affect to the machining productivity and quality of the product. In the design and development stages, then it is very necessary and important to know the behaviors of spindle dynamics for avoiding resonance of critical speeds due to machining operations. To obtain dynamic analysis of spindle system analytically in the early design stage, the Finite Element Method (FEM) has been frequently adopted for modeling of spindle dynamic behaviors. Basically, the FEM model and dynamic analysis for the spindle systems of machine tools is similiar to those developed in rotordynamic. However, the spindle shafts used in machine tools usually have smaller shaft diameters and bearings and possess disk-like in turbomachine components. Thus in this study, we attempt to review some researches relating to the field of rotordynamics. Lin et al. (2013) in their paper stated that the most popular approch for modeling the spindle dynamic behavior is the Finite Element Model (FEM), because of its capability to manage complex geometry and boundary condition and the calculation approaches save time and money while solving the finite element system equation. Lin (2014) developed a Genetic Algorithm (GA) optimization approach to search the optimal location of bearings on the motorized spindle shaft. The goal is to maximize its First-Mode Natural Frequency (FMNF). In order to achieve the results, a spindlebearing system dynamic model is formulated using Finite Element Method (FEM) that was developed in the rotordynamics. Nelson and McVaugh (1976;Nelson, 1980) applied the theory of Timoshenko's beam to build matrix of systems in order to analyze the rotor systems dynamics including the influences of gyroscopic moments, rotational inertia, axial load and shear deformation. Zorzi and Nelson (1977) presented the influences of damping on the rotating systems dynamics. Cao and Altintas (2004) proposed a general method that can be used for modeling the spindle-bearing system. The spindle and spindle housing are considered as element of Timoshenko's beam including the effects of gyroscopic moment and centrifugal force. The stiffness matrix of the bearing, the contact angle, preload and deflection of spindle and spindle housing are all coupled using finite element model of the spindle assembly. Erturk et al. (2006) proposed method of analytical about using receptance coupling and modification of structural in the modeling of spindle-holder-tool assembly. These component assemblies are considered as multi-segment model of Timoshenko's beam and Euler-Bernoulli's beam and they compared its results with those formulation method. They found that the accuracy of Euler-Bernoulli model may give the low accuracy results at higher frequencies. Chatelet et al. (2005) presented a modeling approach of analysis modal method to calculate the rotor system dynamics (the frequencies of natural and mode shapes) of turbo-machinery. Two methods were compared using the model approaches based on finite element model of three-dimensional and one-dimensional. The results show that low accuracy was obtained when the dynamic behavior of system was formulated using one-dimensional beam. Whalley and Abdul-Ameer (2009), by using simple harmonic response methods, they used contoured shaft profiles in order to obtain the natural frequencies and critical speeds of rotor shaft systems. Taplak and Parlak (2012) used the Dynrot program based on finite element method for performing the rotordynamics analysis of gas turbine with certain mechanical and geometrical properties. This program was used to obtain Campbell diagram, determining critical speeds and investigating imbalance response of the rotor due to mass unbalance of compressor. Jalali et al. (2014) investigated full analysis of rotordynamics on high speed turbine-impeller system using finite element model of ANSYS program and onedimensional beam and experiment of modal testing as performed by (Taplak and Parlak, 2012), the diagram of Campbell, critical frequencies and imbalance response analysis due to mass unbalance were obtained using the two models of ANSYS and one-dimensional finite element. Based on the comparing results of theoretical and experiment show that the good accuracy of the behavior of such systems can be achieved using two of these models. Villa et al. (2008) have been investigated the behavior of non linear flexible unbalanced rotor which supported by roller bearings. In this case, the method of harmonic balance was used to find a periodic response of this non-linear system. In frequency form, stability of the system was identified based on a perturbation applied method. They stated that the method of harmonic balance with the AFT application can be used to obtain harmonic solutions. Bai et al. (2012) investigated the dynamic behaviors of a hydroturbine main shaft system by using program named ANSYS. They developed ANSYS Parametric Design Language to generate the geometry model of 3D, analysis of modal and obtaining critical speed at the spin speed. By using ANSYS Parametric Design Language, critical speed determination and analysis of imbalance response for a multi segmen rotor have been presented by Gurudatt et al. (2010;Jagannath, 2012). They showed the advantage of using this method is that with typing one of input script as example shaft diameter, rotor segmen length, loads experienced by the rotor, command of "ANTYPE, MODAL" and "HARMIC" command, all of these commands can be read as variable input and execution command of the program. These scripts are executed by ANSYS Programming Design Language. Results of these programming language are validated with results of theoretical and measurement, which are good agreement of the acceptable limits.
The main purpose of this study is implementation of an ANSYS Parametric Design Language (APDL) program based on finite element method for obtaining full analysis of rotor dynamic in order to investigate the spindle dynamic behaviors. The programs efficiently performed the full analysis by determining the Campbell diagrams, critical speeds and response of imbalance due to mass unbalance at the cutting tool. The imbalance response is carried out to verify the critical speeds obtained from the modal analysis and practically investigating its behavior. In this study, a grinding spindle system with certain mechanical and geometrical properties was used as the analysis model. The results show that the critical speeds were calculated previously are far enough from the spindle speeds operating range, then the spindle would not experience resonance and the maximum response of imbalance occurs on operating speed is also in acceptable limit.

Theoretical Formulation
In this study, a Leadwell STD V-30 spindle system of grinding machine tool is shown in Fig. 1. The spindle is designed to operate at up 8,000 rev/min with a 5.6 kW motor connected to the shaft with a pulley-belt system. In this model, cutting tool, tool-holder, spindle shaft and bearings were included. All components of cutting tool, holder, spindle-bearings assembly are modeled as elements of discrete disk and bearings and the multisegment beams with elasticity and distributed mass. Based on procedure of the finite element discretization in many literatures (Lalanne and Ferraris, 1998;Friswell et al., 2001;Yamamoto and Ishida, 2001), the detail of equations will not be derived here and only the equation of motions are shown below.
Here, q e is the nodal displacement vector, containing the eight degrees-of-freedom for one of shaft element (two rotations and two translations in each node). By combining the individual matrices of each shaft element, one can obtain the global matrices that represent the whole shaft, thus resulting to the following equation of motion: Here, (M G T ) is a matrix of global translational mass, (M G R ) is a matrix of global rotational mass, (K G ) is a matrix of global stiffness, (G G ) is a matrix of global gyroscopic, (F G ) is a global force vector matrix acting on the shaft element and (q) is the displacement vector containing all 4(n e + 1) degrees of freedom of the shaft elements that represent the physical shaft (n e is the number of shaft elements).
Mass element is considered as a rigid disk. The rigid disk is placed at a certain nodal point of finite element. Here, q d is matrix of the nodal displacement vector of the disk center mass. For assuming a constant operating speed Ω then the equation can be expressed as: Where: Coefficients of stiffness and damping represent the dynamic characteristics of bearing. The motion equation for dynamic characteristics of the bearings on shaft element can be written as: br br br br br Here, (K br ) and (C br ) are called matrices the bearing stiffness and the damping matrices and (F br ) is called a bearing force matrix.

Equation of Global System and Analysis of Eigenvalue
Based on the element Equation 2, 3 and 4 then a certain global element equation can be established and other global equations also can be generated for the other elements. These elements are constructed to form the general equation, which represents behavior of the whole system. Then here, the motion equation of the damped system for coordinate of global is expressed as: Here, C G = C br -ΩG e -ΩG d , K G = K e + K br . In order to obtain the natural frequency of system, then eigenvalue must be solved and expressed by Equation 5, the system equation can be written in state space vector form: where, the matrices of A s , B s and displacement x consist of element matrices given as: For assuming harmonic solution x= x 0 e λt of Equation 6, the solution of an eigen-value problem is: For the nontrivial solution, the determinant of Equation 7 equal to zero: where, C a = A s -1 B s and λ is an eigen-value. The eigenvalues are usually as the complex number and conjugate roots: Here, α j and ω j are factor of growth and the damped natural frequencies of j th mode, respectively.

Response of Imbalance
The forces of mass unbalance (F) which is shown in Equation 5 can be expressed as: The steady state response of mass unbalance is considered to be as the form: Substituting Equation 10 and 11 into (5), the equation can be expressed as: By solving equation (12), the response of steady state can be obtained.

ANSYS Parametric Design Language (APDL)
In this study, a practical APDL macro scripting language has been developed to generate all the required results, containing amplitude plots and frequencies plots at all the nodes of the model, with minimal effort of the user. The algorithm incorporated in the macro is as: • Setup the model. Impose boundary conditions and apply excitation force • Performing the analysis of modal for obtaining the natural frequencies and the critical speeds. Set solution using "ANTYPE,MODAL" command. Retrieve mode frequency and critical speed frequency using '*GET' command and store in using '*VFILL' command • Perform harmonic analysis for obtaining imbalance response and provides validation for the frequency found by modal analysis through harmonic analysis. Set solution using "HARMIC" command and set the range of excitation frequencies to increment from 0 to maximum operating speed in number of step (using " NSUBST" command) • Solve for unbalance response. Plot results to get unbalance response at nodal point 'n' • Increment parameter 'n' by 1. If n > 18 (since the spindle-bearing system model here contains 18 nodal points), if ok then go to next step. Otherwise, go back to step 3 • End of program Table 1 shows the mechanical properties and geometric of the element. In this study, the spindle shaft is modeled into 17 beam elements with nodal points at each end of the elements. Two-mass the cutting tool and pulley-belt components can be considered as 1 and 3 elements of rigid disk, respectively. These elements of rigid disk are located at the nodal points number 1, 14, 15 and 16. The parameters of the mass element are tabulated in Table 2. In addition the two set of bearings, located at the nodal points 5 and 12. These bearings are modeled as symmetric isotropic bearings and stiff elastic constrains. Table 3 show model of bearing elements. A schematic of spindle's finite element model is shown in Fig. 2.  As can be shown from Fig. 3, a three dimensional geometric model of the spindle shaft is established using the APDL (ANSYS Parametric Design Language) program. The spindle shaft is considered as the elements of BEAM188 with an internal node and the function of quadratic shape to increase the element accuracy. The characteristic of BEAM188 has two nodal points and having twelve degree of freedom at each element, the motions are translation in the x, y and z axis direction and rotation about x, y and z axis. The element of MASS21 is used for modeling of disk element (mass of rigid disk) and element of COMBIN14 is used for modeling the symmetry bearings. The nodal points, elements, material properties, real constants, boundary conditions and other physical system-defining features that constitute the model have been created by using APDL commands such as RO, PEX, PGXY, MP, ET, MAT, K, N, LSTR, R, RMORE, LATT, LESIZE and E. The effect of shear cannot be ignored in the spindle shaft. The constraints are applied to the element motions of displacement in the x axis direction and rotation about x axis, thus the spindle shaft would not experience any displacements of translation and twist motion about the x axis direction. Parameters for the material and element properties of this spindle shaft model are the same as in beam finite element model.

The Model of Finite Element
By QRDAMP method, a modal analysis on spindle shaft system was performed to obtain speeds of the critical whirl and value of Campbell frequencies. In a stationary frame, "Coriolis Effect" can be included in the rotating structure by activating the CORIOLIS command. The intersection between natural frequencies and spin speeds (speeds of the critical whirl) for gradient 1x is determined. Harmonic analysis also performed with SYNCHRO command to determine amplitude response values.

Results
In this study, the spindle dynamic behaviors at operating speed, the Campbell diagram, critical speed and imbalance response were obtained. Analyzes of the numeric simulation were performed with considering the speed ranges from 0 to 27000 rpm. Table 4 shows the damped natural frequencies of the spindle assembly (at operating speed 8000 rpm) obtained by the APDL FE model. Figure 6 shows the Campbell diagrams obtained by the Finite Element (FE) model constructed in APDL program. In Fig. 4 and 5, mode shapes of the natural frequencies were determined at operating speed (8000 rpm).
As can be seen from Table 5, These critical speeds are obtained from the APDL finite element models which intersecting the frequency modes with spin speed line. There are three Backward Whirls (BW) and two Forward Whirls (FW) modes were considered. They are 1st backward whirl, 1st forward whirl, 2nd backward whirl, 2nd forward whirl and 3rd backward whirl, respectively critical speeds. We also obtained the deflection of mode shapes relating with these five critical speeds. The shape of deflection at critical speed relating to 1st backward whirl, 1st forward whirl, 2nd backward whirl, 2nd forward whirl and 3rd backward whirl are shown in Fig. 7-11.
Analysis of Imbalance response is carried out for investigating the behaviors of spindle dynamic which provides validation for the frequency found by modal analysis through harmonic analysis. An unbalance of 9.981×10 −5 kg.m for center mass of the cutting tool is considered in APDL model. The nodal solutions of imbalance responses have been obtained using the APDL FE model which are tabulated in Table 6. As can be seen from the table, the maximum amplitudes occur near at the 1st forward whirling and 2nd forward whirling critical speeds points which were calculated in the Table 5. It's mean that if the system has damping, the system will resonance when it approaches at these critical points. It is looked at the Table 7, the maximum amplitudes are not coincide 195.9 Hz and 350.3 Hz which were calculated previously in modal analysis. The percentage difference between both model analysis are in very good agreement and the maximum difference is about 0.56%.   Based on this imbalance response analysis, it is easy to understand that the critical speed of the spindle is the speed corresponding to the intersection of the natural frequency (Hz) equal to spindle spin's (rpm) line with only the forward whirl mode.  Figure 12 shows the response of mass imbalance that was evaluated at nodal points 1,5 and 12, which stated disk 1 (cutting tool), bearing set 1 and bearing set 2, respectively. Figure 13 and 15 show the mode shapes at two first critical speeds and Fig. 14 and 16 show the stress contour for two first critical speeds.

Conclusion
A program named ANSYS Parametric Design Language (APDL) which based on finite element method has been implemented to investigate full analysis of spindle dynamic and evaluating the results. A grinding spindle system with certain mechanical properties and geometrical is modeled as beam APDL model technique. The deflection shapes of the spindle at operating speeds, the Campbell diagrams, critical speeds and imbalance response are obtained. Based on this imbalance response analysis that the critical speeds of the spindle are the 1st forward whirling and 2nd forward whirling, which the speed corresponding to the intersection of the natural frequency (Hz) equal to spindle spin's (rpm) line with only the forward whirl mode. These critical speeds are still far from the operating speed range of the spindle, thus, the spindle would not experience resonance and the maximum imbalance response at operating speed is still with acceptable limit. Thus, APDL program can be used by spindle designer as tools in order to increase the product quality, reducing cost and time consuming in the design and development stages.

Ethics
This article is original and contains unpublished material. The corresponding author confirms that all of the other authors have read and approved the manuscript and no ethical issues involved.