# Nonlinear Dynamics of Cutting Process considering Higher-Order Deformation of Composite Cutting Tool

Authors: Donghui Yao; Yongsheng Ren; Yuhuan Zhang; Bole Ma
Journal: Shock and Vibration
 (2021)

---

## Abstract

In this paper, the nonlinear dynamic analysis of the cutting process of composite cutting tool is performed. The cutting tool is simplified to a nonplanar bending rotating shaft. The higher-order bending deformation, structural damping, and gyroscopic effect of cutting tool are considered. It is assumed that cutting tool is subjected to a regenerative two-dimensional cutting force containing the first and second harmonic components. Based on the Hamilton principle, the motion equation of nonlinear chatter of the cutting system is derived. The nonlinear ordinary differential equations in the generalized coordinates are obtained by Galerkin method. In order to analyze the nonlinear dynamic response of cutting process, the multiscale method is used to derive the analytical approximate solution of the forced response for the cutting system under periodic cutting forces. In the forced response analysis, four cases including primary resonance and superharmonic resonance, i.e.,Ω¯=ω1, Ω¯=ω2, 2Ω¯=ω1, and 2Ω¯=ω2, are considered. The influences of ratio of length to diameter, structural damping, cutting force, and ply angle on primary resonance and superharmonic resonance are investigated. The results show that nonlinearity due to higher-order bending deformation significantly affects the dynamic behavior of the milling process and that the effective nonlinearity of the cutting system is of hard type. Multivalued resonance curves and jump phenomenon are presented. The influences of various factors, such as ratio of length to diameter, ply angle, structural damping, cutting force, and rotating speed, are thoroughly discussed.


---

## Body

## 1. Introduction

As a high-efficiency, high-quality, low-cost machining method, high-speed cutting technology has been widely used in aerospace and mold manufacturing. However, chatter reduces the cutting stability during machining operations, causes a decrease in the machining quality and cutting efficiency of the workpiece, damages the tool, and shortens the service life of the machine.The passive chatter control methods are mainly based on various types of dynamic vibration absorbers [1, 2] and impact dampers [3]. However, composite materials are known to have higher static stiffness, damping, and specific stiffness. It has been found that the dynamic stiffness and fundamental natural frequency of cutter bar may be improved simultaneously if composite is employed [4–6]. This is very beneficial for the stability of high rotational speed machining for deep holes.Within the framework of linear chatter theories and based on different theories, several of chatter phenomena as well as their stability boundary in different cutting processes were discovered [7–14]. Tobias [8] introduced the time-delay instability in the cutting force and proposed the regeneration theory, which is considered to provide a complete explanation of the chatter phenomenon so far.However, because linear theory cannot predict some important dynamic behaviors of cutting process, nonlinear modeling of cutting systems has received more attention. The nonlinearity in cutting systems is mainly caused by structural nonlinearity, delayed nonlinearity of cutting force, etc. [15–18]. Hanna and Tobias [15] first proposed a time-delay nonlinear model with quadratic and cubic structural stiffness and cutting forces, which has inspired great interest in analyzing the global dynamics of the cutting system. The effective mathematical methods for the nonlinear theory of cutting systems include center manifold theory, bifurcation theory, perturbation analysis, phase portraits, and Poincaré sections.Pratt [19] used the multiscale method, harmonic balance method, and Floquet theory to study the models of Hanna and Tobias [18]. The results showed that subcritical Hopf bifurcation may occur due to the existence of cubic structural nonlinearity. Moradi et al. [20] studied the existence of different types of bifurcation in the cutting process considering the tool wear and process damping. They developed a 2-dof linear model of the tool and a polynomial nonlinear model of the cutting force. In addition, the multiscale method was used to obtain the approximation analysis solution of primary resonance. Moradi et al. [21] used a model similar to that in literature [20] to analyze the forced vibration of the milling process. In the study, not only the primary resonance but also the superharmonic resonance and the internal resonance were investigated. Moradi et al. [22] investigated the internal resonance and regenerative chatter of the milling process considering both the cutting force and the structural nonlinearity.However, in the above studies, the tool structure was simulated with 1-dof or 2-dof vibration system, in which the concentration parameters such as mass and stiffness were obtained through experiments or empirical method. There is no obvious correlation between these simplified models and the continuous system equations of the tool.In order to investigate the dynamic characteristics of cutting system and the stability mechanism of machining process, it is necessary to conduct a comprehensive analysis on various influencing factors. In this case, if the empirical method is used to establish the model of the tool structure, a large number of repetitive tests are needed to obtain the dynamic parameters of the tool structure of different sizes, geometries, and materials, which is very time consuming and of low effectiveness. Therefore, as a more effective modeling method for cutting tool, the continuous parameter dynamics modeling of cutting tool based on analytical method, has received great attention [23–33].However, in most of the existing continuous system dynamics models, the influence of the nonlinearity of the tool structure has not been considered. Therefore, the existing dynamic model is only suitable for analyzing the linear dynamics of the natural frequency, chatter frequency, dynamic deformation, and chatter stability of the cutting system. In addition, the current continuous system dynamic studies are mainly focused on the tool structure made of isotropic metal materials. Nevertheless, there is no research about the nonlinear dynamics of the process with anisotropy composite tool structures.The structural nonlinearity is essentially arisen by the flexible nature of the cutting tools with small diameter or long extension part. The structural nonlinearity of the cutting tools may be described using higher-order large deformations [34] or geometric nonlinearity [35].In this paper, the nonlinear dynamics of the cutting process with a composite cutting tool are investigated. The composite tool structure is simplified to a nonplanar cantilever rotating shaft. The structural nonlinearity is introduced by considering the higher-order deformation of cutting tool. It is assumed that the cutting tool is subjected to a regenerative cutting force containing harmonic components. Based on the Hamilton principle, the nonlinear chatter equation of the cutting system is derived. The Galerkin method is used to obtain the nonlinear ordinary differential equations in the generalized coordinates. The multiscale method is used to obtain the approximate solution of the forced vibration response of the cutting process subjected to periodic cutting forces. Nonlinear dynamics of the system are studied for four cases of primary and superharmonic resonances; i.e.,Ω¯=ω1, Ω¯=ω2, 2Ω¯=ω1, and 2Ω¯=ω2 are studied. The numerical calculation is conducted to investigate the effect of various parameters on the frequency response of the cutting system.


## 2. Mathematical Model and Solution

### 2.1. Kinetic Energy and Potential Energy

The structural sketch of composite cutting tool with nonplanar bending is shown in Figure1. The kinetic energy of the rotating cutting tool without considering torsional deformation is as follows:(1)T=12∫0lρAv˙2+w˙2+ρIψ˙y2+ψ˙z2−4ΩρIψzψ˙ydx,where v and w represent the displacement at any point on the neutral axis along the y and z directions, respectively. ψz and ψy represent bending rotation angles of the cross section around the y and z axes, respectively. ρA represents the mass per unit length, and ρI represents the mass moment of inertia of the cross section. The superposed dots represent derivatives with respect to the time t.Figure 1
Composite cutting tool with nonplanar bending.Using the expression of the displacement fields for Euler–Bernoulli beam and the stress-strain displacement relations of the composite cutting tool, as shown in Appendix A, one can obtain the following expression for the potential energy of the rotating composite cutting tool with the higher-order deformation:(2)U=12∫0L∬AQ¯11z2∂ψy∂x2+y2∂ψz∂x2−2yz∂ψy∂x∂ψz∂x+14ψy4+14ψz4+12ψy2ψz2−2z∂ψy∂x+y∂ψz∂x12ψy2+12ψz22dAdx,where Q¯11 refers to the off-axis stiffness coefficient of the kth layer for the composite cutting tool.Regardless of the influence of shear deformation, based on the Euler–Bernoulli beam theory, the bending rotation angles and the displacements have the following relationship:(3)ψz=−∂w∂x,ψy=∂v∂x.


### 2.2. Damping Dissipative Energy Function and Virtual Work of Cutting Force

Rayleigh dissipative energy function of the cutting tool can be expressed as follows:(4)Wc=12∫0Lcv˙2+w˙2dx,where c is structural damping coefficient.The virtual work of the cutting forceδW can be expressed as(5)δW=∫0LLvδv+Lwδwdx,where Lv=FyδDx−L, Lw=FzδDx−L, and δDx−L is the delta function.The linear regenerative cutting forcesFy and Fz can be obtained by [20](6)Fy=−N2πα0v−vτ+β0w−wτ+γ0+cf2ς1cos2Ωt−π/2−η1cos2Ωt+η1+η2cosΩt−π/2+ς2cosΩt,Fz=N2πα0′v−vτ+β0′w−wτ+γ0′+cf2η1cos2Ωt−π/2+ς1cos2Ωt−ς1+ς2cosΩt−π/2−η2cosΩt,where(7)α0=0.5ς1+0.25πη1,β0=0.5η1+0.25πς1,α0′=−0.5η1+0.25πς1,β0′=0.5ς1−0.25πη1,γ0=η2+ς2,γ0′=ς2−η2,ς1=Ktca,ς2=Ktea,η1=Krca,η2=Krea,v=vx,t,vτ=vx,t−τ,w=wx,t,wτ=wx,t−τ,in which τ=2π/NΩ is the delay time of milling process; N is the number of teeth; (Ktc, Kte) and (Krc, Kre) are the cutting force coefficients in the tangential and radial directions, respectively; a is the axial cutting depth; and cf is the feed per tooth per revolution.


### 2.3. Dynamic Model of Milling Process

In order to obtain the equation of motion of the cutting system, the principle of Hamilton is used as follows:(8)∫t1t2δU+δWc−δT−δWdt=0.Substituting (1), (2), (4), and (5) into (8), the equations of motion in both y and z directions can be obtained:(9)mequivv¨+cv˙−2ρIequivΩw˙″−ρIequivv¨″+EIequiv∂4v∂x4+EAequiv2v′2+w′2v′′=FyδDx−L,mequivw¨+cw˙+2ρIequivΩv˙″−ρIequivw¨″+EIequiv∂4w∂x4+EAequiv2v′2+w′2w′′=FzδDx−L,where(10)mequiv=ρAequiv=π∑k=1Nρkrk+12−rk2,ρIequiv=π4∑k=1Nρkrk+14−rk4,EIequiv=π4∑k=1NQ¯11rk+14−rk4,EAequiv=π∑k=1NQ¯11rk+12−rk2.Here, (EI)equiv and (EA)equiv represent the equivalent bending and tensile stiffness of the cutting tool, respectively; (ρI)equiv and (m)equiv represent the equivalent diametrical mass moment of inertia and the equivalent mass per unit length, respectively; ρk represents the density of the kth layer; and rk and rk+1 represent the inner diameter and outer diameter of the kth layer, respectively. The detailed expression of Q¯11 can be found in Appendix B. In (9), the primes denote differentiation with respect to x.The following nondimensional quantities are defined:(11)x¯=xL,v¯=vL,w¯=wL,t¯=EIequivmequivL4t=γt,τ¯=γτ,c¯=cL2mequivEIequiv,ρ¯I¯=ρIequivmequivL2,Ω¯=mequivL4EIequivΩ,λ=EAequivL22EIequiv,α¯0=L3α0EIequiv,β¯0=L3β0EIequiv,γ¯0=L2γ0EIequiv,cfς1=L2cfEIequivς1,cfη1=L2cfEIequivη1,α¯0′=L3α0′EIequiv,β¯0′=L3β0′EIequiv,γ¯0′=L2γ0′EIequiv,ς¯2=L2ς2EIequiv,η¯2=L2η2EIequiv.Using the above dimensionless variables, (9) is rewritten as(12)v¨+c¯v˙−2ρ¯IΩ¯w˙″−ρ¯Iv¨″+∂4v∂x4+λv′2+w′2w′′=F¯yδDx−1,w¨+c¯w˙+2ρ¯IΩ¯v˙″−ρ¯Iw¨″+∂4w∂x4+λv′2+w′2w′′=F¯zδDx−1,where(13)F¯y=−N2πα¯0v−vτ+β¯0w−wτ+γ¯0+12cfς1cos2Ω¯t¯−π/2−cfη1cos2Ω¯t¯+cfη1+η¯2cos2Ω¯t¯−π/2+ς¯2cosΩ¯t¯,F¯z=N2πα¯0′v−vτ+β¯′w−wτ+γ0′−12cfη1cos2Ω¯t¯−π/2+cfς1cos2Ω¯t¯−cfς1+ς¯2cos2Ω¯t¯−π/2−η¯2cosΩ¯t¯.In (12), the primes denote differentiation with respect to x¯, and the superposed dots denote derivatives with respect to the time t¯.The solution of (12) can be written as(14)vx¯,t¯=ϕ1x¯Vt¯,wx¯,t¯=ϕ1x¯Wt¯.For cantilever beams,ϕ1x¯ has to meet the following boundary conditions:(15)ϕ10=ϕ1′0=ϕ1″1=ϕ1‴1=0.The mode function that satisfies the boundary condition in (12) can be expressed as follows:(16)ϕ1x¯=cosβ1Lx¯−coshβ1Lx¯−cosβ1L+coshβ1Lsinβ1L+sinhβ1Lsinβ1Lx¯−sinhβ1Lx¯,x¯∈0,1,where β1L=1.875.Substituting (14) into (9), simplifying the equation by the Galerkin method, and dropping the constant cutting force terms, the following ordinary differential equations can be obtained.(17)A¯1V¨+A¯2V˙−A¯3W˙+A¯4V+A¯5VV2+W2=−N2πα¯0V−Vτ+β¯0W−Wτϕ121+F¯y2ϕ11,A¯1W¨+A¯2W˙+A¯3V˙+A¯4W+A¯5WV2+W2=N2πα¯0′V−Vτ+β¯0′W−Wτϕ121+F¯z2ϕ11,where(18)A¯1=∫01ϕ12x¯dx¯−ρ¯I¯∫01ϕ1vϕ1″x¯dx¯,A¯2=c¯∫01ϕ12x¯dx¯,A¯3=2ρ¯I¯Ω¯∫01ϕ1x¯ϕ1″x¯dx¯,A¯4=∫01ϕ1x¯ϕ14x¯dx¯,A¯5=λ∫01ϕ1′x¯3′ϕ1dx¯.Here, the superscript (4) means the fourth-order partial derivative ofx¯.


### 2.4. Perturbation Analysis of Milling Process Using Multiple Scales Method

In order to solve (17) using the multiscale method, the structural damping and nonlinear and exciting force terms are scaled using small parameters ε as follows:(19)A¯2⇒ε2A¯2,N2πα¯0⇒ε2N2πα¯0,N2πβ¯0⇒ε2N2πβ¯0,N2πα¯0′⇒ε2N2πα¯0′,N2πβ¯0′⇒ε2N2πβ¯0′,F¯¯y2⇒ε3F¯¯y2,F¯¯z2⇒ε3F¯¯z2.Vt¯ and Wt¯ are expanded in the form(20)Vt¯=εV1T0,T2+ε3V3T0,T2+⋯,Wt¯=εW1T0,T2+ε3W3T0,T2+⋯,where T0=t¯,T2=ε2t¯.Taking the derivative of (20), one obtains the following equations:(21)V˙=ε∂V1∂T0+ε3∂V3∂T0+∂V1∂T2+⋯,V¨=ε∂2V1∂T02+ε3∂2V3∂T02+2∂2V1∂T0∂T2+⋯,W¨=ε∂W1∂T0+ε3∂W3∂T0+∂W1∂T2+⋯,W¨=ε∂2W1∂T02+ε3∂2W3∂T02+2∂2W1∂T0∂T2+⋯.Substituting (20) and (21) into (17) and considering (19), one can obtain the following equations by comparing the coefficients of ε and ε3:O(ε):(22)A¯1∂2V1∂T02−A¯3∂W1∂T0+A¯4V1=0,A¯1∂2W1∂T02+A¯3∂V1∂T0+A¯4W1=0.O(ε3):(23)A¯1∂2V3∂T02−A¯3∂W3∂T0+A¯4V3=−2A¯1∂2V1∂T0∂T2−A¯2∂V1∂T0+A¯3∂W1∂T2−A¯5V1V12+W12−N2πα¯0V1−V1τ+β¯0W1−W1τϕ121+λV+F¯y2ϕ11,A¯1∂2W3∂T02+A¯3∂V3∂T0+A¯4W3=−2A¯1∂2W1∂T0∂T2−A¯2∂W1∂T0−A¯3∂V1∂T2−A¯5W1V12+W12+N2πα¯0′V1−V1τ+β¯0′W1−W1τϕ121+λW+F¯z2ϕ11.Assume the solution of (22) is as follows:(24)V1T0,T2=F1T2eiω1T0+F2T2eiω2T0+F¯1T2e−iω1T0+F¯2T2e−iω2T0,W1T0,T2=−iF1T2eiω1T0+iF2T2eiω2T0+iF¯1T2e−iω1T0−iF¯2T2e−iω2T0,V1τT0,T2=F1T2eiω1T0−τ¯+F2T2eiω2T0−τ¯+F¯1T2e−iω1T0−τ¯+F¯2T2e−iω2T0−τ¯,W1τT0,T2=−iF1T2eiω1T0−τ¯+iF2T2eiω2T0−τ¯+iF¯1T2e−iω1T0−τ¯−iF¯2T2e−iω2T0−τ¯,where i=−1 represents the imaginary unit, F1T2 and F2T2 are the complex functions to be determined, and F¯1T2 and F¯2T2 represent the complex conjugate. ω1 and ω2 refer to the forward and backward whirling frequencies, respectively, which are expressed as follows:(25)ω1=A¯3+A¯32+4A¯1A¯42A¯1,ω2=−A¯3+A¯32+4A¯1A¯42A¯1.Substituting (24) into (23) and dropping the constant terms λV and λW, one obtains the following equations:(26)A¯1∂2V3∂T02−A¯3∂W3∂T0+A¯4V3=P3eiω1T0+Q3eiω2T0+F¯¯y2,A¯1∂2W3∂T02+A¯3∂V3∂T0+A¯4W3=P4eiω1T0+Q4eiω2T0+F¯¯z2,where(27)P3=Γ1−iΛ1F1′−iA¯2ω1F1−ϕ121N2πα¯¯0−iβ¯¯01−e−iω1τ¯F1,Q3=Γ2−iΛ2F2′−iA¯2ω2F2−ϕ121N2πα¯¯0+iβ¯¯01−e−iω2τ¯F2,P4=−iΓ1−Λ1F1′−A¯2ω1F1+ϕ121N2πα¯¯0′−iβ¯¯0′1−e−iω1τ¯F1,Q4=iΓ2+Λ2F2′+A¯2ω2F2+ϕ121N2πα¯¯0′+iβ¯¯0′1−e−iω2τ¯F2,Λ1=2A¯1ω1+A¯3,Γ1=−4A¯5F12F¯1+2F1F2F¯2,Λ2=2A¯1ω2−A¯3,Γ2=−4A¯5F22F¯2+2F1F2F¯1,F¯¯y2=F¯y2ϕ11,F¯¯z2=F¯z2ϕ11.In this paper, four cases of primary and superharmonic resonances, i.e., case I:Ω¯=ω1, case II: Ω¯=ω2, case III: 2Ω¯=ω1, and case IV: 2Ω¯=ω2, are studied.


#### 2.4.1. Case I: Primary Resonance (Ω¯=ω1)

Assume(28)Ω¯=ω1+ε2σ,where σ represents the detuning parameter.Substituting (24) into (26), one can obtain the following equations:(29)A¯1∂2V3∂T02−A¯3∂W3∂T0+A¯4V3=P¯3eiω1T0+Q¯3eiω2T0,A¯1∂2W3∂T02+A¯3∂V3∂T0+A¯4W3=P¯4eiω1T0+Q¯4eiω2T0,where(30)P¯3=P3+ϕ11q¯1eiσT2,P¯4=P4−iϕ11q¯1eiσT2,Q¯3=Q3,Q¯4=Q4,where q¯1=ς¯2−iη¯2/2.The particular solutions of (29) are(31)V3T0,T2=F11T2eiω1T0+F12T2eiω2T0,W3T0,T2=F21T2eiω1T0+F22T2eiω2T0.Substituting (31) into (29) and equating the coefficient of eiω1T0 and eiω2T0 in both sides of (29), one has(32)A¯4−A¯1ω12F11−iA¯3ω1F21=P¯3,iA¯3ω1F11+A¯4−A¯1ω12F21=P¯4.(33)A¯4−A¯1ω22F12−iA¯3ω2F22=Q¯3,iA¯3ω2F12+A¯4−A¯1ω22F22=Q¯4.Equations (32) and (33) constitute systems of two inhomogeneous algebraic equations for F11, F21 and F12, F22, respectively. Their homogeneous parts have a nontrivial solution. Then their solvability conditions can be written as [36](34)P¯3−iA¯3ω1P¯4A¯4−A¯1ω12=0,A¯4−A¯1ω12P¯3iA¯3ω1P¯4=0.(35)Q¯3−iA¯3ω2Q¯4A¯4−A¯1ω22=0,A¯4−A¯1ω22Q¯3iA¯3ω2Q¯4=0.After simplification, the solvability conditions are reduced to two independent governing equations in the following form:(36)4A¯5F12F¯1+2F1F2F¯2−iΛ1F1′−iA¯2ω1F1+ϕ121N4πβ¯0′−α¯0+iα¯0′+β¯01−e−iω1τ¯F1+ϕ11q¯1eiσT2=0,(37)4A¯5F22F¯2+2F1F2F¯1−iΛ2F2′−iA¯2ω2F2−12ϕ121N2πα¯0+β¯0′+iα¯0′−β¯01−e−iω1τ¯F2=0.Assume that the solutions of (36) and (37) are(38)F1=12a1T2eiθ1T2,F2=12a2T2eiθ2T2,where ajT2 and θjT2 (j = 1,2) refer to the amplitudes and phase angles of the response, respectively.Substituting (38) into (36) and (37), separating the real and the imaginary parts, one can obtain the following equations:(39)Λ1a1′=−A¯2ω1a1+ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯a1+ϕ11ς¯2sinψ1−η¯2cosψ1,Λ1a1ψ1′=Λ1σa1−A¯5a13+2a1a22+ϕ121N4πβ¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯a1+ϕ11η¯2sinψ1+ς¯2cosψ1.(40)Λ2a2′=−A¯2ω2a2−ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯a2,Λ2ψ2′a2=−A¯5a23+2a12a2−ϕ121N4πα¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯a2,where σT2−θj=ψj (j = 1,2).To determine the steady-state forced response, the time derivatives in (39) and (40) are set to zero. It can be immediately concluded from (40) that only solution for a2 is zero. This shows that under the primary resonance Ω≈ω1, only the first mode, i.e., the forward mode, can be excited, while the second mode, i.e., the backward mode, does not participate in the primary resonance and remains stationary.Substitutinga2 = 0 into (39) and solving for σ, one obtains(41)σ=−ϕ121N4πΛ1β¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯+A¯5Λ1a12±1Λ1ϕ121ς¯22+η¯22a12−A¯2ω1−ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯2.To study the stability of the steady-state response of case I, the nature of steady-state response can be investigated numerically by linearizing (39) (with a2 = 0) around (a0, ψ0).


#### 2.4.2. Case II: Primary Resonance (Ω¯=ω2)

According to (26), under such conditions, one introduces detuning parameters as(42)Ω¯=ω2+ε2σ.Using their solvability conditions (which are not shown) and (38) leads to(43)Λ1a1′=−A¯2ω1a1+ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯a1,Λ1a1ψ1′=Λ1σa1−A¯5a13+2a1a22+ϕ121N4πβ¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯a1,Λ2a2′=−A¯2ω2a2−ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯a2−ϕ11ς¯2sinψ2−η¯2cosψ2,Λ2ψ2′a2=−A¯5a23+2a12a2−ϕ121N4πα¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯a2−ϕ11ς¯2cosψ2+η¯2sinψ2.By using steady-state condition (i.e., seta1′ = 0, ψ1′ = 0, a2′ = 0, and ψ2′ = 0), the following steady-state response of case II can be obtained:(44)σ=ϕ121N4πΛ2α¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯+A¯5Λ2a22±1Λ2ϕ121ς¯22+η¯22a22−A¯2ω2+ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯2.From (44), it can be seen that a1=0 and a2≠0, indicating that only the forced vibration response of backward mode is excited, while forward mode remains stationary.


#### 2.4.3. Case III: Superharmonic Resonance (2Ω¯=ω1)

For this case, the formulation of steady-state response is the same as case I, but in (41), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


#### 2.4.4. Case IV: Superharmonic Resonance (2Ω¯=ω2)

For this case, the formulation of steady-state response is the same as case II, but in (43), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


## 2.1. Kinetic Energy and Potential Energy

The structural sketch of composite cutting tool with nonplanar bending is shown in Figure1. The kinetic energy of the rotating cutting tool without considering torsional deformation is as follows:(1)T=12∫0lρAv˙2+w˙2+ρIψ˙y2+ψ˙z2−4ΩρIψzψ˙ydx,where v and w represent the displacement at any point on the neutral axis along the y and z directions, respectively. ψz and ψy represent bending rotation angles of the cross section around the y and z axes, respectively. ρA represents the mass per unit length, and ρI represents the mass moment of inertia of the cross section. The superposed dots represent derivatives with respect to the time t.Figure 1
Composite cutting tool with nonplanar bending.Using the expression of the displacement fields for Euler–Bernoulli beam and the stress-strain displacement relations of the composite cutting tool, as shown in Appendix A, one can obtain the following expression for the potential energy of the rotating composite cutting tool with the higher-order deformation:(2)U=12∫0L∬AQ¯11z2∂ψy∂x2+y2∂ψz∂x2−2yz∂ψy∂x∂ψz∂x+14ψy4+14ψz4+12ψy2ψz2−2z∂ψy∂x+y∂ψz∂x12ψy2+12ψz22dAdx,where Q¯11 refers to the off-axis stiffness coefficient of the kth layer for the composite cutting tool.Regardless of the influence of shear deformation, based on the Euler–Bernoulli beam theory, the bending rotation angles and the displacements have the following relationship:(3)ψz=−∂w∂x,ψy=∂v∂x.


## 2.2. Damping Dissipative Energy Function and Virtual Work of Cutting Force

Rayleigh dissipative energy function of the cutting tool can be expressed as follows:(4)Wc=12∫0Lcv˙2+w˙2dx,where c is structural damping coefficient.The virtual work of the cutting forceδW can be expressed as(5)δW=∫0LLvδv+Lwδwdx,where Lv=FyδDx−L, Lw=FzδDx−L, and δDx−L is the delta function.The linear regenerative cutting forcesFy and Fz can be obtained by [20](6)Fy=−N2πα0v−vτ+β0w−wτ+γ0+cf2ς1cos2Ωt−π/2−η1cos2Ωt+η1+η2cosΩt−π/2+ς2cosΩt,Fz=N2πα0′v−vτ+β0′w−wτ+γ0′+cf2η1cos2Ωt−π/2+ς1cos2Ωt−ς1+ς2cosΩt−π/2−η2cosΩt,where(7)α0=0.5ς1+0.25πη1,β0=0.5η1+0.25πς1,α0′=−0.5η1+0.25πς1,β0′=0.5ς1−0.25πη1,γ0=η2+ς2,γ0′=ς2−η2,ς1=Ktca,ς2=Ktea,η1=Krca,η2=Krea,v=vx,t,vτ=vx,t−τ,w=wx,t,wτ=wx,t−τ,in which τ=2π/NΩ is the delay time of milling process; N is the number of teeth; (Ktc, Kte) and (Krc, Kre) are the cutting force coefficients in the tangential and radial directions, respectively; a is the axial cutting depth; and cf is the feed per tooth per revolution.


## 2.3. Dynamic Model of Milling Process

In order to obtain the equation of motion of the cutting system, the principle of Hamilton is used as follows:(8)∫t1t2δU+δWc−δT−δWdt=0.Substituting (1), (2), (4), and (5) into (8), the equations of motion in both y and z directions can be obtained:(9)mequivv¨+cv˙−2ρIequivΩw˙″−ρIequivv¨″+EIequiv∂4v∂x4+EAequiv2v′2+w′2v′′=FyδDx−L,mequivw¨+cw˙+2ρIequivΩv˙″−ρIequivw¨″+EIequiv∂4w∂x4+EAequiv2v′2+w′2w′′=FzδDx−L,where(10)mequiv=ρAequiv=π∑k=1Nρkrk+12−rk2,ρIequiv=π4∑k=1Nρkrk+14−rk4,EIequiv=π4∑k=1NQ¯11rk+14−rk4,EAequiv=π∑k=1NQ¯11rk+12−rk2.Here, (EI)equiv and (EA)equiv represent the equivalent bending and tensile stiffness of the cutting tool, respectively; (ρI)equiv and (m)equiv represent the equivalent diametrical mass moment of inertia and the equivalent mass per unit length, respectively; ρk represents the density of the kth layer; and rk and rk+1 represent the inner diameter and outer diameter of the kth layer, respectively. The detailed expression of Q¯11 can be found in Appendix B. In (9), the primes denote differentiation with respect to x.The following nondimensional quantities are defined:(11)x¯=xL,v¯=vL,w¯=wL,t¯=EIequivmequivL4t=γt,τ¯=γτ,c¯=cL2mequivEIequiv,ρ¯I¯=ρIequivmequivL2,Ω¯=mequivL4EIequivΩ,λ=EAequivL22EIequiv,α¯0=L3α0EIequiv,β¯0=L3β0EIequiv,γ¯0=L2γ0EIequiv,cfς1=L2cfEIequivς1,cfη1=L2cfEIequivη1,α¯0′=L3α0′EIequiv,β¯0′=L3β0′EIequiv,γ¯0′=L2γ0′EIequiv,ς¯2=L2ς2EIequiv,η¯2=L2η2EIequiv.Using the above dimensionless variables, (9) is rewritten as(12)v¨+c¯v˙−2ρ¯IΩ¯w˙″−ρ¯Iv¨″+∂4v∂x4+λv′2+w′2w′′=F¯yδDx−1,w¨+c¯w˙+2ρ¯IΩ¯v˙″−ρ¯Iw¨″+∂4w∂x4+λv′2+w′2w′′=F¯zδDx−1,where(13)F¯y=−N2πα¯0v−vτ+β¯0w−wτ+γ¯0+12cfς1cos2Ω¯t¯−π/2−cfη1cos2Ω¯t¯+cfη1+η¯2cos2Ω¯t¯−π/2+ς¯2cosΩ¯t¯,F¯z=N2πα¯0′v−vτ+β¯′w−wτ+γ0′−12cfη1cos2Ω¯t¯−π/2+cfς1cos2Ω¯t¯−cfς1+ς¯2cos2Ω¯t¯−π/2−η¯2cosΩ¯t¯.In (12), the primes denote differentiation with respect to x¯, and the superposed dots denote derivatives with respect to the time t¯.The solution of (12) can be written as(14)vx¯,t¯=ϕ1x¯Vt¯,wx¯,t¯=ϕ1x¯Wt¯.For cantilever beams,ϕ1x¯ has to meet the following boundary conditions:(15)ϕ10=ϕ1′0=ϕ1″1=ϕ1‴1=0.The mode function that satisfies the boundary condition in (12) can be expressed as follows:(16)ϕ1x¯=cosβ1Lx¯−coshβ1Lx¯−cosβ1L+coshβ1Lsinβ1L+sinhβ1Lsinβ1Lx¯−sinhβ1Lx¯,x¯∈0,1,where β1L=1.875.Substituting (14) into (9), simplifying the equation by the Galerkin method, and dropping the constant cutting force terms, the following ordinary differential equations can be obtained.(17)A¯1V¨+A¯2V˙−A¯3W˙+A¯4V+A¯5VV2+W2=−N2πα¯0V−Vτ+β¯0W−Wτϕ121+F¯y2ϕ11,A¯1W¨+A¯2W˙+A¯3V˙+A¯4W+A¯5WV2+W2=N2πα¯0′V−Vτ+β¯0′W−Wτϕ121+F¯z2ϕ11,where(18)A¯1=∫01ϕ12x¯dx¯−ρ¯I¯∫01ϕ1vϕ1″x¯dx¯,A¯2=c¯∫01ϕ12x¯dx¯,A¯3=2ρ¯I¯Ω¯∫01ϕ1x¯ϕ1″x¯dx¯,A¯4=∫01ϕ1x¯ϕ14x¯dx¯,A¯5=λ∫01ϕ1′x¯3′ϕ1dx¯.Here, the superscript (4) means the fourth-order partial derivative ofx¯.


## 2.4. Perturbation Analysis of Milling Process Using Multiple Scales Method

In order to solve (17) using the multiscale method, the structural damping and nonlinear and exciting force terms are scaled using small parameters ε as follows:(19)A¯2⇒ε2A¯2,N2πα¯0⇒ε2N2πα¯0,N2πβ¯0⇒ε2N2πβ¯0,N2πα¯0′⇒ε2N2πα¯0′,N2πβ¯0′⇒ε2N2πβ¯0′,F¯¯y2⇒ε3F¯¯y2,F¯¯z2⇒ε3F¯¯z2.Vt¯ and Wt¯ are expanded in the form(20)Vt¯=εV1T0,T2+ε3V3T0,T2+⋯,Wt¯=εW1T0,T2+ε3W3T0,T2+⋯,where T0=t¯,T2=ε2t¯.Taking the derivative of (20), one obtains the following equations:(21)V˙=ε∂V1∂T0+ε3∂V3∂T0+∂V1∂T2+⋯,V¨=ε∂2V1∂T02+ε3∂2V3∂T02+2∂2V1∂T0∂T2+⋯,W¨=ε∂W1∂T0+ε3∂W3∂T0+∂W1∂T2+⋯,W¨=ε∂2W1∂T02+ε3∂2W3∂T02+2∂2W1∂T0∂T2+⋯.Substituting (20) and (21) into (17) and considering (19), one can obtain the following equations by comparing the coefficients of ε and ε3:O(ε):(22)A¯1∂2V1∂T02−A¯3∂W1∂T0+A¯4V1=0,A¯1∂2W1∂T02+A¯3∂V1∂T0+A¯4W1=0.O(ε3):(23)A¯1∂2V3∂T02−A¯3∂W3∂T0+A¯4V3=−2A¯1∂2V1∂T0∂T2−A¯2∂V1∂T0+A¯3∂W1∂T2−A¯5V1V12+W12−N2πα¯0V1−V1τ+β¯0W1−W1τϕ121+λV+F¯y2ϕ11,A¯1∂2W3∂T02+A¯3∂V3∂T0+A¯4W3=−2A¯1∂2W1∂T0∂T2−A¯2∂W1∂T0−A¯3∂V1∂T2−A¯5W1V12+W12+N2πα¯0′V1−V1τ+β¯0′W1−W1τϕ121+λW+F¯z2ϕ11.Assume the solution of (22) is as follows:(24)V1T0,T2=F1T2eiω1T0+F2T2eiω2T0+F¯1T2e−iω1T0+F¯2T2e−iω2T0,W1T0,T2=−iF1T2eiω1T0+iF2T2eiω2T0+iF¯1T2e−iω1T0−iF¯2T2e−iω2T0,V1τT0,T2=F1T2eiω1T0−τ¯+F2T2eiω2T0−τ¯+F¯1T2e−iω1T0−τ¯+F¯2T2e−iω2T0−τ¯,W1τT0,T2=−iF1T2eiω1T0−τ¯+iF2T2eiω2T0−τ¯+iF¯1T2e−iω1T0−τ¯−iF¯2T2e−iω2T0−τ¯,where i=−1 represents the imaginary unit, F1T2 and F2T2 are the complex functions to be determined, and F¯1T2 and F¯2T2 represent the complex conjugate. ω1 and ω2 refer to the forward and backward whirling frequencies, respectively, which are expressed as follows:(25)ω1=A¯3+A¯32+4A¯1A¯42A¯1,ω2=−A¯3+A¯32+4A¯1A¯42A¯1.Substituting (24) into (23) and dropping the constant terms λV and λW, one obtains the following equations:(26)A¯1∂2V3∂T02−A¯3∂W3∂T0+A¯4V3=P3eiω1T0+Q3eiω2T0+F¯¯y2,A¯1∂2W3∂T02+A¯3∂V3∂T0+A¯4W3=P4eiω1T0+Q4eiω2T0+F¯¯z2,where(27)P3=Γ1−iΛ1F1′−iA¯2ω1F1−ϕ121N2πα¯¯0−iβ¯¯01−e−iω1τ¯F1,Q3=Γ2−iΛ2F2′−iA¯2ω2F2−ϕ121N2πα¯¯0+iβ¯¯01−e−iω2τ¯F2,P4=−iΓ1−Λ1F1′−A¯2ω1F1+ϕ121N2πα¯¯0′−iβ¯¯0′1−e−iω1τ¯F1,Q4=iΓ2+Λ2F2′+A¯2ω2F2+ϕ121N2πα¯¯0′+iβ¯¯0′1−e−iω2τ¯F2,Λ1=2A¯1ω1+A¯3,Γ1=−4A¯5F12F¯1+2F1F2F¯2,Λ2=2A¯1ω2−A¯3,Γ2=−4A¯5F22F¯2+2F1F2F¯1,F¯¯y2=F¯y2ϕ11,F¯¯z2=F¯z2ϕ11.In this paper, four cases of primary and superharmonic resonances, i.e., case I:Ω¯=ω1, case II: Ω¯=ω2, case III: 2Ω¯=ω1, and case IV: 2Ω¯=ω2, are studied.


### 2.4.1. Case I: Primary Resonance (Ω¯=ω1)

Assume(28)Ω¯=ω1+ε2σ,where σ represents the detuning parameter.Substituting (24) into (26), one can obtain the following equations:(29)A¯1∂2V3∂T02−A¯3∂W3∂T0+A¯4V3=P¯3eiω1T0+Q¯3eiω2T0,A¯1∂2W3∂T02+A¯3∂V3∂T0+A¯4W3=P¯4eiω1T0+Q¯4eiω2T0,where(30)P¯3=P3+ϕ11q¯1eiσT2,P¯4=P4−iϕ11q¯1eiσT2,Q¯3=Q3,Q¯4=Q4,where q¯1=ς¯2−iη¯2/2.The particular solutions of (29) are(31)V3T0,T2=F11T2eiω1T0+F12T2eiω2T0,W3T0,T2=F21T2eiω1T0+F22T2eiω2T0.Substituting (31) into (29) and equating the coefficient of eiω1T0 and eiω2T0 in both sides of (29), one has(32)A¯4−A¯1ω12F11−iA¯3ω1F21=P¯3,iA¯3ω1F11+A¯4−A¯1ω12F21=P¯4.(33)A¯4−A¯1ω22F12−iA¯3ω2F22=Q¯3,iA¯3ω2F12+A¯4−A¯1ω22F22=Q¯4.Equations (32) and (33) constitute systems of two inhomogeneous algebraic equations for F11, F21 and F12, F22, respectively. Their homogeneous parts have a nontrivial solution. Then their solvability conditions can be written as [36](34)P¯3−iA¯3ω1P¯4A¯4−A¯1ω12=0,A¯4−A¯1ω12P¯3iA¯3ω1P¯4=0.(35)Q¯3−iA¯3ω2Q¯4A¯4−A¯1ω22=0,A¯4−A¯1ω22Q¯3iA¯3ω2Q¯4=0.After simplification, the solvability conditions are reduced to two independent governing equations in the following form:(36)4A¯5F12F¯1+2F1F2F¯2−iΛ1F1′−iA¯2ω1F1+ϕ121N4πβ¯0′−α¯0+iα¯0′+β¯01−e−iω1τ¯F1+ϕ11q¯1eiσT2=0,(37)4A¯5F22F¯2+2F1F2F¯1−iΛ2F2′−iA¯2ω2F2−12ϕ121N2πα¯0+β¯0′+iα¯0′−β¯01−e−iω1τ¯F2=0.Assume that the solutions of (36) and (37) are(38)F1=12a1T2eiθ1T2,F2=12a2T2eiθ2T2,where ajT2 and θjT2 (j = 1,2) refer to the amplitudes and phase angles of the response, respectively.Substituting (38) into (36) and (37), separating the real and the imaginary parts, one can obtain the following equations:(39)Λ1a1′=−A¯2ω1a1+ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯a1+ϕ11ς¯2sinψ1−η¯2cosψ1,Λ1a1ψ1′=Λ1σa1−A¯5a13+2a1a22+ϕ121N4πβ¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯a1+ϕ11η¯2sinψ1+ς¯2cosψ1.(40)Λ2a2′=−A¯2ω2a2−ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯a2,Λ2ψ2′a2=−A¯5a23+2a12a2−ϕ121N4πα¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯a2,where σT2−θj=ψj (j = 1,2).To determine the steady-state forced response, the time derivatives in (39) and (40) are set to zero. It can be immediately concluded from (40) that only solution for a2 is zero. This shows that under the primary resonance Ω≈ω1, only the first mode, i.e., the forward mode, can be excited, while the second mode, i.e., the backward mode, does not participate in the primary resonance and remains stationary.Substitutinga2 = 0 into (39) and solving for σ, one obtains(41)σ=−ϕ121N4πΛ1β¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯+A¯5Λ1a12±1Λ1ϕ121ς¯22+η¯22a12−A¯2ω1−ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯2.To study the stability of the steady-state response of case I, the nature of steady-state response can be investigated numerically by linearizing (39) (with a2 = 0) around (a0, ψ0).


### 2.4.2. Case II: Primary Resonance (Ω¯=ω2)

According to (26), under such conditions, one introduces detuning parameters as(42)Ω¯=ω2+ε2σ.Using their solvability conditions (which are not shown) and (38) leads to(43)Λ1a1′=−A¯2ω1a1+ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯a1,Λ1a1ψ1′=Λ1σa1−A¯5a13+2a1a22+ϕ121N4πβ¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯a1,Λ2a2′=−A¯2ω2a2−ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯a2−ϕ11ς¯2sinψ2−η¯2cosψ2,Λ2ψ2′a2=−A¯5a23+2a12a2−ϕ121N4πα¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯a2−ϕ11ς¯2cosψ2+η¯2sinψ2.By using steady-state condition (i.e., seta1′ = 0, ψ1′ = 0, a2′ = 0, and ψ2′ = 0), the following steady-state response of case II can be obtained:(44)σ=ϕ121N4πΛ2α¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯+A¯5Λ2a22±1Λ2ϕ121ς¯22+η¯22a22−A¯2ω2+ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯2.From (44), it can be seen that a1=0 and a2≠0, indicating that only the forced vibration response of backward mode is excited, while forward mode remains stationary.


### 2.4.3. Case III: Superharmonic Resonance (2Ω¯=ω1)

For this case, the formulation of steady-state response is the same as case I, but in (41), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


### 2.4.4. Case IV: Superharmonic Resonance (2Ω¯=ω2)

For this case, the formulation of steady-state response is the same as case II, but in (43), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


## 2.4.1. Case I: Primary Resonance (Ω¯=ω1)

Assume(28)Ω¯=ω1+ε2σ,where σ represents the detuning parameter.Substituting (24) into (26), one can obtain the following equations:(29)A¯1∂2V3∂T02−A¯3∂W3∂T0+A¯4V3=P¯3eiω1T0+Q¯3eiω2T0,A¯1∂2W3∂T02+A¯3∂V3∂T0+A¯4W3=P¯4eiω1T0+Q¯4eiω2T0,where(30)P¯3=P3+ϕ11q¯1eiσT2,P¯4=P4−iϕ11q¯1eiσT2,Q¯3=Q3,Q¯4=Q4,where q¯1=ς¯2−iη¯2/2.The particular solutions of (29) are(31)V3T0,T2=F11T2eiω1T0+F12T2eiω2T0,W3T0,T2=F21T2eiω1T0+F22T2eiω2T0.Substituting (31) into (29) and equating the coefficient of eiω1T0 and eiω2T0 in both sides of (29), one has(32)A¯4−A¯1ω12F11−iA¯3ω1F21=P¯3,iA¯3ω1F11+A¯4−A¯1ω12F21=P¯4.(33)A¯4−A¯1ω22F12−iA¯3ω2F22=Q¯3,iA¯3ω2F12+A¯4−A¯1ω22F22=Q¯4.Equations (32) and (33) constitute systems of two inhomogeneous algebraic equations for F11, F21 and F12, F22, respectively. Their homogeneous parts have a nontrivial solution. Then their solvability conditions can be written as [36](34)P¯3−iA¯3ω1P¯4A¯4−A¯1ω12=0,A¯4−A¯1ω12P¯3iA¯3ω1P¯4=0.(35)Q¯3−iA¯3ω2Q¯4A¯4−A¯1ω22=0,A¯4−A¯1ω22Q¯3iA¯3ω2Q¯4=0.After simplification, the solvability conditions are reduced to two independent governing equations in the following form:(36)4A¯5F12F¯1+2F1F2F¯2−iΛ1F1′−iA¯2ω1F1+ϕ121N4πβ¯0′−α¯0+iα¯0′+β¯01−e−iω1τ¯F1+ϕ11q¯1eiσT2=0,(37)4A¯5F22F¯2+2F1F2F¯1−iΛ2F2′−iA¯2ω2F2−12ϕ121N2πα¯0+β¯0′+iα¯0′−β¯01−e−iω1τ¯F2=0.Assume that the solutions of (36) and (37) are(38)F1=12a1T2eiθ1T2,F2=12a2T2eiθ2T2,where ajT2 and θjT2 (j = 1,2) refer to the amplitudes and phase angles of the response, respectively.Substituting (38) into (36) and (37), separating the real and the imaginary parts, one can obtain the following equations:(39)Λ1a1′=−A¯2ω1a1+ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯a1+ϕ11ς¯2sinψ1−η¯2cosψ1,Λ1a1ψ1′=Λ1σa1−A¯5a13+2a1a22+ϕ121N4πβ¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯a1+ϕ11η¯2sinψ1+ς¯2cosψ1.(40)Λ2a2′=−A¯2ω2a2−ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯a2,Λ2ψ2′a2=−A¯5a23+2a12a2−ϕ121N4πα¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯a2,where σT2−θj=ψj (j = 1,2).To determine the steady-state forced response, the time derivatives in (39) and (40) are set to zero. It can be immediately concluded from (40) that only solution for a2 is zero. This shows that under the primary resonance Ω≈ω1, only the first mode, i.e., the forward mode, can be excited, while the second mode, i.e., the backward mode, does not participate in the primary resonance and remains stationary.Substitutinga2 = 0 into (39) and solving for σ, one obtains(41)σ=−ϕ121N4πΛ1β¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯+A¯5Λ1a12±1Λ1ϕ121ς¯22+η¯22a12−A¯2ω1−ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯2.To study the stability of the steady-state response of case I, the nature of steady-state response can be investigated numerically by linearizing (39) (with a2 = 0) around (a0, ψ0).


## 2.4.2. Case II: Primary Resonance (Ω¯=ω2)

According to (26), under such conditions, one introduces detuning parameters as(42)Ω¯=ω2+ε2σ.Using their solvability conditions (which are not shown) and (38) leads to(43)Λ1a1′=−A¯2ω1a1+ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯a1,Λ1a1ψ1′=Λ1σa1−A¯5a13+2a1a22+ϕ121N4πβ¯0′−α¯01−cosω1τ¯−α¯0′+β¯0sinω1τ¯a1,Λ2a2′=−A¯2ω2a2−ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯a2−ϕ11ς¯2sinψ2−η¯2cosψ2,Λ2ψ2′a2=−A¯5a23+2a12a2−ϕ121N4πα¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯a2−ϕ11ς¯2cosψ2+η¯2sinψ2.By using steady-state condition (i.e., seta1′ = 0, ψ1′ = 0, a2′ = 0, and ψ2′ = 0), the following steady-state response of case II can be obtained:(44)σ=ϕ121N4πΛ2α¯0−β¯0′1−cosω2τ¯−α¯0′+β¯0sinω2τ¯+A¯5Λ2a22±1Λ2ϕ121ς¯22+η¯22a22−A¯2ω2+ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯2.From (44), it can be seen that a1=0 and a2≠0, indicating that only the forced vibration response of backward mode is excited, while forward mode remains stationary.


## 2.4.3. Case III: Superharmonic Resonance (2Ω¯=ω1)

For this case, the formulation of steady-state response is the same as case I, but in (41), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


## 2.4.4. Case IV: Superharmonic Resonance (2Ω¯=ω2)

For this case, the formulation of steady-state response is the same as case II, but in (43), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


## 3. Numerical Results and Discussions

In this study, the composite material of carbon fiber/epoxy resin is selected as the material of the cutting tool. The mechanical properties of the material are shown in Table1. Coefficients of cutting forces for simulation are given in Table 2. The cutting tool has a hollow structure, the outer diameter of the cross section is D = 8 mm, the inner diameter is D = 4 mm, the thickness of the section is h = 2 mm, and the length L is determined by the given ratio of length to diameter. The composite cutting tool has 16 layers with identical thickness, and the stack sequence is ±θ8. In all cases, Ω¯=200, except where other values are mentioned.Table 1
Mechanical properties of carbon fiber/epoxy composite [37].
ρ(kg/m3)E11 (GPa)E22 (GPa)G12 (GPa)G23 (GPa)ν12167225.88.73.53.50.34Table 2
Coefﬁcients of cutting forces [21].
Cf (mm/rev-tooth)ς¯¯1 (N/mm)η¯¯1 (N)ς¯¯2 (N/mm)η¯¯2 (N)ac (mm)N0.2620208435234Figure2 shows the natural frequency versus rotating speed, which is generally known as a Campbell diagram. In vibration of gyroscopic systems, there are two natural frequencies associated with forward and backward whirling motions. In forward natural frequency, the natural frequency is measured when the rotating cutting tool whirls in direction of the rotation. However, in backward natural frequency, the natural frequency is measured when the cutting tool whirls in the opposite direction of the rotation. The forward natural frequency (black solid line) increases with the increase of the rotating speed, while the backward natural frequency (blue dashed line) decreases with the increase of the rotating speed. The intersections of the curves related to the natural frequencies with the straight line ωn=Ω¯ determine the critical rotating speeds of the rotating cutting tool.Figure 2
Natural frequency versus rotational speed (L/d = 10,θ = 0°).


### 3.1. Stability Lobe Diagram

By removing the nonlinear term and the harmonic cutting force in the right-side term of (17), the stability of the cutting system can be investigated. Figure 3 shows a stability lobe diagram where the Ω versus alim curve separates the space into two regions. Any (Ω, alim) pair that appears above the collective boundary indicates unstable milling process where regenerative chatter or self-excited vibration occurs, while any pair below the boundary is a stable milling process.Figure 3
Stability lobe diagram of the cutting system with composite cutting tool (L/d = 10,c¯ = 0.3, θ = 0°).


### 3.2. Primary Resonance Response

#### 3.2.1. Case I

The numerical solutions for the forced vibration responses of the cutting system with composite cutting tool are presented in Figure4. As shown in Figure 4, by increasing detuning parameters at point A, the amplitude a1 gradually increases until reaching point B. a1 jumps downward from point B to point C. Afterward, a1 gradually drops to point D while continuing increasing the detuning parameter. On the other hand, as the detuning parameter decreases from point D to point E, a1 increases and jumps upward from point E to point F. As the detuning parameter further decreases, a1 drops until arriving at point A. By considering the nonlinearity of higher-order bending deformation, the forced response curve of the cutting system deviates toward the right, suggesting hard spring vibration behavior of Duffing type oscillator. As the detuning parameter varies between σF and σB, the cutting system has three steady-state vibration amplitudes, p1, p2, and p3, among which p1 and p3 are stable and p2 is unstable. This means that the cutting system works in an unstable cutting state. Therefore, cutting conditions and consequently initial conditions should be adjusted such that the stable branch with less vibration amplitudes is realized in practice.Figure 4
Frequency response curve in case I.Figures5–9 show the effects of ratio of length to diameter, structural damping, cutting force, ply angle, and rotating speed on frequency response curve, respectively, for the case Ω¯=ω1. As shown in Figure 5, when ratio of length to diameter increases, vibration amplitudes decrease and frequency response curve bends more strongly toward left. This is physically expected, because the nonlinear stiffness coefficient λ=EAequivL2/2EIequiv is proportional to ratio of length to diameter according to (12). As is observed and physically expected from Figure 6, the increase in the structural damping leads to the decrease in vibration amplitudes. In order to study the effect of cutting force coefficients, we define ς2=Kfς¯¯2, η2=Kfη¯¯2, with the parameters ς¯¯2 and η¯¯2 given in Table 2. As is shown in Figure 7, by increasing cutting forces, vibration amplitudes increase. Figure 8 shows that the primary resonance amplitude increases with the increase of the ply angle because the elastic modulus E11 along the longitudinal direction of the fiber is significantly larger than the transverse elastic modulus E22 (as shown in Table 1). Therefore, as the ply angle is greater, the equivalent bending stiffness (EI)equiv of the composite cutting tool is smaller; thus, the nondimensional cutting force coefficients are greater (as shown in (11)). As a result, the amplitude of the primary resonance response is larger.Figure 5
The effect of ratio of length to diameter on frequency response curve (c¯ = 0.01, θ = 0°, case I).Figure 6
The effect of structural damping on frequency response curve (L/d = 10,θ = 0°, case I).Figure 7
The effect of cutting force coefficient on frequency response curve (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure 8
The effect of ply angle on frequency response curve (L/d = 10,c¯ = 0.01, case I).Figure 9
The effect of rotating speed on frequency response curve (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure9 shows that vibration amplitudes decrease with the decrease of rotating speed. This is physically expected, because the equivalent damping of the cutting system (which will be introduced in the subsequent section) increases with the decrease of rotating speed due to the damping effect from regenerative chatter mechanism.Figures10–14 show the amplitude versus damping coefficient with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, and rotating speeds, respectively, for the case Ω¯=ω1. From these figures, it can be seen that for some values of L/d, σ, Kf, or θ, there are multivalued curves. For example, when L/d = 10 and c¯ ＜0.45, as is shown in Figure 10, the system has two stable and one unstable branches, but for L/d = 10 and c¯ ＞0.45, there exists only one stable branch. For large values of c¯ and large values of L/d, curves are always single-valued.Figure 10
Amplitude versus damping coefficient with different ratios of length to diameter (σ = 2, c¯ = 0.01, θ = 0°, case I).Figure 11
Amplitude versus damping coefficient with different detuning parameter values (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure 12
Amplitude versus damping coefficient with different cutting force coefficients (L/d = 10,σ = 2, θ = 0°, case I).Figure 13
Amplitude versus damping coefficient with different ply angles (L/d = 10,σ = 2, case I).Figure 14
Amplitude versus damping coefficient with different rotating speed (L/d = 10,σ = 2, case I).Figures15–19 show the amplitude versus ply angle with different ratios of length to diameter, detuning parameter values, cutting force coefficients, damping coefficients, and rotating speeds, respectively (Ω¯=ω1). Similar to the cases in Figure 10～14, again, for some values of L/d, σ, Kf, or c¯, multivalued curves can be observed.Figure 15
Amplitude versus ply angle with different ratios of length to diameter (c¯ = 0.01, σ = -1.5, case I).Figure 16
Amplitude versus ply angle with different detuning parameter values (L/d = 10,c¯ = 0.01, case I).Figure 17
Amplitude versus ply angle with different cutting force coefficients (L/d = 10,c¯ = 0.01, σ = 2, case I).Figure 18
Amplitude versus ply angle with different damping coefficients (L/d = 10,σ = 2, case I).Figure 19
Amplitude versus ply angle with different rotating speeds (L/d = 10,c¯ = 0.01, σ = 2, case I).Figures20–24 show the amplitude versus cutting force coefficient with different ratios of length to diameter, detuning parameter values, ply angle, damping coefficients, and rotating speeds, respectively. Again, when cutting force coefficient is less than the specified value, multivalued curves can be obtained for some values of L/d, σ, θ, c¯, or Ω¯.Figure 20
Amplitude versus cutting force coefficient with different ratios of length to diameter (θ = 0°, c¯ = 0.01, σ = 2, case I).Figure 21
Amplitude versus cutting force coefficient with different detuning parameter values (θ = 0°, c¯ = 0.01, L/d = 10, case I).Figure 22
Amplitude versus cutting force coefficient with different ply angles (σ = 2, c¯ = 0.01, L/d = 10, case I).Figure 23
Amplitude versus cutting force coefficient with different damping coefficients (σ = 2, θ = 0°, L/d = 10, case I).Figure 24
Amplitude versus cutting force coefficient with different rotating speeds (σ = 2, c¯ = 0.01, θ = 0°,L/d = 10, case I).Figures25–27 show the effects of ply sequences (see Table 3) on the frequency response curve, the amplitude versus dumping coefficient, and the amplitude versus cutting force coefficient, respectively. It can be seen from Figure 25 that the configuration A leads to the largest vibration amplitude, while the configuration C leads to the smallest vibration amplitude. This phenomenon can also be found in Figure 8 where the vibration amplitude increases with the ply angle when using the same configuration. Again, it can be seen from Figures 26 and 27 that multivalued curves can be observed for some values of L/d, σ, Kf, or c¯.Figure 25
The frequency response curve with different ply sequences (c¯ = 0.05, L/d = 10, case I).Figure 26
Amplitude versus dumping coefficient with different ply sequences (σ = 2, L/d = 10, case I).Figure 27
Amplitude versus cutting force coefficient with different ply sequences (c¯ = 0.01, σ = 2, L/d = 10, case I).Table 3
Configuration of composite.
ConfigurationSequenceA[90°/90°/0°]sB[0°/90°/0°/0°/90°/0°]C[0°]6


#### 3.2.2. Case II

Figures28–33 show the effects of ratio of length to diameter, structural damping, cutting force, ply angle, ply sequence, and rotating speed on frequency response curve, respectively, for the case Ω¯=ω2. As is observed, the frequency responses and effect of various parameters for the cutting system are similar to those in case I; as ratio of length to diameter and structural damping increase, vibration amplitudes decrease. Furthermore, as cutting force coefficient, ply angle, and rotating speed increase, vibration amplitudes increase. However, for this case, vibration amplitudes are less than those of case I.Figure 28
The effect of ratio of length to diameter on frequency response curve (θ = 0°, c¯ = 0.01, case II).Figure 29
The effect of structural damping on frequency response curve (L/d = 10,θ = 0°, case II).Figure 30
The effect of cutting force coefficient on frequency response curve (L/d = 10,θ = 0°, c¯ = 0.01, case II).Figure 31
The effect of ply angle on frequency response curve (L/d = 10,c¯ = 0.01, case II).Figure 32
The effect of rotating speed on frequency response curve (L/d = 10,c¯ = 0.01, case II).Figure 33
Frequency response curve with different ply configurations (c¯ = 0.05, L/d = 10, case II).It should be noted that the equivalent damping of the cutting system is composed of the structural damping and the damping from regenerative chatter mechanism. The equivalent damping coefficient of case I can be expressed as(45)cequiv1=c¯−ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯/ω1∫01ϕ12xdx.In addition, the equivalent damping coefficient of case II can be expressed as(46)cequiv2=c¯+ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯/ω2∫01ϕ12xdx.Figures34–38 show the effect of different parameters on the equivalent damping coefficients. As is observed, the equivalent damping coefficients of case I are always larger than those of case II. This explains the reason why vibration amplitudes of case I are larger than those of case II. It can be also seen that the equivalent damping coefficients increase with ratio of length to diameter, ply angle, cutting force, and structural damping but decrease with rotating speed. Therefore, the increase of rotating speed leads to large vibration amplitudes, as shown in Figures 9 and 32. When rotating speed approaches infinity, the equivalent damping coefficient approaches the structural damping c¯ (see Figure 38).Figure 34
The equivalent damping coefficients versus ratio of length to diameter for cases I and II.Figure 35
The equivalent damping coefficients versus ply angle for cases I and II.Figure 36
The equivalent damping coefficients versus cutting force for cases I and II.Figure 37
The equivalent damping coefficients versus structural damping for cases I and II.Figure 38
The equivalent damping coefficients versus rotating speed for cases I and II.Moreover, for some values ofL/d, σ, Kf, c¯, θ, or Ω¯, multivalued solutions can be found from the amplitude versus damping coefficient with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, and rotating speeds, respectively (not shown).Figure39 shows the frequency response curve for the four cases. As is observed, the curves of cases I and III are very similar, except that the amplitude in case I is larger than that in case III. Likewise, the curves of cases II and IV are very similar, except that the amplitude in case II is larger than that in case IV.Figure 39
The frequency response curve for the four cases (L/d = 10,c¯ = 0.01, θ = 0°, Ω¯ = 100).


#### 3.2.3. Case III

For this case, the equivalent coefficient is identical with that in case I, but the amplitude of excitation forcecfς12+cfη12/4 in this case is lower than ς¯22+η¯22 in case I (e.g., using the same parameters L/d = 10, θ = 90°, and Kf = 3, the amplitudes of excitation force are 0.2695 and 0.2532 for cases I and III, respectively). Therefore, generally, the behavior of nonlinear forced vibration in this case (see Figure 39) is similar to case I, despite having less vibration amplitudes compared with case I.


#### 3.2.4. Case IV

For this case, the equivalent coefficient is identical with that in case II, and the amplitude of excitation force in this case is lower than that in in case II. Thus, under the same conditions, vibration amplitudes are generally lower than those in case II. The results of case IV are not shown to save space.In order to validate the calculated approximate solution using multiscale method, the numerical simulation results according to (17) are also displayed in Figures 40 and 41. These two calculated results using different methods are in good consistency.Figure 40
Comparison of results between the method of multiple scales and the numerical simulation (case I, L/d = 10,c¯ = 0.03, θ = 0°, Ω¯ = 50, Kf = 1).Figure 41
Comparison of results between the method of multiple scales and the numerical simulation (case II, L/d = 10,c¯ = 0.03, θ = 0°, Ω¯ = 100, Kf = 1).


## 3.1. Stability Lobe Diagram

By removing the nonlinear term and the harmonic cutting force in the right-side term of (17), the stability of the cutting system can be investigated. Figure 3 shows a stability lobe diagram where the Ω versus alim curve separates the space into two regions. Any (Ω, alim) pair that appears above the collective boundary indicates unstable milling process where regenerative chatter or self-excited vibration occurs, while any pair below the boundary is a stable milling process.Figure 3
Stability lobe diagram of the cutting system with composite cutting tool (L/d = 10,c¯ = 0.3, θ = 0°).


## 3.2. Primary Resonance Response

### 3.2.1. Case I

The numerical solutions for the forced vibration responses of the cutting system with composite cutting tool are presented in Figure4. As shown in Figure 4, by increasing detuning parameters at point A, the amplitude a1 gradually increases until reaching point B. a1 jumps downward from point B to point C. Afterward, a1 gradually drops to point D while continuing increasing the detuning parameter. On the other hand, as the detuning parameter decreases from point D to point E, a1 increases and jumps upward from point E to point F. As the detuning parameter further decreases, a1 drops until arriving at point A. By considering the nonlinearity of higher-order bending deformation, the forced response curve of the cutting system deviates toward the right, suggesting hard spring vibration behavior of Duffing type oscillator. As the detuning parameter varies between σF and σB, the cutting system has three steady-state vibration amplitudes, p1, p2, and p3, among which p1 and p3 are stable and p2 is unstable. This means that the cutting system works in an unstable cutting state. Therefore, cutting conditions and consequently initial conditions should be adjusted such that the stable branch with less vibration amplitudes is realized in practice.Figure 4
Frequency response curve in case I.Figures5–9 show the effects of ratio of length to diameter, structural damping, cutting force, ply angle, and rotating speed on frequency response curve, respectively, for the case Ω¯=ω1. As shown in Figure 5, when ratio of length to diameter increases, vibration amplitudes decrease and frequency response curve bends more strongly toward left. This is physically expected, because the nonlinear stiffness coefficient λ=EAequivL2/2EIequiv is proportional to ratio of length to diameter according to (12). As is observed and physically expected from Figure 6, the increase in the structural damping leads to the decrease in vibration amplitudes. In order to study the effect of cutting force coefficients, we define ς2=Kfς¯¯2, η2=Kfη¯¯2, with the parameters ς¯¯2 and η¯¯2 given in Table 2. As is shown in Figure 7, by increasing cutting forces, vibration amplitudes increase. Figure 8 shows that the primary resonance amplitude increases with the increase of the ply angle because the elastic modulus E11 along the longitudinal direction of the fiber is significantly larger than the transverse elastic modulus E22 (as shown in Table 1). Therefore, as the ply angle is greater, the equivalent bending stiffness (EI)equiv of the composite cutting tool is smaller; thus, the nondimensional cutting force coefficients are greater (as shown in (11)). As a result, the amplitude of the primary resonance response is larger.Figure 5
The effect of ratio of length to diameter on frequency response curve (c¯ = 0.01, θ = 0°, case I).Figure 6
The effect of structural damping on frequency response curve (L/d = 10,θ = 0°, case I).Figure 7
The effect of cutting force coefficient on frequency response curve (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure 8
The effect of ply angle on frequency response curve (L/d = 10,c¯ = 0.01, case I).Figure 9
The effect of rotating speed on frequency response curve (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure9 shows that vibration amplitudes decrease with the decrease of rotating speed. This is physically expected, because the equivalent damping of the cutting system (which will be introduced in the subsequent section) increases with the decrease of rotating speed due to the damping effect from regenerative chatter mechanism.Figures10–14 show the amplitude versus damping coefficient with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, and rotating speeds, respectively, for the case Ω¯=ω1. From these figures, it can be seen that for some values of L/d, σ, Kf, or θ, there are multivalued curves. For example, when L/d = 10 and c¯ ＜0.45, as is shown in Figure 10, the system has two stable and one unstable branches, but for L/d = 10 and c¯ ＞0.45, there exists only one stable branch. For large values of c¯ and large values of L/d, curves are always single-valued.Figure 10
Amplitude versus damping coefficient with different ratios of length to diameter (σ = 2, c¯ = 0.01, θ = 0°, case I).Figure 11
Amplitude versus damping coefficient with different detuning parameter values (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure 12
Amplitude versus damping coefficient with different cutting force coefficients (L/d = 10,σ = 2, θ = 0°, case I).Figure 13
Amplitude versus damping coefficient with different ply angles (L/d = 10,σ = 2, case I).Figure 14
Amplitude versus damping coefficient with different rotating speed (L/d = 10,σ = 2, case I).Figures15–19 show the amplitude versus ply angle with different ratios of length to diameter, detuning parameter values, cutting force coefficients, damping coefficients, and rotating speeds, respectively (Ω¯=ω1). Similar to the cases in Figure 10～14, again, for some values of L/d, σ, Kf, or c¯, multivalued curves can be observed.Figure 15
Amplitude versus ply angle with different ratios of length to diameter (c¯ = 0.01, σ = -1.5, case I).Figure 16
Amplitude versus ply angle with different detuning parameter values (L/d = 10,c¯ = 0.01, case I).Figure 17
Amplitude versus ply angle with different cutting force coefficients (L/d = 10,c¯ = 0.01, σ = 2, case I).Figure 18
Amplitude versus ply angle with different damping coefficients (L/d = 10,σ = 2, case I).Figure 19
Amplitude versus ply angle with different rotating speeds (L/d = 10,c¯ = 0.01, σ = 2, case I).Figures20–24 show the amplitude versus cutting force coefficient with different ratios of length to diameter, detuning parameter values, ply angle, damping coefficients, and rotating speeds, respectively. Again, when cutting force coefficient is less than the specified value, multivalued curves can be obtained for some values of L/d, σ, θ, c¯, or Ω¯.Figure 20
Amplitude versus cutting force coefficient with different ratios of length to diameter (θ = 0°, c¯ = 0.01, σ = 2, case I).Figure 21
Amplitude versus cutting force coefficient with different detuning parameter values (θ = 0°, c¯ = 0.01, L/d = 10, case I).Figure 22
Amplitude versus cutting force coefficient with different ply angles (σ = 2, c¯ = 0.01, L/d = 10, case I).Figure 23
Amplitude versus cutting force coefficient with different damping coefficients (σ = 2, θ = 0°, L/d = 10, case I).Figure 24
Amplitude versus cutting force coefficient with different rotating speeds (σ = 2, c¯ = 0.01, θ = 0°,L/d = 10, case I).Figures25–27 show the effects of ply sequences (see Table 3) on the frequency response curve, the amplitude versus dumping coefficient, and the amplitude versus cutting force coefficient, respectively. It can be seen from Figure 25 that the configuration A leads to the largest vibration amplitude, while the configuration C leads to the smallest vibration amplitude. This phenomenon can also be found in Figure 8 where the vibration amplitude increases with the ply angle when using the same configuration. Again, it can be seen from Figures 26 and 27 that multivalued curves can be observed for some values of L/d, σ, Kf, or c¯.Figure 25
The frequency response curve with different ply sequences (c¯ = 0.05, L/d = 10, case I).Figure 26
Amplitude versus dumping coefficient with different ply sequences (σ = 2, L/d = 10, case I).Figure 27
Amplitude versus cutting force coefficient with different ply sequences (c¯ = 0.01, σ = 2, L/d = 10, case I).Table 3
Configuration of composite.
ConfigurationSequenceA[90°/90°/0°]sB[0°/90°/0°/0°/90°/0°]C[0°]6


### 3.2.2. Case II

Figures28–33 show the effects of ratio of length to diameter, structural damping, cutting force, ply angle, ply sequence, and rotating speed on frequency response curve, respectively, for the case Ω¯=ω2. As is observed, the frequency responses and effect of various parameters for the cutting system are similar to those in case I; as ratio of length to diameter and structural damping increase, vibration amplitudes decrease. Furthermore, as cutting force coefficient, ply angle, and rotating speed increase, vibration amplitudes increase. However, for this case, vibration amplitudes are less than those of case I.Figure 28
The effect of ratio of length to diameter on frequency response curve (θ = 0°, c¯ = 0.01, case II).Figure 29
The effect of structural damping on frequency response curve (L/d = 10,θ = 0°, case II).Figure 30
The effect of cutting force coefficient on frequency response curve (L/d = 10,θ = 0°, c¯ = 0.01, case II).Figure 31
The effect of ply angle on frequency response curve (L/d = 10,c¯ = 0.01, case II).Figure 32
The effect of rotating speed on frequency response curve (L/d = 10,c¯ = 0.01, case II).Figure 33
Frequency response curve with different ply configurations (c¯ = 0.05, L/d = 10, case II).It should be noted that the equivalent damping of the cutting system is composed of the structural damping and the damping from regenerative chatter mechanism. The equivalent damping coefficient of case I can be expressed as(45)cequiv1=c¯−ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯/ω1∫01ϕ12xdx.In addition, the equivalent damping coefficient of case II can be expressed as(46)cequiv2=c¯+ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯/ω2∫01ϕ12xdx.Figures34–38 show the effect of different parameters on the equivalent damping coefficients. As is observed, the equivalent damping coefficients of case I are always larger than those of case II. This explains the reason why vibration amplitudes of case I are larger than those of case II. It can be also seen that the equivalent damping coefficients increase with ratio of length to diameter, ply angle, cutting force, and structural damping but decrease with rotating speed. Therefore, the increase of rotating speed leads to large vibration amplitudes, as shown in Figures 9 and 32. When rotating speed approaches infinity, the equivalent damping coefficient approaches the structural damping c¯ (see Figure 38).Figure 34
The equivalent damping coefficients versus ratio of length to diameter for cases I and II.Figure 35
The equivalent damping coefficients versus ply angle for cases I and II.Figure 36
The equivalent damping coefficients versus cutting force for cases I and II.Figure 37
The equivalent damping coefficients versus structural damping for cases I and II.Figure 38
The equivalent damping coefficients versus rotating speed for cases I and II.Moreover, for some values ofL/d, σ, Kf, c¯, θ, or Ω¯, multivalued solutions can be found from the amplitude versus damping coefficient with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, and rotating speeds, respectively (not shown).Figure39 shows the frequency response curve for the four cases. As is observed, the curves of cases I and III are very similar, except that the amplitude in case I is larger than that in case III. Likewise, the curves of cases II and IV are very similar, except that the amplitude in case II is larger than that in case IV.Figure 39
The frequency response curve for the four cases (L/d = 10,c¯ = 0.01, θ = 0°, Ω¯ = 100).


### 3.2.3. Case III

For this case, the equivalent coefficient is identical with that in case I, but the amplitude of excitation forcecfς12+cfη12/4 in this case is lower than ς¯22+η¯22 in case I (e.g., using the same parameters L/d = 10, θ = 90°, and Kf = 3, the amplitudes of excitation force are 0.2695 and 0.2532 for cases I and III, respectively). Therefore, generally, the behavior of nonlinear forced vibration in this case (see Figure 39) is similar to case I, despite having less vibration amplitudes compared with case I.


### 3.2.4. Case IV

For this case, the equivalent coefficient is identical with that in case II, and the amplitude of excitation force in this case is lower than that in in case II. Thus, under the same conditions, vibration amplitudes are generally lower than those in case II. The results of case IV are not shown to save space.In order to validate the calculated approximate solution using multiscale method, the numerical simulation results according to (17) are also displayed in Figures 40 and 41. These two calculated results using different methods are in good consistency.Figure 40
Comparison of results between the method of multiple scales and the numerical simulation (case I, L/d = 10,c¯ = 0.03, θ = 0°, Ω¯ = 50, Kf = 1).Figure 41
Comparison of results between the method of multiple scales and the numerical simulation (case II, L/d = 10,c¯ = 0.03, θ = 0°, Ω¯ = 100, Kf = 1).


## 3.2.1. Case I

The numerical solutions for the forced vibration responses of the cutting system with composite cutting tool are presented in Figure4. As shown in Figure 4, by increasing detuning parameters at point A, the amplitude a1 gradually increases until reaching point B. a1 jumps downward from point B to point C. Afterward, a1 gradually drops to point D while continuing increasing the detuning parameter. On the other hand, as the detuning parameter decreases from point D to point E, a1 increases and jumps upward from point E to point F. As the detuning parameter further decreases, a1 drops until arriving at point A. By considering the nonlinearity of higher-order bending deformation, the forced response curve of the cutting system deviates toward the right, suggesting hard spring vibration behavior of Duffing type oscillator. As the detuning parameter varies between σF and σB, the cutting system has three steady-state vibration amplitudes, p1, p2, and p3, among which p1 and p3 are stable and p2 is unstable. This means that the cutting system works in an unstable cutting state. Therefore, cutting conditions and consequently initial conditions should be adjusted such that the stable branch with less vibration amplitudes is realized in practice.Figure 4
Frequency response curve in case I.Figures5–9 show the effects of ratio of length to diameter, structural damping, cutting force, ply angle, and rotating speed on frequency response curve, respectively, for the case Ω¯=ω1. As shown in Figure 5, when ratio of length to diameter increases, vibration amplitudes decrease and frequency response curve bends more strongly toward left. This is physically expected, because the nonlinear stiffness coefficient λ=EAequivL2/2EIequiv is proportional to ratio of length to diameter according to (12). As is observed and physically expected from Figure 6, the increase in the structural damping leads to the decrease in vibration amplitudes. In order to study the effect of cutting force coefficients, we define ς2=Kfς¯¯2, η2=Kfη¯¯2, with the parameters ς¯¯2 and η¯¯2 given in Table 2. As is shown in Figure 7, by increasing cutting forces, vibration amplitudes increase. Figure 8 shows that the primary resonance amplitude increases with the increase of the ply angle because the elastic modulus E11 along the longitudinal direction of the fiber is significantly larger than the transverse elastic modulus E22 (as shown in Table 1). Therefore, as the ply angle is greater, the equivalent bending stiffness (EI)equiv of the composite cutting tool is smaller; thus, the nondimensional cutting force coefficients are greater (as shown in (11)). As a result, the amplitude of the primary resonance response is larger.Figure 5
The effect of ratio of length to diameter on frequency response curve (c¯ = 0.01, θ = 0°, case I).Figure 6
The effect of structural damping on frequency response curve (L/d = 10,θ = 0°, case I).Figure 7
The effect of cutting force coefficient on frequency response curve (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure 8
The effect of ply angle on frequency response curve (L/d = 10,c¯ = 0.01, case I).Figure 9
The effect of rotating speed on frequency response curve (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure9 shows that vibration amplitudes decrease with the decrease of rotating speed. This is physically expected, because the equivalent damping of the cutting system (which will be introduced in the subsequent section) increases with the decrease of rotating speed due to the damping effect from regenerative chatter mechanism.Figures10–14 show the amplitude versus damping coefficient with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, and rotating speeds, respectively, for the case Ω¯=ω1. From these figures, it can be seen that for some values of L/d, σ, Kf, or θ, there are multivalued curves. For example, when L/d = 10 and c¯ ＜0.45, as is shown in Figure 10, the system has two stable and one unstable branches, but for L/d = 10 and c¯ ＞0.45, there exists only one stable branch. For large values of c¯ and large values of L/d, curves are always single-valued.Figure 10
Amplitude versus damping coefficient with different ratios of length to diameter (σ = 2, c¯ = 0.01, θ = 0°, case I).Figure 11
Amplitude versus damping coefficient with different detuning parameter values (L/d = 10,c¯ = 0.01, θ = 0°, case I).Figure 12
Amplitude versus damping coefficient with different cutting force coefficients (L/d = 10,σ = 2, θ = 0°, case I).Figure 13
Amplitude versus damping coefficient with different ply angles (L/d = 10,σ = 2, case I).Figure 14
Amplitude versus damping coefficient with different rotating speed (L/d = 10,σ = 2, case I).Figures15–19 show the amplitude versus ply angle with different ratios of length to diameter, detuning parameter values, cutting force coefficients, damping coefficients, and rotating speeds, respectively (Ω¯=ω1). Similar to the cases in Figure 10～14, again, for some values of L/d, σ, Kf, or c¯, multivalued curves can be observed.Figure 15
Amplitude versus ply angle with different ratios of length to diameter (c¯ = 0.01, σ = -1.5, case I).Figure 16
Amplitude versus ply angle with different detuning parameter values (L/d = 10,c¯ = 0.01, case I).Figure 17
Amplitude versus ply angle with different cutting force coefficients (L/d = 10,c¯ = 0.01, σ = 2, case I).Figure 18
Amplitude versus ply angle with different damping coefficients (L/d = 10,σ = 2, case I).Figure 19
Amplitude versus ply angle with different rotating speeds (L/d = 10,c¯ = 0.01, σ = 2, case I).Figures20–24 show the amplitude versus cutting force coefficient with different ratios of length to diameter, detuning parameter values, ply angle, damping coefficients, and rotating speeds, respectively. Again, when cutting force coefficient is less than the specified value, multivalued curves can be obtained for some values of L/d, σ, θ, c¯, or Ω¯.Figure 20
Amplitude versus cutting force coefficient with different ratios of length to diameter (θ = 0°, c¯ = 0.01, σ = 2, case I).Figure 21
Amplitude versus cutting force coefficient with different detuning parameter values (θ = 0°, c¯ = 0.01, L/d = 10, case I).Figure 22
Amplitude versus cutting force coefficient with different ply angles (σ = 2, c¯ = 0.01, L/d = 10, case I).Figure 23
Amplitude versus cutting force coefficient with different damping coefficients (σ = 2, θ = 0°, L/d = 10, case I).Figure 24
Amplitude versus cutting force coefficient with different rotating speeds (σ = 2, c¯ = 0.01, θ = 0°,L/d = 10, case I).Figures25–27 show the effects of ply sequences (see Table 3) on the frequency response curve, the amplitude versus dumping coefficient, and the amplitude versus cutting force coefficient, respectively. It can be seen from Figure 25 that the configuration A leads to the largest vibration amplitude, while the configuration C leads to the smallest vibration amplitude. This phenomenon can also be found in Figure 8 where the vibration amplitude increases with the ply angle when using the same configuration. Again, it can be seen from Figures 26 and 27 that multivalued curves can be observed for some values of L/d, σ, Kf, or c¯.Figure 25
The frequency response curve with different ply sequences (c¯ = 0.05, L/d = 10, case I).Figure 26
Amplitude versus dumping coefficient with different ply sequences (σ = 2, L/d = 10, case I).Figure 27
Amplitude versus cutting force coefficient with different ply sequences (c¯ = 0.01, σ = 2, L/d = 10, case I).Table 3
Configuration of composite.
ConfigurationSequenceA[90°/90°/0°]sB[0°/90°/0°/0°/90°/0°]C[0°]6


## 3.2.2. Case II

Figures28–33 show the effects of ratio of length to diameter, structural damping, cutting force, ply angle, ply sequence, and rotating speed on frequency response curve, respectively, for the case Ω¯=ω2. As is observed, the frequency responses and effect of various parameters for the cutting system are similar to those in case I; as ratio of length to diameter and structural damping increase, vibration amplitudes decrease. Furthermore, as cutting force coefficient, ply angle, and rotating speed increase, vibration amplitudes increase. However, for this case, vibration amplitudes are less than those of case I.Figure 28
The effect of ratio of length to diameter on frequency response curve (θ = 0°, c¯ = 0.01, case II).Figure 29
The effect of structural damping on frequency response curve (L/d = 10,θ = 0°, case II).Figure 30
The effect of cutting force coefficient on frequency response curve (L/d = 10,θ = 0°, c¯ = 0.01, case II).Figure 31
The effect of ply angle on frequency response curve (L/d = 10,c¯ = 0.01, case II).Figure 32
The effect of rotating speed on frequency response curve (L/d = 10,c¯ = 0.01, case II).Figure 33
Frequency response curve with different ply configurations (c¯ = 0.05, L/d = 10, case II).It should be noted that the equivalent damping of the cutting system is composed of the structural damping and the damping from regenerative chatter mechanism. The equivalent damping coefficient of case I can be expressed as(45)cequiv1=c¯−ϕ121N4πα¯0′+β¯01−cosω1τ¯+β¯0′−α¯0sinω1τ¯/ω1∫01ϕ12xdx.In addition, the equivalent damping coefficient of case II can be expressed as(46)cequiv2=c¯+ϕ121N4πα¯0′+β¯01−cosω2τ¯+α¯0−β¯0′sinω2τ¯/ω2∫01ϕ12xdx.Figures34–38 show the effect of different parameters on the equivalent damping coefficients. As is observed, the equivalent damping coefficients of case I are always larger than those of case II. This explains the reason why vibration amplitudes of case I are larger than those of case II. It can be also seen that the equivalent damping coefficients increase with ratio of length to diameter, ply angle, cutting force, and structural damping but decrease with rotating speed. Therefore, the increase of rotating speed leads to large vibration amplitudes, as shown in Figures 9 and 32. When rotating speed approaches infinity, the equivalent damping coefficient approaches the structural damping c¯ (see Figure 38).Figure 34
The equivalent damping coefficients versus ratio of length to diameter for cases I and II.Figure 35
The equivalent damping coefficients versus ply angle for cases I and II.Figure 36
The equivalent damping coefficients versus cutting force for cases I and II.Figure 37
The equivalent damping coefficients versus structural damping for cases I and II.Figure 38
The equivalent damping coefficients versus rotating speed for cases I and II.Moreover, for some values ofL/d, σ, Kf, c¯, θ, or Ω¯, multivalued solutions can be found from the amplitude versus damping coefficient with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, and rotating speeds, respectively (not shown).Figure39 shows the frequency response curve for the four cases. As is observed, the curves of cases I and III are very similar, except that the amplitude in case I is larger than that in case III. Likewise, the curves of cases II and IV are very similar, except that the amplitude in case II is larger than that in case IV.Figure 39
The frequency response curve for the four cases (L/d = 10,c¯ = 0.01, θ = 0°, Ω¯ = 100).


## 3.2.3. Case III

For this case, the equivalent coefficient is identical with that in case I, but the amplitude of excitation forcecfς12+cfη12/4 in this case is lower than ς¯22+η¯22 in case I (e.g., using the same parameters L/d = 10, θ = 90°, and Kf = 3, the amplitudes of excitation force are 0.2695 and 0.2532 for cases I and III, respectively). Therefore, generally, the behavior of nonlinear forced vibration in this case (see Figure 39) is similar to case I, despite having less vibration amplitudes compared with case I.


## 3.2.4. Case IV

For this case, the equivalent coefficient is identical with that in case II, and the amplitude of excitation force in this case is lower than that in in case II. Thus, under the same conditions, vibration amplitudes are generally lower than those in case II. The results of case IV are not shown to save space.In order to validate the calculated approximate solution using multiscale method, the numerical simulation results according to (17) are also displayed in Figures 40 and 41. These two calculated results using different methods are in good consistency.Figure 40
Comparison of results between the method of multiple scales and the numerical simulation (case I, L/d = 10,c¯ = 0.03, θ = 0°, Ω¯ = 50, Kf = 1).Figure 41
Comparison of results between the method of multiple scales and the numerical simulation (case II, L/d = 10,c¯ = 0.03, θ = 0°, Ω¯ = 100, Kf = 1).


## 4. Conclusions

In this paper, a nonlinear dynamic model of the milling process with composite rotational cutting tool is presented. The cutting tool is simplified to a nonplanar bending Euler–Bernoulli beam, and the structural nonlinearity is attributed to the higher-order bending deformation of the cutting tool. A linear cutting force model considering regenerative mechanism, which includes time-delay terms and the first and second harmonic terms, is used. In addition, structural damping and gyroscopic effect are also considered. Nonlinear ordinary differential equations in the generalized coordinates are derived using the Hamilton principle and the Galerkin method. The lobe diagram of the cutting system is obtained. The multiscale method is used to construct the analytical approximate solutions of the forced vibration frequency response. It is found that, for all cases of primary resonance and superharmonic resonance, excitation of the forward (backward) mode does not produce the vibration responses in the backward (forward) mode because of not taking into account internal resonance in the proposed model. For primary and superharmonic resonance conditions, the effective nonlinearity of the cutting system with higher-order bending deformation rotating composite cutting tool is of a hard type. Jump phenomenon and multivalued solutions can be observed.The influence of higher-order bending deformation causes the frequency response curve to bend more strongly toward right when ratio of length to diameter increases. The study also finds that the vibration amplitudes of the cutting tool increase with cutting force or ply angle. It can be seen that the damping of the cutting system has a significant influence on the vibration amplitude of the composite cutting tool. The damping capacities of the cutting system can be measured by the two different equivalent damping coefficientscequiv1 and cequiv2 which include the structural damping and the damping from regenerative chatter mechanism. Cases I and III have the same equivalent damping coefficient cequiv1, and cases II and IV have the same equivalent damping coefficient cequiv2. It has been found that cequiv2 is larger than cequiv1, so the vibration amplitudes of case I are larger than those of case II. In addition, for all resonant cases, as rotating speed decreases or structural damping increases, the equivalent damping coefficient is strengthened, and consequently less vibration amplitudes are observed. Therefore, from this point of view, low rotating speeds are preferable to keep vibration amplitudes small. [38].



---

*Source: 8699218-2021-12-02.xml*