Analysis on Free Vibration and Stability of Rotating Composite Milling Bar with Large Aspect Ratio

: The free vibration behavior and stability of composite milling bar with large aspect ratio are analyzed. Speciﬁcally, the free vibration equations are derived based on Euler Bernoulli beam theory and Hamilton’s principle and solved by Galerkin method. In addition, to investigate the stability of the cutting system with a rotating composite milling bar, this study develops an analytical model for regeneration and cutting force ﬂuttering, in which the internal damping, external damping, gyroscopic e ﬀ ect, and inertial e ﬀ ect are considered. The subsequent stability lobes are obtained by the time domain method. Then, the stability of composite and conventional metal milling bar is compared. The e ﬀ ects of internal damping, external viscous damping, ply angle, gyroscopic e ﬀ ect, and inertial e ﬀ ect on cutting stability are analyzed.


Introduction
Milling is one of the most extensively used machining methods, where the milling bar for deep groove milling is usually designed as a long and slender cantilever structure. Due to the low dynamic stiffness and low natural frequency of metal milling bar, high speed milling can easily trigger cutting chatter, which can initiate fatigue failure, noise, reduced machining accuracy of the work piece surface, and low machining efficiency.
In high-speed milling applications, in order to avoid chatter and improve stability of milling, the rotating milling bar should have high dynamic stiffness (Resistance to deformation under dynamic load) and high specific stiffness (Ratio of elastic modulus to density of material). The dynamic stiffness is the product of the modal stiffness and the modal damping ratio. Some unique devices, such as tuned bar, impact damper, Lanchester damper, viscoelastic damper [1][2][3][4][5], and constrained layer damping [6], have been proposed for machine systems to enhance the dynamic stiffness. However, due to geometric diversity and installation difficulties of the milling bar, the unique devices are not reliable in suppressing the chatters of boring bar in deep-hole machining. In addition, the metal bars have low specific stiffness and low natural frequencies, which limit the rotating speed of the milling bar during metal cutting. In fact, it is challenging to simultaneously improve the dynamic stiffness and machining speed of metal bar because traditional metals with high stiffness usually have low damping and low specific stiffness [7].
Composite materials are composed of two or more materials with different properties. These composite materials have new properties by physical or chemical methods. Various materials learn from each other in performance and produce synergistic effect, thus the comprehensive performance of composite materials is better than the raw materials and can meet various requirements. In addition, the composite has the advantages of designable mechanical properties. Owing to their high specific strength, specific stiffness, lightweight and anti-fatigue properties, composite materials have been widely applied in aerospace, military industry and mechanical engineering. The composite material offers the possibility to simultaneously enhance the dynamic stiffness and fundamental natural frequency of milling bars with large aspect ratio.
Chatter can be divided into two categories based on their self-excited vibration mechanisms, i.e., regenerative chatter and mode-coupling chatter. Regenerative chatter occurs due to the regeneration of surface corrugation generated by the previous cutting [8]. Tlusty and Polacek [9] analyzed the stability limit of regenerative chatter in an orthogonal cutting. Altintas et al. [10][11][12] obtained the stability lobe diagram of a milling system in frequency domain by the frequency analysis method. This method has a high computational efficiency, and the gyroscopic effects activated by the rotation of milling bar is included in the analysis. On the other hand, in other literature [13][14][15][16], the stability lobe diagrams of the chatter system were obtained in time domain. Urbikain et al. [17][18][19] conducted a series of research on cutting issues, including the milling stability analysis considering the flexibility of cutter and workpiece, the influence of the tilting machining of circular segment end mills on the surface roughness of workpiece, the influence of the tilting machining of the barrel milling cutter on the stability, etc. In Literature [20], the first-order modal frequency domain method and the Chebyshev collocation method were used to analyze the cutting stability of large horizontal lathe. Based on comparison, the Chebyshev polynomial method was found to have higher calculation accuracy, while the first-order modal method had faster calculation speed. Mossadegh et al. [21] used finite element model to analyze milling chatter under gyroscopic effect. The iteration method was used to calculate the stability lobe diagram. Firstly, an initial spindle speed was set to solve the characteristic equation. Then the solved equation was used to determine the flutter frequency, critical cutting width and spindle speed. The next iteration of updating the equation and parameters was determined according to the convergence of spindle speed. Arvajeh et al. [22] and Roukema et al. [23,24] established a cutting process model of rotary two-degree-of-freedom drill pipe and used the frequency domain method to study the stability of the cutting process. Due to the particularity of the drilling force model, it is not necessary to solve the characteristic equation. Instead, the cutting depth expression can be directly derived; this is the stability analysis of isotropic cutting bar, including rotation speed. In addition, the design and manufacture schemes of composite boring bar have been proposed in previous reports [25,26], in which the researchers attempted to improve the static and dynamic rigidity of the boring bar. So far, most of the work has been directed to the structural design and subsequential experimental verification, while only little work has focused on a dynamic analytical model of a rotating composite milling bar. Kim et al. [27] studied the response and stability of a high-speed rotating composite shaft subjected to the typical cutting forces of milling, and Timoshenko beam theory was used for the structural dynamic model of rotating composite shaft. Meanwhile, a deflection-dependent cutting forces model with regeneration delay effects was established and the dynamic milling model was solved by frequency domains method. The influence of inertial force (Square term of rotating speed) was not considered in the model. This paper aims at providing a new dynamic model of composite milling bar, which is a cylindrical tubular composite shaft modeled as a continuous distributed Euler Bernoulli beam without considering the thickness of the bar. For shaft made of isotropic material, when the aspect ratio is larger than 10, the Euler Bernoulli beam theory is usually used. This is a simple model, but accurate results can be obtained [28]. When the Euler beam theory is used, the number of unknown functions of the model is greatly reduced because it does not consider the shear deformation, which greatly promotes the modeling of the cutter bar structure and the stability of the cutting system, and ensures the accuracy of the solution Hamilton's principle was used to derive free vibration equations of the bar, then Galerkin approach was used to obtain the characteristics of the free vibration of the bar. To measure the stability of the bar in a cutting process, the cutting loads including regeneration delay effects were incorporated into the free vibration equations. The stability lobes were obtained in the discrete time domain [14], and the variation of stability lobe diagrams as a function of different parameters including external damping, internal damping, ply angle and gyroscopic and inertial effects was investigated. In addition, the methods in the frequency domains based on zero-order solution [14] were also used in the analysis. The chatter frequency and angular displacement curves corresponding to the stable lobes indicate the physical significance of the regenerated chatter of the milling bar. It was found that the application of composite expanded the stable regions in milling process to a certain extent (15.87% with the given material and aspect ratio). For the first time, the effect of whirl caused by rotation on the cutting stability of composite milling bar was analyzed. The results show that the backward whirl has a negative damping effect, which is not conductive to the cutting stability.

Governing Equations
The bar is modeled as a rotating cantilever shaft as shown in Figure 1. It rotates at a constant angular velocity around its axis. Ω, b, t, r and L represent the outer diameter, wall thickness, the radius of midline of cross-section profile and the length of the milling bar, respectively. (X, Y, Z) denote the inertial coordinate system, and the corresponding unit vectors are (I, J, K); (x, y, z) denote the coordinate systems rotating with shaft axis and the corresponding unit vectors are (i, j, k). (ξ, s, x) is a local coordinate system, s refers to the tangential direction of midline of cross-section profile, ξ designates the normal line, and η is the ply angle.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 18 chatter frequency and angular displacement curves corresponding to the stable lobes indicate the physical significance of the regenerated chatter of the milling bar. It was found that the application of composite expanded the stable regions in milling process to a certain extent (15.87% with the given material and aspect ratio). For the first time, the effect of whirl caused by rotation on the cutting stability of composite milling bar was analyzed. The results show that the backward whirl has a negative damping effect, which is not conductive to the cutting stability.

Governing Equations
The bar is modeled as a rotating cantilever shaft as shown in Figure 1. It rotates at a constant angular velocity around its axis. Ω, b, t, r and L represent the outer diameter, wall thickness, the radius of midline of cross-section profile and the length of the milling bar, respectively. (X, Y, Z) denote the inertial coordinate system, and the corresponding unit vectors are (I, J, K); (x, y, z) denote the coordinate systems rotating with shaft axis and the corresponding unit vectors are (i, j, k). (ξ, s, x) is a local coordinate system, s refers to the tangential direction of midline of cross-section profile, ξ designates the normal line, and η is the ply angle.
u v w denote the X, Y, and Z components of the extension and bending displacements of any point on the bar's axis, as shown in Figure 2, and φ is the rotating angle of the cross-section around x-axis. Based on Euler Bernoulli beam theory which does not consider transvers shear, the linear strain-displacement relationship can be expressed as where u, v, w denote the X, Y, and Z components of the extension and bending displacements of any point on the bar's axis, as shown in Figure 2, and ϕ is the rotating angle of the cross-section around x-axis.
Since the geometry of the bar is a hollow cylinder, it is convenient to use the cylindrical coordinate to express the stress-strain relationship of the bar. The strain components in cylindrical coordinate can be expressed as u v w denote the X, Y, and Z components of the extension and bending displacements of any point on the bar's axis, as shown in Figure 2, and φ is the rotating angle of the cross-section around x-axis. Based on the stress-strain relationship of the composite, the stresses in the kth layer of the bar can be expressed as ij is listed in reference [29]. The strain energy can be expressed as Taking variation of the above strain energy expression and considering Equations (2) and (3), the following equation can be obtained The next step is to derive the kinetic energy of the bar. In order to simplify the calculation, the size and shape of the bar's cross-section are assumed to be unchanged. As shown in Figure 2, the consequent position vector of any point on the bar can be expressed as [30] R P/O = (x + u 1 )i + (y + u 2 )j + (z + u 3 )k (8) where u 1 , u 2 , u 3 are the displacements of any point on the bar along three coordinate axes.

of 17
Considering that the cutter rod rotates around the rotating axis, the speed of any point in the inertial coordinate system can be expressed as where Ω is the spinning speed of the bar. The kinetic energy T of the bar can be expressed as By taking the first order variation of the kinetic energy and considering the circle cross-section, the following can be obtained. where By combining the energy Equation (5) and Equation (11) and adopting the Hamilton's principle, the following equation can be obtained.
The equations of the motion are expressed as follows: The above equations show that the bending motions are decoupled from the axial and torsional motions. The single underlined terms in Equations (15) and (16) represent gyroscopic effect, while the double underlined terms represent the inertial forces. It can be seen from Equations (14)-(17) that the vibration of the boring bar has tension-torsion coupling and bending-bending coupling.
In milling chatter, bending vibration is the main vibration mode, while tension-torsional vibration and bending-bending vibration are independent and have no coupling effect. Thus, only bending vibration is analyzed below.
Damping is important to restrain or enhance instability. The external viscous damping is always beneficial to the stability of the system, and the internal viscous damping may destroy the stability (depending on the rotation speed). So, it is necessary to add damping to the equation of motion for the system.

Considering influence of internal and external damping and the cutting forces, the bending vibration equation can be modified as
in terms of underline in Equations (18) and (19), they indicate the newly added terms. c e is an external viscous damping constant per unit length and c i is an internal viscous damping constant per unit length. Both c e and c i are proportional to the rotational speed. However, c e is related to the speed in the inertial system, while c i is related to the speed in the rotating system [31].
The cutting forces are modeled as concentrated loads at z = L. Then f v and f w in Equation (18) and Equation (19) can be specified by Dirac-delta function as The expression of tangential force F v and radial force F w are based on the report by Tlusty et al. They are deflection dependent and involve delay terms and periodic coefficients [32]. The specific expression of cutting force can be found in Reference 23.

Approximate Solution
Since the exact solutions of Equations (18) and (19) cannot be obtained, Galerkin method is used to simplify the bending deformation expression as follows: where, ξ j (x) denotes the mode shape of bending motion, which is the mode shapes of a uniform, non-rotating, isotropic, and fixed-free Euler-Bernoulli beam, namely Substituting Equation (22) into Equation (18) and Equation (19), the differential equation can be discretized by Galerkin method as Appl. Sci. 2020, 10, 3557 where M is the mass matrix, K is the stiffness matrix, and G is the damping matrix. N is the number of mode shapes.
In the time domain analysis method, the constant force is neglected [27]. Substituting Equation (20) into Equation (23) and eliminating the time delay terms on the right side of the equation, Equation (23) can be re-written as ..

W(t)
where The matrices K f and K p are related to the cutting force. K s is specific cutting force per unit chip area, θ t denotes the angular position of the engaged tooth (0 ≤ θ t < 0.5π). β 0 is the constant angle between the normal of the cut surface and the direction of the cutting force. It is assumed that the cut surface does not fluctuate. The off-diagonal terms in K f have opposite sign, which may cause instability such as chatter. The elements of K p are periodic in time, so it can cause parameter instability.

W(t)
Equation (27) (periodic differential equation with time delay) can be solved by semi discretization method. In the ith semi-discretization interval, Equation (30) can be approximated as [22] By Cauchy transformation, Equation (31) can be written as Appl. Sci. 2020, 10, 3557 w a and w b are weight coefficients. Since the time delay is equal to the time period for the milling system, for any integer of j, w a = w b = 1/2. The diagonal elements of matrix I is 1, and the off-diagonal elements are 0. d i is the cutting depth.
For the initial condition q(t i ) = q i , q i+1 is determined as where From Equation (31), So, the 3rd and 4th columns of matrices B i and R i are zeroes. Therefore, The induced discrete mapping is where the coefficient matrix is where P i,hj and R i,hj are the elements of matrices P i and R i in the hth row and jth column, respectively. The final transition matrix Φ is determined as If the modulus of eigenvalues of Φ is less than one, the system is stable. When the eigenvalues equal to one, the milling system is in critical state. In other conditions, the system is unstable.

Free Vibration Analysis
Using the motion equations derived in this paper, the rotating bar's free vibration characteristics were analyzed, as shown in Figure 3. The material parameters are listed in Table 1. Meanwhile, two milling bars made of different materials, i.e., the isotropy (steel) and anisotropy (boron-epoxy, (0 6 /90 3 / ± 20 3 /90 3 /0 6 )) material were analyzed. The fundamental natural frequencies, which were critical to the stability, were calculated and compared. Figure 4 shows the variation of fundamental natural frequencies of bars with aspect ratio (L/2b). The results were also compared to the reports by Kim in 1999 based on Timoshenko beam theory [29].
The final transition matrix Φ is determined as If the modulus of eigenvalues of Φ is less than one, the system is stable. When the eigenvalues equal to one, the milling system is in critical state. In other conditions, the system is unstable.

Free Vibration Analysis
Using the motion equations derived in this paper, the rotating bar's free vibration characteristics were analyzed, as shown in Figure 3. The material parameters are listed in Table 1. Meanwhile, two milling bars made of different materials, i.e., the isotropy (steel) and anisotropy (boron-epoxy, (06/903/ ± 203/903/06)) material were analyzed. The fundamental natural frequencies, which were critical to the stability, were calculated and compared. Figure 4 shows the variation of fundamental natural frequencies of bars with aspect ratio (L/2b). The results were also compared to the reports by Kim in 1999 based on Timoshenko beam theory [29].      As shown in Figure 4, the fundamental natural frequencies of steel bar and composite bar gradually decreased with the increase of L/2b ratio. In addition, the natural frequencies of composite bar were always higher than the steel bar. For isotropic steel bar, when L/2b is larger than 5, the calculation results by Euler Bernoulli beam theory have high accuracy. However, for the anisotropic composite bar, in order to ensure the calculation accuracy, L/2b should be larger. As shown in Table  2, for composite material, L/2b should be larger than 8 to ensure the error less than 5%, which is the acceptable error of calculation results in many mechanical engineering fields. When the value of L/2b is small, the results of the two beam models are not in a close agreement, which is due to the significant effect of the transverse shear at small aspect ratio. Therefore, Euler beam theory should be used to analyze the composite milling bars with a larger aspect ratio.  As shown in Figure 4, the fundamental natural frequencies of steel bar and composite bar gradually decreased with the increase of L/2b ratio. In addition, the natural frequencies of composite bar were always higher than the steel bar. For isotropic steel bar, when L/2b is larger than 5, the calculation results by Euler Bernoulli beam theory have high accuracy. However, for the anisotropic composite bar, in order to ensure the calculation accuracy, L/2b should be larger. As shown in Table 2, for composite material, L/2b should be larger than 8 to ensure the error less than 5%, which is the acceptable error of calculation results in many mechanical engineering fields. When the value of L/2b is small, the results of the two beam models are not in a close agreement, which is due to the significant effect of the transverse shear at small aspect ratio. Therefore, Euler beam theory should be used to analyze the composite milling bars with a larger aspect ratio. When the effects of transverse shear are not considered, the established model by Galerkin method in this paper has good convergence. The natural frequencies of the bar (L/2b = 15) with different number of modes were calculated, as listed in Table 3. The results of the first three modes are very close, indicating that only one mode shape is enough for convergence. Due to the influence of rotating inertial force and gyroscopic force, the natural vibration frequency can be divided into forward whirl frequency and backward whirl frequency. When the backward whirl frequency is equal to 0, the rotating speed is the critical speed of the cutter bar. When the [η] 36 (η is shown in Figure 1) laying mode is adopted, the natural frequency of the cutter bars with different ply angles changes with the rotation speed, as shown in Figure 5. When the ply angle is 0 • , the critical speed of the cutter bar has the largest value. Due to the influence of rotating inertial force and gyroscopic force, the natural vibration frequency can be divided into forward whirl frequency and backward whirl frequency. When the backward whirl frequency is equal to 0, the rotating speed is the critical speed of the cutter bar. When the [η ]36 (η is shown in Figure 1) laying mode is adopted, the natural frequency of the cutter bars with different ply angles changes with the rotation speed, as shown in Figure 5. When the ply angle is 0°, the critical speed of the cutter bar has the largest value.

Stability Analysis
The bars made of composite and steel as shown in Figure 3 with the ratio of L/2b = 10 were analyzed. Ks is set to a typical value in milling process, i.e., Ks = 556.4Mpa [33].

Stability Analysis
The bars made of composite and steel as shown in Figure 3 with the ratio of L/2b = 10 were analyzed. K s is set to a typical value in milling process, i.e., K s = 556.4 Mpa [33]. Meanwhile, other parameters Appl. Sci. 2020, 10, 3557 11 of 17 are set as follows: β 0 = 73.3 • , ς steel = 0.01, ς composite = 0.015. In addition, the damping values were assigned based on concepts from an oscillator with single-degree of freedom, that is When c i /(c e + c i ) = 0.5, the stability lobes of both steel bar and composite milling bar are shown in Figure 6. Rotating speed (rad/s) Figure 5. The effects of ply angle on the natural frequency when L/2b = 10.

Stability Analysis
The bars made of composite and steel as shown in Figure 3 with the ratio of L/2b = 10 were analyzed. Ks is set to a typical value in milling process, i.e., Ks = 556.4Mpa [33]. Meanwhile, other parameters are set as follows:  Figure 6 shows that after composite material is used, the cutting depth limit is increased due to its high stiffness. At the lowest point of the rightmost lobe of the steel bar, the corresponding cutting depth and rotation speed are 0.063 mm and 647 rad/s, respectively. For the composite material, the corresponding cutting depth and rotation speed are 0.073 mm and 1124 rad/s, respectively. It can be  Figure 6 shows that after composite material is used, the cutting depth limit is increased due to its high stiffness. At the lowest point of the rightmost lobe of the steel bar, the corresponding cutting depth and rotation speed are 0.063 mm and 647 rad/s, respectively. For the composite material, the corresponding cutting depth and rotation speed are 0.073 mm and 1124 rad/s, respectively. It can be seen from the results in Figure 6 that at a large aspect ratio, the cutting depth of the cutter bar made of composite material (material properties are listed in Table 1) is larger than the cutter bar made of steel by 15.87%. The increase in the cutting depth by the use of composite material is not significant. However, when the composite material is used, the increase in the rotating speed is up to 73.7% (from 647 rad/s to 1124 rad/s), which can significantly improve the cutting efficiency. The results in frequency domain shown in Figure 6 are calculated based on the theory in Reference [27] proposed by Kim. The calculation results using the methods in this paper and in the reference are consistent, which demonstrates the reliability of the calculation results in this paper. Figures 7 and 8 show that as the value of ω c T jumps, chatter frequencies have abrupt changes. Physically speaking, ω c T represents the total angular displacement of the vibrating milling bar. During one rotation period, the milling bar vibrates at the chatter frequency. Meanwhile, k is the number of complete waves traversed by the tool. Therefore, ω c T is directly related to the phase difference between successive undulations on the work piece surface. Due to its larger natural frequency, the composite material has higher chatter frequency than steel, as shown in Figure 7. Meanwhile, at the same phase difference, the composite bar has larger rotating speed than the steel bar. Figure 9 shows the effects of the internal damping of materials. From the figure, when the internal damping is not considered, each lobe in the lobe diagram has the same lowest cutting depth. When the internal damping is considered, the lowest cutting depth increases with the rotating speed of the milling bar, as shown in the dot line. Furthermore, at low rotating speed, the influence of internal damping on stability is negligible. As the rotating speed increases, the influence of internal damping progressively becomes noticeable. At an extremely high rotating speed, the internal damping can induce an unstable zone. On the contrary, at a lower rotating speed, the stable zone can be obtained by larger internal damping. Therefore, the internal damping can improve the stability of the bar at a low cutting speed but cause instability at a high cutting speed. This conclusion is consistent with the results in Reference [34], which indicates that the internal damping maybe unstable, depending on rotating speed. Figure 10 shows the displacement response of the points of A, B and C in Figure 9 at different internal damping. The divergence and convergence results of the displacement responses indicate the existence of the unstable zone, which can be initiated by internal damping at an extremely high cutting speed. If rotation is not considered, this phenomenon will not occur. Therefore, the instability of internal damping at high speed is caused by the interaction of rotation and internal damping.
However, when the composite material is used, the increase in the rotating speed is up to 73.7% (from 647 rad/s to 1124 rad/s), which can significantly improve the cutting efficiency. The results in frequency domain shown in Figure 6 are calculated based on the theory in Reference [27] proposed by Kim. The calculation results using the methods in this paper and in the reference are consistent, which demonstrates the reliability of the calculation results in this paper. Figures 7 and 8 show that as the value of c T  jumps, chatter frequencies have abrupt changes.
Physically speaking, c T  represents the total angular displacement of the vibrating milling bar. During one rotation period, the milling bar vibrates at the chatter frequency. Meanwhile, k is the number of complete waves traversed by the tool. Therefore, c T  is directly related to the phase difference between successive undulations on the work piece surface. Due to its larger natural frequency, the composite material has higher chatter frequency than steel, as shown in Figure 7. Meanwhile, at the same phase difference, the composite bar has larger rotating speed than the steel bar.  Figure 9 shows the effects of the internal damping of materials. From the figure, when the internal damping is not considered, each lobe in the lobe diagram has the same lowest cutting depth. When the internal damping is considered, the lowest cutting depth increases with the rotating speed of the milling bar, as shown in the dot line. Furthermore, at low rotating speed, the influence of internal damping on stability is negligible. As the rotating speed increases, the influence of internal damping progressively becomes noticeable. At an extremely high rotating speed, the internal damping can induce an 647 rad/s to 1124 rad/s), which can significantly improve the cutting efficiency. The results in frequency domain shown in Figure 6 are calculated based on the theory in Reference [27] proposed by Kim. The calculation results using the methods in this paper and in the reference are consistent, which demonstrates the reliability of the calculation results in this paper. Figures 7 and 8 show that as the value of c T  jumps, chatter frequencies have abrupt changes.
Physically speaking, c T  represents the total angular displacement of the vibrating milling bar. During one rotation period, the milling bar vibrates at the chatter frequency. Meanwhile, k is the number of complete waves traversed by the tool. Therefore, c T  is directly related to the phase difference between successive undulations on the work piece surface. Due to its larger natural frequency, the composite material has higher chatter frequency than steel, as shown in Figure 7. Meanwhile, at the same phase difference, the composite bar has larger rotating speed than the steel bar.  Figure 9 shows the effects of the internal damping of materials. From the figure, when the internal damping is not considered, each lobe in the lobe diagram has the same lowest cutting depth. When the internal damping is considered, the lowest cutting depth increases with the rotating speed of the milling bar, as shown in the dot line. Furthermore, at low rotating speed, the influence of internal damping on stability is negligible. As the rotating speed increases, the influence of internal damping progressively becomes noticeable. At an extremely high rotating speed, the internal damping can induce an unstable zone. On the contrary, at a lower rotating speed, the stable zone can be obtained by larger internal damping. Therefore, the internal damping can improve the stability of the bar at a low cutting speed but cause instability at a high cutting speed. This conclusion is consistent with the results in Reference [34], which indicates that the internal damping maybe unstable, depending on rotating speed. Figure 10 shows the displacement response of the points of A, B and C in Figure 9 at different internal damping. The divergence and convergence results of the displacement responses indicate the existence of the unstable zone, which can be initiated by internal damping at an extremely high cutting speed. If rotation is not considered, this phenomenon will not occur. Therefore, the instability of internal damping at high speed is caused by the interaction of rotation and internal damping.   It is worth noting that some peculiar peaks appeared near 450 rad/s in the cutting depth curve of the composite milling bar with the damping ratio 0.025 in Figure 11a and near 750 rad/s in the cutting depth curve of steel milling bar in Figure 11b. A similar phenomenon can be seen in Figure 6. This phenomenon might be caused by the following two factors: (1) the accuracy of calculation method; (2) the stability change at different values of damping parameters. The influence of damping parameters on the stability is similar to the influence of radial immersion [35]. The exact reason needs  Figure 11 shows the variation of instability zones with the values of external viscous damping when c i = 0. The results indicate that both the isotropic milling bar and the anisotropy milling bar have the same trend. In addition, as the damping ratio increases, the stable cutting depth is larger. At all rotating speeds, when the external viscous damping intensified from ζ = 0.015 to ζ = 0.02, the stable cutting depth increased by approximately 25%. experimental verification. If the second factor is the reason for these peculiar peaks, the accuracy of the time domain method is higher than that of the frequency domain method.  It can be seen that when  = 0 , the composite bar has the largest stable zone. This result can be explained by the fact that when the fiber orientation is closer to 0 , the elastic modulus of composite is higher (as shown in Table 1), which results in greater stability due to the intensification of dynamic stiffness of the composite bar. However, composite milling bar with the ply angle of 0 might fail during cutting. The analysis results in this paper only focused on the bending deformation, but torsion should also be considered in practical work. The milling bar with the ply angle of 0 has the weakest strength of torsion deformation, which may cause damage in machining. The above stability analysis is based on a fixed natural frequency value at the rotating speed of zero. When the influence of tool bar rotation is considered, the natural frequency value is variable. Figure 13 shows the stability analysis results when the natural frequency changes without considering the internal damping. It can be seen that when the natural frequency is used as the backward whirl frequency, the cutting depth is significantly reduced. Therefore, the rotation effect plays a negative damping role in the stability analysis. It is worth noting that some peculiar peaks appeared near 450 rad/s in the cutting depth curve of the composite milling bar with the damping ratio 0.025 in Figure 11a and near 750 rad/s in the cutting depth curve of steel milling bar in Figure 11b. A similar phenomenon can be seen in Figure 6. This phenomenon might be caused by the following two factors: (1) the accuracy of calculation method; (2) the stability change at different values of damping parameters. The influence of damping parameters on the stability is similar to the influence of radial immersion [35]. The exact reason needs experimental verification. If the second factor is the reason for these peculiar peaks, the accuracy of the time domain method is higher than that of the frequency domain method. Figure 12 shows the effects of ply angle on the composite bar with L/2b = 10, c i /(c e + c i ) = 0.5 and lamination schemes [η] 36 . It can be seen that when η = 0, the composite bar has the largest stable zone. This result can be explained by the fact that when the fiber orientation is closer to 0, the elastic modulus of composite is higher (as shown in Table 1), which results in greater stability due to the intensification of dynamic stiffness of the composite bar. However, composite milling bar with the ply angle of 0 might fail during cutting. The analysis results in this paper only focused on the bending deformation, but torsion should also be considered in practical work. The milling bar with the ply angle of 0 has the weakest strength of torsion deformation, which may cause damage in machining.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 18 experimental verification. If the second factor is the reason for these peculiar peaks, the accuracy of the time domain method is higher than that of the frequency domain method.  It can be seen that when  = 0 , the composite bar has the largest stable zone. This result can be explained by the fact that when the fiber orientation is closer to 0 , the elastic modulus of composite is higher (as shown in Table 1), which results in greater stability due to the intensification of dynamic stiffness of the composite bar. However, composite milling bar with the ply angle of 0 might fail during cutting. The analysis results in this paper only focused on the bending deformation, but torsion should also be considered in practical work. The milling bar with the ply angle of 0 has the weakest strength of torsion deformation, which may cause damage in machining. The above stability analysis is based on a fixed natural frequency value at the rotating speed of zero. When the influence of tool bar rotation is considered, the natural frequency value is variable. Figure 13 shows the stability analysis results when the natural frequency changes without considering the internal damping. It can be seen that when the natural frequency is used as the backward whirl frequency, the cutting depth is significantly reduced. Therefore, the rotation effect plays a negative damping role in the stability analysis. The above stability analysis is based on a fixed natural frequency value at the rotating speed of zero. When the influence of tool bar rotation is considered, the natural frequency value is variable. Figure 13 shows the stability analysis results when the natural frequency changes without considering the internal damping. It can be seen that when the natural frequency is used as the backward whirl frequency, the cutting depth is significantly reduced. Therefore, the rotation effect plays a negative damping role in the stability analysis.

Conclusions
A continuous dynamic model was developed for a rotating composite milling bar. First, the natural frequencies of free vibration were calculated by Galerkin method. Then the stabilities of milling bar subject to deflection-dependent cutting end loads were analyzed in time domain. Furthermore, the model included delay effects due to regenerative, external and internal damping, rotating gyroscopic and inertial force. The results obtained in time domain were compared to the results calculated in frequency domain, and the comparison showed that both results in time and frequency domains were basically consistent.
The following conclusions can be obtained from the analyses: (1) Compared with Kim's model, the model in this study considers not only the gyroscopic effect, internal damping and external damping of the rotating milling bar, but also the effect of the rotary inertia. Moreover, the calculation of this model is simple and convenient, and the results can meet certain accuracy requirements. These advantages are the innovation of this paper. (2) Based on the analysis and calculation results in this paper, the use of composite materials instead of steel for the structural design of milling cutter bar can improve the stability (by 15%) and machining efficiency (increasing cutting speed) to a certain extent. (3) The internal damping affects the stabilities, and the effect increases with the increase of rotating speed. At a high cutting speed, an unstable zone is triggered by internal damping. (4) The amplification of external damping leads to the expansion of the stable zone. (5) Ply angle has a great influence on the stable zone and the largest stable zone is reached at the ply angle of 0 . (6) The rotation effect plays a negative damping role in the stability analysis.
In brief, it has demonstrated that application of composite milling bar can partially improve the chatter stability.

Conclusions
A continuous dynamic model was developed for a rotating composite milling bar. First, the natural frequencies of free vibration were calculated by Galerkin method. Then the stabilities of milling bar subject to deflection-dependent cutting end loads were analyzed in time domain. Furthermore, the model included delay effects due to regenerative, external and internal damping, rotating gyroscopic and inertial force. The results obtained in time domain were compared to the results calculated in frequency domain, and the comparison showed that both results in time and frequency domains were basically consistent.
The following conclusions can be obtained from the analyses: (1) Compared with Kim's model, the model in this study considers not only the gyroscopic effect, internal damping and external damping of the rotating milling bar, but also the effect of the rotary inertia. Moreover, the calculation of this model is simple and convenient, and the results can meet certain accuracy requirements. These advantages are the innovation of this paper. (2) Based on the analysis and calculation results in this paper, the use of composite materials instead of steel for the structural design of milling cutter bar can improve the stability (by 15%) and machining efficiency (increasing cutting speed) to a certain extent. (3) The internal damping affects the stabilities, and the effect increases with the increase of rotating speed. At a high cutting speed, an unstable zone is triggered by internal damping. (4) The amplification of external damping leads to the expansion of the stable zone. (5) Ply angle has a great influence on the stable zone and the largest stable zone is reached at the ply angle of 0. (6) The rotation effect plays a negative damping role in the stability analysis.
In brief, it has demonstrated that application of composite milling bar can partially improve the chatter stability.