# Chatter Stability Prediction Method of the Spindle-Tool Holder-Tool System with Interface Contact Characteristics

Authors: Li Cui; Yin Su
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

To predict chatter stability and suppress chatter vibration, a chatter stability prediction method for the spindle-tool holder-tool system with interface contact characteristics is constructed. A five-DOF model is constructed to determine the spindle-bearing interface dynamic contact stiffness considering the coupling effect of spindle and bearing. A fractal multiscale tool holder-spindle interface dynamic stiffness model is proposed considering time-varying cutting force. The fractal dimensions and cutting force coefficient parameters are identified from the power spectrum experiments and cutting force tests. The cutting force is solved according to the milling stability model. Dynamic model of the spindle-tool holder-tool system is found by the finite element method. Based on extended Floquet theory, chatter stability of the spindle system is studied. Effect of interface parameters, radial cutting depth, and feed rate on milling stability is researched. Milling force tests and milling stability tests are performed in order to verify the reliability of the method. Results find that the increase of front bearing preload and spindle-tool holder’s interference fit are effective to improve the milling stability. The optimal feed rate and the critical radial cutting depth are found. The model proposed in this paper can be used as an instruction for predicting and suppressing the chatter vibration and optimizing cutting parameters and also is helpful for designing the spindle-tool holder-tool system.


---

## Body

## 1. Introduction

Stable cutting is the most important condition for high-speed cutting. However, in actual machine processing, chatter vibration often occurs, which indicates that severe vibration exists. Chatter not only limits the improvement of machining efficiency but also causes great damage to machine tools and cutting tools. Therefore, controlling the generation of chatter and reducing the impact of chatter has become very important. The contact interfaces of the system are the essential factor to influence chatter stability.Dynamic characteristics of spindle-bearing interface had been done by many researchers. A general spindle-bearing model was found according to Timoshenko’s beam and Jones’ bearing model [1]. Bearing stiffness and displacement response are evaluated according to the model. A ball bearing dynamic model considering three degrees of freedom of the balls and surface defects is proposed [2]. Bearing stiffness matrix was solved by Jacobian matrix. Then, stability of the system could be analyzed based on the bearing stiffness model [3]. However, the coupling effect between the spindle and the supported bearing is not considered in the above research studies.Contact stiffness between the spindle and the tool holder was experimentally studied by many researchers. A stiffness model was given based on the spindle system’s modal [4]. However, the recognition result is influenced by test noise seriously, and it is difficult to reflect accurate relationship of stiffness and the surface state. A frequency response calculation model for spindle system is proposed, but the model relied on an impact device [5]. These methods need to rely on experimental results, and the unevenness of contact force cannot be considered, which is not suitable for the design process of the system. Recently, the fractal model has been presented to calculate dynamic characteristics of the system. The fractal contact model was further improved for the spindle-tool holder-tool system [6]. Research found that the single rough region depends on the ratio and the rough plastic-to-elastic mode transition. The tool holder interface contact stiffness was identified by an approximate model, in which only elastic deformation of rough asperities was taken into consideration. The stiffness was solved from every asperity’s stiffness [7]. The results found that cutting force had great effect on stiffness of spindle-tool holder. A spindle-tool holder stiffness model was established by fractal theory. The elastic, elastic-plastic, and plastic deformation were calculated [8]. The cutting force was considered in the model; however, the cutting force was a constant value. In the above model, the real cutting load cannot be taken into account. To calculate the contact properties of tool holder interface, analytical model of cutting force should be constructed. Determination of cutting force coefficients is critical for the prediction of cutting forces. Cutting tests should be done to calculate instantaneous cutting force coefficients [9, 10].Chatter stability is one of the biggest threats to surface quality during the cutting process. Frequency response function (FRF) of tool tip is often used to demonstrate and suppress the chatter stability [11]. A prediction model of tool tip’s frequency response functions was presented and verified by experiment [12]. A model to suppress chatter stability through implementing spindle-tool holder vibration modes was presented by FRF modification [13]. However, accurate FRF depends on experiment, which is not suitable for design process of the spindle system. Stability lobe diagram can predict the cutting stability boundary and divide the stable cutting area and the unstable cutting area. The suitable cutting parameters could be chosen by the stability lobes. In order to obtain the stability lobe diagram of machining center accurately, the dynamic model of the spindle-tool holder-tool system should be constructed. The chatter stability prediction model for spindle-tool holder system was proposed, in which dynamic parameters of the system were achieved by modal identification test [14]. Then, dynamic stability figure was established considering the influence of speed on dynamic behavior. Experiments were done to validate chatter boundary. The spindle-bearing systems dynamic model considering cutting forces was found. A quasi-dynamic model of angular contact ball bearing was used to derive the bearing stiffness. The milling chatter stability of spindle system was studied [15]. The prediction model was experimentally verified by measured displacement response. Performance of the spindle was investigated using the model. Results showed that a rigid preload mechanism is beneficial for the spindle system in high-speed condition. A chatter vibration calculation model for spindle system was established, in which ball bearings and floating bearings dynamic model were considered [16]. The spindle system was modeled by using Timoshenko beam. Stability lobe and chatter stability were researched under different cutting parameters. However, the interface model between tool holder and spindle is not given. Above all, chatter stability prediction model of tool-tool holder-spindle system based on theoretical calculation has not been performed. Deep research should be done to predict chatter stability according to the interface contact characteristics.This paper proposes a new chatter stability prediction method for tool-tool holder-spindle system. Contact model of spindle-bearing interface and spindle-tool holder were constructed. Effect of contact interface parameters, radial cutting depth, and feed rate on milling stability is researched; the results were experimentally verified by comparison with the cutting experiment. The method in this paper is helpful for predicting, suppressing chatter vibration of spindle system in the design process.


## 2. Chatter Stability Prediction Method

### 2.1. Dynamic Contact Stiffness of Spindle-Ball Bearing Interface

To study contact characteristics of ball bearing-spindle interface, a five-DOF model is constructed to determine the contact stiffness. Figure1 illustrates geometry of an angular ball bearing, and the initial contact angle is α, as shown in Figures 1(a) and 1(b), which gives the curvature center position of raceways and the jth ball. Oj, Oi, Oe are the initial position of the jth ball, curvature center of inner, and outer raceway, respectively. Oj′, Oi′, Oe′ are the loading position of the jth ball and curvature center of inner and outer raceway. δa is preload displacement of the outer ring.Figure 1
High-speed ball bearing geometric relationship. (a) Balling bearing under preload. (b) Position of raceways and thejth ball.
(a)(b)It is assumed that the inner ring displacement isXs,Ys,Zs,ϕXs,ϕYs, which are equal to the displacements at the same location on the spindle. The jth ball angular location is θj, and displacement of the jth ball is xj,yj,zj. Then, contact deformation δ1j, δ2j and contact angles α1j, α2j can be deduced from Figure 1.(1)δ1j=Asinα0+δa+zj2+Acosα0+xjsinθj+yjcosθj2−A−h1j,δ2j=Bsinα0+Zs−zj+ϕYssinθj+ϕXscosθj2+Bcosα0+Xs−xjsinθj+Ys−yjcosθj2−B−h2j,α1j=arctanAsinα0+δa+zjAcosα0+xjsinθj+yjcosθj,α2j=arctanBsinα0+Zs−zj+ϕYssinθj+ϕZscosθjBcosα0+X−xjsinθj+Y−yjcosθj,where A=f1−0.5Dw and B=f2−0.5Dw. Dw is the diameter of the ball. f1 and f2are the curvature coefficient of outer raceway and inner raceway. h1j and h2j are the oil film thickness between outer raceway, inner raceway, and the jth ball [17].Then, contact loadsQ1j and Q2j between the outer ring, inner ring, and the jth ball can be solved [17]. Assuming bearing load is FX,FY,FZ,MX,MY, inner ring equations are written as(2)FX−∑j=1zbQ2jcosα2jsinθj=0,FY−∑j=1zbQ2jcosα2jcosθj=0,FZ−∑j=1zbQ2jsinα2j=0,MX−∑j=1zbQ2jsinα2jcosθj=0,MY−∑j=1zbQ2jcosα2jsinθj=0,where zb is the ball number. θj=2π/zbj−1+1/2ω1−Dw/Dpwcosα×t; Dpw is the pitch diameter. ω is the rotating speed.The above equations are solved by Newton–Raphson method. The Jacobian matrix obtained by solving the equations can be defined as the stiffness matrix of the bearing:(3)Kb=∂Fi∂Xj,where Fi=FX,FY,FZ,MX,MY, Xj=Xs,XjYs,Zs,ϕXs,ϕYs, and Kb is bearing stiffness matrix.Because the displacementsXj are varied when the spindle-bearing system is running, dynamic contact stiffness between spindle and ball bearing interface can be solved by equation (3).


### 2.2. Dynamic Contact Stiffness Model of Tool Holder-Spindle Interface

Contact between tool holder and spindle can be seen as a plane stress problem. Figure2 is cross-sectional sketch. Tool holder is divided into nl sections. Elastic modulus, Poisson’s ratio, and density of them are considered to be the same.Figure 2
Sketch of spindle-tool holder interface.Based on classic elasticity theory [18], radial displacements of spindle and the kth slide of tool holder are given:(4)uSk=3+υ4Eρω2rkrs2+1−ν3+νrk2,uTk=1−υ4Eρω2rk3,where uSk is the radial displacement of the spindle at the kth slide location. uTk is the radial displacement of the kth slide of the tool holder. E is the modulus. ρ is the density. υ is Poisson’s ratio. rs is the spindle’s radius. rk is the radius of the kth slide of the tool holder, which is varied from r1 to r2.The actual interference between the spindle and thekth slide can be written as(5)δk=δ0+uTk−uSk−δ′,where δ0 is the initial interference from drawbar load. δ′ is the displacement from the cutting force.According to elastic mechanics theory, the contact stress is(6)σk=Eδkrr2−rk22rkrs2.For the tool holder, the equilibrium equations are expressed as(7)Fy=∑k=1nlσkrklTnlcosβ,(8)Fz=∑k=1nlσkrklTnlsinβ,where Fy is the cutting force in the y-direction.Fz is the drawbar force. β is the taper angle. Solving the above equations by the Newton iterative method, δ′ and σk are obtained.According to fractal theory, a 2D multiscale surface profile is given based on the W-M function [19].(9)Hx=GD−1∑nminnmaxcos2πxγnγ2−Dn1<D<2,where Hx is the height of the profile. G is a characteristic length scale. D is the fractal parameter. γ=1.5 and γnmin=1/L [3]. L is the sample length. nmin and nmax are the minimum and maximum frequency.The critical elastic contact area is given [20].(10)aec=G229−2Dπ3−DEKH2lnγ1/D−1,where K stands for the hardness coefficient. K=0.454+0.41υ. H stands for hardness.The maximum contact area of the asperities is assumed asaL. If aL>aec, elastic contact load is given for an asperity with contact area a.(11)Fce=43×25/2−D×πD−3/2×E×GD−1×lnγ1/2×a3−D/2.The normal contact stiffness of a single asperity can be written as(12)kne=dFcedδ=4E3−D32π2−Da1/2.Elastoplastic deformation of asperities can be regarded as two phases [21]. If aepc≤aL<aec, aepcis the minimum contact area of the first phase.(13)aepc=G229−2D6π3−DlnγEKH21/D−1.Contact load and normal contact stiffness of first phase elastoplastic deformation are given as(14)Fcep1=0.34×23.825−0.85DE0.85G0.85D−1KH0.15π0.425D−1.275lnγ0.425a1.425−0.425D,knep1=dFcep1dδ=1.03×21.825+0.15DKH0.15E0.85a0.425+0.075D1.425−0.425D3π0.075D+0.275G0.15D−1lnγ0.0752−D.Ifapc≤aL<aepc, apcis the minimum contact area of the second phase.(15)apc=G229−2D110π3−DEKH2lnγ1/D−1.Contact load and normal contact stiffness of second phase elastoplastic deformation are given as(16)Fcep2=1.43×22.367−0.526DE0.526G0.526D−1KH0.474π0.263D−0.789lnγ0.263a1.263−0.263D,(17)knep2=dFcep2dδ=1.4×20.367+0.474DKH0.474E0.526a0.237D+0.2631.263−0.263D3π0.237D−0.211G0.474D−1lnγ0.2372−D.The statistical distribution of the asperities can be expressed as [22](18)na=D2η2−D/2aLD/2a−D/2+1.The real contact area is solved as(19)Ar=∫0aLnaada,where η is the domain extension factor.IfaL>aec, interface contact load of the kth slide is(20)Fk=∫apcaepcnaFcep2da+∫aepcaecnaFcep1da+∫aecaLnaFceda.Ifaepc≤aL<aec, contact load of the kth slide is(21)Fk=∫apcaepcnaFcep2da+∫aepcaLnaFcep1da.Ifapc≤aL<aepc, contact load of the kth slide is(22)Fk=∫apcaLnaFcep2da.The following equation can be obtained according to the contact load balance of thekth slide of tool holder-spindle interface:(23)Fk=Ar·σk.ParameteraL can be solved from the above equations. The normal interface contact stiffness of the kth slide can be solved by integral of every asperity. If aL>aec, normal contact stiffness of the kth slide is(24)KkN=∫apcaepcnaknep2da+∫aepcaecnaknep1da+∫aecaLnakneda.Ifaepc≤aL<aec, normal contact stiffness of the kth slide can be written as(25)KkN=∫apcaepcnaknep2da+∫aepcaLnaknep1da.Ifapc≤aL<aepc, the following expression is obtained:(26)KKN=∫apcaLnaknep2da.The total normal interface contact stiffnessKn is solved by every slide’s contact stiffness. Similarly, the total shear interface stiffness Kt can also be solved. Then, radial stiffness Kν and torsional stiffness Kϑ can be determined.(27)Kν=∑k=1nlKncosβ+Ktsinβ,Kϑ=∑k=1nlKnsinβ+Ktcosβrk.Similarly, the stiffness between the tool holder and the tool can also be solved by the same method.From equations (7) and (8), because cutting force is varied when the spindle-bearing system is running, radial stiffness Kν and torsional stiffness Kϑ of the tool holder are dynamic contact stiffness.


### 2.3. Cutting Force Model

According to the milling stability model [23], the instantaneous dynamic cutting thickness can be written as(28)hφj=xsinφj+ycosφjgφj,where x and y stand for the dynamic displacement of two adjacent cutter teeth in X- and Y-directions.(29)gφj=1,φst<φj<φex,0,φj<φstorφj>φex,where φst is cutting angle; φex is cut-out angle. Cutting forces of the jth tooth are(30)Ftj=Ktcaphφj,Frj=KrcFtj,where Ftjis the tangential force. Frj is the radial force. ap is the axial cutting depth. Ktc and Krc are tangential and radial cutting force coefficients. Assuming periodic average of cutting forces Fx¯ and Fy¯ are measured by milling test, then the relationship of the cutting force coefficients and average cutting force are expressed as(31)Fx¯=ap4ωKrcVf+NapπKre,(32)Fy¯=ap4ωKtcVf+NapπKte,where N is the cutter teeth number. Vf is the feed rate. Kre and Kte are radial and tangential edge force coefficients.Cutting forces inX- and Y-directions are written as(33)FxFy=−aphxx−aphxy−aphyx−aphyyxt−xt−Tyt−yt−T,hxx=∑j=1NgφjsinφjKctcosφj+Kcrsinφj,hxy=∑j=1NgφjcosφjKctcosφj+Kcrsinφj,hyx=∑j=1Ngφjsinφj−Kctsinφj+Kcrcosφj,hyy=∑j=1Ngφjcosφj−Kctsinφj+Kcrcosφj,where T is the tooth cutting cycle.


### 2.4. Chatter Stability Prediction Method of the Spindle-Tool Holder-Tool System

Figure3 shows the distribution spring-damping model of spindle-tool holder interface. Assuming the interface is divided into n elements, there are n+1 nodes, i,i+1,…i+n, on the spindle and j,j+1,…j+n on the tool holder, respectively.Figure 3
Distribution spring-damping model.Stiffness coupling matrix of the interface can be expressed as(34)where KST stands for the stiffness matrix of one single spring-damping unit.(35)KST=Kv000000Kv000000000000000000000000000Kϑ.Damping matrix of the interface includes radial structural damping and torsional structural damping and can be determined according to [24].The spindle system consists of CFV Xi motorized spindle, BT-40 holder, and HSS-Co8 square 4 flute end milling cutter. Figure4 gives the FEM model of the system. There are 38 Euler–Bernoulli elements and 41 nodes in the FEM model. The cutting force acts on the 41st node. Front bearings are located at 14th and 15th nodes. Back bearings are located at 3rd and 4th nodes. Spindle-tool holder interface is located from 17th to 21st nodes and from 23rd to 27th nodes. Tool holder-tool interface is located from 30th to 33rd nodes and from 35th to 38th nodes.Figure 4
FEM model of the spindle-tool holder-tool bearing system.Dynamic equation of the system can be given based on the finite element method as follows:(36)MX¨+C+JX˙+KX=Pt+G,X=χ1,χ2,…χ41T,χi=xi,yi,zi,θxi,θyi,θziT,i=1,2,…,41,where M, K, C, and J are mass, stiffness, damping, and gyroscopic matrix, respectively. P, G, and X are cutting force, gravity, and displacement vector, respectively.These matrices are assembled by the integration method. Taking stiffness matrix as an example, Figure5 gives the flow chart of dynamic stiffness assembly process. First, dynamic stiffness matrix of spindle, tool holder, and tool are solved through the Euler–Bernoulli element model. Then, the node stiffness matrix is assembled into the global matrix starting from the first node. If there is bearing at this node, then the bearing stiffness matrix is assembled to the corresponding node in the global stiffness matrix. If there is spindle-holder interface at this node, by using the same method, the holder-tool coupling stiffness matrix can be assembled into the global stiffness matrix. Then, the spindle-holder coupling stiffness matrix is assembled to the corresponding node in the global stiffness matrix.Figure 5
Dynamic stiffness matrix assembly of spindle-tool holder-tool system.Solving the FEM dynamic equations by using the Runge–Kutta method, nonlinear responses at tool tip can be obtained. Periodic motion bifurcation is solved by continuation-shooting method [25]. According to the extended Floquet theory, chatter stability of the response is analyzed [26]. By increasing the cutting depth, the critical stability condition can be found when one of the eigenvalues lies out of the unit circle or one of Floquet multipliers absolute value is larger than 1. Repeating solving the solution and judging the eigenvalues, the stability lobes diagram can be obtained.


## 2.1. Dynamic Contact Stiffness of Spindle-Ball Bearing Interface

To study contact characteristics of ball bearing-spindle interface, a five-DOF model is constructed to determine the contact stiffness. Figure1 illustrates geometry of an angular ball bearing, and the initial contact angle is α, as shown in Figures 1(a) and 1(b), which gives the curvature center position of raceways and the jth ball. Oj, Oi, Oe are the initial position of the jth ball, curvature center of inner, and outer raceway, respectively. Oj′, Oi′, Oe′ are the loading position of the jth ball and curvature center of inner and outer raceway. δa is preload displacement of the outer ring.Figure 1
High-speed ball bearing geometric relationship. (a) Balling bearing under preload. (b) Position of raceways and thejth ball.
(a)(b)It is assumed that the inner ring displacement isXs,Ys,Zs,ϕXs,ϕYs, which are equal to the displacements at the same location on the spindle. The jth ball angular location is θj, and displacement of the jth ball is xj,yj,zj. Then, contact deformation δ1j, δ2j and contact angles α1j, α2j can be deduced from Figure 1.(1)δ1j=Asinα0+δa+zj2+Acosα0+xjsinθj+yjcosθj2−A−h1j,δ2j=Bsinα0+Zs−zj+ϕYssinθj+ϕXscosθj2+Bcosα0+Xs−xjsinθj+Ys−yjcosθj2−B−h2j,α1j=arctanAsinα0+δa+zjAcosα0+xjsinθj+yjcosθj,α2j=arctanBsinα0+Zs−zj+ϕYssinθj+ϕZscosθjBcosα0+X−xjsinθj+Y−yjcosθj,where A=f1−0.5Dw and B=f2−0.5Dw. Dw is the diameter of the ball. f1 and f2are the curvature coefficient of outer raceway and inner raceway. h1j and h2j are the oil film thickness between outer raceway, inner raceway, and the jth ball [17].Then, contact loadsQ1j and Q2j between the outer ring, inner ring, and the jth ball can be solved [17]. Assuming bearing load is FX,FY,FZ,MX,MY, inner ring equations are written as(2)FX−∑j=1zbQ2jcosα2jsinθj=0,FY−∑j=1zbQ2jcosα2jcosθj=0,FZ−∑j=1zbQ2jsinα2j=0,MX−∑j=1zbQ2jsinα2jcosθj=0,MY−∑j=1zbQ2jcosα2jsinθj=0,where zb is the ball number. θj=2π/zbj−1+1/2ω1−Dw/Dpwcosα×t; Dpw is the pitch diameter. ω is the rotating speed.The above equations are solved by Newton–Raphson method. The Jacobian matrix obtained by solving the equations can be defined as the stiffness matrix of the bearing:(3)Kb=∂Fi∂Xj,where Fi=FX,FY,FZ,MX,MY, Xj=Xs,XjYs,Zs,ϕXs,ϕYs, and Kb is bearing stiffness matrix.Because the displacementsXj are varied when the spindle-bearing system is running, dynamic contact stiffness between spindle and ball bearing interface can be solved by equation (3).


## 2.2. Dynamic Contact Stiffness Model of Tool Holder-Spindle Interface

Contact between tool holder and spindle can be seen as a plane stress problem. Figure2 is cross-sectional sketch. Tool holder is divided into nl sections. Elastic modulus, Poisson’s ratio, and density of them are considered to be the same.Figure 2
Sketch of spindle-tool holder interface.Based on classic elasticity theory [18], radial displacements of spindle and the kth slide of tool holder are given:(4)uSk=3+υ4Eρω2rkrs2+1−ν3+νrk2,uTk=1−υ4Eρω2rk3,where uSk is the radial displacement of the spindle at the kth slide location. uTk is the radial displacement of the kth slide of the tool holder. E is the modulus. ρ is the density. υ is Poisson’s ratio. rs is the spindle’s radius. rk is the radius of the kth slide of the tool holder, which is varied from r1 to r2.The actual interference between the spindle and thekth slide can be written as(5)δk=δ0+uTk−uSk−δ′,where δ0 is the initial interference from drawbar load. δ′ is the displacement from the cutting force.According to elastic mechanics theory, the contact stress is(6)σk=Eδkrr2−rk22rkrs2.For the tool holder, the equilibrium equations are expressed as(7)Fy=∑k=1nlσkrklTnlcosβ,(8)Fz=∑k=1nlσkrklTnlsinβ,where Fy is the cutting force in the y-direction.Fz is the drawbar force. β is the taper angle. Solving the above equations by the Newton iterative method, δ′ and σk are obtained.According to fractal theory, a 2D multiscale surface profile is given based on the W-M function [19].(9)Hx=GD−1∑nminnmaxcos2πxγnγ2−Dn1<D<2,where Hx is the height of the profile. G is a characteristic length scale. D is the fractal parameter. γ=1.5 and γnmin=1/L [3]. L is the sample length. nmin and nmax are the minimum and maximum frequency.The critical elastic contact area is given [20].(10)aec=G229−2Dπ3−DEKH2lnγ1/D−1,where K stands for the hardness coefficient. K=0.454+0.41υ. H stands for hardness.The maximum contact area of the asperities is assumed asaL. If aL>aec, elastic contact load is given for an asperity with contact area a.(11)Fce=43×25/2−D×πD−3/2×E×GD−1×lnγ1/2×a3−D/2.The normal contact stiffness of a single asperity can be written as(12)kne=dFcedδ=4E3−D32π2−Da1/2.Elastoplastic deformation of asperities can be regarded as two phases [21]. If aepc≤aL<aec, aepcis the minimum contact area of the first phase.(13)aepc=G229−2D6π3−DlnγEKH21/D−1.Contact load and normal contact stiffness of first phase elastoplastic deformation are given as(14)Fcep1=0.34×23.825−0.85DE0.85G0.85D−1KH0.15π0.425D−1.275lnγ0.425a1.425−0.425D,knep1=dFcep1dδ=1.03×21.825+0.15DKH0.15E0.85a0.425+0.075D1.425−0.425D3π0.075D+0.275G0.15D−1lnγ0.0752−D.Ifapc≤aL<aepc, apcis the minimum contact area of the second phase.(15)apc=G229−2D110π3−DEKH2lnγ1/D−1.Contact load and normal contact stiffness of second phase elastoplastic deformation are given as(16)Fcep2=1.43×22.367−0.526DE0.526G0.526D−1KH0.474π0.263D−0.789lnγ0.263a1.263−0.263D,(17)knep2=dFcep2dδ=1.4×20.367+0.474DKH0.474E0.526a0.237D+0.2631.263−0.263D3π0.237D−0.211G0.474D−1lnγ0.2372−D.The statistical distribution of the asperities can be expressed as [22](18)na=D2η2−D/2aLD/2a−D/2+1.The real contact area is solved as(19)Ar=∫0aLnaada,where η is the domain extension factor.IfaL>aec, interface contact load of the kth slide is(20)Fk=∫apcaepcnaFcep2da+∫aepcaecnaFcep1da+∫aecaLnaFceda.Ifaepc≤aL<aec, contact load of the kth slide is(21)Fk=∫apcaepcnaFcep2da+∫aepcaLnaFcep1da.Ifapc≤aL<aepc, contact load of the kth slide is(22)Fk=∫apcaLnaFcep2da.The following equation can be obtained according to the contact load balance of thekth slide of tool holder-spindle interface:(23)Fk=Ar·σk.ParameteraL can be solved from the above equations. The normal interface contact stiffness of the kth slide can be solved by integral of every asperity. If aL>aec, normal contact stiffness of the kth slide is(24)KkN=∫apcaepcnaknep2da+∫aepcaecnaknep1da+∫aecaLnakneda.Ifaepc≤aL<aec, normal contact stiffness of the kth slide can be written as(25)KkN=∫apcaepcnaknep2da+∫aepcaLnaknep1da.Ifapc≤aL<aepc, the following expression is obtained:(26)KKN=∫apcaLnaknep2da.The total normal interface contact stiffnessKn is solved by every slide’s contact stiffness. Similarly, the total shear interface stiffness Kt can also be solved. Then, radial stiffness Kν and torsional stiffness Kϑ can be determined.(27)Kν=∑k=1nlKncosβ+Ktsinβ,Kϑ=∑k=1nlKnsinβ+Ktcosβrk.Similarly, the stiffness between the tool holder and the tool can also be solved by the same method.From equations (7) and (8), because cutting force is varied when the spindle-bearing system is running, radial stiffness Kν and torsional stiffness Kϑ of the tool holder are dynamic contact stiffness.


## 2.3. Cutting Force Model

According to the milling stability model [23], the instantaneous dynamic cutting thickness can be written as(28)hφj=xsinφj+ycosφjgφj,where x and y stand for the dynamic displacement of two adjacent cutter teeth in X- and Y-directions.(29)gφj=1,φst<φj<φex,0,φj<φstorφj>φex,where φst is cutting angle; φex is cut-out angle. Cutting forces of the jth tooth are(30)Ftj=Ktcaphφj,Frj=KrcFtj,where Ftjis the tangential force. Frj is the radial force. ap is the axial cutting depth. Ktc and Krc are tangential and radial cutting force coefficients. Assuming periodic average of cutting forces Fx¯ and Fy¯ are measured by milling test, then the relationship of the cutting force coefficients and average cutting force are expressed as(31)Fx¯=ap4ωKrcVf+NapπKre,(32)Fy¯=ap4ωKtcVf+NapπKte,where N is the cutter teeth number. Vf is the feed rate. Kre and Kte are radial and tangential edge force coefficients.Cutting forces inX- and Y-directions are written as(33)FxFy=−aphxx−aphxy−aphyx−aphyyxt−xt−Tyt−yt−T,hxx=∑j=1NgφjsinφjKctcosφj+Kcrsinφj,hxy=∑j=1NgφjcosφjKctcosφj+Kcrsinφj,hyx=∑j=1Ngφjsinφj−Kctsinφj+Kcrcosφj,hyy=∑j=1Ngφjcosφj−Kctsinφj+Kcrcosφj,where T is the tooth cutting cycle.


## 2.4. Chatter Stability Prediction Method of the Spindle-Tool Holder-Tool System

Figure3 shows the distribution spring-damping model of spindle-tool holder interface. Assuming the interface is divided into n elements, there are n+1 nodes, i,i+1,…i+n, on the spindle and j,j+1,…j+n on the tool holder, respectively.Figure 3
Distribution spring-damping model.Stiffness coupling matrix of the interface can be expressed as(34)where KST stands for the stiffness matrix of one single spring-damping unit.(35)KST=Kv000000Kv000000000000000000000000000Kϑ.Damping matrix of the interface includes radial structural damping and torsional structural damping and can be determined according to [24].The spindle system consists of CFV Xi motorized spindle, BT-40 holder, and HSS-Co8 square 4 flute end milling cutter. Figure4 gives the FEM model of the system. There are 38 Euler–Bernoulli elements and 41 nodes in the FEM model. The cutting force acts on the 41st node. Front bearings are located at 14th and 15th nodes. Back bearings are located at 3rd and 4th nodes. Spindle-tool holder interface is located from 17th to 21st nodes and from 23rd to 27th nodes. Tool holder-tool interface is located from 30th to 33rd nodes and from 35th to 38th nodes.Figure 4
FEM model of the spindle-tool holder-tool bearing system.Dynamic equation of the system can be given based on the finite element method as follows:(36)MX¨+C+JX˙+KX=Pt+G,X=χ1,χ2,…χ41T,χi=xi,yi,zi,θxi,θyi,θziT,i=1,2,…,41,where M, K, C, and J are mass, stiffness, damping, and gyroscopic matrix, respectively. P, G, and X are cutting force, gravity, and displacement vector, respectively.These matrices are assembled by the integration method. Taking stiffness matrix as an example, Figure5 gives the flow chart of dynamic stiffness assembly process. First, dynamic stiffness matrix of spindle, tool holder, and tool are solved through the Euler–Bernoulli element model. Then, the node stiffness matrix is assembled into the global matrix starting from the first node. If there is bearing at this node, then the bearing stiffness matrix is assembled to the corresponding node in the global stiffness matrix. If there is spindle-holder interface at this node, by using the same method, the holder-tool coupling stiffness matrix can be assembled into the global stiffness matrix. Then, the spindle-holder coupling stiffness matrix is assembled to the corresponding node in the global stiffness matrix.Figure 5
Dynamic stiffness matrix assembly of spindle-tool holder-tool system.Solving the FEM dynamic equations by using the Runge–Kutta method, nonlinear responses at tool tip can be obtained. Periodic motion bifurcation is solved by continuation-shooting method [25]. According to the extended Floquet theory, chatter stability of the response is analyzed [26]. By increasing the cutting depth, the critical stability condition can be found when one of the eigenvalues lies out of the unit circle or one of Floquet multipliers absolute value is larger than 1. Repeating solving the solution and judging the eigenvalues, the stability lobes diagram can be obtained.


## 3. Analysis and Discussion

The spindle system consists of CFV Xi motorized spindle, BT-40 holder, and HSS-Co8 square 4 flute end milling cutter. The workpiece material is aluminum alloy 6061-T6. The dynamometer is installed on the machine table. The cutting force acting on the workpiece can be measured. YOLO-YDC01 three-way dynamometer and a milling force measurement system are used. Figure6(a) shows the setup of the test system. Milling force coefficients can be measured by milling experiment. When spindle speed is 4000 r/min, axial cutting depth is 1 mm, radial cutting depth is 6 mm, and the average milling forces are measured by the test system. Figure 6(b) shows the linear regression diagram of average milling forces in X- and Y-directions. Then, tangential and radial cutting force coefficients are solved according to equations (31) and (32).(37)Ktc=0.073×4×40001=1168N/mm2,Krc=0.051×4×40001=816N/mm2.Figure 6
Setup of milling forces and stability test system. (a) Setup of the test system. (b) Test results of average milling forces.
(a)(b)Similarly, different tangential and radial cutting force coefficients can be obtained by changing the spindle speed, axial cutting depth, radial cutting depth, and feed rate.The fractal parameters of BT-40 holder can be obtained by power spectrum method [8]. The micromorphology of rough surfaces was measured by a laser microscope; then the data are used to identify the fractal parameters. Power spectrum of the measured morphology can be expressed as [27](38)Sχ=G2D−12lnγ1χ5−2D,where Sχ is the power of the spectrum and χ is the reciprocal of the wavelength of roughness. Expressing the previous equation by a function of log-log plot, equation (38) can be written as(39)lgSχ=2D−1lgG−lg2lnγ−5−2Dlgχ.It can be seen that the power law behavior is a straight line, as shown in Figure7. The power spectrum can be written as a linear function.(40)lgSχ=−1.53lgχ−12.627.Figure 7
Power spectrum of the tool holder rough surfaces.Then the fractal dimension parameters can be identified asD = 1.735, G=2.23×10−9m.


### 3.1. Influence of Interface Parameters on Milling Stability

Spindle rotating speed changes from 2000 r/min to 10,000 r/min. Radial cutting depth is set as 6 mm, and the feed rate remains 800 mm/min. The preloads of front bearings and back bearings are set as 8μm. The initial interference fit between spindle and tool holder interface is 8 μm. Figure 8 shows that Floquet multipliers vary with axial cutting depth. Figure 8(a) gives changes of tool tip Floquet multipliers at 3000 r/min when axial cutting depth varies from 2.2 mm to 2.6 mm. The maximum Floquet multiplier comes out of the unit circle by a pair of conjugate complex number at 2.6 mm, which explains that the system becomes unstable. It can be demonstrated that the critical axial cutting depth is 2.6 mm at 3000 r/min. Similarly, Figure 8(b) demonstrates that the critical axial cutting depth is 3.1 mm at 5000 r/min.Figure 8
Floquet multipliers vary with axial cutting depth: (a) 3000 r/min; (b) 5000 r/min.
(a)(b)Solving the Floquet multipliers under different rotating speeds, the chatter stability lobes diagram can be obtained, as shown in Figure9(a). Figure 9(b) gives the influence of front bearings’ preload on the critical axial cutting depth. As front bearings’ preload increases, the stable milling depth increases. The increase of preload of the front bearings is beneficial to improve the milling stability. However, Figure 9(c) shows stable milling depth decrease as back bearings’ preload increases. The increase of preload of the back bearings is not conducive to improving the milling stability. Figure 9(d) illustrates that stable milling depth increases as initial interference fit increases. The increase of spindle-tool holder’s interference fit is significantly beneficial to improve the milling stability.Figure 9
Influence of interface parameters on chatter stability. (a) Original interface parameter. (b) Front bearing interface parameter changes. (c) Back bearing interface parameter changes. (d) Tool holder interface parameter changes.
(a)(b)(c)(d)


### 3.2. Influence of Feed Rate on Milling Stability

Spindle rotating speed varies from 2000 r/min to 10,000 r/min. Radial cutting depth is set as 6 mm, and feed rate is set as 400 mm/min, 800 mm/min, and 1200 mm/min, respectively. Figure10 shows the critical axial cutting depth varies with feed rate. As feed rate increases, the stable milling depth decreases first and then increases with the increase of spindle speed. In order to improve milling stability, feed rate should change with the spindle speed.Figure 10
Influence of feed rate on chatter stability.Figure11 gives the influence of feed rate on critical axial cutting depth. Figure 11 demonstrates that the critical axial cutting depth increases first and then decreases with the increase of feed rate. In order to maximize the critical axial cutting depth, there exists an optimal feed rate corresponding to a spindle speed. When spindle speed is 3000 r/min, 5000 r/min, and 8000 r/min, the optimal feed rate is 600 mm/min, 1200 mm/min, and 1800 mm/min, respectively. It can be concluded that feed rate should match the spindle speed and change with the spindle speed.Figure 11
Influence of feed rate on critical axial cutting depth.


### 3.3. Influence of Radial Cutting Depth on Milling Stability

Spindle rotating speed changes from 2000 r/min to 10,000 r/min, and the feed rate is set as 800 mm/min. Radial cutting depth is set as 2 mm, 6 mm, and 10 mm, respectively. Figure12 shows the influence of radial cutting depth on chatter stability. As radial cutting depth changes from 2 mm to 10 mm, the stable milling depth decreases. The increase of radial cutting depth is not conducive to improving milling stability.Figure 12
Influence of radial cutting depth on chatter stability.Figure13 gives the influence of radial cutting depth on critical axial cutting depth. It illustrates that the critical axial cutting depth decreases with the increase of radial cutting depth. Particularly, the critical axial cutting depth decreases sharply when radial cutting depth exceeds 6 mm. To increase the critical axial cutting depth and, simultaneously, ensure cutting efficiency, the radial cutting depth should not exceed the critical value.Figure 13
Influence of radial cutting depth on critical axial cutting depth.


## 3.1. Influence of Interface Parameters on Milling Stability

Spindle rotating speed changes from 2000 r/min to 10,000 r/min. Radial cutting depth is set as 6 mm, and the feed rate remains 800 mm/min. The preloads of front bearings and back bearings are set as 8μm. The initial interference fit between spindle and tool holder interface is 8 μm. Figure 8 shows that Floquet multipliers vary with axial cutting depth. Figure 8(a) gives changes of tool tip Floquet multipliers at 3000 r/min when axial cutting depth varies from 2.2 mm to 2.6 mm. The maximum Floquet multiplier comes out of the unit circle by a pair of conjugate complex number at 2.6 mm, which explains that the system becomes unstable. It can be demonstrated that the critical axial cutting depth is 2.6 mm at 3000 r/min. Similarly, Figure 8(b) demonstrates that the critical axial cutting depth is 3.1 mm at 5000 r/min.Figure 8
Floquet multipliers vary with axial cutting depth: (a) 3000 r/min; (b) 5000 r/min.
(a)(b)Solving the Floquet multipliers under different rotating speeds, the chatter stability lobes diagram can be obtained, as shown in Figure9(a). Figure 9(b) gives the influence of front bearings’ preload on the critical axial cutting depth. As front bearings’ preload increases, the stable milling depth increases. The increase of preload of the front bearings is beneficial to improve the milling stability. However, Figure 9(c) shows stable milling depth decrease as back bearings’ preload increases. The increase of preload of the back bearings is not conducive to improving the milling stability. Figure 9(d) illustrates that stable milling depth increases as initial interference fit increases. The increase of spindle-tool holder’s interference fit is significantly beneficial to improve the milling stability.Figure 9
Influence of interface parameters on chatter stability. (a) Original interface parameter. (b) Front bearing interface parameter changes. (c) Back bearing interface parameter changes. (d) Tool holder interface parameter changes.
(a)(b)(c)(d)


## 3.2. Influence of Feed Rate on Milling Stability

Spindle rotating speed varies from 2000 r/min to 10,000 r/min. Radial cutting depth is set as 6 mm, and feed rate is set as 400 mm/min, 800 mm/min, and 1200 mm/min, respectively. Figure10 shows the critical axial cutting depth varies with feed rate. As feed rate increases, the stable milling depth decreases first and then increases with the increase of spindle speed. In order to improve milling stability, feed rate should change with the spindle speed.Figure 10
Influence of feed rate on chatter stability.Figure11 gives the influence of feed rate on critical axial cutting depth. Figure 11 demonstrates that the critical axial cutting depth increases first and then decreases with the increase of feed rate. In order to maximize the critical axial cutting depth, there exists an optimal feed rate corresponding to a spindle speed. When spindle speed is 3000 r/min, 5000 r/min, and 8000 r/min, the optimal feed rate is 600 mm/min, 1200 mm/min, and 1800 mm/min, respectively. It can be concluded that feed rate should match the spindle speed and change with the spindle speed.Figure 11
Influence of feed rate on critical axial cutting depth.


## 3.3. Influence of Radial Cutting Depth on Milling Stability

Spindle rotating speed changes from 2000 r/min to 10,000 r/min, and the feed rate is set as 800 mm/min. Radial cutting depth is set as 2 mm, 6 mm, and 10 mm, respectively. Figure12 shows the influence of radial cutting depth on chatter stability. As radial cutting depth changes from 2 mm to 10 mm, the stable milling depth decreases. The increase of radial cutting depth is not conducive to improving milling stability.Figure 12
Influence of radial cutting depth on chatter stability.Figure13 gives the influence of radial cutting depth on critical axial cutting depth. It illustrates that the critical axial cutting depth decreases with the increase of radial cutting depth. Particularly, the critical axial cutting depth decreases sharply when radial cutting depth exceeds 6 mm. To increase the critical axial cutting depth and, simultaneously, ensure cutting efficiency, the radial cutting depth should not exceed the critical value.Figure 13
Influence of radial cutting depth on critical axial cutting depth.


## 4. Chatter Stability Experiment

The milling experiments were performed. The workpiece material is 6061-T6 aluminum alloy. Figure14 shows the surface of unprocessed workpiece. The cutting forces acting on the workpiece are measured by dynamometer.Figure 14
Surface of the unprocessed workpiece.The occurrence of chatter in the experiment can be judged by analyzing the milling force and its frequency spectrum. In the milling force spectrum, if a certain frequency is not an integer multiple of the tooth passing frequency (TPF), chatter vibration occurs and this frequency is the chatter frequency (CF). However, if all frequencies are integer multiples of the cutting frequency of the teeth, no chatter vibration occurs. In addition, it can also provide evidence for chattering judgment results by the roughness of machined workpiece surface.The spindle rotating speed is set as 4000 r/min. Radial cutting depth is set as 6 mm, and feed rate remains 800 mm/min. When axial cutting depth is 1 mm, Figure15 gives the comparison of calculated and measured cutting forces in X- and Y-directions. It shows that the calculated results are in good agreement with the measured milling forces and the maximum relative error of the average milling force between the calculated and measured data is 5.7%.Figure 15
Comparison of calculated and measured milling forces. (a) Milling force in theX-direction. (b) Milling force in the Y-direction.
(a)(b)Figure16 gives test results of the frequency spectrum of the cutting force in the Y-direction and workpiece surface under different axial cutting depth. Figure 16(a) illustrates that when axial cutting depth is 1 mm, the system is stable. Figure 16(b) shows a smooth machined workpiece surface. When axial cutting depth increases to 2 mm, Figures 16(c) and 16(d) demonstrate that chatter frequency occurs in the frequency spectrum and the machined workpiece surface becomes rough, which indicates that the cutting system causes slight chatter vibration. When axial cutting depth increases to 3 mm, Figures 16(e) and 16(f) show that more chatter frequencies occur in the frequency spectrum and the machined workpiece surface becomes rougher, which illustrate that the cutting system exists severe chatter vibration.Figure 16
Frequency spectrum of cutting force in theY-direction and workpiece surface. (a) Frequency spectrum when ap = 1 mm. (b) Machined surface when ap = 1 mm. (c) Frequency spectrum when ap = 2 mm. (d) Machined surface when ap = 2 mm. (e) Frequency spectrum when ap = 3 mm. (f) Machined surface when ap = 3 mm.
(a)(b)(c)(d)(e)(f)Figure17 gives the comparison between calculation results and test results of milling stability. When the milling depth is in the area below the calculation curve, the milling status is stable; otherwise, the milling is in chatter status. It can be concluded that calculation results agree well with experimental results.Figure 17
Comparison between calculation results and test results of milling stability.


## 5. Conclusions

(1)
A chatter stability prediction method of spindle-tool holder-tool system is constructed considering interface contact characteristics. The fractal parameters and cutting force coefficients are identified from the tests. The reliability of the method was proved by milling force tests and milling stability experiments.(2)
The increase of front bearings preload is beneficial while increase of back bearings stiffness is not conducive to improving the milling stability. The increase of spindle-tool holder’s interference fit is significantly beneficial to improve the milling stability.(3)
In order to improve the milling stability, feed rate should vary with spindle speed. There exists an optimal feed rate corresponding to a spindle speed. The optimal feed rate can be determined by the milling stability prediction proposed in this paper. Results find that the optimal feed rate should change with the spindle speed. For the spindle studied in this paper, the optimal feed rate should be 600 mm/min, 1200 mm/min, and 1800 mm/min when the spindle speed changes from 3000 r/min to 5000 r/min and 8000 r/min, respectively.(4)
The critical axial cutting depth decreases with increase of radial cutting depth. Results find that the critical axial cutting depth decreases sharply when radial cutting depth exceeds a critical value. In order to increase the critical axial cutting depth and, simultaneously, ensure cutting efficiency, the radial cutting depth should not exceed the critical value.(5)
The method proposed in this paper can be used as an instruction for predicting and suppressing chatter vibration and optimizing cutting parameters and also is helpful for designing the spindle-tool holder-tool system.



---

*Source: 7121328-2020-12-07.xml*