You have access

Check for updates on crossmark

Research article

On the analysis of the double Hopf bifurcation in machining processes via centre manifold reduction

T. G. Molnar

http://orcid.org/0000-0002-9379-7121

Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

[email protected]

Find this author on PubMed

Search for more papers by this author

,

Z. Dombovari

Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

Find this author on PubMed

Search for more papers by this author

,

T. Insperger

MTA-BME Lendület Human Balancing Research Group, Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

Find this author on PubMed

Search for more papers by this author

and

G. Stepan

Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

Find this author on PubMed

Search for more papers by this author

T. G. Molnar

http://orcid.org/0000-0002-9379-7121

Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

[email protected]

Find this author on PubMed

Search for more papers by this author

,

Z. Dombovari

Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

Find this author on PubMed

Search for more papers by this author

,

T. Insperger

MTA-BME Lendület Human Balancing Research Group, Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

Find this author on PubMed

Search for more papers by this author

and

G. Stepan

Department of Applied Mechanics, Budapest University of Technology and Economics, 1111 Budapest, Hungary

Find this author on PubMed

Search for more papers by this author

Published:15 November 2017https://doi.org/10.1098/rspa.2017.0502

Abstract

The single-degree-of-freedom model of orthogonal cutting is investigated to study machine tool vibrations in the vicinity of a double Hopf bifurcation point. Centre manifold reduction and normal form calculations are performed to investigate the long-term dynamics of the cutting process. The normal form of the four-dimensional centre subsystem is derived analytically, and the possible topologies in the infinite-dimensional phase space of the system are revealed. It is shown that bistable parameter regions exist where unstable periodic and, in certain cases, unstable quasi-periodic motions coexist with the equilibrium. Taking into account the non-smoothness caused by loss of contact between the tool and the workpiece, the boundary of the bistable region is also derived analytically. The results are verified by numerical continuation. The possibility of (transient) chaotic motions in the global non-smooth dynamics is shown.

1. Introduction

Improving the productivity of metal cutting operations is highly important for the manufacturing society. The occurrence of harmful vibrations during machining—known as machine tool chatter—limits the achievable productivity. Therefore, studying the dynamics of machine tool chatter is an active field of research.

The aim of analysing machine tool vibrations is to identify those regions in the space of the technological parameters, which are associated with chatter-free manufacturing processes. These parameter regions are typically illustrated in the plane of the spindle speed and the depth of cut; this picture is called stability lobe diagram or stability chart. The first successful efforts to describe chatter and to derive stability charts are presented in the works of Tobias [1] and Tlusty [2].

Machine tool vibrations are described by delay-differential equations (DDEs), since the chip thickness, which determines the cutting force, is affected both by actual and delayed tool positions due to the surface regeneration effect. Since DDEs have infinite-dimensional phase space representation [3], the dynamics of machine tool chatter is rich and intricate. Furthermore, the modelling DDEs are typically nonlinear, since the cutting force is a nonlinear function of the chip thickness [4].

Here, we restrict ourselves to the analysis of autonomous DDEs, which are the typical models of turning operations. Along the stability boundaries (or stability lobes) of turning, Hopf bifurcation occurs, which gives rise to a periodic solution of the nonlinear system. This periodic solution plays an important role in determining the global dynamics of the cutting process that may involve the high-frequency self-excited vibration called chatter. The effect of nonlinearities on the occurrence of machine tool chatter was investigated, e.g. in [5–8], whereas the bifurcation analysis of machining processes can be found in [9–14] for turning and in [15] for milling operations. The theory of Hopf bifurcation in DDEs is covered by [3,16–19]. The most popular approaches for bifurcation analysis are the rigorous centre manifold reduction [20] and the method of multiple scales [12,13]. Note that other methods also exist for the computation of periodic solutions [21], e.g. the method of small parameters or the theory of averaging [5,22–26].

In this paper, we analyse the intersection of stability boundaries where double Hopf bifurcation takes place. The double Hopf bifurcation complicates the dynamics of metal cutting by giving rise to a quasi-periodic solution [27]. Examples for the analysis of double Hopf bifurcation can be found in [28–38] for various dynamical systems. The double Hopf bifurcation in metal cutting was analysed in [39] using the method of multiple scales. Now we extend this work by analysing a slightly different model using centre manifold reduction. In particular, we investigate the possible topologies in the phase portraits of the system near the double Hopf bifurcation and we also discuss the phenomenon of bistability. We verify the results by numerical continuation. Continuation of periodic solutions of nonlinear DDEs can be carried out by the software DDE-Biftool [40], whereas the arising quasi-periodic solution can be computed by the software Knut [41,42] or by the algorithm used in [15]. In this paper, we use the latter algorithm to verify the results of centre manifold reduction.

The rest of the paper is organized as follows. Section 2 presents the mechanical model of turning. Section 3 discusses the linear stability analysis, whereas §4 shows the bifurcation analysis using centre manifold reduction. The resulting stability charts and bifurcation diagrams are presented together with implications on the global dynamic behaviour in §5. Conclusions are drawn in §6.

2. Equation of motion

We consider the single-degree-of-freedom model of orthogonal cutting shown in figure 1a. The equation of motion for the tool can be given in dimensionless form as follows (see [11] for details):

\begin{aligned} x^{″} (t) + 2 ζ x^{'} (t) + x (t) & = w ((x (t - τ) - x (t)) + η_{2} (x (t - τ) - x {(t))}^{2} + η_{3} (x (t - τ) - x {(t))}^{3}), \end{aligned}

2.1

where t is the dimensionless time measured by the time period of the free undamped natural oscillation of the system, x indicates the dimensionless position of the tool with respect to the constant theoretical chip thickness set to have unit length in figure 1, and ζ is the damping ratio of the dominant vibration mode of the system. The right-hand side of equation (2.1) is proportional to the variation of the dimensionless, x-directional cutting force component, which is assumed to be a cubic polynomial function of the dimensionless chip thickness variation x(t−τ)−x(t). Parameter τ is the dimensionless regenerative delay or, equivalently, the period of workpiece rotation related to the dimensionless angular velocity Ω of the workpiece: τ=2π/Ω. Parameter w is the dimensionless chip width, whereas η₂ and η₃ denote dimensionless quadratic and cubic cutting-force coefficients related to the prescribed feed per revolution, which is just the unit theoretical chip thickness.

Figure 1. — Figure 1. (a) Single-degree-of-freedom mechanical model of orthogonal cutting, (b) stability lobe diagram of the cutting process and (c) a typical bifurcation scenario. (Online version in colour.)

Equation (2.1) is a DDE with cubic nonlinearity. The cubic cutting force model was introduced in [5]. Note that other cutting force characteristics also exist, see [4] and the references therein, but these characteristics can also be approximated by cubic polynomials using Taylor series expansion. Approximation of the nonlinearity by a cubic expression is a necessary step for the subsequent bifurcation analysis in §§4 and 5.

Equation (2.1) is valid only for positive dimensionless chip thickness, h(t)=1+x(t−τ)−x(t)>0. For large-amplitude vibrations, it might occur that the dimensionless chip thickness h(t) drops to zero or to negative values, which indicates that the tool gets out of the workpiece and loses contact with the material. In such cases, the cutting force on the right-hand side of (2.1) becomes zero irrespective of the tool’s position. At this point, we do not investigate loss of contact between the tool and the workpiece and assume that h(t)>0 and (2.1) is valid during the entire cutting process.

In the rest of the paper, we deal with the stability and bifurcation analysis of (2.1). For convenience, the equation is transformed into the first-order form

y^{'} (t) = L y (t) + R y (t - τ) + g (y_{t}),

2.2

where y is the vector of state variables, L and R are the coefficient matrices of the linear and the retarded terms, respectively, and g contains all nonlinear terms:

\begin{aligned} y (t) & = [\begin{matrix} x (t) \\ x^{'} (t) \end{matrix}], L = [\begin{matrix} 0 & 1 \\ - (1 + w) & - 2 ζ \end{matrix}], R = [\begin{matrix} 0 & 0 \\ w & 0 \end{matrix}] \\ and g (y_{t}) & = [\begin{matrix} 0 \\ w (η_{2} (y_{1} (t - τ) - y_{1} {(t))}^{2} + η_{3} (y_{1} (t - τ) - y_{1} {(t))}^{3}) \end{matrix}] . \end{aligned}}

2.3

Here subscript 1 refers to the first component of a vector, i.e. y₁(t)=x(t).

Since DDEs have infinite-dimensional phase space [3], the state of the system is not determined solely by the vector y(t). Instead, we use a function defined in the Banach space $B$ of continuously differentiable vector-valued functions to describe the evolution of the system: we introduce the shift $y_{t} \in B : [- τ, 0] \to R^{2}$ , y_t(θ)=y(t+θ). Using y_t, system (2.1) can be represented in operator differential equation (OpDE) form as:

y_{t}^{'} = A y_{t} + F (y_{t}),

2.4

where

A, F : B \to B

are the linear and the nonlinear operators, respectively, defined by

(A u) (θ) = {\begin{cases} \frac{d}{d θ} u (θ) & if θ \in [- τ, 0), \\ L u (0) + R u (- τ) & if θ = 0 \end{cases}

2.5

and

(F (u)) (θ) = {\begin{cases} 0 & if θ \in [- τ, 0), \\ g (u) & if θ = 0. \end{cases}

2.6

3. Linear stability analysis

Equation (2.1) has a single equilibrium x(t)≡0, which corresponds to stationary cutting. Machine tool chatter is associated with the loss of stability of this equilibrium. Substitution of the exponential trial solution x(t)=Ce^λt, $λ, C \in C$ into (2.1) gives the characteristic function

D (λ) := det (λ I - L - R e^{- λ τ}) = λ^{2} + 2 ζ λ + 1 + w (1 - e^{- λ τ}) .

3.1

Its infinitely many zeros—known as characteristic exponents—are the same as the eigenvalues of the operator

A

in (2.5). The trivial solution is exponentially asymptotically stable if and only if all the characteristic exponents lie in the open left half of the complex plane. The system is at the boundary of stability if a single real root λ=0 or a pair of complex conjugate roots λ=±iω (i²=−1, ω>0) exists, while no characteristic exponents have positive real part. For (3.1), only the latter case with critical roots λ=±iω is possible, which corresponds to a Hopf bifurcation in the nonlinear system. The analysis of the Hopf bifurcation can be found in [9–11] in details. The Hopf bifurcation points can be located after the separation of D(±iω)=0 into real and imaginary parts:

\begin{aligned} R (ω) = & - ω^{2} + 1 + w (1 - \cos (ω τ)) \\ and S (ω) = & 2 ζ ω + w \sin (ω τ) . \end{aligned}}

3.2

Solving R(ω)=0 and S(ω)=0 gives the linear stability boundaries in the form

\begin{aligned} w_{H} (ω) & = \frac{{(ω^{2} - 1)}^{2} + 4 ζ^{2} ω^{2}}{2 (ω^{2} - 1)}, τ_{H} (j, ω) = \frac{2}{ω} (j π - \arctan (\frac{ω^{2} - 1}{2 ζ ω})) \\ and Ω_{H} (j, ω) & = \frac{2 π}{τ_{H} (j, ω)}, \end{aligned}}

3.3

where subscript H indicates that (3.3) gives the location of Hopf bifurcation. Note that (3.3) defines a family of curves called stability lobes where

j \in Z^{+}

is the lobe number. Depicting the stability boundaries in the plane (Ω,w) of the technological parameters leads to so-called stability lobe diagrams (see the example in figure 1b for ζ=0.02 where grey shading indicates the linearly stable region associated with chatter-free cutting process).

As it was shown in [9–11], subcritical Hopf bifurcation occurs along the stability lobes. In the vicinity of the stability boundaries, the bifurcation gives rise to an unstable periodic orbit around the linearly stable equilibrium with approximate period 2π/ω. At the intersection points of the stability lobes, a codimension two bifurcation—often referred to as double Hopf bifurcation—takes place. At these points, two pairs of characteristic exponents λ=±iω_1,2 lie on the imaginary axis (ω_1,2>0). Let us consider the intersection of the lobes with lobe numbers j₁ and j₂. The location of this point (the so-called double Hopf point) is given by w_HH:=w_H(ω₁)=w_H(ω₂), τ_HH:=τ_H(j₁,ω₁)=τ_H(j₂,ω₂) and Ω_HH:=Ω_H(j₁,ω₁)=Ω_H(j₂,ω₂), which can be obtained by solving R(ω₁)=0, S(ω₁)=0, R(ω₂)=0 and S(ω₂)=0 numerically for ω₁, ω₂, w and τ. Subscript HH refers to Hopf–Hopf (double Hopf) bifurcation.

The rest of the paper deals with the analysis of the double Hopf bifurcation via centre manifold reduction and the derivation of the normal form of the resulting four-dimensional centre subsystem.

4. Centre manifold reduction and normal form calculations

In this section, centre manifold reduction is used to analyse the long-term dynamics of the time-delay system (2.1) in the vicinity of the double Hopf bifurcation point. The theory of centre manifold reduction is discussed in [3]. The calculation follows the steps of the single Hopf bifurcation analysis, which is described in [9–11] in details.

At a double Hopf point on the boundary of the linearly stable region, the system has four characteristic exponents located on the imaginary axis, while all the other characteristic roots are located in the open left half-plane. Therefore, the long-term dynamics of the system is determined by the flow on a four-dimensional centre manifold, which is embedded in the infinite-dimensional phase space of (2.1). This flow can be analysed by decomposing the four-dimensional centre subsystem (see (4.24)). In what follows, we perform this decomposition using the theorem given by (3.10) and (3.11) in [3], ch. 7.

The decomposition theorem of [3] introduces the operator $A^{*}$ , which is formally adjoint to operator $A$ relative to a certain bilinear form. The formal adjoint $A^{*} : B^{*} \to B^{*}$ must satisfy $(v, A u) = (A^{*} v, u)$ for any pair of $u \in B : [- τ, 0] \to R^{2}$ and $v \in B^{*} : [0, τ] \to R^{2}$ , where $B^{*}$ is the adjoint space and operation $(,) : B^{*} \times B \to R$ indicates the bilinear form. The definition of the formal adjoint and the bilinear form can be found in [3], see (3.1) and (3.3) in ch. 7, and here they become

(A^{*} v) (φ) = {\begin{cases} - \frac{d}{d φ} v (φ) & if φ \in (0, τ], \\ L^{*} v (0) + R^{*} v (τ) & if φ = 0 \end{cases}

4.1

and

(u, v) = u^{*} (0) v (0) + \int_{0}^{τ} u^{*} (φ) R v (φ - τ) d φ,

4.2

where the superscript * of R, L and u refers to conjugate transpose.

(a) Right and left eigenvectors

The four-dimensional centre subspace of the associated linear system is tangent to the plane spanned by those infinite-dimensional right eigenvectors (eigenfunctions) of operator $A$ , which correspond to the four critical characteristic exponents λ=±iω_k, k=1,2. From this point on, index k=1,2 is used to indicate that an expression is related to one of the angular frequencies, i.e. either to ω₁ or to ω₂. For τ=τ_HH and w=w_HH, the four eigenvectors s_k(θ) and ${\bar{s}}_{k} (θ)$ satisfy

(A s_{k}) (θ) = i ω_{k} s_{k} (θ) and (A {\bar{s}}_{k}) (θ) = - i ω_{k} {\bar{s}}_{k} (θ),

4.3

where overbar indicates complex conjugate. Decomposing the eigenvectors into real and imaginary parts as s_k(θ)=s_k,R(θ)+is_k,I(θ), substituting τ=τ_HH, w=w_HH and the definition (2.5) of operator

A

, we get the boundary-value problem

\frac{d}{d θ} s (θ) = B_{8 \times 8} s (θ), θ \in [- τ_{HH}, 0)

4.4

and

L_{8 \times 8} s (0) + R_{8 \times 8} s (- τ_{HH}) = B_{8 \times 8} s (0) .

4.5

The coefficient matrices and s(θ) are defined as

\begin{aligned} s (θ) & = [\begin{matrix} s_{1, R} (θ) \\ s_{1, I} (θ) \\ s_{2, R} (θ) \\ s_{2, I} (θ) \end{matrix}], B_{8 \times 8} = [\begin{matrix} 0 & - ω_{1} I & 0 & 0 \\ ω_{1} I & 0 & 0 & 0 \\ 0 & 0 & 0 & - ω_{2} I \\ 0 & 0 & ω_{2} I & 0 \end{matrix}] \\ and L_{8 \times 8} & = diag [L, L, L, L] |_{w = w_{HH}}, R_{8 \times 8} = diag [R, R, R, R] |_{w = w_{HH}}, \end{aligned}}

4.6

where I and 0 denote the 2×2 identity and zero matrices, respectively, and diag refers to block-diagonal matrices. The solution of (4.4) can be written in exponential form as

s (θ) = e^{B_{8 \times 8} θ} c

. The constant c can be determined from (4.5). In order to select the norm of the right eigenvectors, four components of c can be chosen arbitrarily. Now we choose c₁=1, c₃=0, c₅=1 and c₇=0, whence we get

s_{k, R} (θ) = [\begin{matrix} \cos (ω_{k} θ) \\ - ω_{k} \sin (ω_{k} θ) \end{matrix}] and s_{k, I} (θ) = [\begin{matrix} \sin (ω_{k} θ) \\ ω_{k} \cos (ω_{k} θ) \end{matrix}] .

4.7

The decomposition theorem of [3] also uses the so-called left eigenvectors: the eigenvectors n_k(φ) and ${\bar{n}}_{k} (φ)$ of the adjoint operator $A^{*}$ . Since the eigenvalues of $A^{*}$ are complex conjugates to those of $A$ , the left eigenvectors n_k(φ) and ${\bar{n}}_{k} (φ)$ satisfy

(A^{*} n_{k}) (φ) = - i ω_{k} n_{k} (φ) and (A^{*} {\bar{n}}_{k}) (φ) = i ω_{k} {\bar{n}}_{k} (φ),

4.8

for τ=τ_HH and w=w_HH. We decompose the eigenvectors into real and imaginary parts as n_k(φ)=n_k,R(φ)+in_k,I(φ), where n_k,R(φ) and n_k,I(φ) can be obtained from the boundary-value problem defined by (4.8) and (4.1) in the same way as we computed s_k,R(θ) and s_k,I(θ). The norm of the left eigenvectors cannot be selected arbitrarily, as they must satisfy the following orthonormality condition in order to apply the decomposition theorem of [3]:

(n_{1, R}, s_{1, R}) = 1, (n_{1, R}, s_{1, I}) = 0, (n_{2, R}, s_{2, R}) = 1 and (n_{2, R}, s_{2, I}) = 0.

4.9

Solving the boundary value problem (4.8) constrained by the orthonormality condition (4.9), we get the left eigenfunctions in the form

\begin{aligned} n_{k, R} (φ) & = \frac{2}{p_{k}^{2} + q_{k}^{2}} [\begin{matrix} (2 ζ p_{k} + ω_{k} q_{k}) \cos (ω_{k} φ) + (ω_{k} p_{k} - 2 ζ q_{k}) \sin (ω_{k} φ) \\ p_{k} \cos (ω_{k} φ) - q_{k} \sin (ω_{k} φ) \end{matrix}] \\ and n_{k, I} (φ) & = \frac{2}{p_{k}^{2} + q_{k}^{2}} [\begin{matrix} (- ω_{k} p_{k} + 2 ζ q_{k}) \cos (ω_{k} φ) + (2 ζ p_{k} + ω_{k} q_{k}) \sin (ω_{k} φ) \\ q_{k} \cos (ω_{k} φ) + p_{k} \sin (ω_{k} φ) \end{matrix}], \end{aligned}}

4.10

where the constants p_k and q_k read

p_{k} = 2 ζ + τ_{HH} (1 + w_{HH} - ω_{k}^{2}) and q_{k} = 2 ω_{k} (1 + ζ τ_{HH}) .

4.11

(b) Decomposition of the solution space

Using the right and left eigenvectors, the four-dimensional centre subspace can now be separated via the decomposition theorem given by (3.10) and (3.11) in [3], ch. 7. We decompose the solution as

y_{t} (θ) = z_{1} (t) s_{1, R} (θ) + z_{2} (t) s_{1, I} (θ) + z_{3} (t) s_{2, R} (θ) + z_{4} (t) s_{2, I} (θ) + y_{tn} (θ),

4.12

where z₁(t), z₂(t), z₃(t) and z₄(t) are local coordinates on the attractive centre manifold, which describe the long-term dynamics of the time-delay system (2.1), whereas y_tn(θ) accounts for the remaining infinite-dimensional stable subsystem with coordinates transverse to the centre manifold. The decomposition theorem gives the formula of these components:

\begin{aligned} z_{1} (t) & = (n_{1, R}, y_{t}), z_{2} (t) = (n_{1, I}, y_{t}), \\ z_{3} (t) & = (n_{2, R}, y_{t}), z_{4} (t) = (n_{2, I}, y_{t}) \\ and y_{tn} (θ) & = y_{t} (θ) - z_{1} (t) s_{1, R} (θ) - z_{2} (t) s_{1, I} (θ) - z_{3} (t) s_{2, R} (θ) - z_{4} (t) s_{2, I} (θ) . \end{aligned}}

4.13

From this point on, we omit the argument t from z₁, z₂, z₃ and z₄ for the sake of simplicity. Differentiating (4.13) with respect to time and using (2.4), (4.12) and (4.3), we get

\begin{aligned} [\begin{matrix} z_{1}^{'} \\ z_{2}^{'} \\ z_{3}^{'} \\ z_{4}^{'} \\ y_{t n}^{'} \end{matrix}] & = [\begin{matrix} 0 & ω_{1} & 0 & 0 & O \\ - ω_{1} & 0 & 0 & 0 & O \\ 0 & 0 & 0 & ω_{2} & O \\ 0 & 0 & - ω_{2} & 0 & O \\ o & o & o & o & A \end{matrix}] [\begin{matrix} z_{1} \\ z_{2} \\ z_{3} \\ z_{4} \\ y_{tn} \end{matrix}] + [\begin{matrix} n_{1, R, 2} (0) F_{2} (0) \\ n_{1, I, 2} (0) F_{2} (0) \\ n_{2, R, 2} (0) F_{2} (0) \\ n_{2, I, 2} (0) F_{2} (0) \\ G (y_{t}) \end{matrix}] \\ and G (y_{t}) & = F (y_{t}) - F_{2} (0) (n_{1, R, 2} (0) s_{1, R} + n_{1, I, 2} (0) s_{1, I} \\ + n_{2, R, 2} (0) s_{2, R} + n_{2, I, 2} (0) s_{2, I}), \end{aligned}}

4.14

where

o : R \to B

is a zero operator,

O : B \to R

is a zero functional, subscript 2 in n_k,R,2 and n_k,I,2 indicates the second component of vectors, and

F_{2} (0)

shortly indicates the second component of

F (y_{t})

at θ=0 with τ=τ_HH and w=w_HH.

(c) Approximation of the centre manifold

Equation (4.14) shows that the four-dimensional centre subsystem is decoupled in the associated linear system, but there is still coupling in the nonlinear terms. In order to fully decouple the four-dimensional subsystem, the dynamics must be restricted to the centre manifold of form $y_{t n} = y_{t n}^{CM} (z_{1}, z_{2}, z_{3}, z_{4})$ . The nonlinear terms in the first four rows of (4.14) must be expanded into Taylor series in terms of z₁, z₂, z₃ and z₄ in order to do normal form analysis—we must expand $F_{2} (0)$ up to third order. The computation of the centre manifold and the Taylor expansion is rather lengthy, but it can be tackled by symbolic algebra. Now we present the milestones of this process.

Substituting the decomposed solution (4.12) into the definition (2.6) of $F$ , we get the Taylor series expansion of $F_{2} (0)$ if the centre manifold $y_{tn} = y_{tn}^{CM} (z_{1}, z_{2}, z_{3}, z_{4})$ is expanded into Taylor series in terms of z_1,2,3,4. The second-order terms in $F_{2} (0)$ are independent of $y_{tn}^{CM}$ , the expansion of the centre manifold is necessary only to obtain the cubic (and higher-order) terms in $F_{2} (0)$ . The quadratic part of $F_{2} (0)$ is of form

\begin{aligned} F_{2}^{2 nd} (0) & = F_{11} z_{1}^{2} + F_{22} z_{2}^{2} + F_{33} z_{3}^{2} + F_{44} z_{4}^{2} + F_{12} z_{1} z_{2} + F_{13} z_{1} z_{3} \\ + F_{14} z_{1} z_{4} + F_{23} z_{2} z_{3} + F_{24} z_{2} z_{4} + F_{34} z_{3} z_{4} \\ and F_{m n} & = \frac{1}{2} {\frac{\partial^{2} F_{2} (0)}{\partial z_{m}^{2}} |}_{0}, if m = n, F_{m n} = {\frac{\partial^{2} F_{2} (0)}{\partial z_{m} \partial z_{n}} |}_{0}, if m < n, \end{aligned}}

4.15

m,n=1,2,3,4, where the subscript 0 stands for the substitution y_t(0)=0. In order to include all cubic terms in

F_{2} (0)

, we need at least the second-order expansion of the centre manifold y^CM_tn in terms of z_1,2,3,4:

\begin{aligned} y_{tn}^{CM} (z_{1}, z_{2}, z_{3}, z_{4}) (θ) & \approx \frac{1}{2} (h_{11} (θ) z_{1}^{2} + h_{22} (θ) z_{2}^{2} + h_{33} (θ) z_{3}^{2} + h_{44} (θ) z_{4}^{2} \\ + 2 h_{12} (θ) z_{1} z_{2} + 2 h_{13} (θ) z_{1} z_{3} + 2 h_{14} (θ) z_{1} z_{4} + 2 h_{23} (θ) z_{2} z_{3} \\ + 2 h_{24} (θ) z_{2} z_{4} + 2 h_{34} (θ) z_{3} z_{4}) . \end{aligned}

4.16

The coefficients h_mn(θ) (m,n=1,2,3,4, m≤n) can be determined as follows. Both sides of (4.16) are differentiated with respect to time, then the first four rows of (4.14) are substituted into the right-hand side and the fifth row of (4.14) is substituted into the left-hand side. The case θ∈[−τ_HH,0) is considered and the definitions (2.5)–(2.6) of $A$ and $F$ are substituted accordingly. Afterwards, the derivative of (4.16) with respect to θ is substituted and the second-order approximation $F_{2} (0) \approx F_{2}^{2 nd} (0)$ is used. Finally, the coefficients of the quadratic terms of z_1,2,3,4 are collected and a polynomial balance is considered. This leads to three decoupled sets of non-homogeneous first-order differential equations:

\frac{d}{d θ} H_{l} (θ) = C_{l} H_{l} (θ) + p_{1, l} \cos (ω_{1} θ) + q_{1, l} \sin (ω_{1} θ) + p_{2, l} \cos (ω_{2} θ) + q_{2, l} \sin (ω_{2} θ),

4.17

l=1,2,3, where

\begin{aligned} H_{1} (θ) = [\begin{matrix} h_{11} (θ) \\ h_{12} (θ) \\ h_{22} (θ) \end{matrix}], H_{2} (θ) = [\begin{matrix} h_{33} (θ) \\ h_{34} (θ) \\ h_{44} (θ) \end{matrix}], H_{3} (θ) = [\begin{matrix} h_{13} (θ) \\ h_{14} (θ) \\ h_{23} (θ) \\ h_{24} (θ) \end{matrix}], \\ Q_{k} = \frac{2}{p_{k}^{2} + q_{k}^{2}} [\begin{matrix} p_{k} & q_{k} \\ q_{k} ω_{k} & - p_{k} ω_{k} \end{matrix}], [p_{k, 1} q_{k, 1}] = [\begin{matrix} 2 F_{11} Q_{k} \\ F_{12} Q_{k} \\ 2 F_{22} Q_{k} \end{matrix}], \\ [p_{k, 2} q_{k, 2}] = [\begin{matrix} 2 F_{33} Q_{k} \\ F_{34} Q_{k} \\ 2 F_{44} Q_{k} \end{matrix}], [p_{k, 3} q_{k, 3}] = [\begin{matrix} F_{13} Q_{k} \\ F_{14} Q_{k} \\ F_{23} Q_{k} \\ F_{24} Q_{k} \end{matrix}] \\ and C_{k} & = [\begin{matrix} 0 & - 2 ω_{k} I & 0 \\ ω_{k} I & 0 & - ω_{k} I \\ 0 & 2 ω_{k} I & 0 \end{matrix}], C_{3} = [\begin{matrix} 0 & - ω_{2} I & - ω_{1} I & 0 \\ ω_{2} I & 0 & 0 & - ω_{1} I \\ ω_{1} I & 0 & 0 & - ω_{2} I \\ 0 & ω_{1} I & ω_{2} I & 0 \end{matrix}] . \end{aligned}}

4.18

The solution of (4.17) takes the form

H_{l} (θ) = e^{C_{l} θ} K_{l} + M_{1, l} \cos (ω_{1} θ) + N_{1, l} \sin (ω_{1} θ) + M_{2, l} \cos (ω_{2} θ) + N_{2, l} \sin (ω_{2} θ),

4.19

l=1,2,3. The coefficients M_k,l and N_k,l can be determined by substituting (4.19) back into (4.17) and considering the harmonic balance of the trigonometric terms, which yields

[\begin{matrix} M_{k, l} \\ N_{k, l} \end{matrix}] = {[\begin{matrix} - C_{l} & ω_{k} I \\ - ω_{k} I & - C_{l} \end{matrix}]}^{- 1} [\begin{matrix} p_{k, l} \\ q_{k, l} \end{matrix}] .

4.20

In order to determine the coefficient K_l in (4.19), a boundary condition corresponding to (4.17) must be derived and satisfied. The boundary condition is formulated similarly as the differential equation (4.17). Both sides of (4.16) are differentiated with respect to time, then the first four rows of (4.14) are substituted into the right-hand side and the fifth row of (4.14) is substituted into the left-hand side as before. Now the case θ=0 is considered when using the definitions (2.5)–(2.6) of $A$ and $F$ . Afterwards, (4.16) is substituted and the second-order approximation $F_{2} (0) \approx F_{2}^{2 nd} (0)$ is used. The coefficients of the quadratic terms of z_1,2,3,4 are collected for a polynomial balance and, by taking τ=τ_HH, w=w_HH, it leads to the boundary conditions for (4.17) in the form

P_{l} H_{l} (0) + R_{l} H_{l} (- τ_{HH}) = p_{1, l} + p_{2, l} + r_{l},

4.21

where

\begin{aligned} L_{k} & = L_{6 \times 6} = diag [L, L, L] |_{w = w_{HH}}, L_{3} = L_{8 \times 8}, P_{l} = L_{l} - C_{l}, \\ R_{k} & = R_{6 \times 6} = diag [R, R, R] |_{w = w_{HH}}, R_{3} = R_{8 \times 8}, \\ r_{1} & = - {[0 2 F_{11} 0 F_{12} 0 2 F_{22}]}^{T}, \\ r_{2} & = - {[0 2 F_{33} 0 F_{34} 0 2 F_{44}]}^{T} \\ and r_{3} & = - {[0 F_{13} 0 F_{14} 0 F_{23} 0 F_{24}]}^{T} . \end{aligned}}

4.22

Finally, substituting the trial solution (4.19) into the boundary condition (4.21), we obtain the coefficient K_l. After some simplifications, we can write K_l in the form

K_{l} = {(P_{l} + R_{l} e^{- C_{l} τ})}^{- 1} r_{l} .

4.23

This way, the solution (4.19) of the above boundary-value problem is constructed and the coefficients h_mn(θ) (m,n=1,2,3,4, m≤n) are available. The second-order approximation (4.16) of the centre manifold is obtained and, with (4.12), it can be used to derive the third-order expansion of $F_{2} (0)$ in terms of z_1,2,3,4. Thus, the nonlinearity in the first four rows of (4.14) is now expressed in terms of the four local coordinates of the centre manifold, and we get a four-dimensional subsystem with cubic nonlinearity in the form

[\begin{matrix} z_{1}^{'} \\ z_{2}^{'} \\ z_{3}^{'} \\ z_{4}^{'} \end{matrix}] = [\begin{matrix} 0 & ω_{1} & 0 & 0 \\ - ω_{1} & 0 & 0 & 0 \\ 0 & 0 & 0 & ω_{2} \\ 0 & 0 & - ω_{2} & 0 \end{matrix}] [\begin{matrix} z_{1} \\ z_{2} \\ z_{3} \\ z_{4} \end{matrix}] + [\begin{matrix} G_{1} (z_{1}, z_{2}, z_{3}, z_{4}) \\ G_{2} (z_{1}, z_{2}, z_{3}, z_{4}) \\ G_{3} (z_{1}, z_{2}, z_{3}, z_{4}) \\ G_{4} (z_{1}, z_{2}, z_{3}, z_{4}) \end{matrix}] .

4.24

The lengthy process of centre manifold reduction has been performed in order to decouple the four-dimensional centre subsystem (4.24). The advantage of this formulation is that the long-term dynamics of the original infinite-dimensional time-delay system (2.1) can be investigated by analysing the finite-dimensional ordinary differential equation (4.24). Hence, from this point on, the bifurcation theorems and normal forms of ordinary differential equations can be used, which are well-known in the literature [17,19]. In the next section, we use normal form theory to analyse the flow on the centre manifold and show the existence of periodic and quasi-periodic motions.

(d) Normal form equations

The four-dimensional system (4.24) can be transformed into polar form with two amplitudes r₁, r₂ and two phase angles θ₁, θ₂ as

\begin{aligned} r_{1}^{'} & = μ_{1} r_{1} + a_{11} r_{1}^{3} + a_{12} r_{1} r_{2}^{2}, θ_{1}^{'} = ω_{1} + c_{11} r_{1}^{2} + c_{12} r_{2}^{2} \\ and r_{2}^{'} & = μ_{2} r_{2} + a_{21} r_{1}^{2} r_{2} + a_{22} r_{2}^{3}, θ_{2}^{'} = ω_{2} + c_{21} r_{1}^{2} + c_{22} r_{2}^{2} . \end{aligned}}

4.25

Now we focus on the amplitudes r₁, r₂ only and investigate the existence of periodic and quasi-periodic solutions by giving analytical formulae for coefficients μ₁, μ₂, and a₁₁, a₁₂, a₂₁, a₂₂.

The coefficients μ₁, μ₂ of the linear terms in (4.25) are unfolding parameters that are functions of the bifurcation parameters. The double Hopf bifurcation is a codimension two bifurcation, hence two bifurcation parameters must be selected. Now we choose the dimensionless chip width w and the dimensionless angular velocity Ω. For convenience, we introduce the shifted parameters $\hat{w} = w - w_{HH}$ and $\hat{Ω} = Ω - Ω_{HH}$ to investigate the system in the vicinity of the double Hopf bifurcation. We can approximate the unfolding parameters by linear functions of the bifurcation parameters as

\begin{aligned} μ_{1} & = γ_{11} \hat{w} + γ_{12} \hat{Ω} \\ and μ_{2} & = γ_{21} \hat{w} + γ_{22} \hat{Ω}, \end{aligned}}

4.26

where the constants γ_mn (m,n=1,2) are called the root tendencies. The constants γ_mn represent the speed by which the four critical characteristic exponents cross the imaginary axis during the double Hopf bifurcation. They can be obtained via implicit differentiation of the characteristic equation D(λ)=0. Using (3.1), this yields

\begin{aligned} γ_{k 1} & = Re ({\frac{\partial λ}{\partial w} |}_{λ = i ω_{k}}) = Re (- {\frac{\partial D / \partial w}{\partial D / \partial λ} |}_{λ = i ω_{k}}) = \frac{p_{k} V_{k} + q_{k} W_{k}}{p_{k}^{2} + q_{k}^{2}} \\ and γ_{k 2} & = Re ({\frac{\partial λ}{\partial Ω} |}_{λ = i ω_{k}}) = Re (- {\frac{(\partial D / \partial τ) (d τ / d Ω)}{\partial D / \partial λ} |}_{λ = i ω_{k}}) \\ = w_{HH} \frac{τ_{HH}^{2}}{2 π} ω_{k} \frac{q_{k} (V_{k} + 1) - p_{k} W_{k}}{p_{k}^{2} + q_{k}^{2}}, \end{aligned}}

4.27

where p_k and q_k are defined by (4.11), whereas V _k and W_k are given by

V_{k} = \cos (ω_{k} τ_{HH}) - 1 and W_{k} = - \sin (ω_{k} τ_{HH}) .

4.28

The coefficients a_mn (m,n=1,2) of the cubic terms in (4.25) can be obtained from the results of centre manifold reduction. These normal form coefficients can directly be calculated from the nonlinear terms of (4.24) by applying the formulae derived in [43], whence we obtain

\begin{aligned} a_{11} & = - \frac{1}{2} w^{2} (V_{1}^{2} + W_{1}^{2}) η_{2}^{2} (\frac{K_{1} R_{1} + L_{1} S_{1}}{K_{1}^{2} + L_{1}^{2}} \frac{p_{1} V_{1} + q_{1} W_{1}}{p_{1}^{2} + q_{1}^{2}} \\ + \frac{K_{1} S_{1} - L_{1} R_{1}}{K_{1}^{2} + L_{1}^{2}} \frac{q_{1} V_{1} - p_{1} W_{1}}{p_{1}^{2} + q_{1}^{2}}) + \frac{3}{4} w (V_{1}^{2} + W_{1}^{2}) η_{3} \frac{p_{1} V_{1} + q_{1} W_{1}}{p_{1}^{2} + q_{1}^{2}} \\ and a_{12} & = - w^{2} (V_{2}^{2} + W_{2}^{2}) η_{2}^{2} ((\frac{K_{3} R_{3} + L_{3} S_{3}}{K_{3}^{2} + L_{3}^{2}} + \frac{K_{4} R_{4} + L_{4} S_{4}}{K_{4}^{2} + L_{4}^{2}}) \frac{p_{1} V_{1} + q_{1} W_{1}}{p_{1}^{2} + q_{1}^{2}} \\ + (\frac{K_{3} S_{3} - L_{3} R_{3}}{K_{3}^{2} + L_{3}^{2}} - \frac{K_{4} S_{4} - L_{4} R_{4}}{K_{4}^{2} + L_{4}^{2}}) \frac{q_{1} V_{1} - p_{1} W_{1}}{p_{1}^{2} + q_{1}^{2}}) \\ + \frac{3}{2} w (V_{2}^{2} + W_{2}^{2}) η_{3} \frac{p_{1} V_{1} + q_{1} W_{1}}{p_{1}^{2} + q_{1}^{2}}, \end{aligned}}

4.29

where η₂ and η₃ are the cutting force coefficients that appear in (2.1), while the auxiliary parameters are

\begin{aligned} R_{k} & = 1 - \cos (2 ω_{k} τ_{HH}), S_{k} = \sin (2 ω_{k} τ_{HH}), \\ R_{3, 4} & = 1 - \cos ((ω_{2} \pm ω_{1}) τ_{HH}), S_{3, 4} = \sin ((ω_{2} \pm ω_{1}) τ_{HH}), \\ K_{k} & = w R_{k} - (4 ω_{k}^{2} - 1), L_{k} = w S_{k} + 4 ζ ω_{k} \\ and K_{3, 4} & = w R_{3, 4} - ({(ω_{2} \pm ω_{1})}^{2} - 1), L_{3, 4} = w S_{3, 4} + 2 ζ (ω_{2} \pm ω_{1}) . \end{aligned}}

4.30

The expressions of a₂₂ and a₂₁ are the same as that of a₁₁ and a₁₂, respectively, but the roles of ω₁ and ω₂ are interchanged. The formula of a₂₂ can be obtained from that of a₁₁ by replacing K₁, L₁, p₁, q₁, V ₁, W₁, R₁ and S₁ with K₂, L₂, p₂, q₂, V ₂, W₂, R₂ and S₂, respectively. Whereas the formula of a₂₁ is the same as that of a₁₂ but with p₂, q₂, V ₂, W₂, V ₁ and W₁ instead of p₁, q₁, V ₁, W₁, V ₂ and W₂, respectively, and replacing L₄, S₄ with −L₄, −S₄ (see the change of the sign when interchanging ω₁ and ω₂ in the definition (4.30) of L₄ and S₄).

The above formulae are valid for non-resonant double Hopf bifurcation, where the angular frequencies ω₁ and ω₂ are not related by small integer numbers. Note that for specific values of the damping ratio ζ, it might be possible to obtain weak resonant double Hopf bifurcations [44], but their analysis is out of scope of this paper. The formula (4.29) of a₁₁ is the same as that obtained for the Poincaré–Lyapunov constant via the analysis of the single Hopf bifurcation in [11]. It is also important to note that the cubic coefficients a_mn (m,n=1,2) are in fact functions of the bifurcation parameters w and Ω, and (4.29) gives only their constant approximation. In §5, we show a linear approximation of the cubic coefficients a_mn in terms of the bifurcation parameters w and Ω following [45]. The linear approximation leads to a higher-order estimation of the amplitude of the arising periodic and quasi-periodic motions.

5. Stability charts and bifurcation diagrams

The analysis of the polar-form system (4.25) can be found in [17]. Accordingly, (4.25) can be further simplified by the transformation $r_{1} = {\bar{r}}_{1} / \sqrt{| a_{11} |}$ , $r_{2} = {\bar{r}}_{2} / \sqrt{| a_{22} |}$ :

\begin{aligned} {\bar{r}}_{1}^{'} & = {\bar{r}}_{1} (μ_{1} + a {\bar{r}}_{1}^{2} + b {\bar{r}}_{2}^{2}) \\ and {\bar{r}}_{2}^{'} & = {\bar{r}}_{2} (μ_{2} + c {\bar{r}}_{1}^{2} + d {\bar{r}}_{2}^{2}), \end{aligned}}

5.1

where a=a₁₁/|a₁₁|, b=a₁₂/|a₂₂|, c=a₂₁/|a₁₁| and d=a₂₂/|a₂₂|. These four parameters are important, since they determine the possible topologies in the phase portraits of the polar-form system (4.25). Simplified phase portraits can be depicted in the plane (r₁,r₂) for different values of the bifurcation parameters w and Ω. As it was shown in [17], the topology of these phase portraits depends on the signs of a, b, c, d and A=ad−bc. Depending on their signs, 12 different cases of unfolding (12 sets of different topologies) can occur ([17], pp. 399–409). For (2.1), the Hopf bifurcation is always subcritical [11], thus a>0 and d>0 hold, which excludes some of these 12 cases. As it is shown below, two of the 12 topologies can be identified.

(a) Phase portraits and topologies

Here we show a numerical case study for damping ratio ζ=0.02 and cutting-force coefficients η₂=1.43059 and η₃=0.738487, which are the dimensionless counterparts of actual measured cutting-force coefficients reported in [5] with feed h₀=250 μm per revolution. We found that topology Case Ia of [17] occurs at the intersection of the first and the second stability lobes, while the intersections of higher-order lobes are all associated with Case Ib of [17]. In what follows, we demonstrate the differences between the two cases.

Consider the intersection of the first (j₁=1) and the second (j₂=2) stability lobes, which is marked by a circle in figure 1b and is also shown by thick lines in more detail in figure 2a. The first row of table 1 presents the location (Ω_HH,w_HH) of the intersection point, the two angular frequencies ω₁, ω₂ and the root tendencies γ_mn (m,n=1,2). The cubic coefficients a_mn (m,n=1,2) and parameters a, b, c, d and A are listed in the first row of table 2. In this example, a, b, c, d and A are all positive, which corresponds to Case Ia in [17]. We also computed coefficients b and c numerically by DDE-Biftool [40] and obtained the same results up to the displayed digits.

Figure 2. — Figure 2. The possible phase portrait topologies in the vicinity of the double Hopf bifurcation at the intersection of the first and the second lobes (a); at the intersection of the fourth and the fifth lobes (b). The intersecting stability lobes are shown by thick blue and cyan lines. Sources, sinks and saddles are indicated by red, green and blue dots, respectively. (Online version in colour.)

Table 1.The parameters of the double Hopf bifurcation points shown in figure 2.
Collapse
j₁	j₂	ω₁	ω₂	Ω_HH	w_HH	γ₁₁	γ₁₂	γ₂₁	γ₂₂
1	2	1.0006	1.52994	1.01018	0.671754	0.005553	0.3623	0.2963	−0.6699
4	5	1.00247	1.15031	0.253092	0.164877	0.02429	1.251	0.3182	−1.530

Table 2.The cubic coefficients of the double Hopf bifurcation points shown in figure 2.
Collapse
j₁	j₂	a₁₁	a₁₂	a₂₁	a₂₂	a	b	c	d	A
1	2	8.090×10⁻⁶	0.008128	0.0002925	0.4317	1	0.01883	36.16	1	0.3193
4	5	1.460×10⁻⁴	0.01560	0.002356	0.1121	1	0.1392	16.14	1	−1.248

The corresponding phase portraits and their topologies are shown in figure 2a. According to Guckenheimer & Holmes [17], six different topologies can be distinguished for Case Ia depending on the unfolding parameters μ₁ and μ₂. Correspondingly, six sectors with different topologies can be separated in the plane (Ω,w) of the bifurcation parameters. The straight lines separating these sectors are

\begin{aligned} Line 1: μ_{1} = 0 & \Rightarrow \hat{w} = - \frac{γ_{12}}{γ_{11}} \hat{Ω}, \\ Line 2: μ_{2} = 0 & \Rightarrow \hat{w} = - \frac{γ_{22}}{γ_{21}} \hat{Ω}, \\ Line 3: μ_{2} = \frac{c}{a} μ_{1} & \Rightarrow \hat{w} = \frac{c γ_{12} - a γ_{22}}{a γ_{21} - c γ_{11}} \hat{Ω} \\ and Line 4: μ_{2} = \frac{d}{b} μ_{1} & \Rightarrow \hat{w} = \frac{d γ_{12} - b γ_{22}}{b γ_{21} - d γ_{11}} \hat{Ω} . \end{aligned}}

5.2

Note that Lines 1 and 2 are the tangents of the intersecting stability lobes at the double Hopf point (the lobes are shown by thick blue and cyan lines in figure 2).

Let us investigate the sectors of plane (μ₁,μ₂) by going around anti-clockwise. In Sector I (μ₁>0, μ₂>0), (4.25) has only a single, trivial equilibrium (r₁,r₂)=(0,0), which corresponds to the trivial equilibrium of the original time-delay system (2.1) and exists in all the six sectors. The trivial equilibrium (0,0) is a source (unstable node) in Sector I. The equilibrium undergoes a pitchfork bifurcation while crossing Line 1, which gives rise to a non-trivial equilibrium $(r_{1}, r_{2}) = (r_{1}^{p}, 0)$ , $r_{1}^{p} \neq 0$ in Sector II (μ₁<0, μ₂>0). This equilibrium corresponds to a periodic solution (P₁) in the time-delay system (2.1), whereas the pitchfork bifurcation in system (4.25) corresponds to a Hopf bifurcation in the DDE (2.1). In Sector II, the equilibrium $(r_{1}^{p}, 0)$ is a source, therefore the corresponding periodic solution P₁ is unstable, while the trivial equilibrium (0,0) becomes a saddle. Another pitchfork bifurcation of (0,0) takes place along Line 2, and a non-trivial equilibrium $(r_{1}, r_{2}) = (0, r_{2}^{p})$ , $r_{2}^{p} \neq 0$ is born when entering Sector III (μ₁<0, cμ₁/a<μ₂<0). This, again, corresponds to a Hopf bifurcation of the equilibrium of the infinite-dimensional time-delay system (2.1), which gives rise to another periodic solution (P₂). In Sector III, the trivial equilibrium becomes a sink (stable node), the equilibrium $(r_{1}^{p}, 0)$ remains a source and the equilibrium $(0, r_{2}^{p})$ is a saddle. When entering Sector IV (μ₁<0, dμ₁/b<μ₂<cμ₁/a), a pitchfork bifurcation of $(r_{1}^{p}, 0)$ gives rise to the non-trivial equilibrium $(r_{1}, r_{2}) = (r_{1}^{qp}, r_{2}^{qp})$ , $r_{1}^{qp} \neq 0$ , $r_{2}^{qp} \neq 0$ . This phenomenon corresponds to a torus bifurcation of the periodic solution P₁ of (2.1), by which a quasi-periodic solution (QP) arises. Thus, in Sector IV, the trivial equilibrium coexists with two periodic solutions and a quasi-periodic one. The trivial equilibrium is still a sink, the equilibrium $(0, r_{2}^{p})$ is still a saddle, while the equilibrium $(r_{1}^{p}, 0)$ becomes a saddle and the equilibrium $(r_{1}^{qp}, r_{2}^{qp})$ is born as a source. The equilibrium $(r_{1}^{qp}, r_{2}^{qp})$ merges with the equilibrium $(0, r_{2}^{p})$ via a pitchfork bifurcation when leaving Sector IV by crossing Line 4 to Sector V (μ₁<0, μ₂<dμ₁/b). Equivalently, a torus bifurcation of P₂ occurs in the infinite-dimensional system (2.1). In Sector V, the equilibrium $(0, r_{2}^{p})$ becomes a source, while the equilibrium $(r_{1}^{p}, 0)$ remains a saddle and the trivial equilibrium remains a sink. Entering Sector VI (μ₁>0, μ₂<0) by crossing Line 1, the equilibrium $(r_{1}^{p}, 0)$ disappears via a pitchfork bifurcation, which corresponds to a Hopf bifurcation of the trivial equilibrium of (2.1). Here, the equilibrium $(0, r_{2}^{p})$ remains a source and the trivial equilibrium becomes a saddle. Finally, the equilibrium $(r_{2}^{p}, 0)$ disappears via a pitchfork bifurcation when returning to Sector I by crossing Line 2, i.e. P₂ vanishes as a result of a Hopf bifurcation in the time-delay system (2.1).

A second example is presented in figure 2b, where the intersection of the fourth (j₁=4) and the fifth (j₂=5) stability lobes is shown, see also the square mark in figure 1b. The parameters of the double Hopf point, the corresponding root tendencies and the cubic coefficients are listed in the second row of tables 1 and 2. In this example, a, b, c and d are positive but now A is negative, thus Case Ib in [17] applies. In this case, the order of Lines 3 and 4 changes when going around the sectors of the plane (μ₁,μ₂) anti-clockwise. Therefore, the quasi-periodic solution QP is born from P₂ (and not P₁) when entering Sector IV along Line 4, and it collapses into P₁ (and not P₂) when leaving Sector IV along Line 3. Interchanging Lines 3 and 4 modifies the topology only in Sector IV: now the equilibria $(r_{1}^{p}, 0)$ and $(0, r_{2}^{p})$ are sources and the equilibrium $(r_{1}^{qp}, r_{2}^{qp})$ is a saddle. Consequently, the periodic solutions P₁ and P₂ behave as sources and the quasi-periodic solution QP behaves as a saddle.

(b) Bifurcation diagrams via higher-order estimation

The amplitudes of the arising periodic and quasi-periodic orbits can be determined by calculating the non-trivial equilibria of (5.1). According to [17], they can be calculated as

\begin{aligned} P_{1} : {\bar{r}}_{1}^{p} = \sqrt{- \frac{μ_{1}}{a}} & \Rightarrow r_{1}^{p} (\hat{w}, \hat{Ω}) = \sqrt{- \frac{γ_{11} \hat{w} + γ_{12} \hat{Ω}}{a_{11}}}, \\ P_{2} : {\bar{r}}_{2}^{p} = \sqrt{- \frac{μ_{2}}{d}} & \Rightarrow r_{2}^{p} (\hat{w}, \hat{Ω}) = \sqrt{- \frac{γ_{21} \hat{w} + γ_{22} \hat{Ω}}{a_{22}}}, \\ and QP : {\bar{r}}_{1}^{qp} = \sqrt{\frac{b μ_{2} - d μ_{1}}{A}} & \Rightarrow r_{1}^{qp} (\hat{w}, \hat{Ω}) = \sqrt{\frac{(a_{12} γ_{21} - a_{22} γ_{11}) \hat{w} + (a_{12} γ_{22} - a_{22} γ_{12}) \hat{Ω}}{a_{11} a_{22} - a_{12} a_{21}}} \\ {\bar{r}}_{2}^{qp} = \sqrt{\frac{c μ_{1} - a μ_{2}}{A}} & \Rightarrow r_{2}^{qp} (\hat{w}, \hat{Ω}) = \sqrt{\frac{(a_{21} γ_{11} - a_{11} γ_{21}) \hat{w} + (a_{21} γ_{12} - a_{11} γ_{22}) \hat{Ω}}{a_{11} a_{22} - a_{12} a_{21}}} . \end{aligned}}

5.3

Here we emphasize that the amplitudes depend on the two bifurcation parameters

\hat{w}

and

\hat{Ω}

. Later on we omit the argument

(\hat{w}, \hat{Ω})

for simplicity.

In [45], it was show that a higher-order estimation of the amplitude of the periodic orbits arising from Hopf bifurcation is possible when special global properties hold in the system. This estimation can also be done for the quasi-periodic orbit arising from the double Hopf bifurcation. The estimation is based on the approximation of the cubic coefficients a_mn (m,n=1,2) by a linear function of the bifurcation parameters Ω, w instead of using the constant values defined by (4.29):

{\hat{a}}_{m n} = a_{m n} + a_{m n}^{Ω} \hat{Ω} + a_{m n}^{w} \hat{w},

5.4

m,n=1,2. Replacing a_mn with

{\hat{a}}_{m n}

in (5.3) yields a higher-order estimation of the amplitude of the periodic and the quasi-periodic orbits as long as the additional coefficients

a_{m n}^{Ω}

and

a_{m n}^{w}

are chosen properly.

The additional coefficients can be calculated based on the following global properties of the system. When the chip width is zero, that is, w=0, the cutting force on the right-hand side of the governing equation (2.1) vanishes and the system reduces to a damped free oscillator. Hence the periodic and the quasi-periodic solutions must disappear at w=0, which happens such that their amplitudes tend to infinity. That is, $lim_{w \to 0} r_{1}^{p} = \infty$ and $lim_{w \to 0} r_{2}^{p} = \infty$ holds for the periodic solutions, whereas $lim_{w \to 0} r_{1}^{qp} = \infty$ and $lim_{w \to 0} r_{2}^{qp} = \infty$ for the quasi-periodic solution. Similarly, when τ=0 or, equivalently, $Ω \to \infty$ , the right-hand side of (2.1) becomes zero. Therefore, the periodic solutions satisfy $lim_{Ω \to \infty} r_{1}^{p} = \infty$ and $lim_{Ω \to \infty} r_{2}^{p} = \infty$ , whereas the quasi-periodic solution satisfies $lim_{Ω \to \infty} r_{1}^{qp} = \infty$ and $lim_{Ω \to \infty} r_{2}^{qp} = \infty$ . It can be shown that these properties are guaranteed with the coefficients

{\hat{a}}_{m n} = a_{m n} + \frac{a_{m n}}{w_{HH}} \hat{w},

5.5

m,n=1,2. Note that approximating the cubic coefficients by such linear functions of the bifurcation parameters instead of using the constants in (4.29) does not change the formulae of Lines 1–4 in (5.2).

The periodic and quasi-periodic solutions of (2.1) can be approximated as

\begin{aligned} y_{t} (θ) & \approx r_{1} (\hat{w}, \hat{Ω}) (\cos (ω_{1} t) s_{1, R} (θ) - \sin (ω_{1} t) s_{1, I} (θ)) \\ + r_{2} (\hat{w}, \hat{Ω}) (\cos (ω_{2} t) s_{2, R} (θ) - \sin (ω_{2} t) s_{2, I} (θ)) \end{aligned}

5.6

and

x (t) = y_{t 1} (0) \approx r_{1} (\hat{w}, \hat{Ω}) \cos (ω_{1} t) + r_{2} (\hat{w}, \hat{Ω}) \cos (ω_{2} t),

5.7

where

(r_{1}, r_{2}) = (r_{1}^{p}, 0)

and

(r_{1}, r_{2}) = (0, r_{2}^{p})

must be substituted for the periodic solutions P₁ and P₂, respectively, whereas

(r_{1}, r_{2}) = (r_{1}^{qp}, r_{2}^{qp})

for the quasi-periodic solution QP with the amplitudes given by (5.3) using

{\hat{a}}_{m n}

instead of a_mn. The bifurcation diagrams presenting the amplitude of the periodic and the quasi-periodic solutions will be presented in figure 3 together with the corresponding global stability charts.

Figure 3. — Figure 3. Stability charts (i) and bifurcation diagrams (ii) close to the double Hopf point (DH) at the intersection of the first and the second lobe (a); the second and the third lobe (b); the fourth and the fifth lobe (c). Analytical results are indicated by solid lines and dashed lines show the numerical results. Blue and cyan lines show the loci of Hopf bifurcation (H) where periodic solutions (P₁ and P₂) are born. Red and purple lines represent the loci of torus bifurcation (T) where a quasi-periodic orbit (QP) emerges from one of the periodic orbits. Green, orange and black lines indicate the loci of Big Bang Bifurcation (B³) where the tool loses contact with the workpiece during periodic or quasi-periodic oscillations. Thin black line indicates period doubling bifurcation (PD). The globally stable region is indicated by light grey shading, whereas dark grey and darker grey shadings show the bistable regions where unstable periodic and quasi-periodic solutions coexist with the linearly stable equilibrium and the large-amplitude stable motion (chatter) that involves loss of contact. (Online version in colour.)

(c) Global dynamics and region of bistability

In the linearly stable region, where the equilibrium is a sink, the coexisting periodic and quasi-periodic solutions are unstable. Owing to these unstable orbits, the basin of attraction of the linearly stable equilibrium is finite, hence the equilibrium is not stable in the global sense. To large enough perturbations, the vibrations of the machine tool amplify in a self-excited manner.

When the vibration amplitude becomes large enough, the tool jumps out of the workpiece and leaves the material during cutting. In such cases, the tool’s motion is governed by non-smooth dynamics: switching occurs between the dynamics of a cutting tool and the dynamics of a ‘flying’ tool (a damped free oscillator). The switching surface is located where the chip thickness is zero and loss of contact takes place. If the tool loses contact with the workpiece and undergoes a free flight, then (2.1) becomes no longer valid, since the cutting force on the right-hand side vanishes. It was shown in [14] for the case of the single Hopf bifurcation that once the unstable periodic motion involves loss of contact and grazes the switching surface, the periodic orbit undergoes a kind of non-smooth fold of limit cycles bifurcation called Big Bang Bifurcation (B³). This bifurcation is illustrated in figure 1c for j=2, ω=1.37 (Ω=0.901) as shown by the dashed line in figure 1b. Via the non-smooth fold, the unstable periodic orbit vanishes by merging with a large-amplitude attractive solution. This solution describes machine tool chatter with intermittent loss of contact, it may be chaotic and it coexists with the unstable periodic orbit and the linearly stable equilibrium. The region of coexistence is also called unsafe zone or region of bistability. Determining the boundary of the bistable region is important, since the cutting process unsafe in this region: it is stable only linearly but not in the global sense.

We can expect similar behaviour in the vicinity of the double Hopf point. Periodic and quasi-periodic solutions exist only up to the point where the tool first loses contact with the workpiece during its large-amplitude motion. Thus, a B³ happens where the dimensionless chip thickness h(t) first drops to zero, that is, when

h (t) = 1 + x (t - τ) - x (t) = 0,

5.8

occurs for any t. Equation (5.8) defines the loci of the B³ and gives the boundary of the bistable region for the periodic and the quasi-periodic solutions. Substituting (5.7) into (5.8), combining the trigonometric terms and approximating τ≈τ_HH, we get

\begin{aligned} r_{1} (\hat{w}, \hat{Ω}) \sqrt{(1 - \cos {(ω_{1} τ_{HH}))}^{2} + \sin^{2} (ω_{1} τ_{HH})} \cos (ω_{1} t + ϕ_{1}) \\ + r_{2} (\hat{w}, \hat{Ω}) \sqrt{(1 - \cos {(ω_{2} τ_{HH}))}^{2} + \sin^{2} (ω_{2} τ_{HH})} \cos (ω_{2} t + ϕ_{2}) = 1, \end{aligned}

5.9

where ϕ₁ and ϕ₂ are certain phase shifts. If there exists any t for which (5.9) holds, then loss of contact takes place and the periodic or the quasi-periodic solutions disappear. The smallest amplitude for which loss of contact might occur is obtained by substituting

\cos (ω_{1} t + ϕ_{1}) = 1

and

\cos (ω_{2} t + ϕ_{2}) = 1

. This motivates the introduction of the scaled amplitude

r := r_{1} (\hat{w}, \hat{Ω}) \sqrt{2 (1 - \cos (ω_{1} τ_{HH}))} + r_{2} (\hat{w}, \hat{Ω}) \sqrt{2 (1 - \cos (ω_{2} τ_{HH}))},

5.10

and r=1 indicates loss of contact. If r>1, then the periodic and the quasi-periodic solutions disappear via the B³. Note that the quasi-periodic solution is located on an invariant torus. Taking

\cos (ω_{1} t + ϕ_{1}) = 1

and

\cos (ω_{2} t + ϕ_{2}) = 1

implies that the solution is assumed to be dense where the torus touches the switching surface (5.8).

Taking r=1, substituting the higher-order estimation of the amplitude of the periodic solutions given by (5.3) with ${\hat{a}}_{m n}$ instead of a_mn, we can obtain the boundary where the periodic solutions disappear via a B³ in the form

{\hat{w}}_{P_{k}}^{B^{3}} (\hat{Ω}) = - \frac{γ_{k 2} \hat{Ω} + a_{k k} / 2 (1 - \cos (ω_{k} τ_{HH}))}{γ_{k 1} + a_{k k} / 2 w_{HH} (1 - \cos (ω_{k} τ_{HH}))} .

5.11

The boundary

{\hat{w}}_{QP}^{B^{3}} (\hat{Ω})

where the quasi-periodic solution vanishes can be determined the same way by symbolic algebra (the resulting formula is too long to be presented here). If any of the unstable periodic or quasi-periodic solutions exists, a large-amplitude attractive solution coexists with the stable equilibrium. Therefore, the boundary of the bistable region is determined by those B³ curves, which are the farthest from the linear stability boundary.

The results are summarized in figure 3 where the double Hopf points (DH) at the intersections of the first and the second, the second and the third, and the fourth and the fifth lobes are considered in rows (a), (b) and (c), respectively. In each panel, the analytical results are indicated by solid lines. The left panels show the stability charts of the system. The right panels show the corresponding bifurcation diagrams with the amplitude r of the periodic orbits P₁, P₂ and the quasi-periodic orbit QP as a function of the bifurcation parameter w. Parameter Ω is fixed to Ω=1.0095, Ω=0.505 and Ω=0.2527 in rows (a), (b) and (c), respectively, which correspond to the dotted vertical lines in the left panels. In the left panels, solid thick blue and cyan lines show the intersecting stability lobes given by (3.3), whereas their tangents (Lines 1 and 2) are indicated by thin blue and cyan lines. Note that in the region depicted, Line 2 overlaps with the solid thick cyan stability lobe. In the right panels, the periodic orbits P₁ and P₂ initiate from the (approximate) Hopf bifurcation points H(P₁) and H(P₂), which correspond to the intersections of the dotted line with Lines 1 and 2 in the left panels of figure 3, respectively. The periodic orbits undergo a torus bifurcation along the red and purple lines (Lines 3 and 4) in the left panels. In the right panel of figure 3a, the quasi-periodic solution QP is born from the periodic solution P₁ through a torus bifurcation at T(P₁), which corresponds to the intersection of the dotted line and Line 3 in the left panel. In figures 3b,c, the quasi-periodic orbit QP arises from the other periodic solution P₂ through a torus bifurcation at T(P₂), see also the intersection of the dotted line and Line 4 in the left panels.

The bifurcation diagrams in the right panels of figure 3 are valid only up to points B³(P₁), B³(P₂) and B³(QP) at r=1, where (2.1) becomes invalid due to loss of contact between the tool and the workpiece. The branches are invalid for r>1, the periodic and the quasi-periodic solutions vanish there and large-amplitude stable solutions are born via a non-smooth fold. In the left panels, solid green and orange lines show the loci of B³ given by (5.11) for the periodic orbits P₁ and P₂, respectively. The B³ curve for the quasi-periodic solution is shown by a black line. Dark grey shading shows the regions where the periodic orbits indeed exist, that is, without losing contact with the workpiece. Similarly, darker grey shading indicates the region where the quasi-periodic solution exists. Light grey shading indicates the region where no periodic or quasi-periodic orbits exist due to loss of contact. Only the stable equilibrium exists here, hence the system is globally stable. In the regions with dark grey and darker grey shading, the system is bistable. Based on the figure, the quasi-periodic orbit occurs in a narrow range of parameters and does not affect the boundary of the globally stable region. Consequently, even close to the intersection of the stability lobes, the globally stable region can be determined by the formulae obtained from the analysis of the single Hopf bifurcation [11].

We calculated the amplitude of solutions by numerical continuation as well, see the dashed lines in figure 3. We used DDE-Biftool [40] for the computation of the periodic orbits, and used the algorithm in [15,41] to compute the quasi-periodic orbit. The loci of torus and B³s have also been continued in two parameters, see the dashed lines in the left panels. Note that the analytical results are only approximations and the numerical results can be considered as exact ones. We can see that the analytical results give only the tangents of the numerical torus bifurcation branches. The agreement between the analytical and numerical results is better in the right panels of figure 3b and 3c than in figure 3a. At the intersection of the first and the second stability lobes in figure 3a, the analytical results for the torus bifurcation curves are valid only in the close vicinity of the double Hopf point, then the numerical branches deviate. Therefore, the analytical bifurcation diagram in the right panel is not accurate even qualitatively: the quasi-periodic solution is born from the other periodic solution as in the numerical case. The numerical counterparts of the B³ branches B³(P₁) and B³(P₂) are also shown by dashed green and orange lines in the left panels of figure 3. Again, the analytical results give only the tangents of the numerical branches. Note that the solid and dashed orange lines overlap. The dashed green and orange lines can also be obtained by analytical formulae using the results of [11,45] on the single Hopf bifurcation, which considers the entire stability lobes not only their intersection. In addition, numerical continuation showed that even a period doubling bifurcation of the periodic solution P₁ can occur, see the dashed thin black line PD(P₁) in the left panel of figure 3a. However, this branch lies beyond the Big Bang Bifurcation B³(P₁), hence the period doubling does not show up when considering loss of contact between the tool and the workpiece.

Finally, figure 4 shows a qualitative picture about the trajectories of the infinite-dimensional non-smooth system including loss of contact between the tool and the workpiece. The infinite-dimensional phase space is illustrated by three axes: x and x′, which are the general coordinates of the corresponding delay-free system, and $x_{\infty}$ , which is introduced only for illustration purposes to represent the other infinitely many dimensions as a cumulative effect of the past. Figure 4a shows the case where only a single periodic orbit (cyan line) exists in the vicinity of the equilibrium. Here, the trajectories which lie within the basin of attraction of the equilibrium tend towards the equilibrium, see the solely green trajectory. However, the basin of attraction is finite and trajectories lying outside spiral outwards, see the red line. In this case, the vibrations amplify until the trajectory hits the switching surface (5.8) indicated by hatches at point A. Here, the tool jumps out of the workpiece. When the tool is out of contact, the infinite-dimensional system reduces to a two-dimensional one as shown by the black line. Temporarily, the system behaves as a damped free oscillator, which is illustrated by the trajectories spiralling towards the origin in the plane (x,x′) (see the green line). The amplitude of the vibrations decreases until the tool gets back into the workpiece at point B. The system becomes infinite-dimensional again (black line) and the trajectory diverges from the origin again (red line).

Figure 4. — Figure 4. Trajectories of the cutting process in the neighbourhood of a single periodic orbit (a); trajectories where two periodic orbits coexist with a quasi-periodic one (b). (Online version in colour.)

Figure 4b shows the case where two periodic solutions (blue and cyan lines) and a quasi-periodic solution (purple line) coexist, cf. Sector IV in figure 2. The green trajectories lie within the basin of attraction of the equilibrium, hence they spiral towards the origin. The red trajectory lies outside and diverges from the origin. Note that this red trajectory is only a possible example, and it can be considered as the representation of the red lines in Sector IV of figure 2b. The solution initiates in the vicinity of one periodic solution, cf. the upper red curve in figure 2b. The trajectory spirals outwards, wraps around the quasi-periodic solution, and runs until it hits the switching surface (5.8) at point A, see the hatches in figure 4b and the dashed line in figure 2b. Here, the tool gets out of the workpiece (black line), the system becomes two-dimensional and restricted to plane (x,x′), where the trajectory spirals towards the origin (green line). Then, the tool enters the workpiece at point B (black line) and the solution continues in the vicinity of the other periodic orbit (red line), cf. the lower red curve in figure 2b. The vibrations start amplifying again until the contact between the tool and the workpiece is lost at point C (black line). A damped motion (green line) follows in the plane (x,x′) until point D is reached. The tool gets back into the workpiece again (black line) and, once more, the infinite-dimensional dynamics prevails. The whole trajectory may even describe chaotic motion. The chaos can be transient, since it can occur that, after a while, the trajectory penetrates the basin of attraction of the equilibrium when the tool gets back into the workpiece. Then, at last, the system becomes globally stable and the trajectory approaches the origin.

6. Conclusion

In this paper, we investigated the dynamics of machine tool vibrations in the vicinity of a double Hopf bifurcation point. Via centre manifold reduction, we derived the normal form coefficients of the four-dimensional centre subsystem, which determines the long-term dynamics in the vicinity of the double Hopf point. Based on the normal form coefficients, we identified the possible topologies in the phase space of the four-dimensional centre subsystem. We found the previously unexplored Case Ia of [17] and also the Case Ib that was predicted in [27], and showed that unstable periodic and quasi-periodic solutions coexist with the linearly stable equilibrium in certain parameter regions. Taking loss of contact between the tool and the workpiece into account, we derived the approximate boundaries of the globally stable region in closed form. We verified the results by numerical continuation using DDE-Biftool [40] and the algorithm in [15,41], which even revealed the existence of period doubling bifurcation of periodic orbits. Finally, we assessed qualitatively the non-smooth global dynamics in the infinite-dimensional phase space and pointed out the possibility of transient chaotic motions.

Data accessibility

This article has no additional data.

Authors' contributions

T.G.M. was responsible for the analytical bifurcation analysis, Z.D. carried out the numerical continuation. G.S. and T.I. supervised the project. All authors contributed to writing the paper and gave final approval for publication.

Competing interests

We declare we have no competing interests.

Funding

This work has been supported by the ÚNKP-16-3-I. New National Excellence Program of the Ministry of Human Capacities. This work has received funding from the János Bolyai Research Scholarship of MTA (BO/00589/13/6). The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC advanced grant agreement no. 340889.

Footnotes

© 2017 The Author(s)

Published by the Royal Society. All rights reserved.

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences cover image

November 2017

Volume 473Issue 2207

Article Information

DOI:https://doi.org/10.1098/rspa.2017.0502
PubMed:29225502
Published by:Royal Society
Online ISSN:1471-2946

History:

Manuscript received22/07/2017
Manuscript accepted17/10/2017
Published online15/11/2017
Published in print30/11/2017

License:

© 2017 The Author(s)

Published by the Royal Society. All rights reserved.

Citations and impact

Keywords

Subjects

mechanical engineering