A Numerical Analysis for Ball End Milling Due to Coupling Effects of a Flexible Rotor-Bearing System Using GPEM

: In this paper, the tool-tip responses for ball end milling, due to the coupling effects of a ﬂexible rotor-bearing system, are investigated numerically. The milling machine tool spindle is modelled as the ﬂexible rotor-bearing system. The critical speeds, natural modes, and unbalance responses of the system are calculated by applying the generalized polynomial expansion method. This generalized polynomial expansion method expresses the displacement as a series formed by the product of generalized coordinates and axial coordinate polynomials. According to the dynamic cutting force obtained by some scholars in the past, combined with the characteristics of the ﬂexible rotor, the dynamic response of the tool-tip for ball end milling is numerically analyzed. The responses, including time histories, orbits, and FFT diagrams, are plotted to analyze the dynamic behaviors of the tool-tip. The coupling effects of the ﬂexible rotor-bearing system on the system for ball end milling are ﬁrst studied using the generalized polynomial expansion method. Unlike previous studies, the natural frequency varies with spindle speed and which of the different modes are included in the tool-tip response depends mainly on the spindle speed. Thanks to the gyroscopic effect, the critical speeds and responses of tool-tips can be discussed with respect to various spindle speed and tool ﬂutes. The natural modes are accurately determined, and will excite critical speeds for certain modes, including forward and backward modes, thereby signiﬁcantly affecting tool-tip response. In addition, the cutting force component associated with the tool-tip response affects the rotor-bearing system parameters, complicating the issue. Milling at higher spindle speed (2160–19,950 rpm), an important new result is found that the tool-tip oscillates with the cutting-force frequency, accompanied by a longer period vibration of the ﬁrst backward mode of the rotor-bearing system. It can also be seen from the frequency spectrum analysis that, as the spindle speed increases, the peak amplitude of the ﬁrst backward mode becomes larger. Milling at lower spindle speed (960, 1320 rpm), the in-plane vibration trajectory of the tool-tip gradually expands outwards clockwise around the origin until a stable loop is reached. This is because only the ﬁrst backward mode of the rotor-bearing system is excited. Considering the coupling effect of the rotor-bearing system to perform the vibration analysis of the milling machine system, the parameters of the system can be designed or the spindle speed can be selected to avoid severe vibration during machining.


Introduction
The process of ball end milling is widely used in modern manufacturing. In many industrial manufacturing processes, large amounts of material must be removed during the machining process. High material removal rates and high-quality machined surfaces become the issues of most concern in the milling process. To generate faster spindle speeds and lower cutting forces, the design trend of rotor systems is to improve their vibration behavior. The main characteristics of the rotor system are critical speed and unbalance response. To analyze the vibration of a milling machine more accurately and efficiently, an industrially used tool spindle can be modeled as a rotor-bearing system.
It is natural to consider that the cutting forces and structural vibrations will be affected by the variation of the depth of cut in machining. In addition, the chatter usually occurs in the interaction of cutting forces and structural vibrations and will severely damage the surface finish as well as the dimensional accuracy of the workpiece. There have been many studies on cutting force models and the models have been successfully analyzed to determine the critical speed by using a variety of methods. However, regarding the cutting force models, there are few studies considering flexible rotor-bearing systems numerically. In this paper, the cutting force model for ball end milling, coupled with the flexible rotorbearing system, is proposed. Furthermore, the system equations are derived by using the generalized polynomial expansion method. A literature review of rotor-bearing dynamics and related studies on ball end milling is given below.

Three Principal Methods for Rotor System
Many researchers have successfully used various methods for determining the critical speeds, modes, and unbalance responses of the rotor-bearing systems. Three principal methods are reviewed. They are the transfer matrix method (TMM), the finite element method (FEM), and the generalized polynomial expansion method (GPEM). These three approaches have been successfully applied to rotordynamics.

Transfer Matrix Method
Lund and Orcutt [1] originally proposed the TMM to calculate critical speed in the linear frequency domain and extended it [2] to determine the unbalance response, the modal response, the stability, and the damped critical speeds. Moreover, the TMM has been applied by Chu and Pilkey [3] to deal with transient analysis in the time domain using the discrete time technique. Gu [4] proposed an improved TMM to study the rotor behaviors.

Finite Element Method
In the published literature [5], the FEM was probably the most popular rotordynamics analysis method. A Rayleigh beam, model that takes translational inertia and bending stiffness into account, was proposed by Ruhl and Booker [6]. Moreover, a series of finite element models was also proposed by Nelson and Mc Vaugh [7], with the effect of different parameters including axial load, rotator inertia, and gyroscopic moments. Subbiah et al. [8] applied a method combining TMM and FEM to solve the transient responses of rotor systems.

Generalized Polynomial Expansion Method
Shiau and Hwang [9,10] proposed the GPEM, which provides a new approach to the dynamic behavior of rotors. This method was first used to represent the response of a flexible rotor-bearing system. They also used the properties of the Rayleigh quotient to calculate the critical speed of the system. Shiau et al. [11] reported on the use of GPEM to study the nonlinear problem. The stability of the steady-state responses of a rotor supported by the nonlinear squeeze film dampers was investigated. Shiau et al. [12,13] used the global assumed mode method to analyze the ball screw system under a moving skew load. The moving skew load consists of an axial driving force and lateral pay load. The transient response is also analyzed using the Runge-Kutta method. Furthermore, the Floquet theorem is employed to determine the system stability.
TMM is not suitable for analyzing nonlinear problems in the frequency domain. FEM requires a lot of computation time, although it can be an effective tool for both linear and nonlinear systems. Compared with other methods, GPEM shows better advantages in Appl. Sci. 2023, 13, 7252 4 of 21 system. The main research object of this study was the parameters of rolling bearings, such as radial clearance, Hertzian contact, and time-varying stiffness.
Li et al. [25] proposed a convenient and non-invasive method for monitoring milling forces. Based on spindle real-time vibration, the milling forces were predicted. This research provided a simple and accurate method for obtaining the milling force of the machine tool, which was of great value to the monitoring of the machining process. A twodegree-of-freedom dynamic model of a spindle in a machine tool is proposed to describe the relationship between the spindle vibration and milling force. The gyroscopic effect was not considered. Dai et al. [26] considered runout to study the stability of a five-axis system of ball end milling under low radial immersion. The system was derived into the state transition matrix by the proposed generalized precise integration method. An experiment was also set up to validate the accuracy of the proposed model. The gyroscopic matrix was not found. Qin et al. [27] designed an algorithm to predict the milling force to obtain the boundary of the cutter workpiece engagement. They developed an improved Z-map method. The effectiveness of the proposed model was verified by the experiment. The spindle was not considered in the analysis. Wojciechowski et al. [28] evaluated the ploughing phenomenon by investigating the ploughing forces at the workpiece interface on the side of the tool during precise ball end milling. The original ploughing force model during ball end milling was developed, involving the effects of the ploughing volume and the minimum thickness of the uncut chip. This study focused on the evaluation of the ploughing phenomenon during precise ball end milling with various machined surface inclinations. The spindle was not considered in this model. The micro-milling method, assisted by a three-dimensional vibration, was proposed by Lv et al. [29]. The resulting shape of the structures was predicted by numerical software and compared with experimental data. A model for predicting surface generation was proposed considering the orthogonal spiral and multi-body kinematics theories.

The Purpose and Novelty of This Work
Several past papers on ball end mills have presented a vibration analysis of the tool-tip during machining. However, some scholars did not include the influence of the entire spindle system, and some scholars regarded the spindle as a rigid rotor for analysis. Although some researchers discretized the spindle into a multi-degree-of-freedom rotor system for research, the gyroscopic effect caused by the rotation of the flexible rotor section was not considered when analyzing the vibration of the ball end mill. Therefore, the natural frequencies of the spindle system were independent of the spindle speed. In practical applications, the operation of the spindle will affect the vibration behavior of the tool during machining process, thereby reducing the machining efficiency.
The aim of this study is to investigate the tool-tip vibration behavior in ball end milling in combination with the coupling effects of a flexible rotor-bearing system. This GPEM is applied to approximate the displacements of the entire flexible shaft using a polynomial multiplied by the time-dependent coefficients. Unlike previous studies, the natural frequency varies with spindle speed and which of the different modes are included in the tool-tip response depends mainly on the spindle speed. Thanks to the gyroscopic effect of the flexible rotor, the relationship between the natural frequencies and critical speeds for various spindle speeds and tool flutes can be discussed. Part of the cutting force associated with the tool-tip response changes the rotor system parameters, complicating the problem. However, the tool-tip responses can be calculated more accurately.

Equation Formulation Using GPEM
In this study, the GPEM is used to solve the dynamic behavior of the flexible rotorbearing system proposed by Shiau and Hwang [9,10]. The lateral displacements are assumed to be very small and all the deflections are assumed to be parallel to the x-y plane. In addition, the shear effect, and the torsional and axial vibrations are neglected. Denote the displacements of two translations and two rotations at any cross section of the shaft as (u, v) and (B, Γ). u and B are along the x-axis. v and Γ are along the y-axis. The deflections of the flexible rotor are described as functions of the time t, the shaft length l, and the axial coordinate z, as follows: where a n (t) and b m (t) represent the corresponding generalized coordinates, that is, the coefficients of the polynomial, and N p represents the total number of terms of the polynomial. The kinetic energy of the shaft and disk, and the strain energy of the shaft and bearing are first composed. Using the Lagrangian approach, the system equation of motion can be expressed as: , Ω denotes the constant rotating speed and M, G, C xx , C xy , C yy , K s , K xx , K xy , and K yy are the N p × N p component matrices shown in Appendix A. The first and the second vector terms on the right side of Equation (2) represent the linear and nonlinear forces. It is noted that G is the gyroscopic matrix caused by the rotation of the shaft section. Therefore, the natural frequencies of the system vary with the rotational speed of the spindle. The whirl map, including forward and backward modes, can be established. The system response at any axial position can be solved by using Equations (1) and (2). In this study, the GPEM is employed to analyze the tool-tip vibration for ball end milling combined with the rotor-bearing system. The left side of Equation (2) will include the mass, damping and stiffness effects, and gyroscopic effects of the overall tool spindle system, while the force terms generated by the ball end mill, discussed in the next section, will be placed on the right side.

Dynamic Responses for Ball End Milling
Shaiu et al. [18] presented the vibration behaviors of a rigid rotor-bearing system for ball end milling under dynamic cutting forces. This rotor system assumes the equivalent system parameters of two degrees of freedom. The gyroscopic effect had not been considered since there was no rotation of the rotor cross section. The aim of this paper is to study the vibration responses of the modeled flexible rotor-bearing system for ball end milling by using the GPEM. Figure 1 shows the ball end milling process and the corresponding model. The milling machine tool spindle is modelled as the flexible rotor-bearing system as shown in Figure 1b. In Figure 1b, the first part (I) is the spindle, the second part (II) is the spindle including the simplified tool holder, and the third part (III) is the overhang length of the tool. Appl. Sci. 2023, 13, x FOR PEER REVIEW 6 of 24 sponding model. The milling machine tool spindle is modelled as the flexible rotor-bearing system as shown in Figure 1b. In Figure 1b, the first part (I) is the spindle, the second part (II) is the spindle including the simplified tool holder, and the third part (III) is the overhang length of the tool.

Dynamic Cutting Forces
Lee and Altintas [14] introduced the linear cutting forces used in this study. The geometry of a helical ball end milling tool is given in Figure 2. Furthermore, the cutting forces, chip thickness, and area of cutting segment for ball end milling are plotted in  acting on the cutter along the tangential, radial, and axial directions, respectively, and are given as: where te re ae ,, K K K represent the edge force coefficients along the tangential, radial, and axial directions and s d represents the differential segment length of the curved cutting edge.

Dynamic Cutting Forces
Lee and Altintas [14] introduced the linear cutting forces used in this study. The geometry of a helical ball end milling tool is given in Figure 2. Furthermore, the cutting forces, chip thickness, and area of cutting segment for ball end milling are plotted in Figure 3. dF t , dF r , and dF a are the differential dynamic cutting force components acting on the cutter along the tangential, radial, and axial directions, respectively, and are given as: dF t = K te ds + K tc t(ψ i ,k) db dF r = K re ds + K rc t(ψ i ,k) db dF a = K ae ds + K ac t(ψ i ,k) db where K te , K re , K ae represent the edge force coefficients along the tangential, radial, and axial directions and d s represents the differential segment length of the curved cutting edge. K tc , K rc , K ac represent the cutting force coefficients (or shear force coefficients) along the tangential, radial, and axial l directionst(ψ i ,k) = h(ψ i ) sink denotes the uncut chip thickness normal to cutting edge, in which ψ i denotes the instantaneous angular immersion in the global coordinate system, andk represents the angle in a vertical plane between the z-axis and a point on the flute, h(ψ i ) is the thickness of the instantaneous chip. d b = d z/ sink represents the differential length of the cutting edge, in which d z is differential axial height. The differential curve length d s, the area of cutting segment d A, and the angular position ψ i are derived as follows: where R 0 is the ball radius of the tool, ϕ = z R 0 tan I 0 represents the lag angle between the tip and a point on the helical flute at height z, i is the number of flute, I 0 is the helix angle, θ is measured clockwise from +y axis, N f is total number of flutes, and z is the Z-coordinate of a point located on the cutting edge. The coefficient where 0 R is the ball radius of the tool,   If the tool position does not vary much during the short cuts [16], the chip thickness becomes:

) Top view and section view of cutting
where f d represents feed per tooth. x(t j ) and y(t j ) represent the current tool-tip deflections.
x(t j−1 ) and y(t j−1 ) represent the tool-tip deflections at the previous time.
The radial, tangential, and axial forces are zero when the tool is not in cut. This can be modeled by multiplying the equations described as the force by a function W(ψ i ), which is defined when a tooth is in or out of cut. Moreover, the tooth is in cut if θ a + 2(n − 1)π ≤ ψ i ≤ θ b + 2(n − 1)π, where n is the number of spindle rotation, and θ a and θ b are the start and exit angles, respectively. This function is given by:  If the tool position does not vary much during the short cuts [16], the chip thickness becomes: The radial, tangential, and axial forces are zero when the tool is not in cut. This can be modeled by multiplying the equations described as the force by a function which is defined when a tooth is in or out of cut. Moreover, the tooth is in cut if where n is the number of spindle rotation, and a  and b  are the start and exit angles, respectively. This function is given by: Substituting Equations (4), (5), and (7) into Equation (3) yields the following equation: where the coefficients Using the coordinate transformation, the differential cutting forces in Equation (9) are expressed in Cartesian coordinate as: Substituting Equations (4), (5) and (7) into Equation (3) yields the following equation: where the coefficients f 2 (ψ i ) and f 3 (ψ i ) are listed in Appendix B.
Using the coordinate transformation, the differential cutting forces in Equation (9) are expressed in Cartesian coordinate as: Integrating the differential cutting forces in Equation (10), the authors [18] gave the cutting force components at the tool rotation angle θ in Cartesian coordinates and expressed in terms of the tool-tip responses as follows: where ϕ a =â R 0 tan I 0 in whichâ represents axial depth of cut. F xi (θ), F yi (θ), and F zi (θ), i = 1, 2, 3 are given in Appendix B.

Responses of System Using GPEM
The dynamic cutting forces expressed in Equation (10) are the external forces of the modeled flexible rotor-bearing system for ball end milling. The GPEM will be applied to solve the vibration problem of the tool-tip for ball end milling combined with the rotor-bearing system. Applying the transformation relation in Equation (1), the tool-tip deflections x(t j ) and y(t j ) can be transformed into the generalized coordinates a n (t) and b m (t) shown below: Substitute Equation (12) into Equation (10) to yield the force terms of ball end mill in generalized coordinates.
Applying the principle of virtual work to the generalized coordinate system and neglecting the z-direction force F z , the external forces F x and F y in Equation (13) can be rewritten as the right side of Equation (2) with generalized coordinates. Therefore, Equation (2) can be modified into the following equation for the whole tool spindle system.
It should be noted that the first part of the milling force, on the right side of Equation (14), can be combined with the stiffness matrix to change the system parameters of the rotor-bearing system. The cutting forces, which depend, in part, on the tool-tip response, are combined with the rotor system to complicate the problem, and the calculated response becomes more accurate. Expressing Equation (14) into the state space form as The X k at time t k can be solved by using the numerical method. The displacement of tooltip in the x direction and the y direction can be determined as a 1 (t) + a 2 (t) + . . . + a N P (t) and b 1 (t) + b 2 (t) + . . . + b N P (t), respectively. Furthermore, the natural frequency can also be determined by solving the eigenvalue problem of Equation (14).

Numerical Results and Discussion
The vibration responses of the tool-tip for ball end milling are solved numerically using the GPEM for the flexible rotor-bearing system in this study. The used cutting conditions and the cutting parameters are given in Table 1 [14,18]. The cutting force coefficients K ac , K tc , K rc are given in Appendix C. Table 2 gives the configuration data for the rotor-bearing system modeled in Figure 1b. E and ρ represent the shaft elastic modulus and density, respectively. d O and d I represent the shaft outer and inner diameters, respectively.
The authors in [18] proposed the approximated critical speeds of the spindle N cr (rpm) as: where f n (rad/s) represents the natural frequency of the tool-tip, N f represents the flute number, and m k represents the vibration cycle number of the tool from one flute to another during cutting. Different combinations of m k and N f lead to different excitation lines in the whirl map of the flexible rotor-bearing system, which will give different relationships between natural frequency and critical speed. Figure 4 is the whirl map which shows that the forward and backward modes are excited in the flexible rotor system under spindle speeds. The natural frequencies of forward mode increase as spindle speeds increase, but those of backward mode decrease as spindle speeds increase. The parameters N f and m k determine the excitation line in the whirl map. In this study, m k = 1, and two cases including N f = 2, 3, are used. Thus, the two whirl ratios are 1/2 and 1/3. The excitation line may intersect the natural mode curves at some corresponding spindle speeds, called the critical speeds. The natural frequencies and critical speeds of various modes can be expressed in Table 3 to show that the critical speeds of forward mode are higher than those of backward mode. Their associated mode shapes of forward and backward modes at critical speeds, with 2 and 3 flutes of the tool, are expressed in Figure 5. It is shown that the eigenvector magnitude for the mode 1 is the largest at tool-tip C, but those for the mode 2 is relatively small. In addition, when the spindle rotates at a higher critical speed, the frequency difference between the forward mode and backward mode becomes larger. This phenomenon is mainly due to the gyroscopic effect of the flexible rotor. Therefore, some lower natural modes will be significantly excited to become the main components of the tool-tip response. This will be verified in this section.            . The mode shapes for the rotor-bearing system at the critical speeds. Figure 5. The mode shapes for the rotor-bearing system at the critical speeds.
N f = 2: The numerical studies for N f = 2 at different spindle speed are shown in Figures 6-8. Figure 6a shows the displacement of the tool-tip with two flutes in the x direction, at the spindle speed of 1320 rpm, which is close to the critical speed of the first backward mode. Figure 6b also shows the x-y plane orbit of the tool-tip. The orbit is an encircled trajectory which expands clockwise around the origin in a stable loop. This mode is verified to be a backward mode. The fast Fourier transformation diagram (FFT diagram) for the tool-tip response in the x direction, at the spindle speed of 1320 rpm, is shown in Figure 6c. The tool is with two flutes. It is found that only one major peak, with the frequency of 276.5 rad/s, dominates the response of the rotor system.  Figure 6a shows the displacement of the tool-tip with two flutes in the x direction, at the spindle speed of 1320 rpm, which is close to the critical speed of the first backward mode. Figure 6b also shows the x-y plane orbit of the tool-tip. The orbit is an encircled trajectory which expands clockwise around the origin in a stable loop. This mode is verified to be a backward mode. The fast Fourier transformation diagram (FFT diagram) for the tool-tip response in the x direction, at the spindle speed of 1320 rpm, is shown in Figure 4c. The tool is with two flutes. It is found that only one major peak, with the frequency of 276.5 rad/s, dominates the response of the rotor system.   Figure 7a shows the displacement in the x direction of the tool-tip with two flutes, at the spindle speed of 6040 rpm, which is close to the critical speed of the second backward mode. It is found that the tool oscillates with the cutting-force frequency combined with a longer cycle vibration. The period of the long cycle is about 0.056 s, which explains how the workpiece surface is often found to have a lot of small waves inside the long cycle wave in the practical milling. From Figure 4, it can be seen that the lower first backward mode with the natural frequency of 112.1 rad/s is excited at the spindle speed of 6040 rpm. Figure 7b plots the tool-tip orbit. It can be seen from the figure that the track of the tip of the tool circles around, and there are many small rings inside. The large ring is excited by the second backward mode, and the small ring corresponds to the first backward mode. This phenomenon can be more clearly observed from the FFT plot shown in Figure 7c. At the spindle speed of 6040 rpm with two flutes, there are two main peaks in the x-response of the tool-tip. The two frequencies governing the behaviors of tool-tip correspond to the first and second backward modes, respectively. The obvious low frequency peak verifies the long cycle wave phenomenon. A very small peak near the right hand of cutting-force frequency is also found. It can be seen, from Figure 4, to coincide with the first forward mode at the spindle speed of 6040 rpm.  Figure 7a shows the displacement in the x direction of the tool-tip with two flutes, at the spindle speed of 6040 rpm, which is close to the critical speed of the second backward mode. It is found that the tool oscillates with the cutting-force frequency combined with a longer cycle vibration. The period of the long cycle is about 0.056 s, which explains how the workpiece surface is often found to have a lot of small waves inside the long cycle wave in the practical milling. From Figure 4, it can be seen that the lower first backward mode with the natural frequency of 112.1 rad/s is excited at the spindle speed of 6040 rpm. Figure 7b plots the tool-tip orbit. It can be seen from the figure that the track of the tip of the tool circles around, and there are many small rings inside. The large ring is excited by the second backward mode, and the small ring corresponds to the first backward mode. This phenomenon can be more clearly observed from the FFT plot shown in Figure 7c. At the spindle speed of 6040 rpm with two flutes, there are two main peaks in the x-response of the tool-tip. The two frequencies governing the behaviors of tool-tip correspond to the first and second backward modes, respectively. The obvious low frequency peak verifies the long cycle wave phenomenon. A very small peak near the right hand of cutting-force frequency is also found. It can be seen, from  Figure 8a shows the displacement of the tool-tip with two flutes in the x direction, at the spindle speed of 19,950 rpm, which is close to the critical speed of the first forward mode. It is found that the tool oscillates with the cutting-force frequency combined with a longer cycle vibration. The period of the long cycle is 0.173 s because the first backward mode is excited and the associated natural frequency is 36.31 rad/s. The period of long cycle at spindle speed 19,950 rpm is larger than that at spindle speed 6040 rpm. Figure 8b plots the orbit of the tool-tip. It shows the orbital trajectory of the tool-tip around a large circle, with many smaller inner circles. The encircled directions of small inner loops are anticlockwise, so this mode is dominated by the forward mode (second forward mode), and the encircled direction of large circle is clockwise, so this mode is dominated by the backward mode (first backward mode). Figure 8c shows the FFT plot of the x-response of the tool-tip with two flutes at the spindle speed of 19,950 rpm. Two frequencies govern the behaviors of the tool-tip vibration. An obvious low-frequency peak is found to verify that the first backward mode is excited. This is the long cycle wave phenomena.    Figure 8a shows the displacement of the tool-tip with two flutes in the x direction, at the spindle speed of 19,950 rpm, which is close to the critical speed of the first forward mode. It is found that the tool oscillates with the cutting-force frequency combined with a longer cycle vibration. The period of the long cycle is 0.173 s because the first backward mode is excited and the associated natural frequency is 36.31 rad/s. The period of long cycle at spindle speed 19,950 rpm is larger than that at spindle speed 6040 rpm. Figure 8b plots the orbit of the tool-tip. It shows the orbital trajectory of the tool-tip around a large circle, with many smaller inner circles. The encircled directions of small inner loops are anticlockwise, so this mode is dominated by the forward mode (second forward mode), and the encircled direction of large circle is clockwise, so this mode is dominated by the backward mode (first backward mode). Figure 8c shows the FFT plot of the x-response of the tool-tip with two flutes at the spindle speed of 19,950 rpm. Two frequencies govern the behaviors of the tool-tip vibration. An obvious low-frequency peak is found to verify that the first backward mode is excited. This is the long cycle wave phenomena.    Figure 9a. The spindle speed is 960 rpm, close to the critical speed of the first backward mode. Figure 9b shows the orbit of the tool-tip. It shows that an encircled trajectory of the tool-tip expands clockwise around the origin in a stable loop. This phenomenon is the same as in Figure 6b, because only the first backward mode of the rotor-bearing system is excited.  Figure 8a shows the displacement of the tool-tip with two flutes in the x direction, at the spindle speed of 19,950 rpm, which is close to the critical speed of the first forward mode. It is found that the tool oscillates with the cutting-force frequency combined with a longer cycle vibration. The period of the long cycle is 0.173 s because the first backward mode is excited and the associated natural frequency is 36.31 rad/s. The period of long cycle at spindle speed 19,950 rpm is larger than that at spindle speed 6040 rpm. Figure 8b plots the orbit of the tool-tip. It shows the orbital trajectory of the tool-tip around a large circle, with many smaller inner circles. The encircled directions of small inner loops are anticlockwise, so this mode is dominated by the forward mode (second forward mode), and the encircled direction of large circle is clockwise, so this mode is dominated by the backward mode (first backward mode). Figure 8c shows the FFT plot of the x-response of the tool-tip with two flutes at the spindle speed of 19,950 rpm. Two frequencies govern the behaviors of the tool-tip vibration. An obvious low-frequency peak is found to verify that the first backward mode is excited. This is the long cycle wave phenomena. N f = 3: Figures 9-12 show the numerical studies for N f = 3 at different spindle speed. The tool-tip displacement in x direction with three flutes is shown in Figure 9a. The spindle speed is 960 rpm, close to the critical speed of the first backward mode. Figure 9b shows the orbit of the tool-tip. It shows that an encircled trajectory of the tool-tip expands clockwise around the origin in a stable loop. This phenomenon is the same as in Figure 6b, because only the first backward mode of the rotor-bearing system is excited. Figure 9c shows the FFT plot of the x-response of the tool-tip with three flutes at the spindle speed of 960 rpm. Only one dominant frequency is found to dominate the response of the tool-tip.    Figure 9a. The spindle speed is 960 rpm, close to the critical speed of the first backward mode. Figure 9b shows the orbit of the tool-tip. It shows that an encircled trajectory of the tool-tip expands clockwise around the origin in a stable loop. This phenomenon is the same as in Figure 6b, because only the first backward mode of the rotor-bearing system is excited. Figure 9c shows the FFT plot of the x-response of the tool-tip with three flutes at the spindle speed of 960 rpm. Only one dominant frequency is found to dominate the response of the tool-tip.   Figure 10a shows the displacement in the x direction of the tool-tip with three flutes, at the spindle speed of 2160 rpm, which is close to the critical speed of the first forward mode. The orbit of tool-tip is plotted in Figure 10b. The result shows an encircled trajectory of the tool-tip expands anticlockwise around the origin in a stable loop, so this mode is dominated by a forward mode. Figure 10c shows the FFT plot of the x-response of the tool-tip with three flutes at the spindle speed of 2160 rpm. A distinct peak with a  Figure 11a shows the displacement of the tool-tip with three flutes in the x direction, at the spindle speed of 4410 rpm, which is close to the critical speed of the second backward mode. The tool oscillates with the cutting-force frequency combined with a longer cycle vibration which is the first backward mode with a natural frequency of 138.54 rad/s. The period of the long cycle is 0.045 s. Figure 11b plots the orbit of the tool-tip. It shows that the trajectory of the tool-tip goes around a large circle with many smaller inner loops. The FFT plot for the x-response of the tool-tip with three flutes at the spindle speed of 4410 rpm is shown in Figure 11c. There are two major peaks that govern the vibration of the tool-tip. The obvious low frequency peak, dominated by the first backward mode, corresponds to the long cycle wave phenomena.   Figure 11a shows the displacement of the tool-tip with three flutes in the x direction, at the spindle speed of 4410 rpm, which is close to the critical speed of the second backward mode. The tool oscillates with the cutting-force frequency combined with a longer cycle vibration which is the first backward mode with a natural frequency of 138.54 rad/s. The period of the long cycle is 0.045 s. Figure 11b plots the orbit of the tool-tip. It shows that the trajectory of the tool-tip goes around a large circle with many smaller inner loops. The FFT plot for the x-response of the tool-tip with three flutes at the spindle speed of 4410 rpm is shown in Figure 11c. There are two major peaks that govern the vibration of the tool-tip. The obvious low frequency peak, dominated by the first backward mode, corresponds to the long cycle wave phenomena.  Figure 12a shows the displacement of the tool-tip with three flutes in the x direction, at the spindle speed of 9850 rpm, which is close to the critical speed of the second forward mode. The tool oscillates with the cutting-force frequency combined with a longer cycle vibration. The long cycle vibration corresponds to the first backward mode with a natural frequency of 71.17 rad/s. The period of the long cycle is 0.088 s. Figure 12b shows the orbit of the tool-tip, which wraps around a large circle with many small inner rings. Figure 12c shows the FFT plot of the x-response of the tool-tip with three flutes at longer cycle vibration. The long cycle vibration corresponds to the first backward mode with a natural frequency of 71.17 rad/s. The period of the long cycle is 0.088 s. Figure 12b shows the orbit of the tool-tip, which wraps around a large circle with many small inner rings. Figure 12c shows the FFT plot of the x-response of the tool-tip with three flutes at the spindle speed of 9850 rpm. It is found that there are two frequencies governing the behaviors of the tool-tip. An obvious low-frequency peak is found. It corresponds to the long cycle wave phenomena, and verifies that the first backward mode is excited.  Figure 10a shows the displacement in the x direction of the tool-tip with three flutes, at the spindle speed of 2160 rpm, which is close to the critical speed of the first forward mode. The orbit of tool-tip is plotted in Figure 10b. The result shows an encircled trajectory of the tool-tip expands anticlockwise around the origin in a stable loop, so this mode is dominated by a forward mode. Figure 10c shows the FFT plot of the x-response of the tool-tip with three flutes at the spindle speed of 2160 rpm. A distinct peak with a frequency of 678.6 rad/s, and a far smaller peak with a frequency of 237.7 rad/s, are found. Figure 11a shows the displacement of the tool-tip with three flutes in the x direction, at the spindle speed of 4410 rpm, which is close to the critical speed of the second backward mode. The tool oscillates with the cutting-force frequency combined with a longer cycle vibration which is the first backward mode with a natural frequency of 138.54 rad/s. The period of the long cycle is 0.045 s. Figure 11b plots the orbit of the tool-tip. It shows that the trajectory of the tool-tip goes around a large circle with many smaller inner loops. The FFT plot for the x-response of the tool-tip with three flutes at the spindle speed of 4410 rpm is shown in Figure 11c. There are two major peaks that govern the vibration of the tool-tip. The obvious low frequency peak, dominated by the first backward mode, corresponds to the long cycle wave phenomena. Figure 12a shows the displacement of the tool-tip with three flutes in the x direction, at the spindle speed of 9850 rpm, which is close to the critical speed of the second forward mode. The tool oscillates with the cutting-force frequency combined with a longer cycle vibration. The long cycle vibration corresponds to the first backward mode with a natural frequency of 71.17 rad/s. The period of the long cycle is 0.088 s. Figure 12b shows the orbit of the tool-tip, which wraps around a large circle with many small inner rings. Figure 12c shows the FFT plot of the x-response of the tool-tip with three flutes at the spindle speed of 9850 rpm. It is found that there are two frequencies governing the behaviors of the tool-tip. An obvious low-frequency peak is found. It corresponds to the long cycle wave phenomena, and verifies that the first backward mode is excited.
From a comparison of Figures 6-8 and comparison of Figures 9-12, it can be seen that the period of the long cycle vibration enlarges as the spindle speed increases. To gain an insight into the phenomenon, it can be seen from Figure 4 that, with the increase in the spindle speed, the natural frequency of the first backward mode decays. Additionally, in the numerical study, it is found that the long cycle vibration is majorly dominated by the first backward mode, thus, the period of the long cycle vibration enlarges as the spindle speed increases. In addition, it can also be seen from the FFT diagram that the peak at the first backward mode (the long-period vibration) becomes larger compared with the peak at the cutting force frequency as the spindle speed becomes larger. The tool-tip response including the forward mode and backward mode is mainly dependent on the spindle speed. The forward and backward modes that vary with the spindle speed are mainly caused by the gyroscopic effect of the flexible rotor. The vibration analysis of a ball end mill is more accurate if the spindle is modelled as a flexible rotor, taking into account the gyroscopic effects caused by the rotation of the rotor section. The cutting forces, which depend, in part, on the tool-tip response, combine with the stiffness matrix of the rotor system to complicate the issue. Therefore, in order to obtain a more accurate vibration analysis for ball end milling, the flexible rotor-bearing system cannot be omitted.

Conclusions
The dynamic behaviors of the flexible rotor-bearing system for ball end milling were first numerically investigated by using the GPEM. The tool-tip responses have been made more realistic due to consideration of the coupling effect of the flexible rotor-bearing system. Thanks to the gyroscopic effect of the flexible rotor, the important whirl map can be drawn. Then, it can be observed from the whirl map that, the critical speeds, which are directly proportional to the natural frequency and inversely proportional to the number of flutes, due to some natural modes, can be excited to significantly affect the tool responses. The flexible isotropic rotor-bearing system with rotating speed can excite the forward and backward modes along the excitation line in the whirl map to form two different critical speeds. In addition, the cutting forces, which depend, in part, on the tool-tip response, combine with the stiffness matrix of the rotor system to complicate the issue. To the best of the authors' knowledge, most previous studies on ball end milling always ignore the rotor system or assume that the rotor is rigid or discretized to a few degrees of freedom, thereby ignoring the gyroscopic effect of the flexible rotor. Unlike the previous literature, the natural frequency in this study varies with the spindle speed. Thanks to the gyroscopic effect of the flexible rotor, which of the different modes are included in the tool-tip response depend mainly on the spindle speed. The model established by considering the coupling effect of the flexible rotor and using GPEM can be more efficient and realistic in analyzing the dynamic behavior of the tool-tip. The important new findings of the present study, by using numerical analysis, are obtained as follows:

1.
In milling at a higher spindle speed (2160-19,950 rpm), the tool-tip oscillates with the cuttingforce frequency combined with a longer cycle vibration of the first backward mode due to the coupling effects of the flexible rotor-bearing system, and the period of the long cycle vibration enlarges as the spindle speed increases. From the FFT plot, it can be seen that as the spindle speed becomes larger, the peak amplitude at the first backward mode (long-period vibration) becomes larger compared to the peak at the cutting force frequency.

2.
Under lower spindle speed milling (960, 1320 rpm), the vibration trajectory of the tool-tip in the co-plane expands clockwise gradually around the origin until arriving at a stable loop. This is because only the first backward mode of the rotor-bearing system is excited.
The authors of [18] calculated the critical speeds of the rotor-bearing system for ball end milling and verified the numerical results with the experimental data. However, the rotor was assumed to be rigid and the lateral rotation was neglected, so the gyroscopic effects of the system could not be included. Therefore, the critical speed was independent of the spindle speed and only one mode was calculated. In this study, the critical speeds of the flexible rotor-bearing system vary with the spindle speed due to the gyroscopic effect. The vibration responses of the tool-tip, including forward and backward modes, are numerically calculated. However, it is also recommended that future work is required to add a comparison by setting up an experiment to validate the presented results.
K s (m, n) = 1 0 (n − 1)(n − 2)(m − 1)(m − 2) E(ẑ)I(ẑ)ẑ n+m−6 dẑ, (A6) where l and z represent the shaft length and axial coordinate, respectively. N p , N d , and N b denote the total number of polynomial, disk, and bearing, respectively. ρ(ẑ) and A(ẑ) represent the shaft density and cross-sectional area, respectively. I D and I P represent the shaft diametral and polar mass moment of inertia, respectively. m d i represents the i-th disk mass. I d Di and I d Pi represent the i-th disk diametral and polar mass moment of inertia, respectively. C b xxj , C b xyj , and C b yyj represent the j-th bearing damping coefficients. E(ẑ) and I(ẑ) represent the shaft elastic modulus and cross-sectional moment of inertia, respectively. K b xxj , K b xyj , and K b yyj represent the j-th bearing elastic constant.