## Abstract

The classical model of regenerative vibration is investigated with new kinds of nonlinear cutting force characteristics. The standard nonlinear characteristics are subjected to a critical review from the nonlinear dynamics viewpoint based on the experimental results available in the literature. The proposed nonlinear model includes finite derivatives at zero chip thickness and has an essential inflexion point. In the case of the one degree-of-freedom model of orthogonal cutting, the existence of unstable self-excited vibrations is proven along the stability limits, which is strongly related to the force characteristic at its inflexion point. An analytical estimate is given for a certain area below the stability limit where stable stationary cutting and a chaotic attractor coexist. It is shown how this domain of bistability depends on the theoretical chip thickness. The comparison of these results with the experimental observations and also with the subcritical *Hopf bifurcation* results obtained for standard nonlinear cutting force characteristics provides relevant information on the nature of the cutting force nonlinearity.

## 1. Introduction

One of the main goals of cutting process optimization methods is to maximize the volume of the chip cut within a certain time. There are several boundaries identified in the operational space of the cutting parameters, namely chip width, chip thickness and the cutting speed. These boundaries are related to the maximum power, maximum cutting force, feed rate, depth of cut, etc. The most difficult boundary to model is the onset of harmful relative vibrations between the tool and the workpiece. The so-called regenerative effect (where the tool interacts with its delayed displacement via the surface profile that is cut one revolution earlier) is considered to be one of the main reasons for these vibrations, which cause poor surface quality or, in extreme cases, damage the machine tool structure.

The regenerative effect for the simple orthogonal cutting model was first introduced and analysed in the middle of the twentieth century (Tlusty & Spacek 1954; Tobias 1965). The central idea of this effect is that the motion of the tool depends on its past motion; in other words, a time delay occurs in the slightly damped oscillator model of the machine tool. This delay is inversely proportional to the cutting speed. The first thorough and detailed experiments on the *nonlinear* regenerative vibrations often showed small domains of attraction around stable stationary cutting (Shi & Tobias 1984). These investigations demonstrated that stable stationary cutting can be quite sensitive to external perturbations and stable, large-amplitude vibrations can appear even in those parameter domains where the stationary cutting is linearly stable.

A rigorous analytical investigation of the nature of the loss of stability of stationary cutting was performed much later only by means of the centre manifold reduction and normal form calculations (Stepan & Kalmár-Nagy 1997; Kalmár-Nagy *et al*. 2001).

The subcritical sense of the cutting process has been shown for several other cutting models too (Nayfeh *et al*. 1997; Kalmár-Nagy & Pratt 1999; Campbell & Stone 2004; Szalai *et al*. 2004; Wahi & Chatterjee 2005). In the substantial situation, the bifurcating local unstable orbit separates two independent attractors, namely stationary cutting and a large-amplitude nonlinear oscillation that is itself stable in a dynamical systems sense. In engineering, this latter vibration is often referred to as *instability*, owing to its harmful nature. The region where this coexistence can occur is called the region of *bistability*, while the terminology *unsafe zone* is used for the same idea in production technology to indicate the possible chatter.

Despite the fact that large-amplitude vibrations are of little interest from a technological standpoint, the location and size of the bistable domain are important. This is because they define the parameter region where the cutting process is more or less sensitive to perturbations caused by, for example, non-homogeneous workpiece material. Some of the existing software packages used for technological parameter optimization are capable of predicting the linearly stable parameter domains for many types of cutting tools, workpieces and machine tool configurations, but they are unable to predict those unsafe zones within the stable domain where unexpected and persistent vibrations can still occur during otherwise stable machining.

In this paper, we give a critical discussion of available nonlinear cutting force characteristics models partly from the point of view of simple geometric characteristics (e.g. derivative at zero chip thickness and the existence of an inflexion point). For the proposed nonlinear model, the subcriticality of the *Hopf bifurcation* is proven, which is related to the cutting force characteristics at a putative inflexion point. With the help of analytical and semi-analytical investigations, the location and the size of the bistable regions are then given. These regions are characterized as a function of the theoretical chip thickness that is strongly related to the existence of an inflexion point in the cutting force. Experimental results (Shi & Tobias 1984; Endres & Loo 2002) in the literature confirm the validity of these conclusions in terms of direct cutting force measurements and the identification of the bistable region.

## 2. Model construction

The description of chip separation during cutting is a challenging task in both modelling and mathematical treatment. Several authors have used finite-element models to try to describe accurately the coupled nonlinear physical effects near the cutting edge (Davies & Burns 2001; Stone *et al*. 2005). While these models simulate realistic cutting processes well, the structure and the dependence on the parameters of the consequent nonlinear vibrations remain hidden. By contrast, physical effects during cutting can be modelled more simply by means of empirical nonlinear cutting force characteristics based on the extensive measurement results of the production technology community. The large-scale geometric and elastic nonlinearities of the machine tool structure have a far smaller influence on the nonlinear vibrations than those of the cutting force characteristics that exist on a small scale. In addition, the linear elasticity of the machine tool is modelled in the direction of the essential mode only, which is related to the lowest natural frequency of the structure and supposed to be perpendicular to the cut surface.

Thus, the mechanical model used here is a simple one degree-of-freedom damped oscillator subjected to a nonlinear cutting force *F* (figure 1). The corresponding governing equation has the form(2.1)where *ω*_{n} and *κ* are, respectively, the natural angular frequency and the damping ratio of the essential vibration mode described by the general coordinate *q* that refers to the tool position. We can express these parameters with the modal mass *m*, the stiffness *k* and the damping factor *b* (figure 1*b*); here, and *κ*=*b*/(2*mω*_{n}), while *F*_{q}(*t*) is the corresponding component of the resultant cutting force *F*(*t*).

### (a) Empirical cutting force characteristics

The most popular and generally applied relation for nonlinear cutting force characteristics is the power law(2.2)where *w* and *h* are, respectively, the chip width and the chip thickness of the orthogonal cutting and *K*_{q} and *ν* are the empirical parameters depending on the cutting conditions (material, tool geometry, etc.). The exponent *ν* may vary from 2/5 (Kalmár-Nagy & Pratt 1999) through to the most popular choice from 3/4 (Kienzle 1957) to 4/5 (Tlusty & Spacek 1954). The power expressions for cutting force were introduced in order to enable linear optimization tasks in the log–log space of the technological parameters (Taylor 1907).

Another approach to modelling cutting force characteristics is to fit a polynomial function to the experimental data. One of these less frequently used expressions is the cubic polynomial approximation (Shi & Tobias 1984) of the cutting force characteristics,(2.3)Clearly, here *ρ*_{1} and *ρ*_{3} are positive because the cutting force has positive gradient both at low and at high values of the chip thickness (figure 2). The potentially negative value for *ρ*_{2} allows the existence of an inflexion point on the cutting force characteristics at(2.4)From a dynamics viewpoint, there are essential differences between the two empirical interpretations (2.2) and (2.3) of the cutting force. One of these differences has a rather theoretical (or philosophical) nature. The traditional power-law characteristics have a vertical tangent at the origin (figure 2*a*), where the tool just touches the surface of the workpiece with zero chip thickness. This feature causes problems in the mathematical treatment of the vibrations at the loss of contact because the uniqueness and regularity of the solutions no longer hold here. Apart from causing uncertainties and unpredictable errors in numerical simulations, this non-uniqueness of solutions in forward time is questionable in physical systems.

The other important difference between the power-law (2.2) and the cubic polynomial force characteristics (2.3) is related to the possible existence of an inflexion point (figure 2*a*).

In the literature, there are several extensive measurement results that provide a basis for cubic polynomial approximations of the cutting force, both for milling and for turning. Shi & Tobias (1984) performed measurements based on full-immersion milling with a face mill of 24 teeth. Owing to the high number of cutting edges, the time periodicity of the parametric tooth–pass excitation could be averaged, and the mean cutting force characteristics with respect to the mean chip thickness were found to exhibit an inflexion point (figure 2*b*).

Another set of experimental data measured by Endres & Loo (2002) during valve-seat turning also shows the existence of an inflexion point in the cutting force characteristics (figure 2*c*). However, the authors fitted an exponential function to the measured data, and hence the inflexion point could not appear in the analytical approximation and its effect was not examined further.

The many combinations of cutting tool and workpiece material require many measurements using different speeds, feeds and different arrangements of the tool–workpiece system. The main goal of these kinds of measurements is to determine the cutting coefficients in different directions, which are essential to predict the stability of certain cutting operations (Altintas & Budak 1995; Altintas 2000; Bayly *et al*. 2003; Warminski *et al*. 2003; Insperger & Stepan 2004; Merdol & Altintas 2004; Zatarain *et al*. 2006). The last set of experimental data in figure 2*d* was collected from a series of orthogonal cutting identification tests (see acknowledgements).

In subsequent sections, we will show that the inflexion of the cutting force characteristics has a central role in the nonlinear vibrations and, consequently, in the size of the unsafe or bistable zone.

### (b) Regenerative effect

Figure 1*b* presents the variation in the instantaneous chip thickness *h*(*t*) as a function of the present position *q*(*t*) and the delayed position *q*(*t*−*τ*) of the tool given by(2.5)where *h*_{0} is the prescribed chip thickness and *τ*=2*π*/*Ω* is the period of the rotating workpiece, where *Ω* is its angular velocity. Since the cubic expression (2.3) of the cutting force depends on the actual chip thickness (2.5), the equation of motion (2.1) has the form(2.6)Its trivial solution provides the equilibrium corresponding to stationary cutting atwhere the spring force is balanced by the cutting force (2.3). Introducing the perturbation *x* around this position by , and rescaling time *t* with the angular natural frequency *ω*_{n}, and the tool position perturbation *x* with the theoretical chip thickness *h*_{0} leads to the dimensionless variables and parameters,(2.7)In the dimensionless form of the equation of motion, the derivatives with respect to dimensionless time are denoted by primes, while all the tildes are henceforth dropped. With the vectorthe new form of the equation of motion (2.6) is(2.8)where the linear part is defined by the non-delayed and delayed coefficient matricesrespectively, and the nonlinear part is given by

### (c) Operator formulation

From the theory of functional differential equations (Hale 1977), the phase space of the delayed systems is infinite dimensional. In order to represent the equation of motion in the infinite-dimensional phase space of the continuous functions defined by the shiftone can transform the delay-differential equation (DDE) (2.8) into the operator differential equation (OpDE)(2.9)where the linear operator and the nonlinear operator are defined byandThe superscript ‘o’ denotes the derivative with respect to the past time *θ*, which equals the derivative with respect to the (dimensionless) timeThe OpDE formulation (2.9) corresponds to the DDE form (2.8) or to the traditional form (2.6) of the equation of motion and makes it possible to carry out a correct nonlinear analytical investigation via linear stability analysis, centre manifold reduction and normal form calculations.

## 3. Linear stability

From an engineering point of view, the linear stability charts are preferred to be represented in a plane formed by the spindle speed and chip width. In the stable domains of these charts, those parameter regions where stationary cutting is possible without chatter are identified. These regions can be determined by the stability investigation of the linear part of the OpDE (2.9). The usual substitution of the exponential trial solution into the linear OpDE leads to the infinite-dimensional eigenvalue/eigenvector problem(3.1)where denotes the unit operator. This provides the characteristic functionfor the complex characteristic exponents *λ* satisfying the characteristic equation(3.2)If the infinitely many complex characteristic roots *λ*_{k} (*k*∈) of the characteristic function (3.2) satisfy Re *λ*_{k}<0 for all *k*, then the trivial solution (i.e. stationary cutting) is locally asymptotically stable; if there exists a *k* with Re *λ*_{k}>0, then it is unstable. At the linear stability boundaries, the characteristic function (3.2) has purely imaginary complex conjugate roots *λ*=±i*ω*, where is the dimensionless angular frequency of the self-excited vibrations arising close to the stability limit.

By substitution of the purely imaginary characteristic roots into (3.2), the stability limits (so-called lobes) can be expressed analytically as a parametric function of the frequency *ω* (e.g. Stepan 2001),(3.3)The minimum value of the stability limits can be expressed in closed form byThe stability boundaries *w*_{stab} will be normalized with respect to this minimal value *w*_{min} in order to show the stability charts independent of the small (approx. 0.01) damping ratio parameter. These expressions provide the stability boundary curves (or lobes) *w*_{stab}(*ω*) in figure 3*a*, while the (dimensionless) vibration frequency *ω* against the cutting speed is given above the lobes at the stability limits in figure 3*b*. Note that the frequencies of the corresponding self-excited vibrations are all larger than the natural frequency of the system, i.e. *ω*>1. The lobe structure in figure 3*a* is typical for delayed oscillators (see Hu & Wang 2002; Kyrychko *et al*. 2006).

## 4. Nonlinear investigation

In this section, an overview of the applied algebraic and numerical techniques is given to clarify the dynamic effects of the nonlinear cutting force characteristics. First, the lengthy Hopf bifurcation calculation is summarized in order to provide analytical estimates for the nonlinear vibrations in the system. Then, the numerical methods are used to check the different kinds of nonlinear behaviours and also to make the results quantitatively more accurate.

### (a) Hopf bifurcation calculation

The fixed point becomes non-hyperbolic at the stability limit (figure 3*a*) with increasing (dimensionless) chip width *w* that will be considered as a bifurcation parameter. At its critical value *w*_{stab}(*ω*) in (3.3), two pure imaginary eigenvalues *λ*_{1,2}=±i*ω* exist. As the bifurcation parameter *w* increases through *w*_{stab}(*ω*), these characteristic roots cross the imaginary axis with non-zero (actually, positive) speed since the implicit differentiation of the characteristic equation (3.2) giveswith the positive expressions(4.1)andfor all *ω*>1 and *κ*>0. This means that the so-called transversality condition is fulfilled (Guckenheimer & Holmes 1983), and a Hopf bifurcation occurs along the linear stability limits. Consequently, a periodic orbit exists in a sufficiently small region of the fixed point. The stability/instability of this periodic motion depends on the super/subcritical sense of the Hopf bifurcation. The Hopf bifurcation calculation results in the *Poincaré–Lyapunov constant* (PLC) Δ(*ω*) (less than or greater than 0) that determines the (super or subcritical, respectively) sense of the bifurcation. For the bifurcation parameter *w* close enough to its critical value *w*_{stab}(*ω*), this also leads to the estimate(4.2)for the amplitude *r* of the corresponding approximately (2*π*/*ω*)-periodic motion that is located in the two-dimensional centre manifold embedded in the infinite-dimensional phase space. At the origin, this centre manifold is tangent to the plane spanned by the real and the imaginary parts of the critical eigenvectors *s*_{1,2} of the operator belonging to the critical eigenvalues *λ*_{1,2}=±i*ω*. Thus, the periodic solution of the OpDE can be approximated as(4.3)and, consequently, the approximation of the periodic solution of the DDE (2.8) assumes the formTo determine the critical eigenvectors, it is easy to solve the boundary-value problem defined by the infinite-dimensional eigenvalue/eigenvector problem (3.1) at the borders of stability,(4.4)However, the calculation of the PLC Δ(*ω*) is rather lengthy (Stepan 1989) and not presented here. The calculation follows the same algebraic steps described for similar delayed oscillator examples in Campbell & Bélair (1995), Nayfeh & Balachandran (1995), Stepan & Kalmár-Nagy (1997), Kalmár-Nagy *et al*. (2001) and Orosz & Stepan (2004), and it can be found in Dombovari (2006) and Dombovari *et al*. (2007) in detail. The result is(4.5)whereandand the always positive *γ*_{1}(*ω*), *τ*(*ω*) and *w*_{stab}(*ω*) are given by (4.1) and (3.3), respectively. Clearly, the sign of Δ(*ω*) depends on the possible negative values of *δ*_{1}(*ω*). A lengthy but straightforward algebraic manipulation proves the following analytical estimate:Actually, numerical investigation shows that the stricter estimate *δ*_{1}(*ω*)/*δ*_{2}(*ω*)>−0.36 holds as well. Also, the cutting force characteristics must always be increasing, which means, for the three parameters *ρ*_{1,2,3} of the cubic polynomial cutting force *F*_{q} in (2.3), thatUsing these estimates and the relations (2.7) between the parameters *η*_{2,3} and the original *ρ*_{1,2,3} parameters of the nonlinear cutting force, lower bounds can be given for the PLC byThus, the Hopf bifurcation is subcritical all along the stability lobes if the gradient of the cubic cutting force function is positive at the inflexion point. A negative derivative at a possible inflexion point on the cutting force characteristic is physically unreasonable and has not been experienced in machine tool dynamics. As an important result of the above Hopf bifurcation calculation, it can be concluded that unstable oscillations exist around the locally stable stationary cutting state when the chip width gets close to its critical value along the linear stability boundaries. This result is a generalization of the same conclusion of Kalmár-Nagy & Pratt (1999) and Wahi & Chatterjee (2005), which was derived for the power-law approximation of the cutting force, where hence the possibility of an inflexion point was excluded.

### (b) Amplitude estimation of unstable oscillation

The approximation of the unstable periodic motion gives useful information about the kind of perturbations that allow stationary cutting. The substitution of the eigenvectors (4.4) into the periodic solution estimate (4.3) gives(4.6)while the substitution of the PLC (4.5) into the amplitude estimate (4.2) results in(4.7)Clearly, unstable oscillations exist for *w*<*w*_{stab}(*ω*), that is, when the chip width is smaller than its critical value. The corresponding dimensionless bifurcation diagrams at three points of the second lobe of the stability chart are presented in figure 4*b* as dashed lines. The curves represent unstable vibration amplitudes for *w*/*w*_{stab}(*ω*)<1 with different curvatures for different (dimensionless) cutting speeds.

### (c) Path-following

The approximation of the unstable limit cycle gets worse as the chip thickness parameter gets far below its critical value at the stability limit. Moreover, below a certain limit, the unstable periodic orbit does not exist at all since the vibration amplitudes become so large that the tool actually loses contact with the workpiece. In order to check the range of validity (and also the correctness) of the analytical treatment, the continuation software DDE-Biftool (Engelborghs *et al*. 2002) was applied to follow the closed unstable orbits in the parameter space. This software uses the collocation method to solve the boundary-value problem of the returned closed orbits of a delayed system.

The corresponding unstable periodic orbits are represented by their amplitude curves denoted by solid lines in the dimensionless bifurcation diagram of figure 4*b*. The results of the path-following method show very good quantitative agreement with the analytical results of the Hopf bifurcation calculation in the region 0.9<*w*/*w*_{stab}(*ω*)≤1.

## 5. Region of bistability

As explained above, the identified unstable periodic motion exists only for that parameter region of the chip thickness where the tool does not leave the workpiece. This region is called bistable or the unsafe zone because the vibrations will settle to the locally stable stationary cutting inside the unstable orbit only, while outside the unstable orbit the vibrations will grow until the motion settles to a large-amplitude complex oscillation that involves the repeated loss of contact between the tool and the workpiece (figure 5).

In order to find the boundaries of the unsafe zone where the unstable periodic orbits exist, one can consider the projection of the infinite-dimensional phase space, that is, the projection to the plane , or equivalently, to the plane of the delayed and the actual dimensionless perturbations. These projections of the unstable orbits (both their analytical and numerical approximations) are presented in figure 4*c* for chip width parameters *w*_{loss}(*ω*), where they just touch the switching line that represents the first loss of contact between the tool and the workpiece.

The definition (2.5) of the actual chip thickness provides a simple algebraic formulation for the condition of the loss of contact in the form(5.1)which determines the switching line of slope 1 shifted by 1 in figure 4*c*. Substituting the analytical estimate (4.6) of the critical unstable periodic motion into (5.1) at the critical chip width parameter *w*_{loss}(*ω*) and using formulae (3.3) for the linear stability limit yieldsWith the help of the amplitude estimate (4.7), the loss of contact condition (5.1) can be reformulated in terms of the chip thickness parameter *w*_{loss} as a function of the frequency parameter *ω* used to parametrize the stability lobes,These values of the chip width are also depicted in figure 4*b* marked with tiny circles and by the text ‘switching point’ for the three different cases of cutting speeds. These values, which characterize the loss of contact parameters, are in the range 0.9–0.92. This range is in agreement with the range measured experimentally by Shi & Tobias (1984). This means that the bistable or unsafe zone is within the rangeand the curves representing the unstable orbits in figure 4*b* are extended hypothetically to the left of the switching points at *w*_{loss}(*ω*), since they do not exist for *w*<*w*_{loss}(*ω*).

Finally, the bistable regions can also be represented in the stability charts. The shaded regions below the stability lobes in figure 6*b* can be determined easily with the help of the ratio of the width of the bistable region relative to the critical chip width parameter at the linear stability limit,(5.2)Note that the parameters *η*_{2,3} of the nonlinear cutting force characteristics depend strongly on the theoretical chip thickness *h*_{0}, as indicated by formulae (2.7). Consequently, the size of the bistable region also depends on *h*_{0} as it can be observed also in figure 6. In the case of the power-law approximation, the width of this bistable region does not vary at all.

Figure 6*a* shows the size of the bistable region in per cent according to formula (5.2) as a function of the chip thickness *h*_{0}, where the nonlinear cutting force parameters *ρ*_{1,2,3} are fixed at the values of figure 2*b* taken from Shi & Tobias (1984). The experimental results from Shi & Tobias (1984) are denoted by filled circles and they match the analytical prediction closely. This is remarkable because the conventional power-law characteristics provide a constant approximately 4 per cent width for the bistable region (Stepan & Kalmár-Nagy 1997). The most critical chip thickness *h*_{cr}, where the size of the bistable region is maximal, can be calculated approximately analytically at a theoretical case when the vibration frequency tends to infinity, i.e.In the meantime, the quantitatively more accurate path-following method determined the switching points more precisely, and the real size of the maximal bistable zone is shown to be about half the size of the analytical estimation (figure 6). This value is nevertheless at approximately 50 per cent, which is approximately 12 times larger than the one predicted by the power-law formulation.

## 6. Conclusions

Stationary cutting force measurements indicate that the nonlinear cutting force characteristics often involve an inflexion point in the otherwise monotonic increasing function of chip thickness. The conventional power-law approximation of the empirical cutting force is not able to describe this inflexion point. We have analysed the effect of this inflexion point from the nonlinear dynamics viewpoint and showed the practical importance of the nonlinear expression of the cutting force.

The conventional power-law approximation has several advantages from the viewpoint of the optimal design of the parameters of the technology, but it has serious disadvantages from the viewpoint of predicting machine tool chatter. On one hand, the application of the power law in a Newtonian equation of machine tool vibration violates the uniqueness of the solution in forward time and also causes unpredictable errors during the numerical simulation of large-amplitude oscillations between the tool and workpiece. On the other hand, the missing inflexion point results in the prediction of a uniform and thin bistable region along the stability limits that does not depend on the mean (or theoretical) chip thickness. This bistable region jeopardizes the chatter-free cutting, even for stable stationary cutting processes designed and optimized by linear theory, and the power-law approximation of the cutting force heavily underestimates the real size of this zone.

With the help of realistic cubic polynomial approximations of the nonlinear cutting force, we have proved the uniform subcriticality of the Hopf bifurcations along all the linear stability limits of regenerative cutting processes. This subcriticality is ensured mathematically by the natural positive slope of the cutting force characteristics at the possible inflexion point.

Furthermore, we have estimated the width of the bistable region analytically. It was shown that its size varies along the stability limits and this bistable region could be much larger than expected before. For both power-law and cubic polynomial cutting forces, the linear stability boundaries increase slightly for increasing values of the theoretical chip thickness. However, the size of the bistable region is much (actually, approx. 12 times) larger in the case of the presence of an inflexion point in the cutting force characteristics than the size of the bistable region predicted by the power-law approximations of the cutting force. The analytical results were checked by continuation software, and the above conclusions were confirmed by the existing experimental results of the specialist literature (Shi & Tobias 1984), which were revisited and re-examined from the nonlinear dynamics points of view.

## Acknowledgments

This research was partially supported by the Hungarian Scientific Research Foundation OTKA grant no. K68910, the Spanish–Hungarian Science and Technology Program grant no. 8/07 and the EU Socrates Action. The orthogonal cutting data for Al7075 were kindly provided by Prof. Y. Altintas, Manufacturing Automation Laboratory, University of British Columbia.

## Footnotes

- Received April 15, 2008.
- Accepted July 18, 2008.

- © 2008 The Royal Society