Nonlinear dynamics of high-speed milling

(1)

Róbert Szalai

Nonlinear dynamics of high-speed milling

Ph.D. dissertation

Advisor: Gábor Stépán, DSc.

September 2005

Budapest University of Technology and Economics

Department of Applied Mechanics

(2)

A dolgozat magyar nyelv¶ címe: Nagysebesség¶ marás nemlineáris rezgései

A dolgozat bírálatai és a védésr®l készült jegyz®könyv a BME Gépészmérnöki kari dékáni hivatalban elérhet®.

Nyilatkozat

Alulírott Szalai Róbert kijelentem, hogy ezt a doktori értekezést magam készítettem és abban csak a megadott forrásokat használtam fel. Minden olyan részt, amelyet szó szerint, vagy azonos tartalomban, de átfogalmazva más forrásból átvettem, egyértelm¶en, a forrás megadásával megjelöltem.

Budapest, 2005. szeptember 23.

Szalai Róbert

(3)

Abstract

In this thesis the local and global dynamics of the high speed milling process is analyzed. We develop dierent analytical and numerical techniques to reveal the dynamical properties of this process. We introduce two types of models, which are a piecewise smooth map and delay dierential equations. The discrete time model (map) is accessible to analytical techniques thus closed form results are deduced about the bifurcations of the stationary cutting and even chaotic motions are demonstrated for certain range a parameters. We prove that one of the delay dierential equation models has unstable islands in its stability chart. With the help of our method of continuing periodic solution bifurcations, global stability of the model is determined, which result clearly shows the sensitivity of the process near the stability boundary. Also, an analytical technique is developed, in order to rigorously prove the subcritical sense of the period doubling bifurcation.

(4)

4 Summary

High speed milling is a common machining process having wide spread appli- cations in the industry. Its application elds include the aerospace industry, the automotive industry, etc., where usually 80%-90% of the workpiece material is re- moved in order to achieve the prescribed sculptured design. One of the the most advantegous feature of this process that it makes possible to use small number of complicated parts instead of the usual large number of simple parts. This feature makes the entire manufacturing process less expensive: although cutting becomes a bit more expensive, in the assembly part great cost reductions can be achieved.

Despite the above advantages, the process has its limitation, as well. During the cutting process appearance of unwanted vibrations is expected, which results in rough surfaces, extended tool wear and can even damage the machine-tool it- self. The common mechanism of these vibrations are described in this thesis, which are found to be of self-excited nature. Unfortunately the unwanted vibrations arise already in the linearly stable regions of the stability chart, which is related to subcritical bifurcations. In the rst and second theses we draw these consequences from a discrete time model, that is, we prove that the period doubling and the Neimark-Sacker bifurcation related to the classical machine-tool chatter are both subcritical. Since the discrete time model is fairly simple as a mathematical ex- pression, the global dynamics can be analyzed in detail. We nd `outer' period two orbits, which are either attractors for some parameters or repellors (unstable) for other parameters. If they are unstable together with the stationary cutting, chaotic motions develop in the system.

The second part of the dissertation deals with delay-dierential equation models, which are capable of quantitative predictions in the case of a rather general class of high-speed milling processes. In contrary to the common belief, we proved that for a certain type of processes there are unstable islands in the stability chart.

A numerical technique was developed to show similar bifurcation results as for the discrete time model. Further, the stability of the arising unwanted vibrations was determined by following the bifurcations of the period-two orbits along two parameters. The parameter region of the unwanted vibrations was also determined by the same numerical technique and it was found that unwanted vibrations overlap with the stable cutting, which makes the process sensitive to perturbations in those critical regions. Existing measurement data in the literature was re-examined and analyzed to underpin the results of our method and show that subcritical and supercritical bifurcations do occur in milling processes. In addition, the numerical results were checked by analytical calculations, that is, the subcritical sense of a period doubling bifurcation was proved.

(5)

5 Összefoglalás

A nagysebesség¶ marás manapság az egyik legkedveltebb és leghatékonyabb forgácsleválasztási eljárás az iparban. Alkalmazzák a repül®gépgyártásban, ahol a könny¶fém munkadarabok térfogatának akár 80-90 százalékát is eltávolítják a kí- vánt szoborfelület elérése érdekében. A folyamat egyik el®nye, hogy sok egyszer¶

alkatrész használata helyett kis számú bonyolult szerkezeti elem alkalmazását teszi lehet®vé, ami egyrészt növeli a megbízhatóságot, másrészt csökkenti a költségeket.

Noha a megmunkálási költségek kis mértékben növekednek nagysebesség¶ marás alkalmazásával, a lecsökkent munkadarabszám miatt a felmerül® raktározási, szál- lítási és szerelési költségek jelent®sen redukálhatók, ezáltal a gyártás gazdaságosabb.

A fenti el®nyök ellenére ezen forgácsolási eljárásnak korlátai is vannak. A for- gácsleválasztás során gyakran lépnek fel olyan (nemkívánatos) rezgések amelyek egyenetlen felületet okoznak a munkadarabon, csökkentik a szerszám élettertamát és akár a szerszámgépet is károsíthatják. Ezek a rezgések rendszerint nem rezo- nanciából, hanem ún. öngerjesztett rezgések, amelyek a szerszám-munkadarab kap- csolatának nemlineáris jellegéb®l származtathatók. A folyamat lineáris stabilitási kérdéseivel a szakirodalom sokat foglalkozott az utóbbi években, de ez nem nyúj- tott megfelelel® magyarázatot a még lineárisan stabil tartományokban kialakuló rezgésekre. Az els® és második tézisben megmutattam, hogy a stabil tartományban fellép® nemkívánatos rezgések a stabilitásvesztéskor bekövetkez® szubkritikus peri- óduskett®z® és Neimark-Sacker bifurkációk következményei. Mivel a diszkrét idej¶

modell matematikaileg kell®en egyszer¶, a globális dinamika is részletesen vizsgál- ható. A modellben küls® dupla periódusú pályákat találtam, amelyek stabilak és instabilak is lehetnek. Abban az esetben ha a küls® pálya instabil a stacionárius vágással együtt, a rendszer kaotikus viselkedést mutat.

A dolgozat második része késleltetett dierenciálegyenlet modellekkel foglalko- zik, amelyek már képesek a fellép® jelenségek kvantitatív leírására a nagysebesség¶

marás több változata esetén is. Megmutattam, hogy a széles körben elterjedt né- zetekkel szemben a nagysebesség¶ marás egy modellje esetében a lineáris perió- duskett®z® stabilitási határok instabil szigeteket formáznak a stabilitási térképen.

Egy új numerikus módszert fejlesztettem ki, és hasonlóan a diszkrét idej¶ modell esetéhez meghatároztam a fellép® bifurkációk értelmét. Kiszámítottam továbbá a nemkívánatos rezgések stabilitását azok bifurkációinak két paraméter menti köve- tésével. Hasonlóképpen meghatároztam ezen rezgések létezési tartományát, és azt találtam, hogy a nemívánt rezgések tartománya átfedi a stabil stacionárius vágás tartományát, ami a folyamat zavarásokra való érzékenységét okozza. Az irodalom- ban megtalálható mérési adatokat újra feldolgozva ellen®riztem numerikus mód- szerem eredményeit, és megmutattam, hogy nagysebeség¶ marás esetén valóban el®fordul szubkritikus és szuperkritikus perióduskett®z® bifurkáció. A numerikus módszer eredményét analitikus módon is ellen®riztem, azaz bebizonyítottam egy perióduskett® bifurkáció szubkritikusságát.

(6)

The idea is to try to give all the information to help others to judge the value of your contribution; not just the information that leads to judgment in one particular direction or another.

Richard Feynman There is a theory which states that if ever anybody discovers exactly what the Universe is for and why it is here, it will instantly disappear and be replaced by something even more bizarre and inexplicable. There is another theory which states that this has already happened.

Douglas Adams

(7)

7 Acknowledgement. I would like to thank to Gábor Stépán for his constant support and invaluable insights. S. John Hogan at the University of Bristol hosted me rst in 2002 when I was writing my Master's thesis and also in the spring of 2004, when I was granted a Hungarian Eötvös Scholarship. This thesis was written at the Massachusetts Institute of Technology where I have spent a year as a Fulbright Fellow. I greatly acknowledge the valuable discussion with George Haller, who was my ocal supervisor at MIT. I would also like to thank Gábor Orosz, James England, Kirk Green from the University of Bristol and Amit Surana, Sabri M. Kilic from MIT. Eusebius Doedel suggested an (undocumented) trick how to obtain starting point data for bifurcation continuation algorithm. I am also indebted to Barnabás Garay and Tamás Insperger who read the draft of this thesis thor- oughly and helped to improve the clearity of the text.

(8)

1 Introduction

Stability and nonlinear vibrations of machining processes are important and non- trivial areas of research, with a history dating back to the early 20th century [81].

Signicant advances however, only started in the 1950s, when it was realized that machining can only be modelled accurately by accounting the so-called regenerative eect. Specically, the tool cuts an already modulated surface produced in the previous cutting period, and hence the chip thickness depends on the previous positions of the tool tip as well as on the current one. As a result, machining must be modeled by delay-dierential equations. The simplest kind of machining, from dynamical point of view, is the turning process, where the slowly moving tool cuts the rotating cylindrical workpiece (see Fig. 1.1). In this case the chip-thickness during steady cutting is constant, waviness on the chip can only be induced by relative vibrations between the tool and the workpiece. Hence, mathematical models of turning are autonomous delay-dierential equations and their stability can be studied analytically [63]. Analyses by Tobias [83], Tlusty [82] and many unnamed authors revealed that a two-parameter stability chart of turning with a lobe-like stability boundaries marks the occurrence of Hopf-bifurcations. This results in self excited vibrations.

Besides the regenerative eect, nonlinearities in the cutting force or in the geom- etry of the process have major contributing eects to the appearance of undesired vibrations. Although, the interaction between the tool and the workpiece is a widely studied area of continuum mechanics (see Shaw [61] and references therein), the dynamics of chip formation is not well understood in the general three dimensional case. However, in case of orthogonal cutting, which is two dimensional in essence, there are some preliminary results on the dynamics of chip formation (see Davies and Burns [17] and Burns and Davies [7]), but these models are dicult to include in machine-tool vibration models. Instead of complicated models, we use the empirical three-quarter formula, which includes non-physical parameters and is suf- ciently precise even for quantitative analysis [82]. This approximation is used by

(12)

1 Introduction 2

Fig. 1.1. Scheme of the turning process

several authors including Davies and Balachandran [16], Paris et al. [57] or Peigne et al. [59] and also widely used in industry.

The strongest nonlinearity in the process is the vanishing of cutting force when the tool leaves the workpiece. This results a non-dierentiable break in the cutting force function and also leads to interruptions in the cutting causing rough surfaces on the workpiece. This self-interruption limits the amplitude of vibration by switching the dynamics to a linear damped oscillator which brings the tool back into cutting. Modeling this phenomena, we include additional delays into the equation of motion to account for cutting of surfaces which left uncut in previous periods. Several simulations in the literature ignore this increase in number and length of delays (see e.g. [53] or [57]). Although ignoring this eect does not inuence stability, it introduces great errors in predicting nonlinear vibrations.

So far, only the major contributory eects to machine tool chatter are mentioned, but there are others, which may change the characteristics of machining greatly. Abele and Fiedler [1] studied the stability of milling processes and found that stability lobes may shift to lower speeds. They argued that the stiness of the spindle is changing with angular speed, which follows from the dynamics of angular contact ball bearings frequently used in spindles. After including the decrease of stiness of the spindle into their model the stability lobes moved to their experimentally predicted positions. Another eect, which appears only at low-speeds, is the process damping related to the ank wear or ploughing of the tool. This phenomenon is explained by Tlusty [82] that the tool at low speeds may produce steep slopes in the workpiece, so the ank surface of the tool contacts the workpiece causing a damping eect (because the force is in counter-phase with the velocity).

Including the process damping into the models results higher stability boundaries at lower speeds. An alternative, but more physical, expalantion of the higher sta-

(13)

1.1 Outline of the thesis 3 bility lobes by Stépán [64] introduces a short time delay in the model. The short distributed delay comes from the force distribution on the rake face of the tool as the chip slides on it. Also, the parameters of the cutting force model can change with speed, which was examined by Faassen et al. [29].

The result of undesired vibrations is a rough surface on the workpiece [57], which has determining inuence on how parts t together and also on the fatigue stiness of the workpiece. Other inaccuracies, such as surface location error occur even when there is no instability [53, 42], but they may be compensated by control. In this thesis we emphasize that unwanted vibrations can arise even when the process is linearly stable, and nd parameter regions where this can happen. We show that this is related to subcritical bifurcations [79] and the persistence of high amplitude orbits in the stable part of the stability chart.

1.1 Outline of the thesis

In the next chapter we introduce a discrete-time nonlinear model of high-speed milling, which is a good approximation in case of low immersion and low number of cutting edges. This approximation makes possible an extensive bifurcation analysis of the process. As a result of the investigation, we conclude that the bifurcations are subcritical.

In chapter 3 we further analyze our discrete-time model to nd hints about its global dynamics. It turns out that in the parameter region of period doubling instabilities there are `global' period-two orbits, which determine the global dynamics. These orbits correspond to vibrations where the tool fails to cut in every second period. The stability and bifurcation analysis of this motion gives that it loses its stability and bifurcates in the same way as the stationary cutting. When both the global period-two orbit and the stationary cutting are unstable, chaos may arise. We demonstrate the proof of the existence of a strange attractor at a special parameter, and conclude that chaotic behavior is characteristic to the system.

Although the discrete-time model possesses most of the dynamical phenomena of high-speed milling, it is not suitable for quantitative predictions. In chapter 4 we take into account that the time instance spent with cutting is not innitesimal. This renement of the model results a time-periodic delay dierential equation, which is used in the stability analysis. It turns out that contrary to the discrete-time model, period-doubling stability boundaries can be closed curves. The method developed in chapter 4 is applied to more general models of up- and down-milling. At the end of the chapter we show how stability charts change by varying the immersion and damping in the models of the above mentioned processes.

(14)

1.1 Outline of the thesis 4 Chapter 5 gives a mathematical result on stability of periodic delay-dierential equations. We prove that for each periodic delay-equation there is a characteristic matrix, which is equivalent to the monodromy operator in the sense of Kaashoek and Lunel [47]. This characteristic matrix construction will be very important for us in chapter 6 in constructing the dening systems of periodic solution bifurcations.

It is not possible to analyze bifurcations of the delay-dierential equation model of high-speed milling fully analytically, therefore we develop a method in chapter 6 for doing it numerically. Our technique is a continuation method which can follow branches of periodic solutions along one parameter. We also introduce a continuation technique for computing bifurcation curves (in most cases stability boundaries) along two parameters, which allows one to produce stability chart for periodic orbits other than the stationary cutting, too. After showing a few examples, we investigate the period-two orbits in our model, which arise at period doubling stability boundaries. We nd sub- and supercritical period doubling bifurcations which are conrmed experimentally.

In chapter 7 we present an analytic method for calculating the sense of period- doubling bifurcation. We use center manifold reduction to reduce our system to one dimension, and then calculate the bifurcation applying the same method as in chapter 2.

(15)

2 Discrete-time model of low immersion high-speed

milling

In this chapter we analyze the bifurcations occurring at the stability loss of high speed milling. We approximate the delay dierential equation model of the milling process by a discrete map in case of innitesimal cutting times. Specically, we consider cutting as an impact. Our model will be a generalization of a model by Davies et al. [18] who gave a linear discrete time model to analyze stability in case of highly interrupted cutting. Our main theoretical tools are the center manifold reduction and the normal form reduction techniques.

2.1 Mechanical model

Fig. 2.1. Scheme of high-speed milling. Feed is provided by the workpiece velocity v0, cutting speed is provided by the (rotating) tool.

We use the simplest possible mechanical model of the process. The tool is modeled as a 1 degree of freedom (DOF) oscillator with undamped natural (angular) frequency ωn =p

k/m, relative damping factorζ =c/(2mωn), and damped natu-

(16)

2.1 Mechanical model 6

Fig. 2.2. Mechanical model. Note the dierence from the model in Fig. 2.1: the feed is provided by the tool while the cutting speed is provided by the motion of the (rotating) workpiece.

ral frequency ω_d =ω_np

1−ζ², where k is the stiness and c denotes the damping factor. The tool rotates with the constant angular speed Ω and has z number of edges. In the case of low immersion high-speed milling, the tool cuts only during a small fragment of the cutting period, so in Fig. 2.1, the workpiece is considered to be thin in the cutting direction and cutting can be approximated as an impact.

Consequently, the motion of the tool can be separated into two parts. As it can be seen in Fig. 2.2, the tool oscillates freely for time period τ1 = τ −τ2 and then cuts the workpiece for time period τ2. The tool starts the free vibration at time instants tj = t0 +jτ, j ∈ Z, enters the workpiece at t⁻_j+1 = tj+1−τ2 and nishes cutting at tj+1 when it starts free-ight, again. For the entire dynamics, we have the equation of motion

¨

x(t) + 2ζωnx(t) +˙ ω_n²x(t) = g(t)

m Fc(h(t)), (2.1)

where

g(t) =

(0 if ∃j ∈Z: tj 6t < t⁻_j+1 1 if ∃j ∈Z: t⁻_j+1 6t < tj+1

.

The cutting force Fc, in general, may depend on many factors of the process like cutting speed, lubrication, etc.. In our model, we consider its dependence merely on geometrical data, namely on the constant chip widthw and on the time-varying chip thickness h(t) in the form of the empirical three-quarter rule [82]

Fc(h(t)) =Kw(h(t))^3/4, (2.2) where K is an experimentally determined parameter. The chip thickness is computed from the current and a delayed tool tip position, i.e.,

(17)

2.1 Mechanical model 7 h(t) =h0+x(t−τ)−x(t),

where h₀ =v₀τ is the feed for a cutting period (see Fig. 2.1).

Introducing the dimensionless time tˆ= ωnt, the governing equation (2.1) becomes

¨

x(ˆt) + 2ζx(ˆ˙ t) +x(ˆt) =g(ˆt)Kw

mω²_n(h0+x(ˆt−ˆτ)−x(ˆt))^3/4.

In what follows by abuse of notation we drop the hat from the dimensionless time variables tˆ, τˆ, ˆτ1, τˆ2. For the free-ight period (g(t)≡0) we have

¨

x(t) + 2ζx(t) +˙ x(t) = 0 t⁻_j 6t < tj, which can be solved and arranged in discrete form:

Ãx(t⁻_j+1)

˙ x(t⁻_j+1)

!

=A

Ãx(t_j)

˙ x(t_j)

!

. (2.3)

When the time period of cutting is very short, i.e., τ2 →0then the time period τ1

of the free-ight can be approximated by the tooth-pass period τ. In this way the coecient matrixA constructed from the linear solution of the free-ight assumes the form

A=



e^−ζτ ³

cos(ˆω_dτ) + _ω_ˆ^ζ

dsin (ˆω_dτ)´

e^−ζτ ˆ

ωd sin (ˆω_dτ)

−^e^−ζτ_ω_ˆ_d sin (ˆωdτ) e^−ζτ

³

cos(ˆωdτ)− _ω_ˆ^ζ_d sin (ˆωdτ)

´



, (2.4)

where the dimensionless damped natural frequency is ˆ

ωd =ωd/ωn =p

1−ζ².

For the cutting period τ2 → 0, we neglect all the forces (spring, damping) except the cutting force, and assume that the position of the tool does not change much during the impact

¨

x(t)≈ Kw

mω_n²(h₀ +x(t_j)−x(t⁻_j+1))^3/4, t⁻_j+1 6t < t_j+1. Integrating the above formula on [t⁻_j+1, tj+1] we nd

˙

x(tj+1) = ˙x(t⁻_j+1) +τ2 Kw

mω_n²(h0+x(tj)−x(t⁻_j+1))^3/4. (2.5) Putting together equations (2.3) and (2.5) yields

(18)

2.3 Stability chart 8 x(tj+1) =A11x(tj) +A12x(t˙ j)

˙

x(tj+1) =A21x(tj) +A22x(t˙ j) + Kwτ2

mω²_n (h0+ (1−A11)x(tj)−A12x(t˙ j))^3/4, where Aij denote the corresponding elements ofA in (2.4).

2.2 Nonlinear discrete map

In the previous section, we constructed the equations mapping the state of system from cutting period to subsequent cutting period. We rephrase it with the state variables xj =x(tj) and vj = ˙x(tj), so the map becomes

Ãxj+1

vj+1

!

=A Ãxj

vj

! +

Ã 0

Kwτ2

mω²_n (h0+ (1−A11)xj −A12vj)^3/4

!

. (2.6)

This map has a xed point Ã

xe

ve

!

= Kwτ2h^3/4₀ mω_n²(1 + detA−trA)

Ã A12

1−A11

!

(2.7) that corresponds to stationary cutting, which is a period-1 motion with period τ. Linearizing around this xed point, we get the local dynamics

Ã xj+1

v_j+1

!

= Ã

0

4 3whˆ ₀

! +B

Ã xj

v_j

!

, B= Ã

A11 A12

A₂₁+ ˆw(1−A₁₁) A₂₂−wAˆ ₁₂

! , (2.8) where the dimensionless chip width

ˆ

w= 3 4h^1/4₀

Kτ2

mω_n²w

is assumed later as the bifurcation parameter. Note, that in the limiting case of τ2 → 0 and τ1 → τ, the original chip width w tends to innity at nite wˆ, which means that the extremely (highly) interrupted cutting could theoretically take place with innite chip width.

2.3 Stability chart

This section summarizes the linear stability results of [18]. In our model (2.6) only ip or Neimark-Sacker bifurcation can occur. For the ip case we have a characteristic multiplier at −1, that is

(19)

2.3 Stability chart 9 det(B( ˆw_cr^f ) +I) = 0.

This equation can be solved for wˆ in the form ˆ

w_cr^f = detA+ trA+ 1 2A12 = ˆωd

cos(ˆω_dτ) + cosh(ζτ)

sin(ˆωdτ) . (2.9) For the Neimark-Sacker case we have

detB( ˆw_cr^ns) = 1,

because B is real and its complex eigenvalues are conjugate pairs. Solving this equation for wˆ yields

ˆ

w^ns_cr = detA−1

A12 =−2ˆωdsinh(ζτ)

sin(ˆωdτ). (2.10)

Fig. 2.3. Stability chart. Grey regions are stable, continuous lines denote period doubling boundaries while dashed lines correspond to Neimark-Sacker bifurcation

The results are shown in Figure 2.3. This so-called stability chart is constructed in the plane of the most important technological parameters. One of these is the cutting speed represented by the product zΩˆ = 2π/τ of the number z of cutting edges and dimensionless angular speedΩˆ, while the other parameter is the dimensionless chip width wˆ.

2.3.1 Lower bounds of the stability lobes

In order to be able to compare this stability chart with the others in chapter 4, we compute the tendency of the minima of lobes. Because sin(ˆω_dτ) ≤ 1 we can estimate the lower bound of wˆ^ns_cr by

(20)

2.3 Stability chart 10 ˆ

w^ns_cr ≥2ˆωdsinh(ζτ)

The lower bound for the period-doubling boundaries can be written as ˆ

w^f_cr≥ωˆdsinh(ζτ).

To prove this, we rewrite the above formula with exponentials in the form using (2.9)

ˆ

ωdcos(ˆωdτ) + cosh(ζτ)

sin(ˆωdτ) ≥sinh(ζτ)

Because in case of period doubling sin(ˆωdτ)>0 we can write that cos(ˆω_dτ) + cosh(ζτ)−sin(ˆω_dτ) sinh(ζτ)>0, which reduces further to

cosh(ζτ) cos(ˆωdτ +ε) + cosh(ζτ)≥0,

whereε is some phase angle. Because cosh(ζτ)>0for all ζτ ∈R, we are left with cos(ˆωdτ +ε) + 1≥0,

which proves our hypothesis. These lower boundary curves can be seen in the stability chart in Fig. 2.4.

Fig. 2.4. Structure of stability chart for highly interrupted cutting with bounds of the stability lobes

(21)

2.4 Local bifurcations 11

2.4 Local bifurcations

In this section, we show that bifurcations occurring at stability boundaries (2.9) and (2.10) are both subcritical. In order to do this, we use standard techniques of local bifurcation theory, which can be found, for example, in [35]. To apply these techniques, we need the Taylor expansion of this map (2.6) at its xed point (2.7).

This requires the expansion of the cutting force function Fc in (2.2) at theoretical chip thickness h0 with respect to the chip thickness variation∆h

Fc(h0+∆h)≈Kwh^3/4₀ + 3Kw

4h^1/4₀ ∆h− 3Kw

32h^5/4₀ (∆h)²+ 5Kw

128h^9/4₀ (∆h)³ or with the dimensionless chip width wˆ

F_c(h₀+∆h) τ₂

mω_n² ≈ 4 ˆw

3 h₀+ ˆw∆h− wˆ 8h0

(∆h)²+ 5 ˆw

96h²₀(∆h)³.

In what follows, we considerwˆas the bifurcation parameter, thus the system around the xed point is now approximated by

Ã xj+1

vj+1

!

= Ã

xe

ve

!

+B( ˆw) Ã

xj

vj

! +

Ã

0 ˆ

wP

q+r=2,3bqrx^q_jv_j^r

!

, (2.11)

where

b₂₀ =−(1−A₁₁)² 8h0

, b₁₁ = (1−A₁₁)A₁₂ 4h0

, b₀₂ =−A²₁₂ 8h0

, b₃₀ = 5(1−A11)³

96h²₀ , b₂₁ =−5(1−A11)²A12

32h²₀ , b₁₂ = 5(1−A11)A²₁₂

32h²₀ , b₀₃ =−5A³₁₂ 96h²₀. 2.4.1 Flip bifurcation

Crossing the stability boundaries at (2.9) a ip bifurcation occurs. At the stability boundaries the critical eigenvalue is

λ1( ˆw^f_cr) =−1

and the linear part of the system (2.11) can be written as

T⁻¹B( ˆw^f_cr)T =

Ã−1 0 0 λ2

! ,

(22)

2.4 Local bifurcations 12 where T is a transformation matrix and contains the eigenvectors s1,2 of B( ˆw^f_cr), i.e.,

T := (s₁,s₂) =

Ã 1 1

−ζ− ^−ˆ^ω^d(^e^ζτ^+cos(ˆ^ωdτ))

sin(ˆωdτ) −ζ+ ^−ˆ^ω_sin(ˆ^d^sinh(ζτ_ω_d_τ) ⁾

!

and

λ₂( ˆw^f_cr) = e^−ζτ(cos(ˆω_dτ) + sinh(ζτ)). (2.12) We do a local analysis in the neighborhood of X^f_cr = ((xe, ve),wˆ_cr^f ) ∈ R² × R. We consider the following perturbation of the system (2.11) around X^f_cr in the coordinate system of the eigenvectors of B( ˆw^f_cr)

Ãξj+1

ηj+1

!

=

Ã−1 +a^f∆wˆ 0 0 C22

! Ãξj

ηj

! +

Ã P

q+r=2,3cqrξ^q_jη_j^r P

q+r=2,3dqrξ_j^qη_j^r

!

, (2.13)

where a^f is the derivative ofλ₁( ˆw) at wˆ^f_cr a^f = ∂λ1

∂wˆ

¯¯

ˆ w^f_cr

=− 2 ˆ ω_d

sin(ˆωdτ)

cos(ˆω_dτ) + cosh(ζτ) + 2 sinh(ζτ) <0

and the nonlinear coecients c_qr, d_qr are computed from b_qr with the help of the transformation matrixT.

From the center manifold theorem we know that there exists an invariant manifold containingX^f_cr that is tangent to the eigenvector corresponding to the critical eigenvalue −1. In the present case, the center manifold is attracting since|λ2|<1 in (2.12) and its graph can be written in the power series [70] (see also Appendix A Section A.2.1)

ξ7→col(ξ, h(ξ)), h(ξ) =h2ξ²+. . . , h2 = d20

1−λ₂.

The projection of (2.13) into the center manifold gives the scalar nonlinear map ξj+1 = (−1 +a^f∆w)ξˆ j+c20ξ_j²+ (c30+ c11d20

1−λ2)ξ_j³+· · · . (2.14) We seek period-2 orbits that arise in the neighborhood of the xed point, therefore consider the second iterate of (2.14) having the form

ξj+2 = (1−2a^f∆w)ξˆ j−2δξ³_j +· · · , (2.15) where

δ =c²₂₀+c30+ c11d20

1−λ2.

(23)

2.4 Local bifurcations 13 The period-2 points on the center manifold can be found as the xed point of (2.15), namely

ξ1,2 =± r

−a^f∆wˆ δ .

Thus, the sense of the ip bifurcation is determined by the sign of δ. A tedious algebraic calculation yields

δ=− 5 12h²₀

cosh(ζτ) + cos(ˆω_dτ)

cosh(ζτ) + 2 sinh(ζτ) + cos(ˆωdτ) <0, (2.16) which shows that the arising bifurcation is always subcritical. Unstable period-2 motions exist around the stable xed point near the stability boundary.

Fig. 2.5. Simulation and analytical results for Eqn. (2.6) (wˆ= 0.827,zΩˆ = 0.78)

The simulation in Fig. 2.5 shows the saddle type period-2 motion represented by (ξ1, h(ξ1)) and(ξ2, h(ξ2)). The basin of attraction of the stable xed point (xe, ve) is approximated by the two straight lines parallel to s2. Starting an iteration at (x0, v0) from between these parallel lines it converges to the xed point. It can be seen that the rate of the convergence in the coordinate direction η is very fast, while in the orientation reversing direction ξ is very slow since the corresponding multiplier is slightly greater than −1.

2.4.2 Neimark-Sacker bifurcation

Along the stability boundaries (2.10), a Neimark-Sacker (also called secondary Hopf) bifurcation occurs, which can be related to the Hopf bifurcation of turning processes. At this type of bifurcation quasiperiodic orbits arise, which are living on an invariant closed curve about the xed point in the case of a discrete system, or on an invariant 2 dimensional torus in the case of a vector eld. The sense of the bifurcation, i.e. whether stable or unstable motion arise, can be calculated using

(24)

2.4 Local bifurcations 14 the normal form reduction [35] (see also appendix A section A.2.2). For the sake of simplicity we use complex transformation here. With the help of the complex eigenvectors s1,2 belonging to the critical eigenvalues of B( ˆw_cr^ns)

λ1,2 = e^±iφ, φ= cos⁻¹ trB 2 ,

we construct the complex transformation matrix T= (s1,s2) in the form

T=

Ã 1 1

−ζ−ωˆdcot(ˆωdτ) + ^ω^ˆ_sin(ˆ^d^e^ζτ+iφ_ω_d_τ₎ −ζ −ωˆdcot(ˆωdτ) + ^ω^ˆ_sin(ˆ^d^e^ζτ−iφ_ω_d_τ)

! .

Similar to the ip case, the following approximation of the system (2.11) can be derived at X^ns_cr = ((x_e, v_e),wˆ^ns_cr) ∈ R²×R in the coordinate system of the above critical eigenvectors of B( ˆw^ns_cr)

Ãzj+1

zj+1

!

= (1 +a^ns∆w)ˆ

Ãe^iφ 0 0 e^−iφ

! Ãzj

zj

! +

Ã P

q+r=2,3cqrz_j^qz^r_j P

q+r=2,3dqrz_j^qz^r_j

!

, (2.17)

where

a^ns= ∂|λ1,2|

∂wˆ

¯¯

ˆ w_cr^ns

= ∂

∂wˆ

pdetB( ˆw_cr^ns) =− 1

2ˆωde^−ζτsin(ˆωdτ)>0.

We can transform out all2^nd degree terms from (2.17), but for the3^rddegree terms we have the resonance related toλ1 =λ²₁λ2 andλ2 =λ1λ²₂ (see appendix A section A.2.2). Thus the normal form yields

Ãzj+1

¯ zj+1

!

=

Ãe^iφ 0 0 e^−iφ

! Ãzj

¯ zj

! +

Ãe21z_j²z¯j

f12zjz¯_j²

!

, (2.18)

where

e21 =f12 = 2 c20c11

1−e^iφ + c11c11

1−e^−iφ + c11c20

e^2iφ−e^iφ + 2 c02c02

e^2iφ−e^−iφ +c21.

By multiplying the two coordinates of (2.18) and using |z|² = zz¯ we transform (2.18) into

|z|² 7→(1 +a^ns∆w)ˆ ²|z|²+ 2(1 +a^ns∆w)δ|zˆ |⁴+e21f12|z|⁶, where

δ = 1

2(e21e^−iφ+f12e^iφ) = Re(e21e^−iφ).

(25)

2.5 New results 15

Fig. 2.6. Simulation showing the computed (dashed line) and true repelling invariant manifold. (wˆ = 0.503,zΩˆ = 0.55)

Neglecting the6^th degree term the solution for the radius|z| of the invariant circle assumes the form

|z|= s

−2a^ns∆wˆ+ (a^ns)²∆wˆ² 2(1 +a^ns∆w)δˆ ≈

r

−a^ns∆wˆ δ .

This circle in the discrete space of(xj, vj)corresponds to a quasiperiodic oscillation of our original mechanical structure during the milling process.

After a long algebraic calculation we obtain δ = e^−5ζτ¡

4e^4ζτ −3e^2ζτ −1¢

(cosh(ζτ)−cos(ˆω_dτ))

32h²₀ >0. (2.19)

Since δ is always positive, we can conclude that the Neimark-Sacker bifurcation is subcritical, too. The consequence for the machining process is similar to the ip case: the basin of attraction of the stationary cutting is inside of an invariant ellipse in the phase space. Thus, if we perturb the system such that we leave the basin of attraction, the system starts large amplitude (initially quasiperiodic) vibrations.

This is called chatter.

The dotted simulations in Fig. 2.6 diverge from the invariant ellipse (continuous line) and its approximated counterpart (dashed line), which represents the unstable quasiperiodic oscillation. We used the inverse map in the simulation to visualize this invariant unstable ellipse.

2.5 New results

Thesis 1 Low immersion high-speed milling processes can lose their stabilty either by period-doubling or by Neimark-Sacker bifurcation. We proved that in the

(26)

2.5 New results 16 nonlinear dicrete time model of this process both bifurcations are subcritical, that is, unstable vibrations coexist with the stable cutting for parameters close to the stability boundary.

(27)

3 Global dynamics of the discrete-time model

3.1 Period two motion with `y-overs'

In the previous chapter we have shown that there is only one xed point (stationary cutting) of the system and if it undergoes a period doubling bifurcation, the resulting orbits may reach the boundary, where the actual chip thicknesshbecomes negative, the tool misses the workpiece and the cutting force vanishes (see Fig. 2.2).

Hence, we suspect that the corresponding orbits continue beyond this boundary. In order to nd these orbits we must modify the governing equations in a way that the tool ies over the workpiece in every second period. The equation for this motion will contain 2 times longer free-ights and impacts with the workpiece with double feed (2h0). Thus, the governing equations become

Ãx_j+2 v_j+2

!

=A⁺ Ãx_j

v_j

! +

Ã 0

Kwτ2

mω²_n

¡2h0+ (1−A⁺₁₁)xj−A⁺₁₂vj

¢_3/4

!

, (3.1)

where A⁺ = A² and A⁺_ij are the elements of A⁺. Note, that A⁺ can be obtained from Ain (2.4) if 2τ is substituted instead of τ. This equation has a unique xed point, again, but this time it physically corresponds to a period-2 motion

Ãx⁺_e v⁺_e

!

= Kwτ₂(2h₀)^3/4 mω_n²(1 + detA⁺−trA⁺)

Ã A⁺₁₂ 1−A⁺₁₁

! .

As explained above, this period-2 motion exists if the tool does not hit the workpiece after the rst period of free-ight. This condition can be checked by the following equation based on the argument of the nonlinear term in (2.6)

0> h(t⁻_j+1) =h0+ (1−A11)x⁺_e −A12v_e⁺ = h0+ Kwτ2

mω²_n (2h0)^3/4

µ(1−A11)A⁺₁₂ −(1−A⁺₁₁)A12

1 + detA⁺−trA⁺

¶

. (3.2)

(28)

3.1 Period two motion with `y-overs' 18 After substituting the matrix elements and the dimensionless parameters we are left with

ˆ

w > 3ˆωd

2^7/4

cos(ˆωdτ) + cosh(ζτ)

sin(ˆωdτ) . (3.3)

Because of the similarity of the governing equations (2.6) and (3.1) the stability boundaries of (x⁺_e, v_e⁺) can be obtained from the same calculation just by changing the matrix A to A⁺ and h0 to 2h0. With the dimensionless parameters this calculation gives the stability boundaries of the period-2 motion in the form

ˆ

w^f+_cr = 2^1/4detA⁺+ trA⁺+ 1

2A⁺₁₂ = 2^1/4ωdcos(2ωdτ1) + cosh(2ζτ1)

sin(2ω_dτ₁) (3.4) and

ˆ

w^ns+_cr = 2^1/4detA⁺−1

A⁺₁₂ =−2^7/4ωdsinh(2ζτ1)

sin(2ωdτ1). (3.5) The same can be said about the bifurcations along these stability boundaries: both bifurcations are subcritical, because replacingτ by2τ andh₀ by2h₀ in Eqns. (2.16) and (2.19) does not change the signs.

Fig. 3.1. Global period-2 motions and their stability

A characteristic region of the corresponding stability chart is shown in Fig.

3.1. The stationary cutting (xe, ve) exists for all the parameters, it is stable below the continuous line representing a ip lobe of the stability chart in Fig. 2.3, and it is unstable above it. The period-2 motion (x⁺_e, v_e⁺) exists above the dashed line according to formula (3.3). In the light gray region only the stationary cutting exists and it is stable. The dash-dotted line `ns' refers to a stability boundary (3.5) of the period-2 motion, where Neimark-Sacker bifurcation occurs, while the dash-dotted line `f' refers to the other stability limit (3.4) where the period-2 motion undergoes another ip bifurcation. Thus, the dark gray region represents parameters, where

(29)

3.1 Period two motion with `y-overs' 19 stable period-2 motion exists, while in the white regions the outer period-2 motion is unstable.

-4 -3 -2 -1 0 1 2 3 4 5 6

0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 x fixed point

stable period-2 point (a)

w^

-0.5 0 0.5 1 1.5 2 2.5 3

0.6 0.8 1.0 1.2 1.4 1.6 1.8

x

fixed point

stable period-2 point chaos

unstable period-2 point (b)

w^{^}

0.5 1 1.5 2 2.5 3

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

x

fixed pointunstable period-2 orbits

chaos (c)

w

^

Fig. 3.2. Bifurcation diagrams for cases a), b), c) in Fig. 3.1. The outer period-2 solution is stable and undergoes a Neimark-Sacker bifurcation (a) or a ip bifurcation (b). The period-2 solution is unstable, which induces chaotic motions around it (c).

To represent the above explained properties of the stability chart in Fig. 3.1, bifurcation diagrams are drawn for the parameter cases a), b), c) in Fig. 3.2. The unstable period doubling branches are computed using AUTO [23] that followed our analytic predictions perfectly. These computations also show that the inner

(30)

3.2 Chaotic oscillation 20 and outer period-2 orbits are connected through a degenerate fold bifurcation, which means that orbits do not change smoothly by varying parameters and their characteristic multipliers are not even continuous.

3.2 Chaotic oscillation

As we have seen in the bifurcation diagrams in Fig. 3.2, there are parameter regions where chaotic motions arise in the system. In this section, we investigate the simplest chaotic case, where both the xed point and the outer period-2 orbit are unstable for case c) and w >ˆ wˆ^f_cr. Both of the corresponding periodic points ((xe, ve), (x⁺_e, v_e⁺)) are saddle-like, each has an unstable orientation reversing subspace as well as a stable orientation preserving subspace. Unfortunately, the local invariant manifolds are not to be extended to global manifolds because of the piecewise structure of the system, consequently, we have just global invariant sets, which are the union of images (preimages) of the local unstable (stable) manifold of the periodic points.

0 0.5 1 1.5 2 2.5 3

0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4

x/h h

0

h F

c

P

₂

P1

W

^s

P1

W

^u

P2

W

^u

P2

W

^s

P

₁

switching line

1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3

0.6 0.8 1 1.2 1.4 1.6 1.8

P2

W

^s

P2

W

^u

P1

W

^s

P1

W

^u

P

₁

P

₂

switching line

x/h

(a) (b)

A B

C

D E

H₁

H₀ V₁

V₀

H₀ H₁

V₁ horseshoe

V₀

0

0 0

v/h

0

Fig. 3.3. The chaotic map (zΩˆ = 0.459,wˆ= 1.64); a) simulation with invariant manifolds and how cutting force varies perpendicular to the switching line; b) horizontal and vertical slabs (see in the text).

(31)

3.2 Chaotic oscillation 21 In Fig. 3.3.a the phase plane is presented for wˆ = 1.64. P1 = (xe, ve) is the xed point representing stationary cutting, P2 = (x⁺_e, v⁺_e) refers to the period-2 motion. The invariant manifoldsW_P^s,u_1,2 of these saddle points were computed using the software dstool module of England et al. [28]. The two parts of the dynamics are separated by the switching line h(t⁻_j+1) = 0 in (3.2), where the cutting force F_c is non-smooth. The cutting force characteristic of Fig. 2.2 is projected onto this plane to visualize the locations of the switching line and the xed pointP₁, where the chip thickness is 0 and h₀, respectively. The gure also shows simulation of a chaotic iteration. Below the switching line, we iterate equation (2.6), above that, (3.1). Since (3.1) describes every second τ-period, we use matrix Ato produce the system's state in the middle of the period-2 free-ight. These additional iterated points show up among other points in the lower part of the gure, below W_P^s₁.

Demonstrating the proof of chaos, we use a geometric method in [84, 56]. In fact, our system is conjugate to a modied version of Smale's horseshoe map. The construction of the map can be seen in Fig. 3.3.b. Let the map f be dened on this phase plane by (3.1) above the switching line and by every second iterate of (2.6) otherwise. Point A is dened as the intersection of W_Pû₂ and the switching line, then B:=f(A), C :=f(B), D:=f(C). Note A is part of both dynamics and as the entire switching line will be mapped by (3.1) and by that second iterate of (2.6) to the same place. Further, letE be the intersection ofW_Pû₁ and the switching line. According to numerical results, the stable manifolds W_P^s_1,2 do not intersect the switching line between A and E. Choose a horizontal slab H1 along the stable manifold W_P^s₂ from P2 until W_Pû₁ to satisfy the following conditions. For a suitable integer k >1, f^k−1(H1) must lie above the switching line such that it contains B and a point between P2 and A on W_Pû₂. Then V1 := f^k(H1) contains C and fully intersects H₁. The thickness of V₁ can be adjusted by an appropriate choice of k together with the vertical size ofH₁. The other slab H₀ is chosen to lie along W_Pû₂ from point A to point C and be thick enough to be fully intersected by V₁. Now we can observe that V0 := f(H0) fully intersects H1, if D is underneath W_P^s₂. In this case we have topological conjugacy with the schematic inset of Fig. 3.3.b.

The calculated manifolds and the iterated sets do indeed satisfy the above conditions, so we can describe the dynamics by means of the horseshoe structure of the inset, that is, by means of the left shift on the space of innite sequences of two symbols equipped with the transition matrix

Ã0 1 1 1

! .

This matrix is irreducible [84, 48], so the motion of the system is indeed chaotic.

(32)

3.4 New results 22

3.3 Discussion

Although our model is very simple the results are qualitatively the same in the case of non-innitesimal cutting time as found by simulations [3] and experiments [54]. In chapter 6 and 7 we will show that period doubling bifurcations in the delay equation model are subcritical and outer period-2 orbits are connected to the unstable orbits through a fold bifurcation (see also [75] and [80]). The demonstration of chaotic motion can be more dicult in the delay equation model because of the innite dimensional phase space. In this case, instead of our purely geometric method of constructing the chaotic attractor one could compute 1 dimensional unstable manifolds of periodic orbits to obtain some information of the structure of chaotic attractor (for a guiding example see Green et al. [34]).

3.4 New results

Thesis 2 A condition for the existence of an additional period-2 orbit is given in the case of the discrete-time model of high-speed milling. This orbit corresponds to a motion where the tool misses the workpiece in every second cutting period. We showed that this solution can lose its stability in the same way as the stationary cutting and these stability losses are similarly subcritical.

If both the period-two solution and the stationary cutting are unstable with characteristic multipliers less than−1, the global dynamics becomes chaotic. We demonstrated the proof of this chaotic motion for a characteristic set of parameters.

Nonlinear dynamics of high-speed milling

Róbert Szalai