An efficient method for calculating damped critical speeds of a build-in motorized spindle

Nowadays, high-speed machining becomes more and more popular for its abilities of increasing productivity and reducing manufacturing cost. Many high-speed machine tools include a build-in motorized spindle which is usually one of the most important parts since its dynamic behaviors directly affect the machining accuracy. There are two traditional methods which can be used to study the dynamic behaviors of the motorized spindle: finite element method (FEM) and transfer matrix method (TMM). The FEM can effectively do the dynamic behavior study, but it also can be time-consuming when dealing with a complex three-dimensional system like the build-in motorized spindle [1-5]. On the other hand, the TMM seems more suitable and efficient for a complex rotor-bearing system since it starts the dynamic behavior study at a certain station and proceeds station by station instead of dealing with the system as a whole. The TMM was first developed by Prohl [6]. Since then, a lot of research work has been done to improve its performances and extend the range of its applications. Lund [7] made a significant advancement in the TMM by using complex variables to express the system eigenvalues and developing a more general formulation of bearing forces when accounting system damping. Bansal and Kirk [8] made a further progress on the basis of Lund's work by taking the effects of the damping and flexibilities of the bearings on the system stability under consideration and using Muller's method to find the eigenvalues of rotor-bearing systems. Horner and Pilkey [9] proposed a modified TMM called Riccati TMM by the use of an existing large catalog of transfer matrices for various structural members, and it can eliminate the numerical instability encountered when calculating high resonant frequencies. Later, some similar work was done by Murphy and Vance [10]. In their work, Bairstow's method was used to find the system eigenvalues. Lu [11] developed an improved Riccati TMM for calculating damped critical speeds and predicting stability of rotor-bearing systems more efficiently. In recent years some research work were done to extend the application of the TMM to more complex rotor-bearing systems. Such as, Zu and Ji [12] proposed an improved TMM for steady-state analysis of nonlinear rotor-bearing systems; Meng et al. [13] used the TMM for dynamic characteristics of a gas turbine rotor; Varney and Green [14] applied the TMM to the rotordynamic analysis of elastomer ring supported rotordynamic system; Jiang and Zheng [15] used the traditional TMM to establish the dynamic model of a motorized spindle supported by angular-contact ball bearings and the vibration behavior is studied. Although the TMM has been successfully applied to many complex rotor-bearing systems, it is rarely used to calculate damped critical speeds of a build-in motorized spindle supported by fluid film bearings. In this paper, an efficient method based on an improved Riccati TMM is suggested for calculating damped critical speeds of a build-in motorized spindle supported by fluid film bearings. Based on the framework of the Riccati TMM, the state vector of a typical station is extended to ten dimensions by taking the electromagnetic force and torque in the motorized spindle under consideration. Then the dynamic model of the spindle is established and the characteristic polynomial is obtained. To calculate damped critical speeds is to find the roots of the characteristic polynomial. In order to improve efficiency, instead of finding all roots of the characteristic polynomial, only the roots in the effective resonance region are found since it is the most cared region for predicting stability. The argument principle is adopted to count the number of roots in the effective resonance region. The micro genetic algorithm is used iteratively to find those roots. At the end of each iteration, one root has been found and the synthetic division is used to reduce the order of characteristic polynomial. The iteration terminates after all roots in the effective resonance region have been found. Finally, this method is successfully applied to calculate the damped critical speeds of a typical build-in motorized spindle.


Introduction
Nowadays, high-speed machining becomes more and more popular for its abilities of increasing productivity and reducing manufacturing cost.Many high-speed machine tools include a build-in motorized spindle which is usually one of the most important parts since its dynamic behaviors directly affect the machining accuracy.There are two traditional methods which can be used to study the dynamic behaviors of the motorized spindle: finite element method (FEM) and transfer matrix method (TMM).The FEM can effectively do the dynamic behavior study, but it also can be time-consuming when dealing with a complex three-dimensional system like the build-in motorized spindle [1][2][3][4][5].On the other hand, the TMM seems more suitable and efficient for a complex rotor-bearing system since it starts the dynamic behavior study at a certain station and proceeds station by station instead of dealing with the system as a whole.
The TMM was first developed by Prohl [6].Since then, a lot of research work has been done to improve its performances and extend the range of its applications.Lund [7] made a significant advancement in the TMM by using complex variables to express the system eigenvalues and developing a more general formulation of bearing forces when accounting system damping.Bansal and Kirk [8] made a further progress on the basis of Lund's work by taking the effects of the damping and flexibilities of the bearings on the system stability under consideration and using Muller's method to find the eigenvalues of rotor-bearing systems.Horner and Pilkey [9] proposed a modified TMM called Riccati TMM by the use of an existing large catalog of transfer matrices for various structural members, and it can eliminate the numerical instability encountered when calculating high resonant frequencies.Later, some similar work was done by Murphy and Vance [10].In their work, Bairstow's method was used to find the system eigenvalues.Lu [11] developed an improved Riccati TMM for calculating damped critical speeds and predicting stability of rotor-bearing systems more efficiently.In recent years some research work were done to extend the application of the TMM to more complex rotor-bearing systems.Such as, Zu and Ji [12] proposed an improved TMM for steady-state analysis of nonlinear rotor-bearing systems; Meng et al. [13] used the TMM for dynamic characteristics of a gas turbine rotor; Varney and Green [14] applied the TMM to the rotordynamic analysis of elastomer ring supported rotordynamic system; Jiang and Zheng [15] used the traditional TMM to establish the dynamic model of a motorized spindle supported by angular-contact ball bearings and the vibration behavior is studied.Although the TMM has been successfully applied to many complex rotor-bearing systems, it is rarely used to calculate damped critical speeds of a build-in motorized spindle supported by fluid film bearings.
In this paper, an efficient method based on an improved Riccati TMM is suggested for calculating damped critical speeds of a build-in motorized spindle supported by fluid film bearings.Based on the framework of the Riccati TMM, the state vector of a typical station is extended to ten dimensions by taking the electromagnetic force and torque in the motorized spindle under consideration.Then the dynamic model of the spindle is established and the characteristic polynomial is obtained.To calculate damped critical speeds is to find the roots of the characteristic polynomial.In order to improve efficiency, instead of finding all roots of the characteristic polynomial, only the roots in the effective resonance region are found since it is the most cared region for predicting stability.The argument principle is adopted to count the number of roots in the effective resonance region.The micro genetic algorithm is used iteratively to find those roots.At the end of each iteration, one root has been found and the synthetic division is used to reduce the order of characteristic polynomial.The iteration terminates after all roots in the effective resonance region have been found.Finally, this method is successfully applied to calculate the damped critical speeds of a typical build-in motorized spindle.

Dynamic modeling by the improved Riccati TMM
A build-in motorized spindle is driven by a built-in motor directly, which eliminates any transmission components such as belts and gears.In order to calculate its damped critical speeds, the dynamic model should be constructed first.In the following part, the dynamic model of a build-in motorized spindle is established based on the framework of the improved Riccati TMM proposed by Lu [11], where the build-in motorized spindle is treated as a rotor-bearing system with multiple disks and two bearings.
Fig. 1 The lumped mass model of the rotor-bearing system with multiple disks and two bearings Fig. 2 The simplified model of the bearings Fig. 1 shows the lumped mass model of this rotorbearing system where the rotor is discretized to N nodes with lumped mass and N-1 elastic shaft segments.The two anisotropic bearings of the rotor-bearing system are represented by eight coefficients as shown in Fig. 2.
For a rotor-bearing system supported by two anisotropic bearings, the state vector of station i is 8-dimensional which consists of force vector {f}i and displacement vector {e}i [16].In the mean time, the electromagnetic force and torque of the build-in motorized spindle should not be ignored.Therefore, the state vector of station i is extended to 10-dimensional, which can be expressed in the XY coordinate axis as follows: Through force analysis of station i, transitive relation of state vector can be built as: where superscripts R and L present the right and left side of the segment respectively, I is the identity matrix, {W}i, {U}i and {R}i are partitioned matrices which can be defined as follows:   where m is the mass of each station, c and K are the dynamic coefficients of the bearings as shown in Fig. 2, Jp is the polar moment of inertia, Jd is the transverse moment of inertia, ω is the rotating speed, M and Q are the electromagnetic torque and force respectively, and Si   is the complex frequency where Ω and λ are the critical speeds and system damping.
  where l is the length of each shaft segment, J is the transverse area moment of inertia of the shaft, E is Young's modulus, and where ks is the transverse area coefficient, which is 0.886 for a Solid circular shaft and 2/3 for a hollow one, G is the shear modulus, and A is the cross-section area of shaft.Expending Eq. (2-3), the following equation can be obtained: If introducing a transform matrix B, then the relationship between displacement vector and force vector on each section can be expressed as: (5) Substituting Eq. ( 5) into Eq.( 4), the following equation can be obtained: .
It will be difficult to use the recurrence relation expressed by Eq. ( 6) to obtain the system characteristic polynomial since it contains an inverse matrix.According to Lu [11], some transformations can be performed to make it easier.First, set

C T T T and
, then substitute them into Eq.( 6) and the following recurrence equation can be obtained: .
According to the left end boundary conditions, The roots of Eq. ( 8) are the complex eigenvalues whose imaginary part is the damped critical speed.

The procedure of damped critical speeds calculation
In this section, the procedure of a damped critical speeds calculation method for a build-in motorized spindle is presented based on the spindle system characteristic polynomial.It is an efficient method which mainly contains three steps.The first step is to retrieve the number of roots of the characteristic polynomial in effective resonance region by argument principle.The second one is to find these roots one by one using the micro genetic algorithm.The third one is to reduce the order of the characteristic polynomial by synthetic division when one root has been found.After all roots in effective resonance region have been found, the damped critical speeds can be obtained since they are the imaginary parts of these roots.In the following parts, all of these steps are discussed in details one by one, and then the procedure of the suggested method is presented.

Retrieving the number of roots in effective resonance region by argument principle
Finding all roots of the characteristic polynomial could be very time-consuming.In the mean time, it may not be necessary for predicting stability in engineering problems.Therefore, in order to improve the efficiency, only the roots in effective resonance region are found.The argument principle is used to retrieve the number of those roots.
The argument principle is a successful method in many science and engineering fields, including in complex analysis [17].For a complex function f(x) such as the characteristic polynomial described in Eq. ( 8), if let it be meromorphic and C be a simple closed contour which does not pass through any poles or zeros of f(x), then according to the argument principle, the winding number is equal to the difference between the number of zeros and poles of f(x) that are inside the contour.Since the winding number of f(x) can be found by moving around to a given circle and in the case of finite length sequences the poles are only trivial ones, the number of zeros of f(x) inside the circle can be retrieved.
In this case, the number of roots of Eq. ( 8) in the effective resonance region equals to the winding number of the determinant value of the characteristic polynomial on the right hand of the original point.

Finding roots by the micro genetic algorithm
The next step is to find the roots of the characteristic polynomial in effective resonance region.There are many methods which have been used to find the roots.The Newton-Raphson approach was the first one which was tried.But it has many difficulties when the rotor-bearing system becomes more complicated.After that more efficient methods have been used, such as Bairstow method, Bairstow-Newton method and Muller's method.However they still have some difficulties when dealing with a complex rotor-bearing system, such as being sensitive to the initial points and finding the roots disorderly.So all roots need to be found when predicting stability, which will be very timeconsuming.As a kind of global optimization methods which are not sensitive to the initial points, genetic algorithms (GAs) are widely used in many fields, especially in function optimizations.Finding the roots of the characteristic polynomial that described by Eq. ( 8) can be transformed into an optimization problem of finding the minimum value of the following function In this paper, the μGA developed by Krishnakumar [18] is used to find the minimum point of function (9), which is an improvement of traditional GAs.It uses a very small population (typically 5～8) and a similar evolutionary strategy of traditional GAs.In order to avoid the generation converging fast to a local optimum, a restart strategy is introduced to preserve the genetic diversity.It performs better than the traditional GAs.

Reducing the order of characteristic polynomial by synthetic division
The next step is to reduce the order of the characteristic polynomial after one root has been found, which can make it easier to find the next root.In this paper, synthetic division is used to reduce the order of Eq. ( 8), which can be described as follows: where Si is the current root that has been found and i S  is its conjugate root, 1, 2, , im  ; m is the number of roots that need to be found.

Procedure of the present method
The procedure of the present method is illustrated in Fig. 3. Fig. 3 The procedure of the present method

Numerical example
In this section, the present method is applied to calculate the damped critical speeds of a typical build-in motorized spindle which is supported by two identical fluid film cylindrical bearings.The operating speed of the spindle is 1000 rad/s, and the motor power is 35.3 kW.According to Zhong et al. [1], if the spindle is operated under the first order critical speed nc1, its operating speed should not be higher than 0.75 × nc1; otherwise, it should be between 1.4 × nc1 and 0.7 × nc2, where nc2 is the second order critical speed.Therefore, the upper bound of the most cared range of the critical speed could be 2000 rad/s.That is, the effective resonance region of this spindle is [0, 2000] rad/s.Fig. 4 shows the dynamic model of the motorized spindle which is treated as a rotor-bearing system with multiple rigid disks and bearings.The rotor is divided into six shaft sections, which has seven nodes with lumped mass and six elastic shaft segments.Parameters of the shaft sections are given in Table 1.The grinding wheel, the motor rotor and the baffle ring are simplified as rigid disks 1, 2 and 3 respectively.Parameters of the rigid disks are summarized in Table 2.
Parameters of the bearings are: the clearance-to-radius ratio ψ = 2.0‰, the length-to-diameter ratio L/D = 0.5 and oil viscosity η = 0.0184 N s/m 2 .The stiffness and damping coefficients with seven different eccentricity ratios are calculated by the method developed by Xiong et al. [19], as shown in Table 3.By using the improved Riccati TMM introduced in Section 2, the characteristic polynomial equation can be obtained.The number of roots in effective resonance region retrieved by argument principle is 4.Then, the μGA is used iteratively to find these four roots.The parameters used by μGA are: size of the population N = 5, the crossover probability pc = 0.9, the mutation probability pm = 0.05, the convergence criterion is . For comparing, Muller's method is also used to find the roots.4 shows the roots found by Muller's method and the present method, between which the differences are very small.In order to obtain the four roots in effective resonance region, Muller's method need to find all roots of the characteristic polynomial and sort them from small to large by the imaginary part, which costs 620 seconds.However, the present method can find four roots directly without finding all roots, which only costs 23 seconds.The fifth root of the sorted roots found by Muller's method is also shown in Table 4, the imaginary part of which is beyond the effective resonance region [0,2000].It verifies that the number of roots in the effective resonance region has been retrieved correctly.
Once the roots of the characteristic polynomial have been found, the damped critical speeds can be obtained.Since the operating speed of the spindle is much larger than the fourth order of the critical speed and much smaller than its fifth order, the spindle can be considered stably at operating speed.The calculation results also can be used for further rotordynamic stability analysis.

Conclusions
In this paper, an efficient method based on an improved Riccati TMM is presented to calculate the damped critical speeds of the build-in motorized spindle.By the using of the improved Riccati TMM, the dynamic model of the build-in motorized spindle treated as a rotor-bearing system with multiple rigid disks and bearings is established and the characteristic polynomial is obtained.Instead of finding all roots of the characteristic polynomial, the suggested method can only find the roots in the effective resonance region, the number of which is retrieved by argument principle.Then the μGA is used to find those roots iteratively.Once a root has been found, the synthetic division is adopted to reduce the order of the characteristic polynomial.Unlike most of the calculation methods it can find the roots in order of their imaginary parts from small to large.Therefore the present method is able to save a lot of calculation time.The simulation results of the damped critical speeds calculation of a typical build-in motorized spindle demonstrate that the present method can efficiently calculate damped critical speeds of a complex rotor-bearing system.
shearing force; x, y are the deflection in XOZ and YOZ-plane; x  , y  are the slope in XOZ and YOZ-plane.

Table 1
Parameters of the shaft segments of the dynamic model The dynamic of the typical build-in motorized spindle

Table 2
Parameters of the rigid disks of the dynamic model

Table 3
The bearing's direct non-dimensional film stiffness and damping coefficients with different eccentricity ratios