Stability and bifurcation analysis of micro-electromechanical nonlinear coupling system with delay

J. Math. Anal. Appl. 461 (2018) 577–590

Contents lists available at ScienceDirect

Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa

Stability and bifurcation analysis of micro-electromechanical nonlinear coupling system with delay Yuting Ding a,∗ , Liyuan Zheng b , Jinli Xu a a b

Department of Mathematics, Northeast Forestry University, Harbin, 150040, PR China Library, Northeast Forestry University, Harbin, 150040, PR China

a r t i c l e

i n f o

Article history: Received 30 January 2017 Available online xxxx Submitted by J. Shi Keywords: Micro-electromechanical coupling system Hopf bifurcation Hopf-pitchfork bifurcation Multiple time scales Normal form

a b s t r a c t In this paper, we study dynamics in delayed micro-electromechanical nonlinear coupling system, with particular attention focused on Hopf and Hopf-pitchfork bifurcations. Based on the distribution of eigenvalues, we prove that a sequence of Hopf and Hopf-pitchfork bifurcations occur at the trivial equilibrium as the delay increases and obtain the critical values of two types of bifurcations. Next, by applying the multiple time scales method, the normal forms near the Hopf and Hopf-pitchfork bifurcations critical points are derived. Finally, bifurcation analysis and numerical simulations are presented to demonstrate the application of the theoretical results. We show the regions near above bifurcation critical points in which the micro-electromechanical nonlinear coupling system exists stable fixed point or stable periodic solution. Detailed numerical analysis using MATLAB extends the local bifurcation analysis to a global picture, and stable windows are observed as we change control parameters. Namely, the stable fixed point and stable periodic solution can exist in large regions of unfolding parameters as the unfolding parameters increase away from the critical value. © 2018 Elsevier Inc. All rights reserved.

1. Introduction Micro-electromechanical systems (MEMS, also written as micro-electro-mechanical, micro-electromechanical or microelectronic and microelectromechanical systems and the related micromechatronics) are the technology of microscopic devices, particularly those with moving parts. They merge at the nano-scale into nanoelectromechanical systems (NEMS) and nanotechnology. MEMS are separate and distinct from the hypothetical vision of molecular nanotechnology or molecular electronics. They usually consist of a central unit that processes data (the microprocessor) and several components that interact with the surroundings such as microsensors. At these size scales, the standard constructs of classical physics are not always sufficient. Because of the large surface area to volume ratio of MEMS, surface effects such as electrostatics and wetting dominate over volume effects such as inertia or thermal mass. MEMS become practical once * Corresponding author. E-mail addresses: [email protected] (Y. Ding), [email protected] (L. Zheng), [email protected] (J. Xu). https://doi.org/10.1016/j.jmaa.2018.01.032 0022-247X/© 2018 Elsevier Inc. All rights reserved.

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

578

Fig. 1. Dynamic model of the electromechanical coupling system driven by an AC motor.

they could be fabricated using modified semiconductor device fabrication technologies, normally used to make electronics. These include molding and plating, wet etching (KOH, TMAH) and dry etching (RIE and DRIE), electro discharge machining (EDM), and other technologies capable of manufacturing small devices. An early example of MEMS device is the resonistor, that is, an electromechanical monolithic resonator [7,16]. The MEMS represent a very important class of systems having applications in all fields. Because of the nonlinearity and the complexity, modeling and dynamics of MEMS have strongly attracted scholar’s attention. MEMS often involve the nonlinear coupling of electrostatic and mechanical physical fields in engineering, so the dynamic characters are complicated. The transmission system driven by an alternating current (AC) motor as a core part of rotating machinery plays an irreplaceable role in electricity, energy, transportation, metallurgy, and national defense fields. The torsional vibration of the rotating machinery will affect the normal work and cause the equipment to be damaged [9,10]. The rolling mill main drive system driven by an AC motor can be abstracted into a two-mass torsional vibration model, as shown in Fig. 1. By considering the energy in the air-gap field of the AC motor, the dynamical equation of the electromechanical coupling transmission system, associated with a two masses relative rotation system driven by an AC motor, is deduced as follows [9,10]: 

J1 ϕ¨1 + K(ϕ1 − ϕ2 ) + C(ϕ˙ 1 − ϕ˙ 2 ) + Ce ϕ˙ 1 − Ke ϕ1 = Me , J2 ϕ¨2 − K(ϕ1 − ϕ2 ) − C(ϕ˙ 1 − ϕ˙ 2 ) = Md ,

(1)

where Ji > 0 (i = 1, 2) is the moment of inertia; ϕi (i = 1, 2) and ϕ˙ i (i = 1, 2) are angle of rotation and rotational speed, respectively; K > 0 is the torsional stiffness of drive shaft; C > 0 is the shafting damping coefficient; Ce > 0 is the electromagnetic damping coefficient (generated by the rotor damper bar); Me is the electromagnetic torque, and Md is the load torque. The detailed deducing can be referred to Refs. [9,10]. Equation (1) is the nonlinear dynamic equation of two masses relative rotation system considering both electrical parameters and mechanical parameters, which is the basis of dynamic behavior analysis of electromechanical coupling relative rotation system. Time delay, as a kind of basic nature of the physical phenomenon, exists widely in industrial systems. It has a great influence on the analysis and the control of a system, such as leading to the instability of the system and the control law failure. However, reasonably introducing and processing the time delay can improve the performance of a system. It is well known that delay differential equations (DDE) may exhibit higher codimension singularities more frequently than that in ordinary differential equations (ODE) [1,2,8, 13–15,17,20]. In order to control the dynamic behaviors of the electromechanical coupling system, the nonlinear time delay state feedback term just like k˜1 ϕ2 (t − τ ) + k˜2 ϕ32 (t − τ ) is introduced, that is, Md = k˜1 ϕ2 (t − τ ) + k˜2 ϕ32 (t − τ ). In the case of no eccentricity, Me = ∂Wm /∂ϕ1 = μ0 + μ1 ϕ1 + μ2 ϕ21 + μ3 ϕ31 , where Wm is air-gap magnetic field energy between the stator and rotor, μ0 = μ ˜0 sin(ψ1 + ψ2 ), μ1 = μ ˜1 cos(ψ1 + ψ2 ), μ2 = μ ˜2 sin(ψ1 + ψ2 ), μ3 = μ ˜3 cos(ψ1 + ψ2 ),

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

579

with the phase angle of the stator synthetic magneto motive force ψ1 + ψ2 . Actually, ψ1 + ψ2 can equal π through proper rotation transformation, then μ0 = 0, μ1 = μ ˜1 , μ2 = 0 and μ3 = μ ˜3 . Let x = ϕ1 , y = ϕ˙ 1 , z = ϕ2 , w = ϕ˙ 2 , a1 = (μ1 + Ke − K)/J1 , a2 = −(C + Ce )/J1 , a3 = K/J1 , a4 = C/J1 , l1 = μ3 /J1 , b1 = K/J2 , b2 = C/J2 , k1 = k˜1 /J2 , k2 = k˜2 /J2 , a model for micro-electromechanical nonlinear coupling system with delay can be described mathematically as follows: ⎧ x˙ = y, ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ y˙ = a1 x + a2 y + a3 z + a4 w + l1 x3 , ⎪ z˙ = w, ⎪ ⎪ ⎪ ⎪ ⎩ w˙ = b1 x + b2 y − b1 z − b2 w + k1 z(t − τ ) + k2 z 3 (t − τ ).

(2)

In this paper, by applying the local stability theory, we investigate the existence of several types of bifurcations for this micro-electromechanical nonlinear coupling system with delay (2). For example, it can bring fold bifurcation, Hopf bifurcation and Hopf-pitchfork bifurcation, thus makes system more complicated. Multiple time scales (MTS) [3,4,11,12] and center manifold reduction (CMR) [4–6,18] are two useful methods for computing normal forms of differential equations. Yu et al. [19] proved the equivalence of the MTS method and the CMR method for computing the normal forms of ordinary differential equations, and shown the conditions under which the normal forms, derived by using the MTS and CMR methods, are identical up to third order for delay differential equations. Actually, the normal forms of system (2) associated with Hopf and Hopf-zero bifurcations, derived by using the MTS and CMR methods, are identical up to third order due to Ref. [19]. Thus, we derive the normal forms near Hopf and Hopf-pitchfork critical points by applying the multiple time scales (MTS) method, and show the bifurcation analysis and numerical simulations associated with the above two bifurcations. The rest of the paper is organized as follows. In Section 2, we consider the stability of the trivial equilibrium and the existence of fold bifurcation, Hopf bifurcation and Hopf-pitchfork bifurcation in microelectromechanical nonlinear coupling system with delay. Then we derive the normal forms associated with Hopf and Hopf-pitchfork bifurcations in Section 3. Some examples, associated with above bifurcation analysis and numerical simulations, are presented in Section 4. Finally, conclusion is drawn in Section 5. 2. Stability of the trivial equilibrium and the existence of several types of bifurcations In this section, system (2) is considered, and we will consider the trivial equilibrium of this system. The characteristic equation of system (2), evaluated at origin, is given by λ4 + c3 λ3 + c2 λ2 + c1 λ + c0 − k1 e−λτ (λ2 − a2 λ − a1 ) = 0,

(3)

where c3 = b2 − a2 , c2 = b1 − a2 b2 − b2 a4 − a1 , c1 = −a2 b1 − b2 a3 − a1 b2 − b1 a4 , c0 = −a1 b1 − b1 a3 . Case 1: Fold bifurcation. Note that when c0 + k1 a1 = 0, that is, k1 a1 − a1 b1 − b1 a3 = 0, the characteristic equation (3) has one zero root λ = 0, and system (2) undergoes a fold bifurcation. When τ = 0 and k1 a1 − a1 b1 − b1 a3 = 0, Eq. (3) becomes λ[λ3 + c3 λ2 + (c2 − k1 )λ + c1 + k1 a2 ] = 0. We give following assumption: (H1) c3 > 0, c3 (c2 − k1 ) − c1 − k1 a2 > 0, c1 + k1 a2 > 0. Then, we have the following conclusions for system (2).

(4)

580

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

Theorem 2.1. When k1 a1 − a1 b1 − b1 a3 = 0, system (2) undergoes a fold bifurcation. Particularly, when (H1) holds, all the roots of Eq. (4) have negative real part expect one zero root. Case 2: Hopf bifurcation. When τ = 0, Eq. (3) becomes λ4 + c3 λ3 + (c2 − k1 )λ2 + (c1 + k1 a2 )λ + c0 + a1 k1 = 0.

(5)

We give following assumption: (H2) c3 > 0, c3 (c2 −k1 )−c1 −k1 a2 > 0, c0 +a1 k1 > 0, c3 (c2 −k1 )(c1 +k1 a2 )−c23 (c0 +a1 k1 )−(c1 +k1 a2 )2 > 0, then all the roots of Eq. (5) have negative real parts due to Routh–Hurwitz criterion. To find possible periodic solutions, which may bifurcate from a Hopf critical point, let λ = iβ (i2 = −1, β > 0) be a root of Eq. (3). Substituting λ = iβ into (3) and separating the real and imaginary parts yields β 4 − c2 β 2 + c0 = −k1 a2 β sin(βτ ) − k1 (a1 + β 2 ) cos(βτ ), c1 β − c3 β 3 = k1 (β 2 + a1 ) sin(βτ ) − k1 a2 β cos(βτ ).

(6)

Let Z = β 2 . Then it follows from (6) that F (Z) := Z 4 + m3 Z 3 + m2 Z 2 + m1 Z + m0 = 0,

(7)

where m3 = c23 − 2c2 , m2 = c22 + 2c0 − 2c1 c3 − k12 , m1 = c21 − k12 a22 − 2k12 a1 − 2c0 c2 , m0 = c20 − k12 a21 . Assume Eq. (7) has positive real root. Without loss of generality, we assume that Eq. (7) has k (k ≤ 4) positive √ roots: Zn (n = 1, · · · , k), thus, βn = Zn . Due to Eq. (6), we obtain Qn := sin(βn τ ) =

βn k1 (a1 + βn2 )(c1 − c3 βn2 ) − k1 a2 βn (βn4 − c2 βn2 + c0 ) , k12 a22 βn2 + k12 (a1 + βn2 )2

k1 (βn4 − c2 βn2 + c0 )(βn2 + a1 ) + k1 a2 βn2 (c1 − c3 βn2 ) Pn := cos(βn τ ) = − . k12 a22 βn2 + k12 (a1 + βn2 )2

(8)

The time delay τ can be determined from (8) as  τn(j)

=

1 βn [arccos(Pn ) + 2jπ], 1 βn [2π − arccos(Pn ) + 2jπ],

Qn ≥ 0, Qn < 0,

(9)

where n = 1, · · · , k; j = 0, 1, 2, · · · . (j) (j) Let λ(τ ) = α(τ ) + iω(τ ) be the root of Eq. (3) satisfying α(τn ) = 0, β(τn ) = βn , n = 1, · · · , k; j = 0, 1, 2, · · · . Then, we have the transversality conditions: Sign[Re(

dλ

)−1 ] = Sign[ (j)

dτn

k12 a22 βn2

F  (Zn ) ] = Sign[F  (Zn )], + k12 (a1 + βn2 )2

(10)

where n = 1, · · · , k; j = 0, 1, 2, · · · . Combining the above results, we have the following theorem. Theorem 2.2. If Eq. (7) has k(k ≤ 4) positive real root Zn , n = 1, · · · , k, system (2) undergoes a Hopf (j) (j) bifurcation at τ = τn (n = 1, · · · , k; j = 0, 1, 2 · · · ), where τn is given by (9). In particular, if (H2) holds, (0) when τ ∈ [0, τ0 ) with τ0 = min{τn }, all the roots of Eq. (3) have negative real part.

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

581

Case 3: Hopf-zero bifurcation. When k1 a1 − a1 b1 − b1 a3 = 0, characteristic equation (3) has one zero root. Then Eq. (7) becomes F (Z) := Z[Z 3 + m3 Z 2 + m2 Z + m1 ] = 0. Let F1 (Z) := Z 3 + m3 Z 2 + m2 Z + m1 and 1 (Z) equation dFdZ = 0 has two real roots: Z1∗

dF1 (Z) dZ

−2m3 − = 6

√

Δ

(11)

:= 3Z 2 + 2m3 Z + m2 . When Δ = 4m23 − 12m2 > 0,

,

Z2∗

−2m3 + = 6

√

Δ

.

Therefore, we give the following assumptions: (H3a) Δ > 0, m1 < 0, Z1∗ < 0. (H3b) Δ > 0, m1 > 0, Z2∗ > 0, F1 (Z2∗ ) < 0. (H3c) Δ > 0, m1 < 0, Z1∗ > 0, F1 (Z1∗ ) < 0, F1 (Z2∗ ) > 0. Under (H3a), Eq. (7) has one positive root Z1 , and F  (Z1 ) > 0. Under (H3b), Eq. (7) has two positive roots, Z1 and Z2 . Suppose Z1 < Z2 , then F  (Z1 ) < 0, F  (Z2 ) > 0. Under (H3c), Eq. (7) has three positive roots, Z1 , Z2 and Z3 . Suppose Z1 < Z2 < Z3 , then F  (Z1 ) > 0, F  (Z2 ) < 0 and F  (Z3 ) > 0. Combining the above results, we have the following theorem. Theorem 2.3. When k1 a1 − a1 b1 − b1 a3 = 0, if (H3a), (H3b) or (H3c) holds, Eq. (7) has 1, 2, 3 positive real (j) roots, respectively, and system (2) undergoes a Hopf-zero bifurcation at τ = τn (n = 1, 2, 3; j = 0, 1, 2 · · · ), (j) (0) where τn is given by (9). In particular, if (H1) also holds, when τ ∈ [0, τ0 ) with τ0 = min{τn }, Eq. (3) has a zero root, and all the other roots have negative real part. 3. Normal forms of Hopf-zero and Hopf bifurcations In this section, we derive the normal forms of Hopf-zero and Hopf bifurcations by using the multiple time scales method. We only show the computation process associated with the normal form of Hopf-zero bifurcation, and the details associated with Hopf bifurcation can be deduced similarly. 3.1. Normal form of Hopf-zero bifurcation (j)

When k1 a1 − a1 b1 − b1 a3 = 0 and τ = τn (n = 1, 2, 3; j = 0, 1, 2, · · · ), the characteristic equation (3) (j) has one zero root λ = 0 and one pair of pure imaginary roots λ = ±iβn , where τn is given by (9). We treat the feedback parameter k1 and the delay τ as two bifurcation parameters. Remark 1. In the process of motor running, complex dynamic behaviors may occur due to the feedback term for electromechanical coupling system. Thus, feedback delay τ and linear feedback strength k1 are selected as two major parameters to analyze the system dynamic characteristics. Suppose system (2) undergoes a Hopf-zero bifurcation from the trivial equilibrium at the critical point: k1 = k1c , τ = τc , and then the characteristic equation (3) has one zero root λ = 0 and one pair of pure imaginary roots λ = ±iβ. Further, we write Eq. (2) as follows: U˙ (t) = AU (t) + BU (t − τ ) + F (U (t), U (t − τ )),  T  T where U (t) = x(t), y(t), z(t), w(t) , U (t − τ ) = x(t − τ ), y(t − τ ), z(t − τ ), w(t − τ ) ,

(12)

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

582

⎛

0 ⎜ a1 A=⎝ 0 b1

1 a2 0 b2

0 a3 0 −b1

⎞ ⎛ 0 0 0 a4 ⎟ ⎜0 0 , B=⎝ 0 0 1 ⎠ 0 0 −b2

0 0 0 k1

⎞ ⎛ ⎞ 0 0 3 0⎟ l1 x ⎜ ⎟ , F (U (t), U (t − τ )) = ⎝ ⎠. 0⎠ 0 0 k2 z 3 (t − τ )

Define the linearized equation of (12) as U˙ (t) = AU (t) + Bc U (t − τc ) := Lc (U (t), U (t − τc )), whose characteristic equation has a pair of purely imaginary roots ±iβ and a zero root, and no other roots with zero real part, where ⎛

0 ⎜0 Bc = ⎝ 0 0

0 0 0 0 0 0 0 k1c

⎞ 0 0⎟ . 0⎠ 0

Let p1 = (p11 , p12 , p13 , p14 )T and p2 = (p21 , p22 , p23 , p24 )T be the two eigenvectors of the linear operator Lc corresponding to the eigenvalues iβ and 0, respectively. Further, let p∗1 = (p∗11 , p∗12 , p∗13 , p∗14 )T and p∗2 = (p∗21 , p∗22 , p∗23 , p∗24 )T be the two normalized eigenvectors of the adjoint operator L∗c of the linear operator Lc corresponding to the eigenvalues −iβ and 0, respectively, satisfying the inner product T

p∗i , pi  = p¯∗i pi = 1, i = 1, 2. By a simple calculation, we have a1 + β 2 + a2 βi , iβp13 )T , a3 + a4 βi a1 p2 = (p21 , p22 , p23 , p24 )T = (1, 0, − , 0)T , a3

p1 = (p11 , p12 , p13 , p14 )T = (1, iβ, −

(13)

p∗1 = d1 (p∗11 , p∗12 , p∗13 , p∗14 )T , p∗2 = d2 (p∗21 , p∗22 , p∗23 , p∗24 )T , where p∗11 = 1, p∗12 = p∗14 =

b1 − iβb2 , a1 b2 − b1 (iβ + a2 )

iβ(iβ + a2 ) − a1 , p∗ = −a4 p∗12 − (iβ − b2 )p∗14 , a1 b2 − b1 (iβ + a2 ) 13

d1 = (1 + p∗12 p¯12 + p∗13 p¯13 + p∗14 p¯14 )−1 , a2 b1 b1 a4 b1 − b2 , p∗22 = − , p∗23 = + b2 , p∗24 = 1, a1 a1 a1 a 1 a3 d2 = . (a2 b1 − b2 a1 )a3 − (a4 b1 + b2 a1 )a1 p∗21 =

Next, we use the multiple time scales (MTS) method to derive the normal form of system (2) associated with Hopf-zero bifurcation. The solution of (12) is assumed to take the form: U (t) = U (T0 , T1 , · · · ) = 1/2 U1 (T0 , T1 , · · · ) + 3/2 U2 (T0 , T1 , · · · ) + · · · ,

(14)

 T where U (T0 , T1 , · · · ) = x(T0 , T1 , · · · ), y(T0 , T1 , · · · ), z(T0 , T1 , · · · ), w(T0 , T1 , · · · ) and Uk (T0 , T1 , · · · ) = T  xk (T0 , T1 , · · · ), yk (T0 , T1 , · · · ), zk (T0 , T1 , · · · ), wk (T0 , T1 , · · · ) , k = 1, 2, · · · . The derivative with respect to t is now transformed into

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

583

∂ ∂ d = + + · · · = D0 + D1 + · · · , dt ∂T0 ∂T1 ∂ where the differential operator Di = ∂T , i = 0, 1, 2, · · · . i T Define Uj = (xj , yj , zj , wj ) = Uj (T0 , T1 , · · · ) and Uj,τc = (xj,τc , yj,τc , zj,τc , wj,τc )T = Uj (T0 − τc , T1 , · · · ) with j = 1, 2, · · · . From (14), we obtain

U˙ (t) = 1/2 D0 U1 + 3/2 D1 U1 + 3/2 D0 U2 + · · · .

(15)

We take perturbations as k1 = k1c + k and τ = τc + τ in system (12), where k and τ are called detuning parameters [12]. To deal with the delayed terms, we expand Uj (t − τ ) at Uj,τc for j = 1, 2, · · · , that is, U (t − τc − τ , (t − τc − τ ), · · · ) = 1/2 U1,τc − 3/2 τ D0 U1,τc − 3/2 τc D1 U1,τc + 3/2 U2,τc + · · · ,

(16)

where Uj,τc = Uj (T0 − τc , T1 , · · · ) with j = 1, 2, · · · . Then, substituting the solutions with the multiple scales (14)–(16) into (12) and balancing the coefficients of j/2 (j = 1, 3, · · · ) yields a set of ordered linear differential equations. First, for the 1/2 -order terms, we have D0 U1 − AU1 − Bc U1,τc = 0.

(17)

Since ±iβ and 0 are the eigenvalues of the linear part of (12), the solution of (17) can be expressed in the form of ¯ 1 (T1 , T2 , · · · )e−iβT0 p¯1 + G2 (T1 , T2 , · · · )p2 , U1 (T1 , T2 , · · · ) = G1 (T1 , T2 , · · · )eiβT0 p1 + G

(18)

where p1 and p2 are given by (13). Next, for the 3/2 -order terms, we obtain D0 U2 − AU2 − Bc U2,τc = −D1 U1 + f,

(19)

 T 3 where f = 0, l1 x31 , 0, k z1,τc − k1c (τ D0 z1τc + τc D1 z1τc ) + k2 z1τ . c Nonhomogeneous equation (19) has a solution if and only if the so-called solvability condition is satisfied [11]. That is, the right-hand side of nonhomogeneous equation (19) is orthogonal to every solution of the adjoint homogeneous problem. Substituting solution (18) into the right expression of (19), we obtain the coefficient of term eiβT0 , denoted as δ1 , and the constant term, denoted as δ2 . As a matter of fact, finding the solvability conditions is equivalent to finding the conditions resulted from eliminating the secular terms. ∂G2 1 Let p∗1 , δ1  = 0 and p∗2 , δ2  = 0, where p∗j (j = 1, 2) is given by (13). Then ∂G ∂T1 and ∂T1 are solved to yield ∂G1 ¯ 1 + Q2 G1 G22 , = α1 G1 + Q1 G21 G ∂T1 ∂G2 ¯ 1 G2 , = α2 G2 + Q3 G32 + Q4 G1 G ∂T1 where α1 =

p23 k d2 p∗24 d¯1 p¯∗14 (p13 e−iβτc k − k1c p13 iβe−iβτc τ ) , α = , 2 1 + k1c τc p23 d2 p∗24 1 + k1c τc d¯1 p¯∗14 e−iβτc p13

(20)

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

584

Q1 =

3d¯1 (¯ 3d¯1 (¯ p∗12 l1 + k2 p¯∗14 p213 p¯13 e−iβτc ) p∗12 l1 + k2 p¯∗14 p13 p223 e−iβτc ) , Q2 = , ∗ −iβτ cp 1 + k1c τc d¯1 p¯14 e 1 + k1c τc d¯1 p¯∗14 e−iβτc p13 13

Q3 =

l1 p∗22 d2 + k2 p323 d2 p∗24 6d2 (l1 p∗22 + k2 p13 p¯13 p23 p∗24 ) , Q4 = . ∗ 1 + k1c τc p23 d2 p24 1 + k1c τc p23 d2 p∗24

Now, let G1 = reiθ and G2 = v. Substituting these expressions into (20), we obtain the following normal form in cylindrical coordinates: ⎧ dr ⎪ ⎪ = Re(α1 )r + Re(Q1 )r3 + Re(Q2 )rv 2 , ⎪ ⎪ dt ⎪ ⎪ ⎨ dv = α2 v + Q3 v 3 + Q4 vr2 , ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ dθ = Im(α1 ) + Im(Q1 )r2 + Im(Q2 )v 2 . dt

(21)

3.2. Normal form of Hopf bifurcation Assume when τ = τc , the characteristic equation (3) has a pair of pure imaginary roots λ = ±iβ. We treat the feedback delay τ as bifurcation parameter, and take perturbation as τ = τc + τ , where τ is called detuning parameters. The normal form associated with Hopf bifurcation derived by multiple time scales is ∂G ¯ = q1 τ G + q2 G2 G, ∂T1

(22)

where q1 = −

3d¯1 (¯ p∗12 l1 + k2 p¯∗14 p213 p¯13 e−iβτc ) d¯1 p¯∗14 k1 p13 iβe−iβτc , q2 = , ∗ −iβτ ¯ cp 1 + k1 τc d1 p¯14 e 1 + k1 τc d¯1 p¯∗14 e−iβτc p13 13

with p13 , p∗12 and p∗14 are given by (13). Substituting G = reiθ into (22), we obtain the following normal form in polar coordinates: ⎧ ⎪ ⎪ dr = Re(q1 )τ r + Re(q2 )r3 , ⎨ dt ⎪ dθ ⎪ ⎩ = Im(q1 ) + Im(q2 )r2 . dt

(23)

4. Bifurcation analysis and numerical simulations In this section, we first give a bifurcation analysis and some numerical simulation results associated with Hopf bifurcation, and then present the bifurcation analysis and numerical simulations for Hopf-pitchfork bifurcation. 4.1. Hopf bifurcation Considering the first equation of the normal form (23), that is, dr = Re(q1 )τ r + Re(q2 )r3 , dt then we have the following theorem.

(24)

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

585

Fig. 2. Simulated solution of system (2) associated with the time history for τ = 3, showing a stable equilibrium.

Theorem 4.1. The trivial equilibrium r0 = 0 of (24) is stable for Re(q1 )τ < 0, and unstable for Re(q1 )τ > 0.

Re(q1 )τ 1 )τ − Re(q Re(q2 ) when Re(q2 ) < 0, which corresponds to the  1 )τ periodic solution of system (2). The non-trivial equilibrium r1 = − Re(q Re(q2 ) is stable for Re(q1 )τ > 0, and

There exists another non-trivial equilibrium r1 =

unstable for Re(q1 )τ < 0. That is, system (2) exist bifurcating periodic solutions when are stable for Re(q1 )τ > 0, and unstable for Re(q1 )τ < 0.

Re(q1 )τ Re(q2 )

< 0, which

We choose a1 = −2.5, a2 = −2, a3 = 1, a4 = 1, b1 = 2, b2 = 2, k1 = 1, k2 = 1, l1 = 1, which satisfy the assumption (H2), namely, all the roots of Eq. (5) have negative real parts. Equation (7) has two positive root: Z1 =0.2504 ˙ and Z2 =0.9219, ˙ thus, from (8)–(10), we obtain  (0) (0) Z1 =0.5004, ˙ τ1 =11.4960, ˙ Re(λ (τ1 )) < 0,  (0) (0) (k) ˙ τ1 =4.7452, ˙ τ2 =11.2891, ˙ Re(λ (τ2 )) > 0 (k = 1, 2). β2 = Z2 =0.9602,

β1 =

(0)

We take τc = τ2 =4.7452. ˙ Thus, the characteristic equation (3) has a pair of purely imaginary eigenvalues ±iβ = ±0.9602i, and all the other eigenvalues have negative real part. Then system (2) undergoes a Hopf bifurcation from the trivial equilibrium. By a simple calculation from (22), we obtain Re(q1 ) = 0.0276 > 0, Re(q2 ) = 1.1715 > 0, (0)

thus, from Theorem 4.1, the trivial equilibrium is stable when τ ∈ [0, τ2 ) = [0, 4.7452), and system (2) (0) undergoes Hopf bifurcation when τ = τ2 , the bifurcating periodic solutions are stable for τ > 0. For example, we choose two groups of perturbation parameter value: τ = 3 and τ = 4.75 (that is, τ = 0.0048 > 0) with initial function being ϕ(θ) = [0.01, 0.05, 0.01, 0.05], θ ∈ [−τ, 0], resulting in a stable fixed point shown in Fig. 2 and stable periodic solution shown in Fig. 3, respectively. It is clear that the numerical simulations agree with the analytical predictions. 4.2. Hopf-pitchfork bifurcation Considering the first two equations of the normal form (21), associated with Hopf-zero bifurcation, that is,

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

586

Fig. 3. Simulated solution of system (2) associated with the time history for τ = 4.75, showing a stable periodic solution.

⎧ dr ⎪ ⎪ = Re(α1 )r + Re(Q1 )r3 + Re(Q2 )rv 2 , ⎨ dt ⎪ dv ⎪ ⎩ = α2 v + Q3 v 3 + Q4 vr2 . dt

(25)

dv Equilibrium solutions of (25) are obtained by simply setting dr dt = dt = 0. Note that E0 = (r, v) = (0, 0) corresponds to the original trivial equilibrium, and the other ones are

 Re(α1 ) Re(α1 ) , 0) for < 0, Re(Q1 ) Re(Q1 )  α2 α2 = (0, ± − ) for < 0, Q3 Q3   α2 Re(Q2 ) − Re(α1 )Q3 Re(α1 )Q4 − α2 Re(Q1 ) =( ,± ) Re(Q1 )Q3 − Re(Q2 )Q4 Re(Q1 )Q3 − Re(Q2 )Q4

E1 = ( − E2± E3±

for

Re(α1 )Q4 − α2 Re(Q1 ) α2 Re(Q2 ) − Re(α1 )Q3 > 0, > 0. Re(Q1 )Q3 − Re(Q2 )Q4 Re(Q1 )Q3 − Re(Q2 )Q4

The semi-trivial equilibria E1 and E2± are bifurcating from the origin on the critical lines L1 : Re(α1 ) = 0 and L2 : α2 = 0, respectively. The pair of nontrivial equilibria E3± collide with the semi-trivial equilibria E1 Re(α1 ) and E2± respectively on the critical lines L3 : Re(α1 )Q4 −α2 Re(Q1 ) = 0 for Re(Q < 0 and L4 : α2 Re(Q2 ) − 1) α2 Re(α1 )Q3 = 0 for Q3 < 0. Near Hopf-zero bifurcation critical points, the solutions on the center manifold determine the local asymptotic behavior of solutions of the original system (2). So, for (25), equilibria on the v-axis remain equilibria, while equilibria away from the r-axis become periodic orbits. Obviously, the system (2) undergoes Hopf(j) pitchfork bifurcation when k1 a1 − a1 b1 − b1 a3 = 0 and τ = τn . In order to give a more clear bifurcation picture, we choose a1 = −5, a2 = −3, a3 = 3, a4 = 2, b1 = 2.5, b2 = 1.5, k1 = 1, l1 = 1, k2 = 1, which satisfy the conditions k1 a1 − a1 b1 − b1 a3 = 0, (H1) and (H3a). Under these parameter values, the characteristic equation (3) with τ = 0 has a zero root, a negative root, and a pair of complex conjugate roots with negative real part. Equation (11) has only one positive root: Z1 =1.3020, ˙ thus, from (8)–(10), we obtain

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

587

Fig. 4. Critical bifurcation lines in the (τε , kε ) parameter plane near (τc , kc ) and the corresponding phase portraits in the (r, v) plane.

β=



Z1 =1.1410, ˙ τ1 =3.3365, ˙ Re(λ (τ1 )) > 0. (0)

(0)

(0)

We take τc = τ1 =3.3365. ˙ Thus, the characteristic equation (3) has a zero root and a pair of purely imaginary eigenvalues ±iβ = ±1.1410i, and all the other eigenvalues have negative real part. Then system (2) undergoes a Hopf-pitchfork bifurcation from the trivial equilibrium. By a simple calculation, we obtain Re(Q1 )=0.6736, ˙ Re(Q2 )=1.1202, ˙ Q3 =0.8022, ˙ Q4 =3.2643, ˙ Re(Q1 )Q3 − Re(Q2 )Q4 < 0, ˙ L1 : Re(α1 )=0.1503k  + 0.1088τ = 0, ˙ L2 : α2 =2607k  = 0, ˙ L3 : Re(α1 )Q4 − α2 Re(Q1 )=0.3150k  + 0.3550τ = 0, for Re(α1 ) < 0, ˙ L4 : α2 Re(Q2 ) − Re(α1 )Q3 =0.1714k  − 0.0872τ = 0, for α2 < 0. For the above chosen parameter values, the critical bifurcation lines become: L1 : k = −0.7237τ , L2 : k = 0, L3 : k = −1.1271τ , for k < −0.7237τ , L4 : b = 0.5090τ , for k < 0, as shown in the bifurcation diagram (see Fig. 4). Fig. 4 shows the critical bifurcation lines in the (τ , k ) parameter plane near the critical point (τc , kc ) and the corresponding phase portraits in the (r, z) plane, whose origin is the Hopf-pitchfork critical point. The bifurcation behaviors of the original system (2) in the neighborhood of (0, 0, 0, 0) can be observed from Fig. 4. Note that the bifurcation boundaries divide the (τ , k ) parameter plane into six regions. Also, it is seen from the phase portraits that the orbits are symmetric with respect to the r coordinate, therefore, we only draw the orbits in the first quadrant.

588

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

Fig. 5. Simulated solution of system (2) associated with the time history for τ = −0.1 and k = −0.1, showing a stable equilibrium.

For system (2), we use Fig. 4 to describe the bifurcations in the clockwise direction, starting from B1 and ending at B1 . First, in region B1 , there is only one trivial equilibrium which is a source. When the parameters are varied across line L2 (τ -axis) from region B1 to B2 , the trivial equilibrium becomes a saddle, and a pair of unstable fixed points, which are source (for convenience, also called E2±) appear from the trivial solution due to a pitchfork bifurcation. When the parameters are changed from region B2 to B3 , an unstable periodic solution, which is saddle (for convenience, also called E1 ) occurs from the trivial solution due to a Hopf bifurcation, and the trivial solution becomes a sink from a saddle. In region B4 , a pair of periodic solutions (for convenience, also called E3± ), which are saddles, occur from E1 due to a pitchfork bifurcation, and E1 becomes a source from a saddle. When the parameters are further changed from region B4 to B5 crossing line L4 , the pair of periodic solutions E3± collide with the semi-trivial fixed point E2± , and then disappear, and E2± become saddles. When the parameters are further varied across line L2 from region B5 to B6 , the semi-trivial fixed points E2± collide with the trivial solution and then disappear, and the trivial solution becomes a saddle from a sink. Finally, when the parameters are varied across line L1 from region B6 to B1 , the semi-trivial fixed point E1 collides with the trivial solution and then disappears, and the trivial solution becomes a source from a saddle. Under above parameters values: a1 = −5, a2 = −3, a3 = 3, a4 = 2, b1 = 2.5, b2 = 1.5, k1 = 1, k2 = 1, l1 = 1, and the bifurcation critical parameters k1c = 1 and τc = 3.3365, system (2) undergoes a Hopf-pitchfork bifurcation from the trivial equilibrium. From Fig. 4, in regions B3 , B4 and B5 , there exists stable fixed point. For example, we choose one group of perturbation parameter value: (τ, k ) = (−0.1, −0.1) with initial value being ϕ(θ) = [0.01, 0.05, 0.01, 0.05], θ ∈ [−τ, 0], which belongs to the region B4 , resulting in a stable fixed point shown in Fig. 5. It is clear that the numerical simulations agree with the analytical predictions. The normal form (20) can only describe the local behaviors near the Hopf-pitchfork bifurcation critical point. We concern whether or not the stable fixed point can exist in large regions of unfolding parameters. Therefore, we give a bifurcation diagram for τ = k ∈ (−3.3365, 0), which belongs to region B4 (see Fig. 6). It is plotted by using Matlab software, which depicts the maximum value of solutions of system (2) in the x-axis as the unfolding parameter increases. From this figure, we have seen that the stable fixed point can exist in large regions of unfolding parameters as the unfolding parameters increase away from the critical value in region B4 . Then, system (2) appears globally existed stable periodic solution, whose amplitude increases to a maximum then decreases. When the amplitude decreases to zero, system (2) appears stable fixed point again.

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

589

Fig. 6. Bifurcation diagram for τ = k ∈ (−3.3365, 0), where the ordinate depicts the maximum value of solutions of system (2) in the x-axis.

Remark 2. Hopf-pitchfork bifurcation parameters: k1 = k1c +k and τ = τc +τ with k1c = 1 and τc = 3.3365, thus, we choose τ ∈ (−3.3365, 0) such that τ ≥ 0. Actually, near above Hopf-pitchfork bifurcation critical point, the stable fixed point can exist in large regions of unfolding parameters as the unfolding parameters decrease away from the critical value in region B4 , then system (2) undergoes Hopf bifurcation. On the other side, under above parameters with τ = k = −3.3365, system (2) has a stable equilibrium, then system (2) undergoes Hopf bifurcation as the delay τ and feedback parameter k1 increase, two families of stable periodic solutions, which are resulted from Hopf bifurcation, exist in a large region of delay, and they merge into a family of stable and globally existed periodic solutions. Remark 3. For micro-electromechanical nonlinear coupling system with delay, we need to control the linear feedback strength k1 and feedback delay τ effectively. x and z are angle of rotation, y and w are rotational speed, and we expect them to vary in stable states (stable equilibrium or stable periodic solutions). Otherwise, the torsional vibration of the rotating machinery will affect the normal work and cause the equipment to be damaged. Actually, in accordance with above theoretical analysis, even for different parameter values the delayed micro-electromechanical nonlinear coupling system can be controlled in new stable states by adjusting control parameters, and we give the regions of parameters under which system may have a stable equalized motion or stable periodic motion. Therefore, according to the above theoretical analysis, we can choose proper controller parameters in delayed micro-electromechanical nonlinear coupling system in order to achieve various applications. The method presented in this paper may provide a very good tool in the study of real problems. 5. Conclusion In this paper, we have considered a micro-electromechanical nonlinear coupling system with delay, and analyzed the stability of the trivial equilibrium and the existence of several types of bifurcations. To study dynamical motion, we have derived the normal forms associated with Hopf and Hopf-pitchfork bifurcations by using the multiple time scales method. Two examples, associated with micro-electromechanical nonlinear coupling system with estimated values of parameters, are presented to demonstrate the application of the theoretical results. We show that system (2) exists stable fixed point or stable periodic solution near the

590

Y. Ding et al. / J. Math. Anal. Appl. 461 (2018) 577–590

above critical points. The stable fixed point and stable periodic solution can exist in large regions of unfolding parameters as the unfolding parameters decrease away from the Hopf-pitchfork critical value. Acknowledgments This work was supported by the National Nature Science Foundation of China (No. 11501091), Fundamental Research Funds for the Central Universities (No. 2572017CB22), Heilongjiang Provincial Postdoctoral Scientific Research Foundation (No. LBH-Q17007). References [1] S.A. Campbell, Y. Yuan, Zero singularities of codimension two and three in delay differential equations, Nonlinearity 21 (2008) 2671–2691. [2] Y. Choi, V.G. LeBlanc, Toroidal normal forms for bifurcations in retarded functional differential equations I: multiple Hopf and transcritical/multiple Hopf interaction, J. Differential Equations 227 (2006) 166–203. [3] S.L. Das, A. Chatterjee, Multiple scales without center manifold reductions for delay differential equations near Hopf bifurcations, Nonlinear Dynam. 30 (2002) 323–335. [4] Y. Ding, Dynamic analysis of nonlinear variable frequency water supply system with time delay, Nonlinear Dynam. 90 (2017) 561–574. [5] T. Faria, L.T. Magalhães, Normal form for retarded functional differential equations with parameters and applications to Hopf bifurcation, J. Differential Equations 122 (1995) 181–200. [6] T. Faria, L.T. Magalhães, Normal form for retarded functional differential equations and applications to Bogdanov–Takens singularity, J. Differential Equations 122 (1995) 201–224. [7] B.A. James, C.T. Stephen, W.B. Phillip, Silicon micromechanical devices, Sci. Am. 248 (1983) 44–55. [8] Z. Jiang, W. Ma, J. Wei, Global Hopf bifurcation and permanence of a delayed SEIRS epidemic model, Math. Comput. Simulation 122 (2016) 35–54. [9] S. Liu, S. Zhao, B. Sun, W. Zhang, Bifurcation and chaos analysis of a nonlinear electromechanical coupling relative rotation system, Chin. Phys. B 23 (2014) 094501. [10] S. Liu, S. Zhao, Z. Wang, H. Li, Stability and Hopf bifurcation of a nonlinear electromechanical coupling system with time delay feedback, Chin. Phys. B 24 (2015) 014501. [11] A.H. Nayfeh, Introduction to Perturbation Techniques, Wiley-Interscience, New York, 1981. [12] A.H. Nayfeh, Order reduction of retarded nonlinear systems—the method of multiple scales versus center-manifold reduction, Nonlinear Dynam. 51 (2008) 483–500. [13] Y. Song, T. Zhang, Y. Peng, Turing–Hopf bifurcation in the reaction–diffusion equations and its applications, Commun. Nonlinear Sci. Numer. Simul. 33 (2016) 229–258. [14] Y. Song, H. Jiang, Q. Liu, Y. Yuan, Spatiotemporal dynamics of the diffusive mussel-algae model near Turing–Hopf bifurcation, SIAM J. Appl. Dyn. Syst. 16 (2017) 2030–2062. [15] Y. Song, X. Tang, Stability, steady-state bifurcations, and Turing patterns in a predator–prey model with herd behavior and prey-taxis, Stud. Appl. Math. 139 (2017) 371–404. [16] J.B. Waldner, Nanocomputers and Swarm Intelligence, ISTE John Wiley & Sons, London, 2010. [17] H. Wang, W. Jiang, Hopf-pitchfork bifurcation in van der Pol’s oscillator with nonlinear delayed feedback, J. Math. Anal. Appl. 368 (2010) 9–18. [18] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer, New York, 1990. [19] P. Yu, Y. Ding, W. Jiang, Equivalence of the MTS method and CMR method for differential equations associated with semisimple singularity, Internat. J. Bifur. Chaos 24 (2014) 1450003. [20] F. Zhang, J. Li, C. Zheng, Dynamics of an HBV/HCV infection model with intracellular delay and cell proliferation, Commun. Nonlinear Sci. Numer. Simul. 42 (2017) 464–476.