- Email: [email protected]
A Nonlinear Delay-Differential Equation with Harmonic Excitation
Proceedings, Proceedings, 14th 14th IFAC IFAC Workshop Workshop on on Time Time Delay Delay Systems Systems Pesti Vigadó, Vigadó, Budapest, Budapest, Hungary, June June 28-30, 2018 Proceedings, 14th IFAC Workshop on Time Delay Systems Pesti Hungary, 28-30, 2018 Available online at www.sciencedirect.com Pesti Vigadó, Budapest, Hungary, June 28-30, 2018 Proceedings, 14th IFAC Workshop on Time Delay Systems Pesti Vigadó, Budapest, Hungary, June 28-30, 2018
ScienceDirect
A A A A
IFAC PapersOnLine 51-14 (2018) 224–229
Nonlinear Nonlinear Delay-Differential Delay-Differential Equation Equation Nonlinear Delay-Differential Equation with Harmonic Excitation with Harmonic Excitation Nonlinear Delay-Differential Equation with Harmonic Excitation ∗ ∗ with Harmonic Excitation J´ a a a J´ anos nos Lelkes Lelkes ∗ ,, Tam´ Tam´ ass Kalm´ Kalm´ ar-Nagy r-Nagy ∗
J´ anos Lelkes ∗ , Tam´ as Kalm´ ar-Nagy ∗ ∗ Department of Fluid Mechanics, Faculty of anos Lelkes , Tam´ as Kalm´ ar-Nagy ∗ Engineering, DepartmentJ´ of Fluid Mechanics, Faculty of Mechanical Mechanical Engineering, Department of Fluid Mechanics, Faculty of Mechanical Engineering, Budapest University of Technology and Economics, 1111, Budapest University of Technology and Economics, Budapest Budapest 1111, ∗ Budapest University Hungary (e-mail: [email protected], of Technology and Economics, Budapest 1111, Department Hungary of Fluid Mechanics, Faculty of Mechanical Engineering, (e-mail: [email protected], [email protected]). Hungary (e-mail: [email protected], Budapest University of Technology and Economics, Budapest 1111, [email protected]). [email protected]). Hungary (e-mail: [email protected], [email protected]). Abstract: Abstract: A A machining machining tool tool can can be be subject subject to to different different kinds kinds of of excitations. excitations. The The forcing forcing may may Abstract: A machining tool can be subject to different kinds of excitations. The forcing have external external sources sources (such (such as as rotating rotating imbalance or misalignment misalignment of the the workpiece) workpiece) or it it may can have imbalance or of or can arise from cutting process itself (e.g. chip formation). We the tool have sources rotating imbalance or misalignment of the workpiece) or it may can Abstract: A machining toolascan be subject to different kinds of investigate excitations. Theclassical forcing arise external from the the cutting(such process itself (e.g. chip formation). We investigate the classical tool vibration model which is a delay-differential equation with a quadratic and cubic nonlinearity arise from the cutting process itself (e.g. chip formation). We investigate the classical tool have external sources (such as rotating imbalance or misalignment of the workpiece) or it can vibration model which is a delay-differential equation with a quadratic and cubic nonlinearity and periodic forcing. The method of multiple scales gave an excellent approximation of vibration model which is a delay-differential equation with a quadratic and cubic nonlinearity arise from the cutting process itself (e.g. chip formation). We investigate the classical tool and periodic forcing. The method of multiple scales gave an excellent approximation of the the and periodic forcing. ofhere multiple scales to gave excellent approximation of the solution. The resonance curves found here areequation similar to those for the Duffing-equation, Duffing-equation, having vibration model whichThe iscurves a method delay-differential with aanquadratic and cubic nonlinearity solution. The resonance found are similar those for the having a We Hopf and saddle-node bifurcations. solution. Thecharacteristic. resonance curves here are similar those the Duffing-equation, and periodic forcing. The method multiple scales gave an for excellent approximation having of the a hardening hardening characteristic. We found foundofsubcritical subcritical Hopf to and saddle-node bifurcations. a hardening We found subcritical Hopf to and saddle-node bifurcations. solution. Thecharacteristic. resonance curves here are similar those for the Duffing-equation, having © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. aKeywords: hardeningdelay-differential characteristic. We found subcritical and scales, saddle-node bifurcations. equation, method multiple harmonic excitation, Keywords: delay-differential equation, method of of Hopf multiple scales, harmonic excitation, Hopf Hopf excitation, Hopf bifurcation, primary resonance Keywords: delay-differential equation, method of multiple scales, harmonic bifurcation, primary resonance bifurcation,delay-differential primary resonance Keywords: equation, method of multiple scales, harmonic excitation, Hopf 2 bifurcation, primary resonance 1. INTRODUCTION (3) + εx + O x (t) = x 0 (t 0 ,, t 1) 1 (t 0 ,, t 1) 1. INTRODUCTION εε2 ,, (3) (t t ) + εx (t t ) + O x (t) = x 0 0 1 1 0 1 1. INTRODUCTION (3) x (t) = x0 (t0 , t1 ) + εx1 (t0 , t1 ) + O ε2 , where the timescales are defined as where the timescales are defined as 2 Delay-differential equations (DDEs) are important in 1. INTRODUCTION , (3) , t1 )defined + εx1 (tas (t) = x0 (t0are Delay-differential equations (DDEs) are important in where thex timescales 0 , t1 ) + O ε Delay-differential equations and (DDEs) areFor important in many science. example t, tt1 = εt. (4) tt0 = many areas areas of of engineering engineering and science. For example in = t, = εt. (4) 0 1 where the timescalestare defined many areas of engineering science. exampleand in t1 =asεt. (4) stability analysis, computational techniques (symbolic Delay-differential equations and (DDEs) areFor important in With 0 = t, stability analysis, computational techniques (symbolic and With the the differential differentialtoperators operators stability analysis, computational techniques numerical), automotive engineering, manufacturing, neu(symbolic and many areas of engineering and science. For example in t1 = εt. (4) 0 = t, numerical), automotive engineering, manufacturing, neu- With the differential operators ∂ numerical), automotive engineering, manufacturing, roscience, and control theory (see Balachandran et al. ∂ , D1 = ∂ ∂ , stability computational techniques (symbolic and roscience,analysis, and control theory (see Balachandran etneual. With the differential = (5) D 0 operators (5) D0 = ∂ , D1 = ∂ 1 , roscience, and control theory (see Balachandran al. [2009]). numerical), automotive engineering, manufacturing,etneu[2009]). ∂t00 , D1 = ∂t ∂t1 , (5) D0 = ∂t [2009]). ∂ ∂ ∂t ∂t roscience, andthe control theory (seeform Balachandran et al. time differentiation 0 written as 1 can be In this paper nondimensional of a regenerative = , D = , (5) D time differentiation can be written as 0 1 In this paper the nondimensional form of a regenerative [2009]). ∂t be0 written as ∂t21 In this paper the nondimensional of a model regenerative d one-degree-of-freedom machine tool toolform vibration model in the the time differentiation d =can one-degree-of-freedom machine vibration in D (6) + εD D00be (6) + written εD11 + +O Oas εε22 ,, one-degree-of-freedom machine toolform vibration in the time differentiation d =can case of orthogonal cutting is derivation In the nondimensional of (for a model regenerative dt casethis of paper orthogonal cutting is investigated investigated (for derivation dt = D , (6) + εD + O ε 0 1 case of orthogonal cutting is investigated (for derivation see Kalm´ a r-Nagy et al. [2001]) d dt one-degree-of-freedom machine tool vibration model in the 2 see Kalm´ ar-Nagy et al. [2001]) and similarly and similarly = D , (6) + εD + O ε 0 1 see Kalm´ ar-Nagy etcutting al. [2001]) 2 (for derivation 3 case of2ζorthogonal is qinvestigated and similarly d22dt 2 − (x − xτ )3 x) + (x − x x ¨ + x ˙˙ + x = p (x τ − τ) − x) + q (x − x ) − (x − x ) x ¨ + 2ζ x + x = p (x τ [2001]) τ 2 τ 3 seex r-Nagy (7) =D D022 + + 2εD 2εD0 D D1 + O ε22 . and similarly d ¨ Kalm´ + 2ζ x˙a+ x = pet(xal. τ − x), + q (x − xτ ) − (x − xτ ) d222 = (7) 0 1 + O ε2 . +A 0 dt 2 +A cos cos (ωt) (ωt) , 2 3 dt (7) = D + 2εD 0 D1 + O ε . 22 0 x ¨ + 2ζ x˙ + x = (xτ (ωt) − x), + q (x − xτ ) − (x − xτ ) (1) +Ap cos d dt 2 2 (1) . (7) = D + 2εD D + O ε 0 1 0 where x x is is the the+A tool displacement (xτ = = x x (t (t − − ττ )) is is (1) its cosdisplacement (ωt) , dt2 where tool (x its τ 2. LINEAR LINEAR STABILITY STABILITY ANALYSIS ANALYSIS delayed ττ > ζζ τ> where x value is thewith tooldelay displacement − τrelative ) is (1) its 2. delayed value with delay > 0), 0), (x > =0 0 xis is(tthe the relative OF THE UNFORCED SYSTEM 2. LINEAR STABILITY ANALYSIS damping factor, p > 0 is the nondimensional cutting force, delayed value with delay τ > 0), ζ > 0 is the relative OF THE UNFORCED SYSTEM where x is the tool displacement (x = x (t − τ ) is its τ damping factor, p > 0 is the nondimensional cutting force, OF THE UNFORCED SYSTEM 2. LINEAR STABILITY ANALYSIS qqdelayed > 0 is the coefficient of nonlinearity, A is the amplitude damping factor, p > 0 is the nondimensional cutting force, value with delay τ > 0), ζ > 0 is the relative > 0 is the coefficient of nonlinearity, A is the amplitude OF THE UNFORCED SYSTEM qof > 0 is the coefficient of nonlinearity, A is the amplitude The stability stability analysis analysis of of the the x x= = 00 solution solution of the the linearized linearized of the forcing forcing and ω 0is isisits itsthe frequency. damping factor, p> nondimensional cutting force, The of the and ω frequency. equation The stability analysis of the x = 0 solution of the linearized of the forcing and ω is its frequency. qPerturbation > 0 is the coefficient of nonlinearity, A is the amplitude equation methods were successfully Perturbation methods were successfully applied applied to to delaydelay- The equation stability analysis of˙ the x =p0(xsolution of the linearized of the forcing and ω is its frequency. x ¨¨ + (8) Perturbation methods were applied toee method delaydifferential equations (such assuccessfully the Linstedt-Poincar´ Linstedt-Poincar´ method x + 2ζ 2ζ x x˙ + +x x= = p (xττ − − x) x) ,, (8) differential equations (such as the equation x ¨ + 2ζ x ˙ + x = p (x − x) , (8) τ (Casal and Freedman [1980]), the method of multiple differential equations (such as the Linstedt-Poincar´ e method was performed performed in, in, for for example, example, Hanna Hanna and and Tobias Tobias [1974], [1974], Perturbation methods were successfully applied delay- was (Casal and Freedman [1980]), the method of to multiple x ¨ +for 2ζ example, x˙ + x = pHanna (xτ − x) , Tobias (8) scales et [1997], [2008], et (Casal andequations Freedman [1980]), methodOztepe of multiple St´ aan stability boundaries of performed in, and [1974], differential asNayfeh thethe Linstedt-Poincar´ e method scales (Nayfeh (Nayfeh et al. al. (such [1997], Nayfeh [2008], Oztepe et al. al. was St´eep´ p´ n [1989]. [1989]. The The stability boundaries of (8) (8) is is obtained obtained [2015]) or the combination of the method of multiple scales scales (Nayfeh et al. [1997], Nayfeh [2008], Oztepe et al. by substituting the trial solution x(t) = C exp(iΩt) into St´ e p´ a n [1989]. The stability boundaries of (8) is obtained (Casal and Freedman [1980]), the method of multiple was performed in, for example, Hanna and Tobias [1974], [2015]) or the combination of the method of multiple scales by substituting the trial solution x(t) = C exp(iΩt) into [2015]) or the combination of the method of multiple scales and Linstedt-Poincar´ Linstedt-Poincar´ method (Pakdemirli and Karahan Karahan The stability diagram in Figure 11 is in byep´ substituting the trial solution exp(iΩt) into scales (Nayfeh et al. ee[1997], Nayfeh [2008], Oztepe et al. (8). St´ an [1989]. The stability ofC(8) is parametobtained and method (Pakdemirli and (8). The stability diagram in boundaries Figurex(t) is=given given in paramet[2010]). (8). The in stability in Figure 1 is=given in parametand Linstedt-Poincar´ e method (Pakdemirli and Karahan Kalm´ aadiagram r-Nagy [2002] as [2015]) by form substituting the trial solution C exp(iΩt) into [2010]).or the combination of the method of multiple scales ric ric form in Kalm´ r-Nagy [2002] as x(t) We approximate solution by the method [2010]). [2002] Kalm´adiagram r-Nagy and Linstedt-Poincar´ method (Pakdemirli (8).form The in stability in2 Figure is given in paramet2 as 1 We will will approximatee the the solution by using usingand theKarahan method ric 2 2 2 Ω of scales. following we will that We will approximate the method 11 − −[2002] Ω22 2 + +as4ζ 4ζ 22 Ω Ω22 , [2010]). ar-Nagy of multiple multiple scales. In In the the solution followingby weusing will assume assume that ric form in Kalm´ p = (9) p = (9) 2 1 − Ω + 4ζ of multiple scales. In the following we will assume that damping is small, the nonlinearity and forcing are weak We will approximate solution by theare method 2 − 1) Ω , damping is small, the the nonlinearity andusing forcing weak (Ω p = 1 −22Ω(Ω (9) 2 22 − 1) 2 2 , + 4ζ Ω (see Kalm´ a r-Nagy [2002]). In particular 2 (Ω − 1) damping is small, the nonlinearity and forcing are weak of scales. [2002]). In the following we will assume that (seemultiple Kalm´ ar-Nagy In particular 2 p = , (9) Ω 22 − (see Kalm´ In ,particular damping isar-Nagy the∼ Ω222(Ω ζ,small, p, q, q, [2002]). A ∼nonlinearity O (ε) (ε) ωand =1 1 forcing + εσ, εσ, are weak (2) −211−1) jπ − arctan ,, jj = 1, 2, .. .. .. ,, (10) ττ = ζ, p, A O , ω = + (2) Ω jπ − arctan = 1, 2, (10) = 2 − 1 (see Kalm´ ar-Nagy In ,particular ζ, p,and q, [2002]). A∼ O (ε) ω = 1 + εσ, We also (2) , j = 1, 2, . . . , where (10) Ω jπ − arctan 2ζΩ 2ζΩ τ=Ω 2 where εε 1 1 and σ σ is is the the detuning detuning frequency. frequency. We also Ω2ζΩ 2 − 1 ‘lobe’ and Ω > 1 because Ω ζ, p, q, A ∼ O (ε) , ω = 1 + εσ, (2) where ε 1 and σ is the detuning frequency. We also where j corresponds to the jth assume that the solution of (1) can be well approximated jπ − arctan , j = 1, 2, . . . , (10) τ = assume that the solution of (1) can be well approximated where j corresponds to the jth ‘lobe’ and Ω > 1 because Ω 2ζΩ p by expansion assume that solution of (1) can befrequency. well approximated where 1the and σ is the detuning We also where p> > 0. 0. j corresponds to the jth ‘lobe’ and Ω > 1 because by the the εtwo-scale two-scale expansion p > 0. j corresponds to the jth ‘lobe’ and Ω > 1 because by the two-scale expansion assume that the solution of (1) can be well approximated where p > 0. by Elsevier Ltd. All rights reserved. by the two-scale expansion Copyright © 2018, 2018 IFAC 224Hosting 2405-8963 © IFAC (International Federation of Automatic Control) ∗ ∗ ∗
Copyright © 2018 IFAC 224 Peer review©under of International Federation of Automatic Copyright 2018 responsibility IFAC 224Control. 10.1016/j.ifacol.2018.07.227 Copyright © 2018 IFAC 224
2018 IFAC TDS Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
225
Now we divide (20) by eiβ(t1 ) and introduce φ (t1 ) = σt1 − β (t1 ) ,
(21)
to get 1 0 = − Aeiφ(t1 ) + i (D1 α (t1 ) + i (σ − D1 φ(t1 ))) 2 1 + α (t1 ) p − e−iτ p + 2iζ . 2
(22) Separating the real and imaginary parts of equation (22) results in two ordinary differential equations describing the evolution of the amplitude and phase A p (23) D1 α = − ζ + sin τ α + sin φ, 2 2 αD1 φ =
Fig. 1. The linear stability chart of (1) At the minima (‘notches’) of the stability lobes, Ω, p, τ assume the particularly simple forms (11) Ωcrit = 1 + 2ζ, (12) pcrit = 2ζ (ζ + 1) , 2 jπ − arctan √ 1 1+2ζ √ , j = 1, 2, . . . (13) τcrit = 1 + 2ζ 3. THE LINEAR DELAY-DIFFERENTIAL EQUATION WITH HARMONIC FORCING Here we study the harmonically forced linear equation (14) x ¨ + 2ζ x˙ + x = p (xτ − x) + A cos (ωt) . With assumption (2), substituting the differential operators (6) and (7) into (14) and equating like powers of ε one obtains (15) ε0 : D02 x0 (t0 , t1 ) + x0 (t0 , t1 ) = 0, D02 x1 (t0 , t1 ) + x1 (t0 , t1 ) = −2D0 D1 x0 (t0 , t1 ) ε : −2ζD0 x0 (t0 , t1 ) + p [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] +A cos (t0 + σt1 ) . (16) Solving (15) for x0 (t0 , t1 ), knowing that the slow temporal variable t1 is implicit in the constants of integration, we get x0 (t0 , t1 ) = a (t1 ) eit0 + a ¯ (t1 ) e−it0 , (17) We substitute this solution into (16) and eliminate the secular terms which would give rise to unbounded terms 1 0 = − Aeiσt1 + 2iD1 a (t1 ) + a (t1 ) p − e−iτ p + 2iζ . 2 (18) Writing a (t1 ) in a polar form (where α (t1 ) is the amplitude and β (t1 ) is the phase) 1 1 a ¯ (t1 ) = α (t1 ) e−iβ(t1 ) , (19) a (t1 ) = α (t1 ) eiβ(t1 ) , 2 2 then substituting into (18), we get 1 0 = − Aeiσt1 + i (D1 α (t1 ) + iα (t1 ) D1 β (t1 )) eiβ(t1 ) 2 1 + α (t1 ) eiβ(t1 ) p − e−iτ p + 2iζ . 2 (20) 1
225
A 1 (p cos τ + 2σ − p) α + cos φ. 2 2
(24)
To get the amplitude and phase of the steady-state primary resonance (i.e. the fixed point of (23) and (24)) we set the left-hand sides of (23) and (24) to zero and solve the resulting algebraic equations (by eliminating the trigonometric terms from them). For the amplitude α∗ and phase φ∗ we get A α∗ = 2 (p (1 − cos τ ) − 2(ω − 1)) + 4 ζ +
p 2
sin τ
α p . φ∗ = arcsin 2 ζ + sin τ 2 A
2 .
(25)
(26)
The stability of the fixed point (α∗ , φ∗ ) is determined by the Jacobian p A ∗ cos φ −ζ − 2 sin τ 2 J = (27) . A A sin φ∗ − ∗ 2 cos φ∗ − 2α 2α The trace and determinant of J are
T = −2ζ − p sin τ,
(28)
2 1 p ∆ = ζ + sin τ + (p(1 − cos τ ) − 2ω + 2)2 . 2 4
(29)
T 2 − 4∆ = −(p(1 − cos τ ) − 2ω + 2)2 ≤ 0,
(30)
The determinant ∆ is always non-negative. The discriminant
therefore the fixed point is a stable spiral when T < 0, a center when T = 0 (at pbif = −2 sinζ τ ), and an unstable spiral when T > 0. We chose p as the bifurcation parameter. The phase portraits of the “slow-flow” (Eqs. (23) and (24)) are shown in Figs. 2, 3, 4 for different p values.
2018 IFAC TDS 226 Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
3.1 Numerical results The direct numerical solution of (14) was determined in Mathematica using the Dormand-Prince method with an initial function α∗ cos(ωt + φ∗ ), t ∈ [−τ, 0], where α∗ and φ∗ were determined from (25) and (26). The following numerical values were chosen for the simulation ζ = 0.01, A = 0.01. (31) We have chosen p values at the first ‘notch’ (j = 1). Using (12) and (13) gives (32) pcrit = 0.0202, τcrit = 4.676. Time series x (t) for (33) p = {0.75pcrit , 0.90pcrit , 0.99pcrit , 1.05pcrit } are shown in Figure 5.
Fig. 2. Stable spiral p < pbif
(a) p = 0.75pcrit , stable
(b) p = 0.90pcrit , stable
(c) p = 0.99pcrit , stable
(d) p = 1.05pcrit , unstable
.
Fig. 3. Center p = pbif
. Fig. 5. Time series of the numerical solution (ω = 1.02) The grey areas indicate the time interval [300τ, 450τ ] where the steady state vibration amplitude of the numerical solution was determined as (34) αnum = max|x(t)|, t ∈ [300τ, 450τ ]. In the case of p = 0.90pcrit and p = 0.99pcrit the vibration becomes modulated, so the amplitude of the vibration also varies between a certain range during the simulation. For a better comparison of the method of multiple scales and numerical results we introduce the method of multiple scales amplification factor G∗ and numerical amplification factor Gnum as α∗ αnum , Gnum := . (35) G∗ := A A The comparison of the amplification factors G∗ and Gnum is shown in Figure 6, where the solid line indicates the method of multiple scales (MMS) results, the dots the numerical results.
Fig. 4. Unstable spiral p > pbif 226
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
227
tial equations describing the evolution of the amplitude and phase τ p A τ D1 α = − ζ + sin τ α − 3q cos sin3 α3 + sin φ, 2 2 2 2 (41) 1 A τ αD1 φ = (p cos τ + 2σ − p) α − 3q sin4 α3 + cos φ. 2 2 2 (42) To get the amplitude of the steady-state primary resonance i.e. the fixed points of (41) and (42) we set the lefthand sides of (41) and (42) to zero and solve the resulting algebraic equations (by eliminating the trigonometric terms from them) τ 2 1 p (1 − cos τ ) − 2σ + 6qα∗ 2 sin4 4 2 2 A2 p τ 3 τ ∗2 α∗ 2 − = 0. + ζ + sin τ + 3qα cos sin 2 2 2 4 (43) p 2 τ 3 τ ∗3 ∗ ∗ ζ + sin τ α + 3q cos sin α φ = arcsin . A 2 2 2 (44) Equation (43) can be solved for α∗ resulting 1, 2 or 3 pair (±) of real roots.
MMS
(a) p = 0.75pcrit MMS
(b) p = 0.90pcrit MMS
Stability of the fixed points (α∗ , φ∗ ) is determined by the Jacobian τ p A 3 τ ∗2 ∗ cos φ −ζ − 2 sin τ − 9qα cos 2 sin 2 2 J = . τ A A −6α∗ q sin4 − ∗ 2 cos φ∗ − ∗ sin φ∗ 2 2α 2α (45) The stability of the fixed points can be investigated by calculating the determinant ∆, the trace T and the discriminant T 2 − 4∆ of J .
(c) p = 0.99pcrit
Fig. 6. The amplification factor G as a function of ω 4. THE HARMONICALLY FORCED NONLINEAR DELAY-DIFFERENTIAL EQUATION Substituting the differential operators (6) and (7) into (1) and equating like powers of ε one obtains ε0 : D02 x0 (t0 , t1 ) + x0 (t0 , t1 ) = 0, (36) D02 x1 (t0 , t1 ) + x1 (t0 , t1 ) = −2D0 D1 x0 (t0 , t1 ) −2ζD0 x0 (t0 , t1 ) + p [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] ε1 : 2 +q [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] 3 +q [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] + A cos (t0 + σt1 ) . (37) Solving (36) for x0 (t0 , t1 ) yields x0 (t0 , t1 ) = a (t1 ) eit0 + a ¯ (t1 ) e−it0 . (38) Then we substitute this solution into (37) and eliminate the secular terms which would give rise to unbounded terms 1 0 = − Aeiσt1 + 2iD1 a (t1 ) 2 3 ¯ (t1 ) . +a (t1 ) p − e−iτ p + 2iζ − 3e−2iτ eiτ − 1 qa (t1 ) a (39) Writing a (t1 ) in a polar form (where α (t1 ) is the amplitude and β (t1 ) is the phase) 1 1 a (t1 ) = α (t1 ) eiβ(t1 ), a ¯ (t1 ) = α (t1 ) e−iβ(t1 ). (40) 2 2 then substituting into (39) using (21) and separating the real and imaginary parts results in two ordinary differen227
The nondimensional cutting force coefficient was taken to be p = 0.5pcrit , the coefficient of nonlinearity q = 0.003, forcing amplitude A = 0.01. The quantity ω was chosen as bifurcation parameter, increased from 0.98 to 1.03. 2.0
1.5
MMS Numerical ÈΑÈ
2018 IFAC TDS Budapest, Hungary, June 28-30, 2018
_ SBH _ < SN
1.0
<
0.5
0.0 0.98
0.99
1.00
1.01
1.02
1.03
Ω
Fig. 7. Bifurcation diagram (or resonance curve), where SBH corresponds to subcritical Hopf bifurcation and SN to saddle-node bifurcation.
2018 IFAC TDS 228 Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
Figure 7 shows the amplitude α∗ as a function of forcing frequency ω (the solid line indicates stable solution, the dashed line the unstable solution, both method of multiple scales solution). The stable numerical equlibrium solutions are shown with filled circles, unstable limit cycles are depicted by empty circles. The thin vertical lines in Fig. 7 indicate the ω values (1.00, 1.005, 1.013, 1.02) for which the phase portraits of Figs. (8-11) were generated.
A saddle and a node connected with a heteroclinic orbit can be observed at ω ≈ 1.013, see Figure 10.
Fig. 11. Stable spiral and unstable limit cycle at ω = 1.02 5. CONCLUSIONS
Fig. 8. Stable spiral at ω = 1, the arrowheads indicate two points of the unstable limit cycle
The method of multiple scales gave an excellent approximation of the solution of the nonlinear delay-differential equation (1). The resonance curves found here are similar to those for the Duffing-equation (Nayfeh and Mook [1979]), having a hardening characteristic. We found subcritical Hopf bifurcations and saddle-node bifurcations. Equation (1) without forcing admits a subcritical Hopf bifurcation (Kalm´ar-Nagy [2002]) and the interaction of forcing and Hopf bifurcation may lead to incredibly complex phenomena (Gambaudo [1985]). Plaut and Hsieh (Plaut and Hsieh [1987]) studied a one-DOF mechanical system with delay and excitation and found periodic, chaotic and unbounded responses. Forcing of a Duffing-type equation without delay can also result in complicated bifurcation structure (Sanchez and Nayfeh [1990]; Zavodney et al. [1990]). REFERENCES
Fig. 9. Unstable spiral at ω = 1.005 Subcritical Hopf bifurcation occurs at ω ≈ 1.003, see Figure 8 and 9.
Fig. 10. Stable spiral at ω = 1.013 and two unstable fixed points 228
Balachandran, B., Kalm´ar-Nagy, T., and Gilsinn, D.E. (2009). Delay differential equations. Springer. Casal, A. and Freedman, M. (1980). A Poincar´e-Lindstedt approach to bifurcation problems for differential-delay equations. IEEE Transactions on Automatic Control, 25(5), 967–973. Gambaudo, J.M. (1985). Perturbation of a Hopf bifurcation by an external time-periodic forcing. Journal of differential equations, 57(2), 172–199. Hanna, N. and Tobias, S. (1974). A theory of nonlinear regenerative chatter. ASME Journal of Engineering for Industry, 96(1974), 247–255. Kalm´ar-Nagy, T. (2002). Delay-differential models of cutting tool dynamics with nonlinear and mode-coupling effects. PhD thesis. Kalm´ar-Nagy, T., St´ep´an, G., and Moon, F.C. (2001). Subcritical Hopf bifurcation in the delay equation model for machine tool vibrations. Nonlinear Dynamics, 26(2), 121–142. Nayfeh, A., Chin, C., and Pratt, J. (1997). Applications of perturbation methods to tool chatter dynamics. Dynamics and chaos in manufacturing processes, 193–213.
2018 IFAC TDS Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
Nayfeh, A.H. (2008). Order reduction of retarded nonlinear systems–the method of multiple scales versus centermanifold reduction. Nonlinear Dynamics, 51(4), 483– 500. Nayfeh, A.H. and Mook, D.T. (1979). Forced oscillations of systems having a single degree of freedom. Nonlinear Oscillations, 161–257. Oztepe, G.S., Choudhury, S.R., and Bhatt, A. (2015). Multiple scales and energy analysis of coupled Rayleighvan der Pol oscillators with time-delayed displacement and velocity feedback: Hopf bifurcations and amplitude death. Far East Journal of Dynamical Systems, 26(1), 31. Pakdemirli, M. and Karahan, M.M.F. (2010). A new perturbation solution for systems with strong quadratic and cubic nonlinearities. Mathematical Methods in the Applied Sciences, 33(6), 704–712. Plaut, R. and Hsieh, J.C. (1987). Chaos in a mechanism with time delays under parametric and external excitation. Journal of Sound and Vibration, 114(1), 73–90. Sanchez, N.E. and Nayfeh, A.H. (1990). Prediction of bifurcations in a parametrically excited Duffing oscillator. International Journal of Non-Linear Mechanics, 25(23), 163–176. St´ep´ an, G. (1989). Retarded dynamical systems: stability and characteristic functions. Longman Scientific & Technical. Zavodney, L.D., Nayfeh, A., and Sanchez, N. (1990). Bifurcations and chaos in parametrically excited singledegree-of-freedom systems. Nonlinear Dynamics, 1(1), 1–21.
229
229
ScienceDirect
A A A A
IFAC PapersOnLine 51-14 (2018) 224–229
Nonlinear Nonlinear Delay-Differential Delay-Differential Equation Equation Nonlinear Delay-Differential Equation with Harmonic Excitation with Harmonic Excitation Nonlinear Delay-Differential Equation with Harmonic Excitation ∗ ∗ with Harmonic Excitation J´ a a a J´ anos nos Lelkes Lelkes ∗ ,, Tam´ Tam´ ass Kalm´ Kalm´ ar-Nagy r-Nagy ∗
J´ anos Lelkes ∗ , Tam´ as Kalm´ ar-Nagy ∗ ∗ Department of Fluid Mechanics, Faculty of anos Lelkes , Tam´ as Kalm´ ar-Nagy ∗ Engineering, DepartmentJ´ of Fluid Mechanics, Faculty of Mechanical Mechanical Engineering, Department of Fluid Mechanics, Faculty of Mechanical Engineering, Budapest University of Technology and Economics, 1111, Budapest University of Technology and Economics, Budapest Budapest 1111, ∗ Budapest University Hungary (e-mail: [email protected], of Technology and Economics, Budapest 1111, Department Hungary of Fluid Mechanics, Faculty of Mechanical Engineering, (e-mail: [email protected], [email protected]). Hungary (e-mail: [email protected], Budapest University of Technology and Economics, Budapest 1111, [email protected]). [email protected]). Hungary (e-mail: [email protected], [email protected]). Abstract: Abstract: A A machining machining tool tool can can be be subject subject to to different different kinds kinds of of excitations. excitations. The The forcing forcing may may Abstract: A machining tool can be subject to different kinds of excitations. The forcing have external external sources sources (such (such as as rotating rotating imbalance or misalignment misalignment of the the workpiece) workpiece) or it it may can have imbalance or of or can arise from cutting process itself (e.g. chip formation). We the tool have sources rotating imbalance or misalignment of the workpiece) or it may can Abstract: A machining toolascan be subject to different kinds of investigate excitations. Theclassical forcing arise external from the the cutting(such process itself (e.g. chip formation). We investigate the classical tool vibration model which is a delay-differential equation with a quadratic and cubic nonlinearity arise from the cutting process itself (e.g. chip formation). We investigate the classical tool have external sources (such as rotating imbalance or misalignment of the workpiece) or it can vibration model which is a delay-differential equation with a quadratic and cubic nonlinearity and periodic forcing. The method of multiple scales gave an excellent approximation of vibration model which is a delay-differential equation with a quadratic and cubic nonlinearity arise from the cutting process itself (e.g. chip formation). We investigate the classical tool and periodic forcing. The method of multiple scales gave an excellent approximation of the the and periodic forcing. ofhere multiple scales to gave excellent approximation of the solution. The resonance curves found here areequation similar to those for the Duffing-equation, Duffing-equation, having vibration model whichThe iscurves a method delay-differential with aanquadratic and cubic nonlinearity solution. The resonance found are similar those for the having a We Hopf and saddle-node bifurcations. solution. Thecharacteristic. resonance curves here are similar those the Duffing-equation, and periodic forcing. The method multiple scales gave an for excellent approximation having of the a hardening hardening characteristic. We found foundofsubcritical subcritical Hopf to and saddle-node bifurcations. a hardening We found subcritical Hopf to and saddle-node bifurcations. solution. Thecharacteristic. resonance curves here are similar those for the Duffing-equation, having © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. aKeywords: hardeningdelay-differential characteristic. We found subcritical and scales, saddle-node bifurcations. equation, method multiple harmonic excitation, Keywords: delay-differential equation, method of of Hopf multiple scales, harmonic excitation, Hopf Hopf excitation, Hopf bifurcation, primary resonance Keywords: delay-differential equation, method of multiple scales, harmonic bifurcation, primary resonance bifurcation,delay-differential primary resonance Keywords: equation, method of multiple scales, harmonic excitation, Hopf 2 bifurcation, primary resonance 1. INTRODUCTION (3) + εx + O x (t) = x 0 (t 0 ,, t 1) 1 (t 0 ,, t 1) 1. INTRODUCTION εε2 ,, (3) (t t ) + εx (t t ) + O x (t) = x 0 0 1 1 0 1 1. INTRODUCTION (3) x (t) = x0 (t0 , t1 ) + εx1 (t0 , t1 ) + O ε2 , where the timescales are defined as where the timescales are defined as 2 Delay-differential equations (DDEs) are important in 1. INTRODUCTION , (3) , t1 )defined + εx1 (tas (t) = x0 (t0are Delay-differential equations (DDEs) are important in where thex timescales 0 , t1 ) + O ε Delay-differential equations and (DDEs) areFor important in many science. example t, tt1 = εt. (4) tt0 = many areas areas of of engineering engineering and science. For example in = t, = εt. (4) 0 1 where the timescalestare defined many areas of engineering science. exampleand in t1 =asεt. (4) stability analysis, computational techniques (symbolic Delay-differential equations and (DDEs) areFor important in With 0 = t, stability analysis, computational techniques (symbolic and With the the differential differentialtoperators operators stability analysis, computational techniques numerical), automotive engineering, manufacturing, neu(symbolic and many areas of engineering and science. For example in t1 = εt. (4) 0 = t, numerical), automotive engineering, manufacturing, neu- With the differential operators ∂ numerical), automotive engineering, manufacturing, roscience, and control theory (see Balachandran et al. ∂ , D1 = ∂ ∂ , stability computational techniques (symbolic and roscience,analysis, and control theory (see Balachandran etneual. With the differential = (5) D 0 operators (5) D0 = ∂ , D1 = ∂ 1 , roscience, and control theory (see Balachandran al. [2009]). numerical), automotive engineering, manufacturing,etneu[2009]). ∂t00 , D1 = ∂t ∂t1 , (5) D0 = ∂t [2009]). ∂ ∂ ∂t ∂t roscience, andthe control theory (seeform Balachandran et al. time differentiation 0 written as 1 can be In this paper nondimensional of a regenerative = , D = , (5) D time differentiation can be written as 0 1 In this paper the nondimensional form of a regenerative [2009]). ∂t be0 written as ∂t21 In this paper the nondimensional of a model regenerative d one-degree-of-freedom machine tool toolform vibration model in the the time differentiation d =can one-degree-of-freedom machine vibration in D (6) + εD D00be (6) + written εD11 + +O Oas εε22 ,, one-degree-of-freedom machine toolform vibration in the time differentiation d =can case of orthogonal cutting is derivation In the nondimensional of (for a model regenerative dt casethis of paper orthogonal cutting is investigated investigated (for derivation dt = D , (6) + εD + O ε 0 1 case of orthogonal cutting is investigated (for derivation see Kalm´ a r-Nagy et al. [2001]) d dt one-degree-of-freedom machine tool vibration model in the 2 see Kalm´ ar-Nagy et al. [2001]) and similarly and similarly = D , (6) + εD + O ε 0 1 see Kalm´ ar-Nagy etcutting al. [2001]) 2 (for derivation 3 case of2ζorthogonal is qinvestigated and similarly d22dt 2 − (x − xτ )3 x) + (x − x x ¨ + x ˙˙ + x = p (x τ − τ) − x) + q (x − x ) − (x − x ) x ¨ + 2ζ x + x = p (x τ [2001]) τ 2 τ 3 seex r-Nagy (7) =D D022 + + 2εD 2εD0 D D1 + O ε22 . and similarly d ¨ Kalm´ + 2ζ x˙a+ x = pet(xal. τ − x), + q (x − xτ ) − (x − xτ ) d222 = (7) 0 1 + O ε2 . +A 0 dt 2 +A cos cos (ωt) (ωt) , 2 3 dt (7) = D + 2εD 0 D1 + O ε . 22 0 x ¨ + 2ζ x˙ + x = (xτ (ωt) − x), + q (x − xτ ) − (x − xτ ) (1) +Ap cos d dt 2 2 (1) . (7) = D + 2εD D + O ε 0 1 0 where x x is is the the+A tool displacement (xτ = = x x (t (t − − ττ )) is is (1) its cosdisplacement (ωt) , dt2 where tool (x its τ 2. LINEAR LINEAR STABILITY STABILITY ANALYSIS ANALYSIS delayed ττ > ζζ τ> where x value is thewith tooldelay displacement − τrelative ) is (1) its 2. delayed value with delay > 0), 0), (x > =0 0 xis is(tthe the relative OF THE UNFORCED SYSTEM 2. LINEAR STABILITY ANALYSIS damping factor, p > 0 is the nondimensional cutting force, delayed value with delay τ > 0), ζ > 0 is the relative OF THE UNFORCED SYSTEM where x is the tool displacement (x = x (t − τ ) is its τ damping factor, p > 0 is the nondimensional cutting force, OF THE UNFORCED SYSTEM 2. LINEAR STABILITY ANALYSIS qqdelayed > 0 is the coefficient of nonlinearity, A is the amplitude damping factor, p > 0 is the nondimensional cutting force, value with delay τ > 0), ζ > 0 is the relative > 0 is the coefficient of nonlinearity, A is the amplitude OF THE UNFORCED SYSTEM qof > 0 is the coefficient of nonlinearity, A is the amplitude The stability stability analysis analysis of of the the x x= = 00 solution solution of the the linearized linearized of the forcing forcing and ω 0is isisits itsthe frequency. damping factor, p> nondimensional cutting force, The of the and ω frequency. equation The stability analysis of the x = 0 solution of the linearized of the forcing and ω is its frequency. qPerturbation > 0 is the coefficient of nonlinearity, A is the amplitude equation methods were successfully Perturbation methods were successfully applied applied to to delaydelay- The equation stability analysis of˙ the x =p0(xsolution of the linearized of the forcing and ω is its frequency. x ¨¨ + (8) Perturbation methods were applied toee method delaydifferential equations (such assuccessfully the Linstedt-Poincar´ Linstedt-Poincar´ method x + 2ζ 2ζ x x˙ + +x x= = p (xττ − − x) x) ,, (8) differential equations (such as the equation x ¨ + 2ζ x ˙ + x = p (x − x) , (8) τ (Casal and Freedman [1980]), the method of multiple differential equations (such as the Linstedt-Poincar´ e method was performed performed in, in, for for example, example, Hanna Hanna and and Tobias Tobias [1974], [1974], Perturbation methods were successfully applied delay- was (Casal and Freedman [1980]), the method of to multiple x ¨ +for 2ζ example, x˙ + x = pHanna (xτ − x) , Tobias (8) scales et [1997], [2008], et (Casal andequations Freedman [1980]), methodOztepe of multiple St´ aan stability boundaries of performed in, and [1974], differential asNayfeh thethe Linstedt-Poincar´ e method scales (Nayfeh (Nayfeh et al. al. (such [1997], Nayfeh [2008], Oztepe et al. al. was St´eep´ p´ n [1989]. [1989]. The The stability boundaries of (8) (8) is is obtained obtained [2015]) or the combination of the method of multiple scales scales (Nayfeh et al. [1997], Nayfeh [2008], Oztepe et al. by substituting the trial solution x(t) = C exp(iΩt) into St´ e p´ a n [1989]. The stability boundaries of (8) is obtained (Casal and Freedman [1980]), the method of multiple was performed in, for example, Hanna and Tobias [1974], [2015]) or the combination of the method of multiple scales by substituting the trial solution x(t) = C exp(iΩt) into [2015]) or the combination of the method of multiple scales and Linstedt-Poincar´ Linstedt-Poincar´ method (Pakdemirli and Karahan Karahan The stability diagram in Figure 11 is in byep´ substituting the trial solution exp(iΩt) into scales (Nayfeh et al. ee[1997], Nayfeh [2008], Oztepe et al. (8). St´ an [1989]. The stability ofC(8) is parametobtained and method (Pakdemirli and (8). The stability diagram in boundaries Figurex(t) is=given given in paramet[2010]). (8). The in stability in Figure 1 is=given in parametand Linstedt-Poincar´ e method (Pakdemirli and Karahan Kalm´ aadiagram r-Nagy [2002] as [2015]) by form substituting the trial solution C exp(iΩt) into [2010]).or the combination of the method of multiple scales ric ric form in Kalm´ r-Nagy [2002] as x(t) We approximate solution by the method [2010]). [2002] Kalm´adiagram r-Nagy and Linstedt-Poincar´ method (Pakdemirli (8).form The in stability in2 Figure is given in paramet2 as 1 We will will approximatee the the solution by using usingand theKarahan method ric 2 2 2 Ω of scales. following we will that We will approximate the method 11 − −[2002] Ω22 2 + +as4ζ 4ζ 22 Ω Ω22 , [2010]). ar-Nagy of multiple multiple scales. In In the the solution followingby weusing will assume assume that ric form in Kalm´ p = (9) p = (9) 2 1 − Ω + 4ζ of multiple scales. In the following we will assume that damping is small, the nonlinearity and forcing are weak We will approximate solution by theare method 2 − 1) Ω , damping is small, the the nonlinearity andusing forcing weak (Ω p = 1 −22Ω(Ω (9) 2 22 − 1) 2 2 , + 4ζ Ω (see Kalm´ a r-Nagy [2002]). In particular 2 (Ω − 1) damping is small, the nonlinearity and forcing are weak of scales. [2002]). In the following we will assume that (seemultiple Kalm´ ar-Nagy In particular 2 p = , (9) Ω 22 − (see Kalm´ In ,particular damping isar-Nagy the∼ Ω222(Ω ζ,small, p, q, q, [2002]). A ∼nonlinearity O (ε) (ε) ωand =1 1 forcing + εσ, εσ, are weak (2) −211−1) jπ − arctan ,, jj = 1, 2, .. .. .. ,, (10) ττ = ζ, p, A O , ω = + (2) Ω jπ − arctan = 1, 2, (10) = 2 − 1 (see Kalm´ ar-Nagy In ,particular ζ, p,and q, [2002]). A∼ O (ε) ω = 1 + εσ, We also (2) , j = 1, 2, . . . , where (10) Ω jπ − arctan 2ζΩ 2ζΩ τ=Ω 2 where εε 1 1 and σ σ is is the the detuning detuning frequency. frequency. We also Ω2ζΩ 2 − 1 ‘lobe’ and Ω > 1 because Ω ζ, p, q, A ∼ O (ε) , ω = 1 + εσ, (2) where ε 1 and σ is the detuning frequency. We also where j corresponds to the jth assume that the solution of (1) can be well approximated jπ − arctan , j = 1, 2, . . . , (10) τ = assume that the solution of (1) can be well approximated where j corresponds to the jth ‘lobe’ and Ω > 1 because Ω 2ζΩ p by expansion assume that solution of (1) can befrequency. well approximated where 1the and σ is the detuning We also where p> > 0. 0. j corresponds to the jth ‘lobe’ and Ω > 1 because by the the εtwo-scale two-scale expansion p > 0. j corresponds to the jth ‘lobe’ and Ω > 1 because by the two-scale expansion assume that the solution of (1) can be well approximated where p > 0. by Elsevier Ltd. All rights reserved. by the two-scale expansion Copyright © 2018, 2018 IFAC 224Hosting 2405-8963 © IFAC (International Federation of Automatic Control) ∗ ∗ ∗
Copyright © 2018 IFAC 224 Peer review©under of International Federation of Automatic Copyright 2018 responsibility IFAC 224Control. 10.1016/j.ifacol.2018.07.227 Copyright © 2018 IFAC 224
2018 IFAC TDS Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
225
Now we divide (20) by eiβ(t1 ) and introduce φ (t1 ) = σt1 − β (t1 ) ,
(21)
to get 1 0 = − Aeiφ(t1 ) + i (D1 α (t1 ) + i (σ − D1 φ(t1 ))) 2 1 + α (t1 ) p − e−iτ p + 2iζ . 2
(22) Separating the real and imaginary parts of equation (22) results in two ordinary differential equations describing the evolution of the amplitude and phase A p (23) D1 α = − ζ + sin τ α + sin φ, 2 2 αD1 φ =
Fig. 1. The linear stability chart of (1) At the minima (‘notches’) of the stability lobes, Ω, p, τ assume the particularly simple forms (11) Ωcrit = 1 + 2ζ, (12) pcrit = 2ζ (ζ + 1) , 2 jπ − arctan √ 1 1+2ζ √ , j = 1, 2, . . . (13) τcrit = 1 + 2ζ 3. THE LINEAR DELAY-DIFFERENTIAL EQUATION WITH HARMONIC FORCING Here we study the harmonically forced linear equation (14) x ¨ + 2ζ x˙ + x = p (xτ − x) + A cos (ωt) . With assumption (2), substituting the differential operators (6) and (7) into (14) and equating like powers of ε one obtains (15) ε0 : D02 x0 (t0 , t1 ) + x0 (t0 , t1 ) = 0, D02 x1 (t0 , t1 ) + x1 (t0 , t1 ) = −2D0 D1 x0 (t0 , t1 ) ε : −2ζD0 x0 (t0 , t1 ) + p [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] +A cos (t0 + σt1 ) . (16) Solving (15) for x0 (t0 , t1 ), knowing that the slow temporal variable t1 is implicit in the constants of integration, we get x0 (t0 , t1 ) = a (t1 ) eit0 + a ¯ (t1 ) e−it0 , (17) We substitute this solution into (16) and eliminate the secular terms which would give rise to unbounded terms 1 0 = − Aeiσt1 + 2iD1 a (t1 ) + a (t1 ) p − e−iτ p + 2iζ . 2 (18) Writing a (t1 ) in a polar form (where α (t1 ) is the amplitude and β (t1 ) is the phase) 1 1 a ¯ (t1 ) = α (t1 ) e−iβ(t1 ) , (19) a (t1 ) = α (t1 ) eiβ(t1 ) , 2 2 then substituting into (18), we get 1 0 = − Aeiσt1 + i (D1 α (t1 ) + iα (t1 ) D1 β (t1 )) eiβ(t1 ) 2 1 + α (t1 ) eiβ(t1 ) p − e−iτ p + 2iζ . 2 (20) 1
225
A 1 (p cos τ + 2σ − p) α + cos φ. 2 2
(24)
To get the amplitude and phase of the steady-state primary resonance (i.e. the fixed point of (23) and (24)) we set the left-hand sides of (23) and (24) to zero and solve the resulting algebraic equations (by eliminating the trigonometric terms from them). For the amplitude α∗ and phase φ∗ we get A α∗ = 2 (p (1 − cos τ ) − 2(ω − 1)) + 4 ζ +
p 2
sin τ
α p . φ∗ = arcsin 2 ζ + sin τ 2 A
2 .
(25)
(26)
The stability of the fixed point (α∗ , φ∗ ) is determined by the Jacobian p A ∗ cos φ −ζ − 2 sin τ 2 J = (27) . A A sin φ∗ − ∗ 2 cos φ∗ − 2α 2α The trace and determinant of J are
T = −2ζ − p sin τ,
(28)
2 1 p ∆ = ζ + sin τ + (p(1 − cos τ ) − 2ω + 2)2 . 2 4
(29)
T 2 − 4∆ = −(p(1 − cos τ ) − 2ω + 2)2 ≤ 0,
(30)
The determinant ∆ is always non-negative. The discriminant
therefore the fixed point is a stable spiral when T < 0, a center when T = 0 (at pbif = −2 sinζ τ ), and an unstable spiral when T > 0. We chose p as the bifurcation parameter. The phase portraits of the “slow-flow” (Eqs. (23) and (24)) are shown in Figs. 2, 3, 4 for different p values.
2018 IFAC TDS 226 Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
3.1 Numerical results The direct numerical solution of (14) was determined in Mathematica using the Dormand-Prince method with an initial function α∗ cos(ωt + φ∗ ), t ∈ [−τ, 0], where α∗ and φ∗ were determined from (25) and (26). The following numerical values were chosen for the simulation ζ = 0.01, A = 0.01. (31) We have chosen p values at the first ‘notch’ (j = 1). Using (12) and (13) gives (32) pcrit = 0.0202, τcrit = 4.676. Time series x (t) for (33) p = {0.75pcrit , 0.90pcrit , 0.99pcrit , 1.05pcrit } are shown in Figure 5.
Fig. 2. Stable spiral p < pbif
(a) p = 0.75pcrit , stable
(b) p = 0.90pcrit , stable
(c) p = 0.99pcrit , stable
(d) p = 1.05pcrit , unstable
.
Fig. 3. Center p = pbif
. Fig. 5. Time series of the numerical solution (ω = 1.02) The grey areas indicate the time interval [300τ, 450τ ] where the steady state vibration amplitude of the numerical solution was determined as (34) αnum = max|x(t)|, t ∈ [300τ, 450τ ]. In the case of p = 0.90pcrit and p = 0.99pcrit the vibration becomes modulated, so the amplitude of the vibration also varies between a certain range during the simulation. For a better comparison of the method of multiple scales and numerical results we introduce the method of multiple scales amplification factor G∗ and numerical amplification factor Gnum as α∗ αnum , Gnum := . (35) G∗ := A A The comparison of the amplification factors G∗ and Gnum is shown in Figure 6, where the solid line indicates the method of multiple scales (MMS) results, the dots the numerical results.
Fig. 4. Unstable spiral p > pbif 226
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
227
tial equations describing the evolution of the amplitude and phase τ p A τ D1 α = − ζ + sin τ α − 3q cos sin3 α3 + sin φ, 2 2 2 2 (41) 1 A τ αD1 φ = (p cos τ + 2σ − p) α − 3q sin4 α3 + cos φ. 2 2 2 (42) To get the amplitude of the steady-state primary resonance i.e. the fixed points of (41) and (42) we set the lefthand sides of (41) and (42) to zero and solve the resulting algebraic equations (by eliminating the trigonometric terms from them) τ 2 1 p (1 − cos τ ) − 2σ + 6qα∗ 2 sin4 4 2 2 A2 p τ 3 τ ∗2 α∗ 2 − = 0. + ζ + sin τ + 3qα cos sin 2 2 2 4 (43) p 2 τ 3 τ ∗3 ∗ ∗ ζ + sin τ α + 3q cos sin α φ = arcsin . A 2 2 2 (44) Equation (43) can be solved for α∗ resulting 1, 2 or 3 pair (±) of real roots.
MMS
(a) p = 0.75pcrit MMS
(b) p = 0.90pcrit MMS
Stability of the fixed points (α∗ , φ∗ ) is determined by the Jacobian τ p A 3 τ ∗2 ∗ cos φ −ζ − 2 sin τ − 9qα cos 2 sin 2 2 J = . τ A A −6α∗ q sin4 − ∗ 2 cos φ∗ − ∗ sin φ∗ 2 2α 2α (45) The stability of the fixed points can be investigated by calculating the determinant ∆, the trace T and the discriminant T 2 − 4∆ of J .
(c) p = 0.99pcrit
Fig. 6. The amplification factor G as a function of ω 4. THE HARMONICALLY FORCED NONLINEAR DELAY-DIFFERENTIAL EQUATION Substituting the differential operators (6) and (7) into (1) and equating like powers of ε one obtains ε0 : D02 x0 (t0 , t1 ) + x0 (t0 , t1 ) = 0, (36) D02 x1 (t0 , t1 ) + x1 (t0 , t1 ) = −2D0 D1 x0 (t0 , t1 ) −2ζD0 x0 (t0 , t1 ) + p [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] ε1 : 2 +q [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] 3 +q [x0 (t0 − τ, t1 ) − x0 (t0 , t1 )] + A cos (t0 + σt1 ) . (37) Solving (36) for x0 (t0 , t1 ) yields x0 (t0 , t1 ) = a (t1 ) eit0 + a ¯ (t1 ) e−it0 . (38) Then we substitute this solution into (37) and eliminate the secular terms which would give rise to unbounded terms 1 0 = − Aeiσt1 + 2iD1 a (t1 ) 2 3 ¯ (t1 ) . +a (t1 ) p − e−iτ p + 2iζ − 3e−2iτ eiτ − 1 qa (t1 ) a (39) Writing a (t1 ) in a polar form (where α (t1 ) is the amplitude and β (t1 ) is the phase) 1 1 a (t1 ) = α (t1 ) eiβ(t1 ), a ¯ (t1 ) = α (t1 ) e−iβ(t1 ). (40) 2 2 then substituting into (39) using (21) and separating the real and imaginary parts results in two ordinary differen227
The nondimensional cutting force coefficient was taken to be p = 0.5pcrit , the coefficient of nonlinearity q = 0.003, forcing amplitude A = 0.01. The quantity ω was chosen as bifurcation parameter, increased from 0.98 to 1.03. 2.0
1.5
MMS Numerical ÈΑÈ
2018 IFAC TDS Budapest, Hungary, June 28-30, 2018
_ SBH _ < SN
1.0
<
0.5
0.0 0.98
0.99
1.00
1.01
1.02
1.03
Ω
Fig. 7. Bifurcation diagram (or resonance curve), where SBH corresponds to subcritical Hopf bifurcation and SN to saddle-node bifurcation.
2018 IFAC TDS 228 Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
Figure 7 shows the amplitude α∗ as a function of forcing frequency ω (the solid line indicates stable solution, the dashed line the unstable solution, both method of multiple scales solution). The stable numerical equlibrium solutions are shown with filled circles, unstable limit cycles are depicted by empty circles. The thin vertical lines in Fig. 7 indicate the ω values (1.00, 1.005, 1.013, 1.02) for which the phase portraits of Figs. (8-11) were generated.
A saddle and a node connected with a heteroclinic orbit can be observed at ω ≈ 1.013, see Figure 10.
Fig. 11. Stable spiral and unstable limit cycle at ω = 1.02 5. CONCLUSIONS
Fig. 8. Stable spiral at ω = 1, the arrowheads indicate two points of the unstable limit cycle
The method of multiple scales gave an excellent approximation of the solution of the nonlinear delay-differential equation (1). The resonance curves found here are similar to those for the Duffing-equation (Nayfeh and Mook [1979]), having a hardening characteristic. We found subcritical Hopf bifurcations and saddle-node bifurcations. Equation (1) without forcing admits a subcritical Hopf bifurcation (Kalm´ar-Nagy [2002]) and the interaction of forcing and Hopf bifurcation may lead to incredibly complex phenomena (Gambaudo [1985]). Plaut and Hsieh (Plaut and Hsieh [1987]) studied a one-DOF mechanical system with delay and excitation and found periodic, chaotic and unbounded responses. Forcing of a Duffing-type equation without delay can also result in complicated bifurcation structure (Sanchez and Nayfeh [1990]; Zavodney et al. [1990]). REFERENCES
Fig. 9. Unstable spiral at ω = 1.005 Subcritical Hopf bifurcation occurs at ω ≈ 1.003, see Figure 8 and 9.
Fig. 10. Stable spiral at ω = 1.013 and two unstable fixed points 228
Balachandran, B., Kalm´ar-Nagy, T., and Gilsinn, D.E. (2009). Delay differential equations. Springer. Casal, A. and Freedman, M. (1980). A Poincar´e-Lindstedt approach to bifurcation problems for differential-delay equations. IEEE Transactions on Automatic Control, 25(5), 967–973. Gambaudo, J.M. (1985). Perturbation of a Hopf bifurcation by an external time-periodic forcing. Journal of differential equations, 57(2), 172–199. Hanna, N. and Tobias, S. (1974). A theory of nonlinear regenerative chatter. ASME Journal of Engineering for Industry, 96(1974), 247–255. Kalm´ar-Nagy, T. (2002). Delay-differential models of cutting tool dynamics with nonlinear and mode-coupling effects. PhD thesis. Kalm´ar-Nagy, T., St´ep´an, G., and Moon, F.C. (2001). Subcritical Hopf bifurcation in the delay equation model for machine tool vibrations. Nonlinear Dynamics, 26(2), 121–142. Nayfeh, A., Chin, C., and Pratt, J. (1997). Applications of perturbation methods to tool chatter dynamics. Dynamics and chaos in manufacturing processes, 193–213.
2018 IFAC TDS Budapest, Hungary, June 28-30, 2018
János Lelkes et al. / IFAC PapersOnLine 51-14 (2018) 224–229
Nayfeh, A.H. (2008). Order reduction of retarded nonlinear systems–the method of multiple scales versus centermanifold reduction. Nonlinear Dynamics, 51(4), 483– 500. Nayfeh, A.H. and Mook, D.T. (1979). Forced oscillations of systems having a single degree of freedom. Nonlinear Oscillations, 161–257. Oztepe, G.S., Choudhury, S.R., and Bhatt, A. (2015). Multiple scales and energy analysis of coupled Rayleighvan der Pol oscillators with time-delayed displacement and velocity feedback: Hopf bifurcations and amplitude death. Far East Journal of Dynamical Systems, 26(1), 31. Pakdemirli, M. and Karahan, M.M.F. (2010). A new perturbation solution for systems with strong quadratic and cubic nonlinearities. Mathematical Methods in the Applied Sciences, 33(6), 704–712. Plaut, R. and Hsieh, J.C. (1987). Chaos in a mechanism with time delays under parametric and external excitation. Journal of Sound and Vibration, 114(1), 73–90. Sanchez, N.E. and Nayfeh, A.H. (1990). Prediction of bifurcations in a parametrically excited Duffing oscillator. International Journal of Non-Linear Mechanics, 25(23), 163–176. St´ep´ an, G. (1989). Retarded dynamical systems: stability and characteristic functions. Longman Scientific & Technical. Zavodney, L.D., Nayfeh, A., and Sanchez, N. (1990). Bifurcations and chaos in parametrically excited singledegree-of-freedom systems. Nonlinear Dynamics, 1(1), 1–21.
229
229