# Forced Vibration in Cutting Process considering the Nonlinear Curvature and Inertia of a Rotating Composite Cutter Bar

Authors: Yongsheng Ren; Donghui Yao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Forced vibration of the cutting system with a three-dimensional composite cutter bar is investigated. The composite cutter bar is simplified as a rotating cantilever shaft which is subjected to a cutting force including regenerative delay effects and harmonic exciting items. The nonlinear curvature and inertia of the cutter bar are taken into account based on inextensible assumption. The effects of the moment of inertia, gyroscopic effect, and internal and external damping are also considered, but shear deformation is neglected. Equation of motion is derived based on Hamiltonʼs extended principle and discretized by the Galerkin method. The analytical solutions of the steady-state response of the cutting system are constructed by using the method of multiple scales. The response of the cutting system is studied for primary and superharmonic resonances. The effects of length-to-diameter ratio, damping ratio, cutting force coefficients, ply angle, rotating speed, and internal and external damping are investigated. The results show that nonlinear curvature and inertia imposed a significant effect on the dynamic behavior of the cutting process. The equivalent nonlinearity of the cutting system shows hard spring characteristics. Multiple solutions and jumping phenomenon of typical Duffing system are found in forced response curves.


---

## Body

## 1. Introduction

In boring or milling machining applications, the chatter problem during the cutting process is increasingly more significant because, for these operations, a flexible cutter bar is used. Chatter is detrimental for surface finish and dimensional accuracy.The existing modeling and analysis of cutting chatter were established within the framework of linear theory. Although the chatter phenomenon was found by Taylor as early as in 1907 [1] and identified as the limit of productivity, the theoretical explanations for chatter generation were proposed after five to six decades. Regenerative chatter occurs as a result of dynamic interaction between previous and current cuts and remains the most comprehensive explanation for chatter. Several efforts on modeling chatter under linear analysis have been performed, particularly for the purpose of prediction of chatter occurrence and the conditions under which such instability occurs in the milling process [2–8].Since linear theory cannot predict some important dynamic behaviors in cutting process, nonlinear modeling of the cutting system has received increasing attention. Hanna and Tobias first proposed a time-delay nonlinear model including quadratic and cubic structural stiffness and cutting force [9]; this study triggered a strong interest from many researchers to analyze the global dynamics of the problem [10–18].Basically, two different approaches for modeling slender cutter bar can be found in the literature [9–18]. One is purely analytical; that is, it is based on physical and dynamic equations governing the dynamic behavior of cutter bar [17, 18]. The other is purely empirical: the simple one-dof or two-dof mass-spring-damper model is used; modal parameters of the cutter bar are measured by performing experimental tests, without explicit connection to the underlying dynamic equations [9–16].In order to model the dynamic behaviors of a cutter bar, it is necessary to comprehensively examine the influence of various factors, such as cutter bar geometry and material, on the modal parameters. For this purpose, a great number of repetitive tests have to be performed if the empirical method is used. It will be quite time-consuming and low-efficiency. In view of this, as a more reliable and practical modeling method, dynamic modeling of the continuous system of cutter bar based on the analytical method is needed.In recent years, fiber-reinforced composite materials have received great interest in the structural design of cutter bars [19–24] because of a much higher specific stiffness and higher damping than conventional metal materials. The superior performance of the composite is beneficial for enhancing the stability against chatter.However, most of the works emphasized the design of composite boring bars [20–22] or analytical prediction of chatter stability in the cutting process [23–25]; until now, no work seems to have been reported towards the study for the nonlinear dynamics of the cutting system with the cutter bar made up of composite materials.The aim of this work is to present a dynamic model of the cutting process by considering the geometrical nonlinearity of a rotating composite cutter bar. The composite cutter bar is modeled as an inextensible continuous cantilever rotating shaft with nonlinear curvature and nonlinear inertia. The dynamic model includes the moment of inertia, gyroscopic effect, and external damping. It also is assumed that the cutter bar is subjected to a cutting force including linear regenerative terms and harmonic components. Neglecting the effect of shear deformation, the equations of motion are derived by using Hamilton principal. The Galerkin method is used to discretize the equations of motion. The approximate response solutions of the nonlinear cutting system under cutting force are derived by means of the method of multiple scales. Primary and superharmonic resonances are considered. Moreover, the effects of different parameters including ratio of length to diameter, damping ratio, cutting coefficient, ply angle, rotating speed, and external and internal damping on the forced response of the cutting system are examined.


## 2. Mathematical Model and Solution

### 2.1. Kinetic Energy and Potential Energy

Figure1 displays a rotating composite cutter bar with a length of L. The cutter bar is fixed in one end while another end is free. The rectangular coordinate system (X, Y, Z) denotes an inertia coordinate system, while (x, y, z) denotes a local coordinate system. The inertial coordinate axis is consistent with the principal axis of the cross section of the cutter bar. The coordinate origin is located at x on the centerline of the deformed cutter bar. In addition, ux,t,vx,t, and wx,t denote displacements of point x on the deformed cutter bar along X, Y, and Z directions, respectively, and ϕx,t denotes the torsional deformation of cross section around x-axis.Figure 1
Schematic of the structure for the rotating composite cutter bar.The kinetic energy of a composite cutter bar can be written as [26](1)T=12∫0Lmu˙2+v˙2+w˙2+2ρIω12+ρIω22+ω32dx,where “.” represents partial derivative with respect to time t. m and ρI denote mass and cross-section inertia moment per unit length, respectively, which can be written as(2)m=π∑k=1Nρkrk+12−rk2,ρI=π4∑k=1Nρkrk+14−rk4,where N denotes ply numbers, ρk denotes the density of the k–th layer, and rk and rk+1 denote both inner and external diameters of the k-layer.The rotating angular velocities of the coordinate system (x, y, z) with respect to (X, Y, Z), ω1, ω2, and ω3 can be written as(3)ω1=φ˙−ψ˙ysinψz,ω2=ψ˙ysinφcosψz+ψ˙zcosφ,ω3=ψ˙ycosφcosψz−ψ˙zsinφ,where ψ=φ+Ωt. ψy and ψz denote the rotation angles of cross section around z-axis and y-axis, respectively, and ψz and ψy can be expressed as [26](4)ψz=sin−1v′1+u′2+v′2,ψy=−sin−1w′1+u′2+v′2+w′2,where “′” represents partial derivative with respect to time x.The elastic potential energy of a composite Euler–Bernoulli cutter bar can be written as(5)U=12∫Vσxεx+τxαγxαdV,where dV=rdrdαdx denotes a differential volume element in the cylindrical coordinate system and α and r are polar angle and polar diameter, respectively.By neglecting shear deformation, the stress-strain relations of thek–th layer expressed in the cylindrical coordinate system can be expressed as [27](6)σx=Q¯11εx+Q¯16γxα,(7)τxα=Q¯16εx+Q¯66γxα,where σx and τxα denote normal stress and shear stress of the point in the cylindrical coordinate x,r,α and Q¯iji,j=1,6 denotes the off-axis stiffness coefficient of a single layer of the composite.The strain-displacement equation can be written as(8)εx=ε+ρ2rsinα−ρ3rcosα,(9)γxα=rρ1,where ρii=1,2,3 denotes the curvature of the cutter bar. By virtue of Kirchhoffʼs kinetic analog [28], we can obtain the cutter barʼs curvatures ρ1, ρ2, and ρ3 by replacing the time derivatives ∂/∂t with the spatial derivatives ∂/∂x in the angular velocity expressions. Therefore, the curvatures are given by [29](10)ρ1=φ′−ψy′sinψz,ρ2=ψy′sinφcosψz+ψz′cosφ,ρ3=ψy′cosφcosψz−ψz′sinφ.The variation of the elastic potential energy can be written as(11)U=12∫Vσxεx+τxαγxαdV.By substituting equations (6)–(9) into equation (11), the following expression can be derived:(12)U=12∫0LA11ε2+2B16ρ1ε+D66ρ12+D66ρ22+ρ32dx,where(13)D11=π4∑k=1NQ¯11rk+14−rk4,D66=π2∑k=1NQ¯66rk+14−rk4,B16=2π3∑k=1NQ¯16rk+13−rk3,A11=π∑k=1NQ¯11rk+12−rk2.The strain along the elastic axis of a differential element dx is defined by(14)ε=1+u′2+v′2+w′2−1.Assuming that the composite cutter bar is inextensible in axial direction, that is, the strainε equals 0, the following expression can be derived:(15)u′=1−v′2−w′2−1≅−12v′2+w′2+⋅⋅⋅.Assuming that lateral displacementsv and w are of order Oe, e≪1, then u=Oe2. Substituting equation (4) into equations (3) and (10), expanding the outcomes in Taylor series, and retaining terms up to Oe3, one can obtain the angular velocities and curvatures up to Oe3. Substituting the resulted curvatures and angular velocities into equations (1) and (12) and using equation (15), the final expressions for kinetic and strain energy are obtained.


### 2.2. Virtual Work of the Cutting Force

The virtual work of the cutting force can be written as(16)δW=∫0LLvδv+Lwδwdx,where Lv=FyδDx−L, Lw=FzδDx−L, and δD denotes Delta function.The regenerative cutting force,Fy and Fz, can be expressed as [15](17)Fy=−Nα0v−vτ+β0w−wτ+γ02π+cfς1cos2Ωt−π/2−η1cos2Ωt+η12+η2cosΩt−π2+ς2cosΩt,Fz=Nα0′v−vτ+β0′w−wτ+γ0′2π+cfη1cos2Ωt−π/2+ς1cos2Ωt−ς12+ς2cosΩt−π2−η2cosΩt,where(18)α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−τ.Here, τ denotes delay time τ=2π/NΩ, N denotes the number of cutting teeth, (Ktc, Kte) and (Krc, Kre) are the cutting force coefficients along tangential and radial directions, respectively, a denotes the axial cutting depth, cf denotes the feed per tooth per revolution and corresponds to the static part of the chip thickness cfsinϕj, and vjt−τ−vjt denotes the dynamic chip thickness produced due to vibrations of the cutter bar at the present vjt and previous vjt−τ tooth periods (see Figure 2).Figure 2
Sketches of the milling process. (a) Cross-section view of the milling process. (b) Dynamic chip thickness.
(a)(b)


### 2.3. Equation of Motion in Milling Process

The Lagrangian is defined asL=T−U. Using equations (1) and (12) and introducing the Lagrange multiplier λ to enforce the inextensionality conditions [30], one has(19)L=12∫0Lmu˙2+v˙2+w˙2+2ρIω12+ρIω22+ω32−A11ε2+2B16ρ1ε+D66ρ12+D66ρ22+ρ32+λ1−1+u′2−v′2−w′2dx.To derive the governing equations of motion, we use Hamilton’s extended principle as(20)∫t1t2δL+δWdt=0.The bending-bending coupled nonlinear motion equation of a rotating composite cutter bar can be obtained as(21)mv¨+v′h′+ce+civ˙−ρI2Ωw˙″+v¨″+D11g1+ciΩw=Fyδx−L,(22)mw¨+w′h′+ce+ciw˙+ρI2Ωv˙″−w¨″+D11g2−ciΩv=Fzδx−L,where(23)g1=4v″v′v‴+w4w′v′+w′v″w‴+3w″v′w‴+w″2v″+v″3+v4v′2+v4,g2=4w″w′w‴+v4v′w′+v′w″v‴+3v″w′v‴+v″2w″+w″3+w4w′2+w4,h=∫Lx∫0xv˙′2+v¨′v′+w˙′2+w¨′w′dxdx.It should be noted that the torsional fundamental frequency of the cutting bar with circular cross section is much greater than the bending frequency; thus, the torsional inertia item can be neglected [29].Using the above assumption, the following expression can be obtained:(24)ϕ=−∫0xv″w′dx+⋯.Since the cutting bar is slender with small rotation inertia, the nonlinear item multiplied with the rotational inertia can also be neglected [26]. Accordingly, bending-bending-torsion coupled nonlinear motion equation can be appropriately simplified into equations (21) and (22).In equation (21), the terms g1 and g2 show the effect of the geometric nonlinearities, and the term h shows the effect of the inertia nonlinearities. In fact, h is given by h=−λ=m/2∫Lx∂2/∂t2∫0xv′2+w′2dxdx. One can interpret λ as an axial force acting at the cutter barʼs tip to prevent it from stretching.The termsce and ci are external and internal viscous damping coefficients, respectively. Damping is an important parameter which influences the dynamic behavior of the cutting process [23]. Two commonly used models in rotating system dynamics are external damping (nonrotating damping) and internal damping (rotating damping). The external damping is generally treated like viscous damping, but internal damping can be described by the viscous damping model or the hysteretic damping model [31, 32]. The damping mechanism is complex, particularly for hysteretic damping. Inclusion of hysteretic damping could be achieved using the complex stiffness; however, this leads to considerably increased complexity of the present dynamic model. It can be noted that hysteretic damping has a minor effect compared to viscous damping for an isotropic cantilevered shaft [31]. In light of above, the classical viscous damping model is employed in our work to describe the external and internal damping. ce is proportional to the speed in the inertial coordinate system and ci is proportional to the speed in the rotating coordinate system [23, 31, 32].Using the dimensionless transformation(25)x¯=xL,v¯=vL,w¯=wL,t¯=D11mL4t=γt,τ¯=γτ,c¯e=ceL2mD11,c¯i=ciL2mD11,ρ¯I¯=ρImL2,Ω¯=mL4D11Ω,nonlinear equation of motion in dimensionless form can be written as(26)v¨+ce+civ˙−2ΩρIw˙″−ρIv¨″+v′h′+g1+ciΩw=F¯yδx−1,(27)w¨+ce+ciw˙+2ΩρIv˙″−ρIw¨″+w′h′+g2−ciΩv=F¯zδx−1,where(28)F¯y=−Nα¯0v−vτ+β¯0w−wτ+γ¯02π+cfς1cos2Ωt−π/2−cfη1cos2Ωt+cfη12+η¯2cosΩt−π2+ς¯2cosΩt,F¯z=Nα¯0′v−vτ+β¯0′w−wτ+γ¯0′2π−cfη1cos2Ωt−π/2+cfς1cos2Ωt−cfς12+ς¯2cosΩt−π2−η¯2cosΩt.The coefficients in equation (26) can be defined as(29)α¯0=L3α0D11,β¯0=L3β0D11,γ¯0=L2γ0D11,cfς1=L2cfς1D11,cfη1=L2cfη1D11,α¯0′=L3α0′D11,β¯0′=L3β0′D11,γ¯0′=L2γ0′D11,ς¯2=L2ς2D11,η¯2=L2η2D11.For simplicity, the symbol “–” above all dimensionless variables in equations (26) and (27) has been removed.Before applying the multiple scales method, the partial differential equations of motion are discretized by using the Galerkin method.Let(30)vx,t=ϕ1xVt,wx,t=ϕ1xWt.For cantilever beam, the mode shapeϕ1x can be written as(31)φ1x=cosβ1Lx−coshβ1Lx−cosβ1L+coshβ1Lsinβ1L+sinhβ1Lsinβ1Lx−sinhβ1Lx,x∈0,1,where β1L=1.875.By substituting equation (31) into equations (26) and (27), taking the inner product of each equation with its corresponding mode shape, the following ordinary differential equations can be obtained:(32)A1V¨+A1ce+ciV˙−A2ρI2ΩW˙+V¨+A3V+A1ciΩW+A4VV˙2+VV¨+W˙2+WW¨+A5VV2+W2=∫01F¯yδx−1dx,A1W¨+A1ce+ciW˙+A2ρI2ΩV˙−W¨+A3W−A1ciΩV+A4WV˙2+VV¨+W˙2+WW¨+A5WV2+W2=∫01F¯zδx−1dx,in which(33)A1=∫01φ12dx,A2=∫01φ1φ1″dx,A3=∫01φ1φ14dx,A4=∫01φ1φ1′F1′dx,A5=∫01φ1F2dx,F1=∫1x∫0xφ1′2dxdx′,F2=4φ1′φ1″φ1‴+φ1″3+φ14φ1′2.To simplify the delayed expressions in the cutting force (28), the first-order Pade approximation as e−sτ≅1−sτ is used for delayed terms (i.e., a truncated Taylor series expansion for e−sτ with one delayed term is used) [15].Substituting the cutting forces in equation (28) into equations (26) and (27), the following expression can be derived:(34)A1−A2ρIV¨+μVV˙−2ΩρIA2W˙+A3V+A1ciΩW+A4VV˙2+VV¨+W˙2+WW¨+A5VV2+W2=−λV+F¯y2φ11,A1−A2ρIW¨+μWW˙+2ΩρIA2V˙+A3W−A1ciΩV+A4WV˙2+VV¨+W˙2+WW¨+A5WV2+W2=λW+F¯z2φ11,where F¯y2 and F¯z2 denote the sum of the second and third items in the expressions of cutting forces (28) and the parameters in equation (34) are defined as(35)μV=A1ce+ci+Nα¯0τφ1212π,μW=A1ce+ci−Nβ¯0′τφ1212π,λV=Nγ¯0φ112π,λW=Nγ¯0′φ112π.According to equation (35), the terms arisen from regenerative chatter mechanism (Nα¯0τφ121/2π and Nβ¯0τφ121/2π) are proportional to delay time τ and consequently inversely proportional to the spindle speed Ω. In this paper, a truncated Taylor series expansion for e−sτ is used to simplify the regenerative chatter mechanism. Therefore, it must be noticed that the results of this research are valid for small-time delays which occur essentially at high-speed machining. Here, for the considered case study, it can be shown that using the first-order truncated Taylor series for e−sτ and modifying only the damping coefficients are sufficient (as done in the works of Moradi et al. [15]).


### 2.4. Method of Multiple Scales

The multiple scales method is used to find an approximate solution for the proposed nonlinear dynamics of the milling process. In order to apply multiple scale method,V and W can be written in the expansion form as(36)Vt=εV1T0,T2+ε3V3T0,T2+⋯,Wt=εW1T0,T2+ε3W3T0,T2+⋯,where ε is a small dimensionless parameter. T0=t and T2=ε2t are fast and slow time scales, respectively. The damping and exciting force terms should be scaled, so that their effects are balanced with nonlinearities. Thus, μV, μW, ci, λV, λW, F¯y2, and F¯z2 are replaced with ε2μV, ε2μW, ε2ci, ε3λV, ε3λW, ε3F¯y2, and ε3F¯z2.Time derivatives become expansions in terms of the partial derivatives with respect toT0 and T2 by using the following chain rule:(37)ddt=∂∂T0+ε2∂∂T2+⋯;d2dt2=∂2∂T02+2ε2∂∂T2∂∂T0+⋯.By substituting equations (36) and (37) into equation (34) and equating the coefficients of ε and ε3, one hasOε(38)A1−A2ρI∂2V1∂T02−2A2ρIΩ∂W1∂T0+A3V1=0,A1−A2ρI∂2W1∂T02+2A2ρIΩ∂V1∂T0+A3W1=0.Oε3(39)A1−A2ρI∂2V3∂T02−2A2ρIΩ∂W3∂T0+A3V3=−μV∂V1∂T0−A1ciΩW1−2A1−A2ρI∂2V1∂T0∂T2+2A2ρIΩ∂W1∂T2−A4V1∂V1∂T02+V1∂2V1∂T02+∂W1∂T02+W1∂2W1∂T02−A5V1V12+W12−λV+F¯y2φ11,A1−A2ρI∂2W3∂T02+2A2ρIΩ∂V3∂T0+A3W3=−μW∂W1∂T0+A1ciΩV1−2A1−A2ρI∂2W1∂T0∂T2−2A2ρIΩ∂V1∂T2−A4W1∂V1∂T02+V1∂2V1∂T02+∂W1∂T02+W1∂2W1∂T02−A5W1V12+W12+λW+F¯z2φ11.The solution to equation (38) can be written as(40)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,where i=−1 denotes the imaginary unit; F1T2 and F2T2 are the complex-valued functions to be determined, respectively; F¯1T2 and F¯2T2 are complex conjugates; ω1 and ω2 denote forward and backward linear natural frequencies, respectively, and can be written as(41)ω1=A2ρIΩ+A2ρIΩ2+A1−A2ρIA3A1−A2ρI,ω2=−A2ρIΩ+A2ρIΩ2+A1−A2ρIA3A1−A2ρI.By substituting (40) into equation (39), one has(42)A1−A2ρI∂2V3∂T02−2A2ρIΩ∂W3∂T0+A3V3=P¯3eiω1T0+Q¯3eiω2T0−λV+F¯y2φ11+cc,A1−A2ρI∂2W3∂T02+2A2ρIΩ∂V3∂T0+A3W3=P¯4eiω1T0+Q¯4eiω2T0+λW+F¯z2φ11+cc,where cc denotes the complex conjugate.The coefficients of the harmonic items on the left end of equation (42) have the following form:(43)P¯3=ΓF1F¯2F22−4A5F¯1F12+iΛ1F1′2−iμVω1F1+iA1ciΩF1,Q¯3=ΓF1F¯1F22−4A5F¯2F22+iΛ2F2′2−iμVω2F2−iA1ciΩF2,P¯4=−iΓF1F¯2F22+i4A5F¯1F12+Λ1F1′2−μWω1F1+A1ciΩF1,Q¯4=iΓF1F¯1F22−i4A5F¯2F22−Λ2F2′2+μWω2F2+A1ciΩF2,where(44)Γ=A4ω12+ω222−16A5,Λ1=4−A1−A2ρIω1+A2ρIΩ,Λ2=−4A1−A2ρIω2+A2ρIΩ.In this paper, both types of primary and superharmonic resonances in the cutting process, namely, (I)Ω≈ω1 and (II) 2Ω≈ω1, are investigated.(I) Primary resonance condition (Ω≈ω1): For this case, the expression Ω=ω1+ε2σ is used, where σ is a detuning parameter for controlling the nearness of Ω to ω1.By substitutingΩ=ω1+ε2σ into equation (42), the following equations can be derived:(45)A1−A2ρI∂2V3∂T02−2A2ρIΩ∂W3∂T0+A3V3=P¯¯3eiω1T0+Q¯¯3eiω2T0,A1−A2ρI∂2W3∂T02+2A2ρIΩ∂V3∂T0+A3W3=P¯¯4eiω1T0+Q¯¯4eiω2T0,where(46)P¯¯3=P¯3+q¯1eiσT2,P¯¯4=P¯4−iq¯1eiσT2,Q¯¯3=Q¯3,Q¯¯4=Q¯4.Here, q¯1=φ11−iη¯2+ς¯2/2.The particular solutions of equation (45) are(47)V3T0,T2=F11T2eiω1T0+F12T2eiω2T0,W3T0,T2=F21T2eiω1T0+F22T2eiω2T0.Substituting equation (47) into (45) and equating the coefficient of eiω1T0 and eiω2T0 in both sides of equation (45), one has(48)A3−A1−A2Iω12F11−i2A2ρIΩω1F21=P¯¯3,i2A2ρIΩω1F11+A3−A1−A2Iω12F21=P¯¯4,(49)A3−A1−A2Iω22F12−i2A2ρIΩω2F22=Q¯¯3,i2A2ρIΩω2F12+A3−A1−A2Iω22F22=Q¯¯4.Equations (48) and (49) 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 [33](50)P¯¯3−i2A2ρIΩω1P¯¯4A3−A1−A2Iω12=0,A3−A1−A2Iω12P¯¯3i2A2ρIΩω1P¯¯4=0,Q¯¯3−i2A2ρIΩω2Q¯¯4A3−A1−A2Iω22=0,A3−A1−A2Iω22Q¯¯3i2A2IΩω2Q¯¯4=0.After simplification, the above four solvability conditions are reduced into two independent equations governingF1 and F2 in the following form:(51)iΛ1∂F1∂T2+i2A1ciΩF1−iμv+μwω1F1+ΓF1F2F¯2−8A¯5F12F¯1+2q¯1eiσT2=0,(52)iΛ2∂F2∂T2−i2A1ciΩF2−iμv+μwω2F2+ΓF1F2F¯1−8A¯5F22F¯2=0.ExpressF1 and F2 in the polar form; that is, Fj=ajT2eiθjT2/2, where ajT2 and θjT2j=1,2 are amplitudes and phase angles of the response, respectively.SubstitutingF1 and F2 into equations (51) and (52) and separating the real and imaginary parts, the following expression can be obtained:(53)Λ1ψ1′a12−σa12+Γa1a228−A5a13+φ11ς¯2cosψ1+η¯2sinψ1=0,Λ1a1′2−μv+μwω1a12+A1ciΩa1+φ11ς¯2sinψ1−η¯2cosψ1=0,Λ2a2ψ2′−Γa12a24−2A5a23=0,Λ2a2′2−μv+μwω2a22−A1ciΩa2=0,where σT2−θj=ψjj=1,2.To determine the steady-state forced response, the time derivatives in equation (53) are equated to zero. It can be immediately concluded that only solution for a2T2 is zero. This shows that, under the primary resonance Ω≈ω1, only the first mode, that is, the forward mode, can be excited, while the second mode, that is, the backward mode, does not participate in the primary resonance and remains stationary.Substitutinga2=0 into equation (51) and solving σ,(54)σ=−2A5a12Λ1±2φ121ς¯22+η¯22/a12−cd2Λ1,where(55)cd=A1ceω1+A1ciω1−Ω+φ121Nα¯0−β¯0′ω1τ4π.According to equation (55), the total damping is arisen by the external damping ce, the internal damping ci, and the induced terms due to the regenerative chatter mechanism.Damping coefficientsce and ci are determined by using concepts from a single-degree-of-freedom vibration system; namely,(56)ce+ci=2ςω1Ω=0A1−A2ρIA1,where ς is the damping ratio. Given a small value of ς, ce+ci can be obtained from equation (56), and the reasonableness value of ce+ci can be checked by inspecting the decay rate in a numerically obtained solution to an initial value problem.To study the stability of the steady-state response, consider the following equation:(57)Λ1ψ1′a12−σa12−A5a13+φ11ς¯2cosψ1+η¯2sinψ1=0,Λ1a1′2−μv+μwω1a12+A1ciΩa1+φ11ς¯2sinψ1−η¯2cosψ1=0.The nature of singular points is investigated by linearization of equation (55) around a0,ψ0 as(58)a1′ψ1′=2Λ1cd−σ2+A5a02a0σ2a0+3A5a0cda1ψ1.The stability of the steady-state response depends on the eigenvalue of the coefficient matrix (Jacobian matrix) on the right-hand side of equation (58):(59)λ1,2=cd±−σ2+3A5a02σ2+A5a02.To have a stable steady-state response, the real part of eigenvalues must be negative. On the other hand, the stable steady-state response will be unstable when the following condition is met:(60)cd2+σ2+3A5a02σ2+A5a02<0.(II) Superharmonic resonance condition (2Ω≈ω1): For this case, the formulation of steady-state response is the same as case I, but in equation (54), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


## 2.1. Kinetic Energy and Potential Energy

Figure1 displays a rotating composite cutter bar with a length of L. The cutter bar is fixed in one end while another end is free. The rectangular coordinate system (X, Y, Z) denotes an inertia coordinate system, while (x, y, z) denotes a local coordinate system. The inertial coordinate axis is consistent with the principal axis of the cross section of the cutter bar. The coordinate origin is located at x on the centerline of the deformed cutter bar. In addition, ux,t,vx,t, and wx,t denote displacements of point x on the deformed cutter bar along X, Y, and Z directions, respectively, and ϕx,t denotes the torsional deformation of cross section around x-axis.Figure 1
Schematic of the structure for the rotating composite cutter bar.The kinetic energy of a composite cutter bar can be written as [26](1)T=12∫0Lmu˙2+v˙2+w˙2+2ρIω12+ρIω22+ω32dx,where “.” represents partial derivative with respect to time t. m and ρI denote mass and cross-section inertia moment per unit length, respectively, which can be written as(2)m=π∑k=1Nρkrk+12−rk2,ρI=π4∑k=1Nρkrk+14−rk4,where N denotes ply numbers, ρk denotes the density of the k–th layer, and rk and rk+1 denote both inner and external diameters of the k-layer.The rotating angular velocities of the coordinate system (x, y, z) with respect to (X, Y, Z), ω1, ω2, and ω3 can be written as(3)ω1=φ˙−ψ˙ysinψz,ω2=ψ˙ysinφcosψz+ψ˙zcosφ,ω3=ψ˙ycosφcosψz−ψ˙zsinφ,where ψ=φ+Ωt. ψy and ψz denote the rotation angles of cross section around z-axis and y-axis, respectively, and ψz and ψy can be expressed as [26](4)ψz=sin−1v′1+u′2+v′2,ψy=−sin−1w′1+u′2+v′2+w′2,where “′” represents partial derivative with respect to time x.The elastic potential energy of a composite Euler–Bernoulli cutter bar can be written as(5)U=12∫Vσxεx+τxαγxαdV,where dV=rdrdαdx denotes a differential volume element in the cylindrical coordinate system and α and r are polar angle and polar diameter, respectively.By neglecting shear deformation, the stress-strain relations of thek–th layer expressed in the cylindrical coordinate system can be expressed as [27](6)σx=Q¯11εx+Q¯16γxα,(7)τxα=Q¯16εx+Q¯66γxα,where σx and τxα denote normal stress and shear stress of the point in the cylindrical coordinate x,r,α and Q¯iji,j=1,6 denotes the off-axis stiffness coefficient of a single layer of the composite.The strain-displacement equation can be written as(8)εx=ε+ρ2rsinα−ρ3rcosα,(9)γxα=rρ1,where ρii=1,2,3 denotes the curvature of the cutter bar. By virtue of Kirchhoffʼs kinetic analog [28], we can obtain the cutter barʼs curvatures ρ1, ρ2, and ρ3 by replacing the time derivatives ∂/∂t with the spatial derivatives ∂/∂x in the angular velocity expressions. Therefore, the curvatures are given by [29](10)ρ1=φ′−ψy′sinψz,ρ2=ψy′sinφcosψz+ψz′cosφ,ρ3=ψy′cosφcosψz−ψz′sinφ.The variation of the elastic potential energy can be written as(11)U=12∫Vσxεx+τxαγxαdV.By substituting equations (6)–(9) into equation (11), the following expression can be derived:(12)U=12∫0LA11ε2+2B16ρ1ε+D66ρ12+D66ρ22+ρ32dx,where(13)D11=π4∑k=1NQ¯11rk+14−rk4,D66=π2∑k=1NQ¯66rk+14−rk4,B16=2π3∑k=1NQ¯16rk+13−rk3,A11=π∑k=1NQ¯11rk+12−rk2.The strain along the elastic axis of a differential element dx is defined by(14)ε=1+u′2+v′2+w′2−1.Assuming that the composite cutter bar is inextensible in axial direction, that is, the strainε equals 0, the following expression can be derived:(15)u′=1−v′2−w′2−1≅−12v′2+w′2+⋅⋅⋅.Assuming that lateral displacementsv and w are of order Oe, e≪1, then u=Oe2. Substituting equation (4) into equations (3) and (10), expanding the outcomes in Taylor series, and retaining terms up to Oe3, one can obtain the angular velocities and curvatures up to Oe3. Substituting the resulted curvatures and angular velocities into equations (1) and (12) and using equation (15), the final expressions for kinetic and strain energy are obtained.


## 2.2. Virtual Work of the Cutting Force

The virtual work of the cutting force can be written as(16)δW=∫0LLvδv+Lwδwdx,where Lv=FyδDx−L, Lw=FzδDx−L, and δD denotes Delta function.The regenerative cutting force,Fy and Fz, can be expressed as [15](17)Fy=−Nα0v−vτ+β0w−wτ+γ02π+cfς1cos2Ωt−π/2−η1cos2Ωt+η12+η2cosΩt−π2+ς2cosΩt,Fz=Nα0′v−vτ+β0′w−wτ+γ0′2π+cfη1cos2Ωt−π/2+ς1cos2Ωt−ς12+ς2cosΩt−π2−η2cosΩt,where(18)α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−τ.Here, τ denotes delay time τ=2π/NΩ, N denotes the number of cutting teeth, (Ktc, Kte) and (Krc, Kre) are the cutting force coefficients along tangential and radial directions, respectively, a denotes the axial cutting depth, cf denotes the feed per tooth per revolution and corresponds to the static part of the chip thickness cfsinϕj, and vjt−τ−vjt denotes the dynamic chip thickness produced due to vibrations of the cutter bar at the present vjt and previous vjt−τ tooth periods (see Figure 2).Figure 2
Sketches of the milling process. (a) Cross-section view of the milling process. (b) Dynamic chip thickness.
(a)(b)


## 2.3. Equation of Motion in Milling Process

The Lagrangian is defined asL=T−U. Using equations (1) and (12) and introducing the Lagrange multiplier λ to enforce the inextensionality conditions [30], one has(19)L=12∫0Lmu˙2+v˙2+w˙2+2ρIω12+ρIω22+ω32−A11ε2+2B16ρ1ε+D66ρ12+D66ρ22+ρ32+λ1−1+u′2−v′2−w′2dx.To derive the governing equations of motion, we use Hamilton’s extended principle as(20)∫t1t2δL+δWdt=0.The bending-bending coupled nonlinear motion equation of a rotating composite cutter bar can be obtained as(21)mv¨+v′h′+ce+civ˙−ρI2Ωw˙″+v¨″+D11g1+ciΩw=Fyδx−L,(22)mw¨+w′h′+ce+ciw˙+ρI2Ωv˙″−w¨″+D11g2−ciΩv=Fzδx−L,where(23)g1=4v″v′v‴+w4w′v′+w′v″w‴+3w″v′w‴+w″2v″+v″3+v4v′2+v4,g2=4w″w′w‴+v4v′w′+v′w″v‴+3v″w′v‴+v″2w″+w″3+w4w′2+w4,h=∫Lx∫0xv˙′2+v¨′v′+w˙′2+w¨′w′dxdx.It should be noted that the torsional fundamental frequency of the cutting bar with circular cross section is much greater than the bending frequency; thus, the torsional inertia item can be neglected [29].Using the above assumption, the following expression can be obtained:(24)ϕ=−∫0xv″w′dx+⋯.Since the cutting bar is slender with small rotation inertia, the nonlinear item multiplied with the rotational inertia can also be neglected [26]. Accordingly, bending-bending-torsion coupled nonlinear motion equation can be appropriately simplified into equations (21) and (22).In equation (21), the terms g1 and g2 show the effect of the geometric nonlinearities, and the term h shows the effect of the inertia nonlinearities. In fact, h is given by h=−λ=m/2∫Lx∂2/∂t2∫0xv′2+w′2dxdx. One can interpret λ as an axial force acting at the cutter barʼs tip to prevent it from stretching.The termsce and ci are external and internal viscous damping coefficients, respectively. Damping is an important parameter which influences the dynamic behavior of the cutting process [23]. Two commonly used models in rotating system dynamics are external damping (nonrotating damping) and internal damping (rotating damping). The external damping is generally treated like viscous damping, but internal damping can be described by the viscous damping model or the hysteretic damping model [31, 32]. The damping mechanism is complex, particularly for hysteretic damping. Inclusion of hysteretic damping could be achieved using the complex stiffness; however, this leads to considerably increased complexity of the present dynamic model. It can be noted that hysteretic damping has a minor effect compared to viscous damping for an isotropic cantilevered shaft [31]. In light of above, the classical viscous damping model is employed in our work to describe the external and internal damping. ce is proportional to the speed in the inertial coordinate system and ci is proportional to the speed in the rotating coordinate system [23, 31, 32].Using the dimensionless transformation(25)x¯=xL,v¯=vL,w¯=wL,t¯=D11mL4t=γt,τ¯=γτ,c¯e=ceL2mD11,c¯i=ciL2mD11,ρ¯I¯=ρImL2,Ω¯=mL4D11Ω,nonlinear equation of motion in dimensionless form can be written as(26)v¨+ce+civ˙−2ΩρIw˙″−ρIv¨″+v′h′+g1+ciΩw=F¯yδx−1,(27)w¨+ce+ciw˙+2ΩρIv˙″−ρIw¨″+w′h′+g2−ciΩv=F¯zδx−1,where(28)F¯y=−Nα¯0v−vτ+β¯0w−wτ+γ¯02π+cfς1cos2Ωt−π/2−cfη1cos2Ωt+cfη12+η¯2cosΩt−π2+ς¯2cosΩt,F¯z=Nα¯0′v−vτ+β¯0′w−wτ+γ¯0′2π−cfη1cos2Ωt−π/2+cfς1cos2Ωt−cfς12+ς¯2cosΩt−π2−η¯2cosΩt.The coefficients in equation (26) can be defined as(29)α¯0=L3α0D11,β¯0=L3β0D11,γ¯0=L2γ0D11,cfς1=L2cfς1D11,cfη1=L2cfη1D11,α¯0′=L3α0′D11,β¯0′=L3β0′D11,γ¯0′=L2γ0′D11,ς¯2=L2ς2D11,η¯2=L2η2D11.For simplicity, the symbol “–” above all dimensionless variables in equations (26) and (27) has been removed.Before applying the multiple scales method, the partial differential equations of motion are discretized by using the Galerkin method.Let(30)vx,t=ϕ1xVt,wx,t=ϕ1xWt.For cantilever beam, the mode shapeϕ1x can be written as(31)φ1x=cosβ1Lx−coshβ1Lx−cosβ1L+coshβ1Lsinβ1L+sinhβ1Lsinβ1Lx−sinhβ1Lx,x∈0,1,where β1L=1.875.By substituting equation (31) into equations (26) and (27), taking the inner product of each equation with its corresponding mode shape, the following ordinary differential equations can be obtained:(32)A1V¨+A1ce+ciV˙−A2ρI2ΩW˙+V¨+A3V+A1ciΩW+A4VV˙2+VV¨+W˙2+WW¨+A5VV2+W2=∫01F¯yδx−1dx,A1W¨+A1ce+ciW˙+A2ρI2ΩV˙−W¨+A3W−A1ciΩV+A4WV˙2+VV¨+W˙2+WW¨+A5WV2+W2=∫01F¯zδx−1dx,in which(33)A1=∫01φ12dx,A2=∫01φ1φ1″dx,A3=∫01φ1φ14dx,A4=∫01φ1φ1′F1′dx,A5=∫01φ1F2dx,F1=∫1x∫0xφ1′2dxdx′,F2=4φ1′φ1″φ1‴+φ1″3+φ14φ1′2.To simplify the delayed expressions in the cutting force (28), the first-order Pade approximation as e−sτ≅1−sτ is used for delayed terms (i.e., a truncated Taylor series expansion for e−sτ with one delayed term is used) [15].Substituting the cutting forces in equation (28) into equations (26) and (27), the following expression can be derived:(34)A1−A2ρIV¨+μVV˙−2ΩρIA2W˙+A3V+A1ciΩW+A4VV˙2+VV¨+W˙2+WW¨+A5VV2+W2=−λV+F¯y2φ11,A1−A2ρIW¨+μWW˙+2ΩρIA2V˙+A3W−A1ciΩV+A4WV˙2+VV¨+W˙2+WW¨+A5WV2+W2=λW+F¯z2φ11,where F¯y2 and F¯z2 denote the sum of the second and third items in the expressions of cutting forces (28) and the parameters in equation (34) are defined as(35)μV=A1ce+ci+Nα¯0τφ1212π,μW=A1ce+ci−Nβ¯0′τφ1212π,λV=Nγ¯0φ112π,λW=Nγ¯0′φ112π.According to equation (35), the terms arisen from regenerative chatter mechanism (Nα¯0τφ121/2π and Nβ¯0τφ121/2π) are proportional to delay time τ and consequently inversely proportional to the spindle speed Ω. In this paper, a truncated Taylor series expansion for e−sτ is used to simplify the regenerative chatter mechanism. Therefore, it must be noticed that the results of this research are valid for small-time delays which occur essentially at high-speed machining. Here, for the considered case study, it can be shown that using the first-order truncated Taylor series for e−sτ and modifying only the damping coefficients are sufficient (as done in the works of Moradi et al. [15]).


## 2.4. Method of Multiple Scales

The multiple scales method is used to find an approximate solution for the proposed nonlinear dynamics of the milling process. In order to apply multiple scale method,V and W can be written in the expansion form as(36)Vt=εV1T0,T2+ε3V3T0,T2+⋯,Wt=εW1T0,T2+ε3W3T0,T2+⋯,where ε is a small dimensionless parameter. T0=t and T2=ε2t are fast and slow time scales, respectively. The damping and exciting force terms should be scaled, so that their effects are balanced with nonlinearities. Thus, μV, μW, ci, λV, λW, F¯y2, and F¯z2 are replaced with ε2μV, ε2μW, ε2ci, ε3λV, ε3λW, ε3F¯y2, and ε3F¯z2.Time derivatives become expansions in terms of the partial derivatives with respect toT0 and T2 by using the following chain rule:(37)ddt=∂∂T0+ε2∂∂T2+⋯;d2dt2=∂2∂T02+2ε2∂∂T2∂∂T0+⋯.By substituting equations (36) and (37) into equation (34) and equating the coefficients of ε and ε3, one hasOε(38)A1−A2ρI∂2V1∂T02−2A2ρIΩ∂W1∂T0+A3V1=0,A1−A2ρI∂2W1∂T02+2A2ρIΩ∂V1∂T0+A3W1=0.Oε3(39)A1−A2ρI∂2V3∂T02−2A2ρIΩ∂W3∂T0+A3V3=−μV∂V1∂T0−A1ciΩW1−2A1−A2ρI∂2V1∂T0∂T2+2A2ρIΩ∂W1∂T2−A4V1∂V1∂T02+V1∂2V1∂T02+∂W1∂T02+W1∂2W1∂T02−A5V1V12+W12−λV+F¯y2φ11,A1−A2ρI∂2W3∂T02+2A2ρIΩ∂V3∂T0+A3W3=−μW∂W1∂T0+A1ciΩV1−2A1−A2ρI∂2W1∂T0∂T2−2A2ρIΩ∂V1∂T2−A4W1∂V1∂T02+V1∂2V1∂T02+∂W1∂T02+W1∂2W1∂T02−A5W1V12+W12+λW+F¯z2φ11.The solution to equation (38) can be written as(40)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,where i=−1 denotes the imaginary unit; F1T2 and F2T2 are the complex-valued functions to be determined, respectively; F¯1T2 and F¯2T2 are complex conjugates; ω1 and ω2 denote forward and backward linear natural frequencies, respectively, and can be written as(41)ω1=A2ρIΩ+A2ρIΩ2+A1−A2ρIA3A1−A2ρI,ω2=−A2ρIΩ+A2ρIΩ2+A1−A2ρIA3A1−A2ρI.By substituting (40) into equation (39), one has(42)A1−A2ρI∂2V3∂T02−2A2ρIΩ∂W3∂T0+A3V3=P¯3eiω1T0+Q¯3eiω2T0−λV+F¯y2φ11+cc,A1−A2ρI∂2W3∂T02+2A2ρIΩ∂V3∂T0+A3W3=P¯4eiω1T0+Q¯4eiω2T0+λW+F¯z2φ11+cc,where cc denotes the complex conjugate.The coefficients of the harmonic items on the left end of equation (42) have the following form:(43)P¯3=ΓF1F¯2F22−4A5F¯1F12+iΛ1F1′2−iμVω1F1+iA1ciΩF1,Q¯3=ΓF1F¯1F22−4A5F¯2F22+iΛ2F2′2−iμVω2F2−iA1ciΩF2,P¯4=−iΓF1F¯2F22+i4A5F¯1F12+Λ1F1′2−μWω1F1+A1ciΩF1,Q¯4=iΓF1F¯1F22−i4A5F¯2F22−Λ2F2′2+μWω2F2+A1ciΩF2,where(44)Γ=A4ω12+ω222−16A5,Λ1=4−A1−A2ρIω1+A2ρIΩ,Λ2=−4A1−A2ρIω2+A2ρIΩ.In this paper, both types of primary and superharmonic resonances in the cutting process, namely, (I)Ω≈ω1 and (II) 2Ω≈ω1, are investigated.(I) Primary resonance condition (Ω≈ω1): For this case, the expression Ω=ω1+ε2σ is used, where σ is a detuning parameter for controlling the nearness of Ω to ω1.By substitutingΩ=ω1+ε2σ into equation (42), the following equations can be derived:(45)A1−A2ρI∂2V3∂T02−2A2ρIΩ∂W3∂T0+A3V3=P¯¯3eiω1T0+Q¯¯3eiω2T0,A1−A2ρI∂2W3∂T02+2A2ρIΩ∂V3∂T0+A3W3=P¯¯4eiω1T0+Q¯¯4eiω2T0,where(46)P¯¯3=P¯3+q¯1eiσT2,P¯¯4=P¯4−iq¯1eiσT2,Q¯¯3=Q¯3,Q¯¯4=Q¯4.Here, q¯1=φ11−iη¯2+ς¯2/2.The particular solutions of equation (45) are(47)V3T0,T2=F11T2eiω1T0+F12T2eiω2T0,W3T0,T2=F21T2eiω1T0+F22T2eiω2T0.Substituting equation (47) into (45) and equating the coefficient of eiω1T0 and eiω2T0 in both sides of equation (45), one has(48)A3−A1−A2Iω12F11−i2A2ρIΩω1F21=P¯¯3,i2A2ρIΩω1F11+A3−A1−A2Iω12F21=P¯¯4,(49)A3−A1−A2Iω22F12−i2A2ρIΩω2F22=Q¯¯3,i2A2ρIΩω2F12+A3−A1−A2Iω22F22=Q¯¯4.Equations (48) and (49) 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 [33](50)P¯¯3−i2A2ρIΩω1P¯¯4A3−A1−A2Iω12=0,A3−A1−A2Iω12P¯¯3i2A2ρIΩω1P¯¯4=0,Q¯¯3−i2A2ρIΩω2Q¯¯4A3−A1−A2Iω22=0,A3−A1−A2Iω22Q¯¯3i2A2IΩω2Q¯¯4=0.After simplification, the above four solvability conditions are reduced into two independent equations governingF1 and F2 in the following form:(51)iΛ1∂F1∂T2+i2A1ciΩF1−iμv+μwω1F1+ΓF1F2F¯2−8A¯5F12F¯1+2q¯1eiσT2=0,(52)iΛ2∂F2∂T2−i2A1ciΩF2−iμv+μwω2F2+ΓF1F2F¯1−8A¯5F22F¯2=0.ExpressF1 and F2 in the polar form; that is, Fj=ajT2eiθjT2/2, where ajT2 and θjT2j=1,2 are amplitudes and phase angles of the response, respectively.SubstitutingF1 and F2 into equations (51) and (52) and separating the real and imaginary parts, the following expression can be obtained:(53)Λ1ψ1′a12−σa12+Γa1a228−A5a13+φ11ς¯2cosψ1+η¯2sinψ1=0,Λ1a1′2−μv+μwω1a12+A1ciΩa1+φ11ς¯2sinψ1−η¯2cosψ1=0,Λ2a2ψ2′−Γa12a24−2A5a23=0,Λ2a2′2−μv+μwω2a22−A1ciΩa2=0,where σT2−θj=ψjj=1,2.To determine the steady-state forced response, the time derivatives in equation (53) are equated to zero. It can be immediately concluded that only solution for a2T2 is zero. This shows that, under the primary resonance Ω≈ω1, only the first mode, that is, the forward mode, can be excited, while the second mode, that is, the backward mode, does not participate in the primary resonance and remains stationary.Substitutinga2=0 into equation (51) and solving σ,(54)σ=−2A5a12Λ1±2φ121ς¯22+η¯22/a12−cd2Λ1,where(55)cd=A1ceω1+A1ciω1−Ω+φ121Nα¯0−β¯0′ω1τ4π.According to equation (55), the total damping is arisen by the external damping ce, the internal damping ci, and the induced terms due to the regenerative chatter mechanism.Damping coefficientsce and ci are determined by using concepts from a single-degree-of-freedom vibration system; namely,(56)ce+ci=2ςω1Ω=0A1−A2ρIA1,where ς is the damping ratio. Given a small value of ς, ce+ci can be obtained from equation (56), and the reasonableness value of ce+ci can be checked by inspecting the decay rate in a numerically obtained solution to an initial value problem.To study the stability of the steady-state response, consider the following equation:(57)Λ1ψ1′a12−σa12−A5a13+φ11ς¯2cosψ1+η¯2sinψ1=0,Λ1a1′2−μv+μwω1a12+A1ciΩa1+φ11ς¯2sinψ1−η¯2cosψ1=0.The nature of singular points is investigated by linearization of equation (55) around a0,ψ0 as(58)a1′ψ1′=2Λ1cd−σ2+A5a02a0σ2a0+3A5a0cda1ψ1.The stability of the steady-state response depends on the eigenvalue of the coefficient matrix (Jacobian matrix) on the right-hand side of equation (58):(59)λ1,2=cd±−σ2+3A5a02σ2+A5a02.To have a stable steady-state response, the real part of eigenvalues must be negative. On the other hand, the stable steady-state response will be unstable when the following condition is met:(60)cd2+σ2+3A5a02σ2+A5a02<0.(II) Superharmonic resonance condition (2Ω≈ω1): For this case, the formulation of steady-state response is the same as case I, but in equation (54), ς¯2 and η¯2 are replaced with cfς1/2 and cfη1/2, respectively.


## 3. Numerical Results and Discussion

In this study, the composite material of carbon fiber/epoxy resin is selected as the material of the cutter bar. The mechanical properties of the material are shown in Table1. Coefficients of cutting forces for simulation are given in Table 2. The cutter bar has a hollow structure, the outer diameter of the cross section is D = 8 mm, the inner diameter is D = 40 mm, the thickness of the section is h = 4 mm, and length L is determined by the given ratio of length to diameter. The composite cutter bar has 16 layers with identical thickness, and the stack sequence is ±θ8.Table 1
Mechanical properties of carbon fiber/epoxy composite [34].
ρ (kg/m3)E11 (GPa)E22 (GPa)G12 (GPa)G23 (GPa)υ12167225.88.73.53.50.34Table 2
Coefficients of cutting forces [15].
Cf (mm/rev-tooth)ξ1 (N/mm)η1 (N)ξ2 (N/mm)η2 (N)ac (mm)N0.2620432085234Figure3 shows the natural frequency versus rotating speed when the gyroscopic effect of the cutter bar is considered (equation (41)). 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 cutter bar whirls in direction of the rotation. But, in backward natural frequency, the natural frequency is measured when the cutter bar whirls in the opposite direction of the rotation. The forward natural frequency (black solid line) increases with the increase of the rotational speed, while the backward natural frequency (blue dashed line) decreases with the increase of the rotational speed.Figure 3
Natural frequency versus rotational speed (L/d = 10; θ = 0°).


### 3.1. Primary ResonanceΩ¯≈ω1

Figure4 shows the effect of ratio of length to diameter on the frequency response of the vibration at the cutter tip. The increase of ratio of length to diameter results in less amplitude vibration. This is mathematically because, according to equation (29) and consequently equation (28), by increasing ratio of length to diameter, total damping increases due to the damping term arisen from the regenerative chatter mechanism. Figure 5 shows the effect of the damping ratio on the frequency response of the vibration at the cutter tip. As it is shown and physically expected, when the damping ratio increases, vibration amplitudes of the cutter tip decrease. In order to study the effect of cutting force coefficients, we let ςi=Kfς¯¯i and ηi=Kfη¯¯ii=1,2, with nominal parameters ς¯¯i and η¯¯i given in Table 2. Figure 6 shows the effect of the cutting force coefficient on the frequency response of the vibration at the cutter tip. As the cutting force coefficient increases, the vibration amplitudes of the cutter tip increase. Figure 7 shows the effect of ply angle on the frequency response of the vibration at the cutter tip. The vibration amplitudes of the cutter tip increase as the ply angle increase. This is physically expected. As listed in Table 1, the longitudinal elastic modulus along fiber E11 is significantly greater than the lateral elastic modulus E22. Therefore, at a larger ply angle, the flexural rigidity of the composite cutter bar decreases, thereby causing a larger vibration response at the cutter tip. Figure 8 shows the effect of the values of β=ci/ci+ce on the frequency response of the vibration at the cutter tip. The cases of β = 0, 0.5, and 1 refer no internal damping, internal damping equal to external damping, and no external damping, respectively. As it is shown, vibration amplitudes of the cutter tip are proportional to the internal damping and inversely proportional to the external damping.Figure 4
The effect of ratios of length to diameter on the frequency response of the cutter tip vibration for the primary resonance (ς = 0.01, Kf = 1, and θ = 00).Figure 5
The effect of damping ratio on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, Kf = 1, and θ = 00).Figure 6
The effect of cutting force coefficient on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, and θ = 00).Figure 7
The effect of ply angle on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, and Kf = 1).Figure 8
The effect of external and internal damping on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, Kf = 1, and θ = 00).The effect of rotating speed on the frequency response of the vibration at the cutter tip is shown in Figure9. As it is observed, decreasing rotating speed values leads to a decrease in vibration amplitudes. This is physically expected, because total damping increases with the decrease of rotating speed due to the damping effect from the regenerative chatter mechanism.Figure 9
The effect of rotating speed on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, Kf = 1, and θ = 00).It is clear from Figures4–9 that the curves bend to the right side because of the existence of nonlinear factors. Therefore, the effective nonlinearity for the cutting system is of hardening type. For a specified value of σ, there are three values for amplitude. Consequently, the phenomena of jumping and bifurcation occur in the system which means that the cutting system works in an unstable cutting condition.Figures10–15 show the amplitude versus damping ratio with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, the values of β, and rotating speeds, respectively. From these figures, it can be seen that, for some values of L/d, σ, Kf, θ, β, or Ω, there are multivalued curves. For example, as shown in Figure 10, when L/d = 10 and ς < 0.018, the system has two stable and one unstable branches, but for L/d = 10 and ς > 0.018, there exists only one stable branch. For large values of ς and arbitrary values of L/d, curves are always single-valued.Figure 10
Amplitude versus damping ratio with different ratios of length to diameter (σ = 0.25, Kf = 1, and θ = 00).Figure 11
Amplitude versus damping ratio with different detuning parameter values (L/d = 10, Kf = 1, and θ = 00).Figure 12
Amplitude versus damping ratio with different cutting force coefficients (L/d = 10, σ = 0.25, and θ = 00).Figure 13
Amplitude versus damping ratio with different ply angles (L/d = 10, Kf = 1, and σ = 0.25).Figure 14
Amplitude versus damping ratio with different values ofβ (L/d = 10, Kf = 1, σ = 0.25, and θ = 00).Figure 15
Amplitude versus damping ratio with different rotating speeds (L/d = 10, Kf = 1, σ = 0.25, and θ = 00).Figures16–21 show the amplitude versus ply angle with different ratios of length to diameter, detuning parameter values, cutting force coefficients, damping ratios, the values of β, and rotating speeds, respectively. Similar to the results in Figures 10–15, for some values of L/d, σ, Kf, ς, β, or Ω, multivalued curves appeared.Figure 16
Amplitude versus ply angle with different ratios of length to diameter (ς = 0.01, Kf = 1, and σ = 0.25).Figure 17
Amplitude versus ply angle with different detuning parameter values (L/d = 10, Kf = 1, and ς = 0.01).Figure 18
Amplitude versus ply angle with different cutting force coefficients (L/d = 10, σ = 0.25, and ς = 0.01).Figure 19
Amplitude versus ply angle with different damping ratios (L/d = 10, σ = 0.25, and Kf = 1).Figure 20
Amplitude versus ply angle with different values ofβ (L/d = 10, σ = 0.25, Kf = 1, and ς = 0.01).Figure 21
Amplitude versus ply angle with different rotating speeds (L/d = 10, σ = 0.25, Kf = 1, and ς = 0.01).Figures22–27 show the amplitude versus cutting force coefficient with different ratios of length to diameter, detuning parameter values, ply angle, damping ratios, the values of β, and rotating speeds, respectively. Again, for some values of L/d, σ, θ, ς, β, or Ω, multivalued curves can be seen.Figure 22
Amplitude versus cutting force coefficient with different ratios of length to diameter (θ = 00, ς = 0.01, and σ = −0.25).Figure 23
Amplitude versus cutting force coefficient with different detuning parameter values (θ = 00, ς = 0.01, and L/d = 10).Figure 24
Amplitude versus cutting force coefficient with different ply angles (σ = 0.25, ς = 0.01, and L/d = 10).Figure 25
Amplitude versus cutting force coefficient with different damping ratios (L/d = 10, θ = 00, σ = 0.25, and Kf = 1).Figure 26
Amplitude versus cutting force coefficient with different values ofβ (L/d = 10, θ = 00, σ = 0.25, and ς = 0.01).Figure 27
Amplitude versus cutting force coefficient with different rotating speeds (L/d = 10, θ = 00, σ = 0.25, and ς = 0.01).


### 3.2. Superharmonic Resonance2Ω¯≈ω1

For the superharmonic resonance case, similar to case I, botha1-σ curves and the amplitude versus different parameter curves show jump phenomenon and consequently multiple solutions at a given rotating speed. The effective nonlinearity of the hardening type can also be observed. The stability of the fixed points can also be determined by investigating the eigenvalues of the Jacobian matrix of the system. Moreover, comparing primary resonance and superharmonic resonance cases indicates that under superharmonic resonance condition vibration amplitudes are generally lower than those under primary resonance. Accordingly, the calculated results in case II are not shown here for the sake of space.To verify the perturbation results, a numerical simulation has been utilized in Figures28 and 29. The results achieved from multiple scales method show a good agreement with those of numerical simulation. In Figures 28 and 29, it is seen that if A5 (the cubic nonlinearity) is removed from equation (54), the nonlinear solution reduces the (single-valued) linear solution. In addition, an interesting result is observed that vibration amplitudes of the linear system are less than those of nonlinear ones. This is because internal damping provides extra negative damping to the peak amplitudes of the nonlinear system, while internal damping has no effect on the peak amplitudes of the linear system when σ = 0. If internal damping is absent, vibration amplitudes of the linear system are larger than those of nonlinear ones, as can be seen in Figures 30 and 31.Figure 28
Comparison of results between the method of multiple scales and the numerical simulation (case I with internal damping).Figure 29
Comparison of results between the method of multiple scales and the numerical simulation (case II with internal damping).Figure 30
Comparison of results between the method of multiple scales and the numerical simulation (case I without internal damping).Figure 31
Comparison of results between the method of multiple scales and the numerical simulation (case II without internal damping).


## 3.1. Primary ResonanceΩ¯≈ω1

Figure4 shows the effect of ratio of length to diameter on the frequency response of the vibration at the cutter tip. The increase of ratio of length to diameter results in less amplitude vibration. This is mathematically because, according to equation (29) and consequently equation (28), by increasing ratio of length to diameter, total damping increases due to the damping term arisen from the regenerative chatter mechanism. Figure 5 shows the effect of the damping ratio on the frequency response of the vibration at the cutter tip. As it is shown and physically expected, when the damping ratio increases, vibration amplitudes of the cutter tip decrease. In order to study the effect of cutting force coefficients, we let ςi=Kfς¯¯i and ηi=Kfη¯¯ii=1,2, with nominal parameters ς¯¯i and η¯¯i given in Table 2. Figure 6 shows the effect of the cutting force coefficient on the frequency response of the vibration at the cutter tip. As the cutting force coefficient increases, the vibration amplitudes of the cutter tip increase. Figure 7 shows the effect of ply angle on the frequency response of the vibration at the cutter tip. The vibration amplitudes of the cutter tip increase as the ply angle increase. This is physically expected. As listed in Table 1, the longitudinal elastic modulus along fiber E11 is significantly greater than the lateral elastic modulus E22. Therefore, at a larger ply angle, the flexural rigidity of the composite cutter bar decreases, thereby causing a larger vibration response at the cutter tip. Figure 8 shows the effect of the values of β=ci/ci+ce on the frequency response of the vibration at the cutter tip. The cases of β = 0, 0.5, and 1 refer no internal damping, internal damping equal to external damping, and no external damping, respectively. As it is shown, vibration amplitudes of the cutter tip are proportional to the internal damping and inversely proportional to the external damping.Figure 4
The effect of ratios of length to diameter on the frequency response of the cutter tip vibration for the primary resonance (ς = 0.01, Kf = 1, and θ = 00).Figure 5
The effect of damping ratio on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, Kf = 1, and θ = 00).Figure 6
The effect of cutting force coefficient on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, and θ = 00).Figure 7
The effect of ply angle on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, and Kf = 1).Figure 8
The effect of external and internal damping on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, Kf = 1, and θ = 00).The effect of rotating speed on the frequency response of the vibration at the cutter tip is shown in Figure9. As it is observed, decreasing rotating speed values leads to a decrease in vibration amplitudes. This is physically expected, because total damping increases with the decrease of rotating speed due to the damping effect from the regenerative chatter mechanism.Figure 9
The effect of rotating speed on the frequency response of the cutter tip vibration for the primary resonance (L/d = 10, ς = 0.01, Kf = 1, and θ = 00).It is clear from Figures4–9 that the curves bend to the right side because of the existence of nonlinear factors. Therefore, the effective nonlinearity for the cutting system is of hardening type. For a specified value of σ, there are three values for amplitude. Consequently, the phenomena of jumping and bifurcation occur in the system which means that the cutting system works in an unstable cutting condition.Figures10–15 show the amplitude versus damping ratio with different ratios of length to diameter, detuning parameter values, cutting force coefficients, ply angles, the values of β, and rotating speeds, respectively. From these figures, it can be seen that, for some values of L/d, σ, Kf, θ, β, or Ω, there are multivalued curves. For example, as shown in Figure 10, when L/d = 10 and ς < 0.018, the system has two stable and one unstable branches, but for L/d = 10 and ς > 0.018, there exists only one stable branch. For large values of ς and arbitrary values of L/d, curves are always single-valued.Figure 10
Amplitude versus damping ratio with different ratios of length to diameter (σ = 0.25, Kf = 1, and θ = 00).Figure 11
Amplitude versus damping ratio with different detuning parameter values (L/d = 10, Kf = 1, and θ = 00).Figure 12
Amplitude versus damping ratio with different cutting force coefficients (L/d = 10, σ = 0.25, and θ = 00).Figure 13
Amplitude versus damping ratio with different ply angles (L/d = 10, Kf = 1, and σ = 0.25).Figure 14
Amplitude versus damping ratio with different values ofβ (L/d = 10, Kf = 1, σ = 0.25, and θ = 00).Figure 15
Amplitude versus damping ratio with different rotating speeds (L/d = 10, Kf = 1, σ = 0.25, and θ = 00).Figures16–21 show the amplitude versus ply angle with different ratios of length to diameter, detuning parameter values, cutting force coefficients, damping ratios, the values of β, and rotating speeds, respectively. Similar to the results in Figures 10–15, for some values of L/d, σ, Kf, ς, β, or Ω, multivalued curves appeared.Figure 16
Amplitude versus ply angle with different ratios of length to diameter (ς = 0.01, Kf = 1, and σ = 0.25).Figure 17
Amplitude versus ply angle with different detuning parameter values (L/d = 10, Kf = 1, and ς = 0.01).Figure 18
Amplitude versus ply angle with different cutting force coefficients (L/d = 10, σ = 0.25, and ς = 0.01).Figure 19
Amplitude versus ply angle with different damping ratios (L/d = 10, σ = 0.25, and Kf = 1).Figure 20
Amplitude versus ply angle with different values ofβ (L/d = 10, σ = 0.25, Kf = 1, and ς = 0.01).Figure 21
Amplitude versus ply angle with different rotating speeds (L/d = 10, σ = 0.25, Kf = 1, and ς = 0.01).Figures22–27 show the amplitude versus cutting force coefficient with different ratios of length to diameter, detuning parameter values, ply angle, damping ratios, the values of β, and rotating speeds, respectively. Again, for some values of L/d, σ, θ, ς, β, or Ω, multivalued curves can be seen.Figure 22
Amplitude versus cutting force coefficient with different ratios of length to diameter (θ = 00, ς = 0.01, and σ = −0.25).Figure 23
Amplitude versus cutting force coefficient with different detuning parameter values (θ = 00, ς = 0.01, and L/d = 10).Figure 24
Amplitude versus cutting force coefficient with different ply angles (σ = 0.25, ς = 0.01, and L/d = 10).Figure 25
Amplitude versus cutting force coefficient with different damping ratios (L/d = 10, θ = 00, σ = 0.25, and Kf = 1).Figure 26
Amplitude versus cutting force coefficient with different values ofβ (L/d = 10, θ = 00, σ = 0.25, and ς = 0.01).Figure 27
Amplitude versus cutting force coefficient with different rotating speeds (L/d = 10, θ = 00, σ = 0.25, and ς = 0.01).


## 3.2. Superharmonic Resonance2Ω¯≈ω1

For the superharmonic resonance case, similar to case I, botha1-σ curves and the amplitude versus different parameter curves show jump phenomenon and consequently multiple solutions at a given rotating speed. The effective nonlinearity of the hardening type can also be observed. The stability of the fixed points can also be determined by investigating the eigenvalues of the Jacobian matrix of the system. Moreover, comparing primary resonance and superharmonic resonance cases indicates that under superharmonic resonance condition vibration amplitudes are generally lower than those under primary resonance. Accordingly, the calculated results in case II are not shown here for the sake of space.To verify the perturbation results, a numerical simulation has been utilized in Figures28 and 29. The results achieved from multiple scales method show a good agreement with those of numerical simulation. In Figures 28 and 29, it is seen that if A5 (the cubic nonlinearity) is removed from equation (54), the nonlinear solution reduces the (single-valued) linear solution. In addition, an interesting result is observed that vibration amplitudes of the linear system are less than those of nonlinear ones. This is because internal damping provides extra negative damping to the peak amplitudes of the nonlinear system, while internal damping has no effect on the peak amplitudes of the linear system when σ = 0. If internal damping is absent, vibration amplitudes of the linear system are larger than those of nonlinear ones, as can be seen in Figures 30 and 31.Figure 28
Comparison of results between the method of multiple scales and the numerical simulation (case I with internal damping).Figure 29
Comparison of results between the method of multiple scales and the numerical simulation (case II with internal damping).Figure 30
Comparison of results between the method of multiple scales and the numerical simulation (case I without internal damping).Figure 31
Comparison of results between the method of multiple scales and the numerical simulation (case II without internal damping).


## 4. Conclusions

This study presents the dynamic model of the cutting process considering the nonlinearity of a rotating composite cutter bar. Based on the inextensible assumption, the composite cutter bar was simplified as Euler–Bernoulli cantilever beam with nonlinear curvature and inertia. The cutter tip is subjected to regenerative cutting force. The cutting force includes harmonic exciting items. Meanwhile, the effect of internal and external damping is also considered. Nonlinear ordinary differential motion equation of the cutting system is derived by combining Hamiltonʼs extended principle and the Galerkin method. The analytical expressions of the nonlinear forced vibration response of the cutting system are derived by means of the method of multiple scales. Furthermore, this study demonstrates that, in primary resonance and superharmonic resonance of the cutting process with nonlinear cutter bar, only forward modes of the cutter bar are excited, the response amplitude of the backward mode equals 0, and the amplitudes of the forced vibration responses alongy-axis and z-axis directions only depended on a1. Due to the effect of nonlinearity, the resonance curves show hard spring vibration behavior of Duffing type oscillator; meanwhile, the jumping phenomenon and multiple solutions exist in the response curves. From the perspective of engineering application, cutting conditions and consequently initial conditions can be adjusted to achieve the stable branch with fewer vibration amplitudes. It is also found that the forced vibration resonance amplitude increases with the increase of cutting force, ply angle, internal damping, or rotating speed of composite cutter bar but decreases with the increase of damping ratio or ratio of length to diameter. External damping can reduce or limit the vibration amplitudes. The regenerative chatter mechanism (inversely proportional to rotating speed) and internal damping provide the total damping of the cutting system with positive and negative damping, respectively.



---

*Source: 2463136-2020-10-10.xml*