#### 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., , , , and , 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.

#### 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., , , , and 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 Figure 1. The kinetic energy of the rotating cutting tool without considering torsional deformation is as follows:where and represent the displacement at any point on the neutral axis along the *y* and *z* directions, respectively. and represent bending rotation angles of the cross section around the *y* and *z* axes, respectively. represents the mass per unit length, and represents the mass moment of inertia of the cross section. The superposed dots represent derivatives with respect to the time *t*.

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:where refers to the off-axis stiffness coefficient of the *k*th 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:

##### 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:where *c* is structural damping coefficient.

The virtual work of the cutting force *δW* can be expressed aswhere , , and is the delta function.

The linear regenerative cutting forces *F*_{y} and *F*_{z} can be obtained by [20]wherein which is the delay time of milling process; *N* is the number of teeth; (*K*_{tc}, *K*_{te}) and (*K*_{rc}, *K*_{re}) are the cutting force coefficients in the tangential and radial directions, respectively; is the axial cutting depth; and *c*_{f} 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:

Substituting (1), (2), (4), and (5) into (8), the equations of motion in both *y* and *z* directions can be obtained:where

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; represents the density of the *k*th layer; and *r*_{k} and *r*_{k+1} represent the inner diameter and outer diameter of the *k*th layer, respectively. The detailed expression of can be found in Appendix B. In (9), the primes denote differentiation with respect to *x*.

The following nondimensional quantities are defined:

Using the above dimensionless variables, (9) is rewritten aswhere

In (12), the primes denote differentiation with respect to , and the superposed dots denote derivatives with respect to the time .

The solution of (12) can be written as

For cantilever beams, has to meet the following boundary conditions:

The mode function that satisfies the boundary condition in (12) can be expressed as follows:where .

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.where

Here, the superscript (4) means the fourth-order partial derivative of .

##### 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:

and are expanded in the formwhere .

Taking the derivative of (20), one obtains the following equations:

Substituting (20) and (21) into (17) and considering (19), one can obtain the following equations by comparing the coefficients of and :

O():

O():

Assume the solution of (22) is as follows:where represents the imaginary unit, and are the complex functions to be determined, and and represent the complex conjugate. and refer to the forward and backward whirling frequencies, respectively, which are expressed as follows:

Substituting (24) into (23) and dropping the constant terms and , one obtains the following equations:where

In this paper, four cases of primary and superharmonic resonances, i.e., case I: , case II: , case III: , and case IV: , are studied.

###### 2.4.1. Case I: Primary Resonance ()

Assumewhere *σ* represents the detuning parameter.

Substituting (24) into (26), one can obtain the following equations:wherewhere .

The particular solutions of (29) are

Substituting (31) into (29) and equating the coefficient of and in both sides of (29), one has

Equations (32) and (33) constitute systems of two inhomogeneous algebraic equations for , and , , respectively. Their homogeneous parts have a nontrivial solution. Then their solvability conditions can be written as [36]

After simplification, the solvability conditions are reduced to two independent governing equations in the following form:

Assume that the solutions of (36) and (37) arewhere and (*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:where (*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 is zero. This shows that under the primary resonance , 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.

Substituting = 0 into (39) and solving for *σ*, one obtains

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 *a*_{2} = 0) around (*a*_{0}, *ψ*_{0}).

###### 2.4.2. Case II: Primary Resonance ()

According to (26), under such conditions, one introduces detuning parameters as

Using their solvability conditions (which are not shown) and (38) leads to

By using steady-state condition (i.e., set = 0, = 0, = 0, and = 0), the following steady-state response of case II can be obtained:

From (44), it can be seen that and , indicating that only the forced vibration response of backward mode is excited, while forward mode remains stationary.

###### 2.4.3. Case III: Superharmonic Resonance ()

For this case, the formulation of steady-state response is the same as case I, but in (41), and are replaced with and , respectively.

###### 2.4.4. Case IV: Superharmonic Resonance ()

For this case, the formulation of steady-state response is the same as case II, but in (43), and are replaced with and , 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 Table 1. 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 . In all cases, , except where other values are mentioned.

Figure 2 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 determine the critical rotating speeds of the rotating cutting tool.

##### 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 *a*_{lim} curve separates the space into two regions. Any (, *a*_{lim}) 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.

##### 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 Figure 4. As shown in Figure 4, by increasing detuning parameters at point A, the amplitude *a*_{1} gradually increases until reaching point B. *a*_{1} jumps downward from point B to point C. Afterward, *a*_{1} 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*, *a*_{1} increases and jumps upward from point *E* to point F. As the detuning parameter further decreases, *a*_{1} 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, *p*_{1}, *p*_{2}, and *p*_{3}, among which *p*_{1} and *p*_{3} are stable and *p*_{2} 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.

Figures 5–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 . 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 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 , , with the parameters and 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 *E*_{11} along the longitudinal direction of the fiber is significantly larger than the transverse elastic modulus *E*_{22} (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 9 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.

Figures 10–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 . From these figures, it can be seen that for some values of *L*/*d*, *σ*, *K*_{f}, or , there are multivalued curves. For example, when *L*/*d* = 10 and ＜0.45, as is shown in Figure 10, the system has two stable and one unstable branches, but for *L*/*d* = 10 and ＞0.45, there exists only one stable branch. For large values of and large values of *L*/*d*, curves are always single-valued.

Figures 15–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 (). Similar to the cases in Figure 10～14, again, for some values of *L*/*d*, *σ*, *K*_{f}, or , multivalued curves can be observed.

Figures 20–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*, *σ*, *θ*, , or .

Figures 25–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*, *σ*, *K*_{f}, or .

###### 3.2.2. Case II

Figures 28–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 . 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.

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

In addition, the equivalent damping coefficient of case II can be expressed as

Figures 34–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 (see Figure 38).

Moreover, for some values of *L*/*d*, *σ*, *K*_{f}, , , 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).

Figure 39 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.

###### 3.2.3. Case III

For this case, the equivalent coefficient is identical with that in case I, but the amplitude of excitation force in this case is lower than in case I (e.g., using the same parameters L/d = 10, *θ* = 90°, and *K*_{f} = 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.

#### 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 coefficients and which include the structural damping and the damping from regenerative chatter mechanism. Cases I and III have the same equivalent damping coefficient , and cases II and IV have the same equivalent damping coefficient . It has been found that is larger than , 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].

#### Appendix

#### A. The Stain Energy of the Composite Cutting Tool

If the torsional deformation is negligible, the bending displacements of the composite cutting tool in *x*, y, and *z* directions are as follows:

The strain in the *x* direction iswherewhere and are the linear strain and the nonlinear strain generated in high-order bending deformation, respectively.

Bending energy for the composite cutting tool can be expressed as

The elastic stress-strain relations of the laminate cutting tool can be expressed as

Employing (A.2), (A.3), and (A.5), (A.4) can be rewritten as

Due to the symmetry of the cross section, the third and seventh terms in the above equation are equal to zero.

Therefore, the strain energy of the composite cutting tool can be reduced as follows:

#### B. The Expression of the Off-Axis Stiffness Coefficient

The off-axis stiffness coefficient of the *k*th layer for the composite cutting tool is defined aswhere *θ*^{(k)} is the layer angle of each layer of the material.

The expressions of *C*_{11}, *C*_{12}, *C*_{22}, and *C*_{66} are as follows:where

#### Abbreviations

: | Kinetic energy |

L: | Length of the cutting tool |

t: | Time |

: | Density |

A: | Cross-sectional area |

I: | Second moment of area of the beam cross section |

v, w: | Transverse deflections of the shaft element |

x, y, z: | Variable coordinates |

: | Rotation angles of the cross section around the y and z axes |

: | Rotating speed of the shaft |

: | Off-axis stiffness coefficient of the kth layer for the composite cutting tool |

: | Potential energy |

: | Rayleigh dissipative energy function of the cutting tool |

c: | Structural damping coefficient |

δ: | Variational operator |

δW: | Virtual work of the cutting force |

F_{y,}F_{z}: | Linear regenerative cutting forces |

: | Delay time of milling process |

N: | The number of teeth |

K_{tc}: | Cutting force coefficients in the tangential directions |

K_{te}: | Edge force coefficients in the tangential directions |

K_{rc}: | Cutting force coefficients in the radial directions |

K_{re}: | Edge force coefficients in the radial directions |

: | Axial cutting depth |

c_{f}: | Feed per tooth per revolution |

(EI)_{equiv}, (EA)_{equiv}: | Equivalent bending and tensile stiffness of the cutting tool |

(ρI)_{equiv}, (m)_{equiv}: | Equivalent diametrical mass moment of inertia and equivalent mass per unit length |

: | Density of the kth layer |

r_{k}, r_{k+1}: | Inner diameter and outer diameter of the kth layer |

: | Mode function |

, : | Forward and backward whirling frequencies |

: | Small parameter |

i: | . |

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The research was funded by the https://doi.org/10.13039/501100001809National Natural Science Foundation of China (Grant No. 11672166).