- Email: [email protected]
Periodic and chaotic motions of a harmonically forced piecewise linear system
ARTICLE IN PRESS
International Journal of Mechanical Sciences 46 (2004) 1807–1825 www.elsevier.com/locate/ijmecsci
Periodic and chaotic motions of a harmonically forced piecewise linear system J.C. Ji, A.Y.T. Leung Department of Building and Construction, City University of Hong Kong, Kowloon, Hong Kong Received 1 July 2003; received in revised form 13 October 2004; accepted 19 October 2004
Abstract The dynamics of a harmonically excited single degree-of-freedom linear system with a feedback control, in which the actuator is subjected to dead zone and saturation constraints, is investigated in detail. The controlled system is mathematically modeled by a set of three piecewise linear equations. It is found that the system may exhibit nine types of symmetric and asymmetric period-one motions, which are characterized by a different number of crossing dead zone and saturation region per cycle. A solution for the symmetric period-one motion with a doubly crossing dead zone and saturation region is analytically constructed and its stability characteristics is examined. Other types of dynamic response such as sub-harmonic periodic motions and chaotic motions, found through numerical simulations, are also presented. r 2004 Elsevier Ltd. All rights reserved. Keywords: Periodic motions; Chaotic motions; Piecewise linear system; Dead zone; Saturation
1. Introduction Various non-linearities of the structural and control elements are frequently encountered in modeling many of controlled systems. Typical nonlinearities in such systems are the non-smooth nature of actuator devices such as dead zone, saturation, relay and hysteresis. These Corresponding author. Department of Mechanical Engineering, The University of Adelaide, Adelaide, SA 5005,
Australia. Tel.: +61 8 8303 6941; fax: +61 8 8303 4367. E-mail addresses: [email protected] (J.C. Ji), [email protected] (A.Y.T. Leung). 0020-7403/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2004.10.005
ARTICLE IN PRESS 1808
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
non-linearities, which were often neglected in designing and modeling of the control system, may deteriorate the performance of a feedback control system and lead to a poor dynamical behavior of the controlled system [1]. Dead zone, which describes the component’s insensitivity to small signals, is a static input–output relationship which for a range of input values gives no output. Once the output appears, the slope between the input and output is constant. Saturation may be caused by the limitation of the physical devices. When an input signal reaches a certain value, further increases in the input signal produce only the same output, which is usually called the saturation value. Investigation of the effect of dead zone and saturation is of great importance for the control systems designers due to the fact that these nonlinearities exist in a large number of actuators [2]. On the other hand, a thorough understanding of the dynamical behavior can be of prime importance to design the system and to provide useful information for the fault diagnosis. The main objective of the present paper is to investigate the dynamic behavior of the controlled system with dead zone and saturation non-linearities. For simplicity, the non-linearities considered in the present paper are assumed to be approximated by piecewise linear characteristics. Due to the presence of dead zone and saturation constraints, equations of motion of the controlled system appear in a strongly nonlinear form with piecewise linear characteristic, which is switched by the displacement and velocity constraints. Though the dynamical behavior of piecewise linear systems with a bilinear or tri-linear characteristic has been an active topic of intensive research e.g. [3–17], the dynamics of a piecewise linear system with quintuple linear elements has not yet been reported in the literature to the authors’ knowledge. It is found that the system considered here may exhibit symmetric and asymmetric period-one motions. Period-one motions may lose their stability and give their places to sub-harmonic or chaotic motions. In addition, sub-harmonic motions may coexist with a period-one motion, and a chaotic motion may coexist with two asymmetric period-one motions. The present paper is organized as follows; in Section 2, the equations of motion are derived. An analytical method is developed to determine period-one motions in Section 3, and their stability characteristics are examined in Section 4. Other types of motions such as sub-harmonic periodic motions and chaotic motions, are discussed in Section 5. Conclusions are presented in Section 6.
2. Equations of motion Consider a simple linear feedback control system in which a periodic excitation (disturbance) is presented, as shown in Fig. 1. The equation governing the motion of the plant (i.e., a singledegree-of-freedom linear oscillator) can be written as Mx00 þ Cx0 þ Kx ¼ F cos OT U;
Fig. 1. Block diagram of a linear system with feedback control.
(1)
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1809
where M is the mass, C is the damping coefficient, K is the stiffness, F and O are the amplitude and frequency of the excitation, U denotes the output of the actuator, and a prime denotes differentiation with respect to the time T. In the presence of dead zone and saturation constraints, the output of the actuator, as shown in Fig. 2, is of the form 8 if V XV s ; K a ðV s V d Þ; > > > > > if V d pVpV s ; > < K a ðV V d Þ; if V d pVpVd ; U ¼ 0; (2) > > > K a ðV þ V d Þ; if V s pV p V d ; > > > : K a ðV s V d Þ; if V p V s ; where K a denotes the actuator gain, V d and V s denote the maximum threshold values of dead zone and saturation. For simplicity, the feedback control is assumed to employ a PD-controller with the form V ¼ K p x þ K d x0 ; where K p and K d are the proportional and derivative gains, respectively. Substituting Eq. (2) into Eq. (1) and then introducing nondimensional variables x ¼ yV d and t ¼ O1 T; yields the equations of motion in non-dimensional form y€ þ 2dy_ þ y ¼ f cos ot;
_ for jkp y þ kd yjp1;
y€ þ 2my_ þ ð1 þ kÞy f d sgnðyÞ ¼ f cos ot; y€ þ 2dy_ þ y þ f s sgnðyÞ ¼ f cos ot;
_ for 1pjkp y þ kd yjpu s;
_ for jkp y þ kd yjXu s;
(3a) (3b) (3c)
where, O21 ¼ K=M; 2d ¼ C=MO1 ; o ¼ O=O1 ; f ¼ F =MV d O21 ; f d ¼ K a =M O21 ; f s ¼ K a ðV s V d Þ=MO21 V d ; ka ¼ K a =M O21 ; kp ¼ K p ; kd ¼ K d O1 ; k ¼ ka kp ; us ¼ ðV s V d Þ=V d ; 2m ¼ 2d þ ka kd ; sgn( ) denotes the sign function, and a superscript dot indicates differentiation with respect to the time t.
Fig. 2. Output of the actuator subjected to dead zone and saturation constraints.
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1810
It is easy to note that system (3) is continuous and satisfies the Lipschitz condition. Therefore, the solution of Eq. (3) exists and is unique. However, except for a small amplitude response, the analytical solutions of Eq. (3) cannot be obtained in closed form, since, as will be seen, the exact times when the input of the actuator crosses the boundaries of dead zone and saturation region cannot be found explicitly. These time instances will be subsequently referred to as the crossing times. In the next section, a method is developed to determine various types of period-one motions.
3. Period-one motions The response of the overall system given by Eq. (3) may exhibit nine types of period-one motions, which are categorized by the characteristics of motions and the number of crossing dead zone and saturation region per cycle, as shown in Fig. 3. The period-one response will be subsequently referred to as small amplitude motion without crossing dead zone, symmetric and asymmetric motions with a doubly crossing dead zone and saturation region, asymmetric motions with a doubly crossing dead zone and singly crossing saturation region, symmetric and asymmetric motions with a doubly crossing dead zone, respectively. 3.1. Small amplitude motion For a small amplitude motion as shown in Fig. 3a, the actuator does not give any output. The response of the oscillator is linear and the steady state solution is of the form yðtÞ ¼ f 1 cosðot yÞ;
(4)
where
2do f1 ¼ ; y ¼ arctan : 1 o2 ½ð1 o2 Þ2 þ 4d2 o2 1=2 f
The critical value of the excitation amplitude at which the actuator gives an output, can be _ which is given by f c ¼ ð½ð1 o2 Þ2 þ obtained by applying the condition jkp y þ kd yjX1; 2 2 1=2 2 2 2 1=2 4d o Þ=ð½kp þ o kd Þ: When the excitation amplitude is larger than f c ; the feedback control will apply a control output to the plant. The motion of the system may cross the boundaries of dead zero and/or saturation regions and then leave and enter dead zone and/or saturation regions at least once over a period of period-one response. 3.2. Motions with a doubly crossing dead zone and saturation region A symmetric or an asymmetric periodic motion with a doubly crossing dead zone and saturation region, as shown in Figs. 3b, c, and d, may appear when the excitation amplitude is much larger than f c : The input signal of the actuator and the time trajectory of the response (as shown in Fig. 3b) are entirely composed of eight segments according to the time intervals: ½t0 ; t1 ; ½t1 ; t2 ; ½t2 ; t3 ; ½t3 ; t4 ; ½t4 ; t5 ; ½t5 ; t6 ; ½t6 ; t7 and ½t7 ; t8 ; where t8 ¼ t0 þ 2p o and ti represents the time instant that the input signal passes the boundaries of dead zero and saturation regions. An
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1811
Fig. 3. Time histories of the system response (left subfigure) and the corresponding input signal of the actuator (right subfigure) for nine types of period-one motions: (a) a small amplitude motion without crossing dead zone; (b) a symmetric motion with a doubly crossing dead zone and saturation region (where 1; 2; 9 denote the points at which the stiffness and damping discontinuities occur); (c) and (d) asymmetric motions with a doubly crossing dead zone and saturation region; (e) and (f) asymmetric motions with a doubly crossing dead zone and singly crossing saturation region; (g) asymmetric motion with a doubly crossing dead zone; (h) and (i) asymmetric motions with a doubly crossing dead zone.
ARTICLE IN PRESS 1812
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Fig. 3. (Continued)
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1813
Fig. 3. (Continued)
analytical construction of the solution for each interval is based on the fact that Eq. (3) is linear in _ _ _ 1pjkp y þ kd yjpu each of the regions jkp y þ kd yjp1; s ; and jk p y þ kd yjXu s ; respectively. In the present paper, emphasis is placed on the analytical determination of the symmetric period-one motion with a doubly crossing dead zone and saturation region. Other types of period-one solutions can be obtained in a similar procedure. Because of symmetry of the motion examined, only four segments of the periodic solution need to be considered. _ y2 ðtÞ from t1 in Let y1 ðtÞ denote the solution starting from t0 in the region jkp y þ kd yjp1; _ _ _ 1pðkp y þ kd yÞpu s ; y3 ðtÞ from t2 in ðk p y þ k d yÞXu s ; and y4 ðtÞ from t3 in 1pðkp y þ k d yÞpu s; respectively, then the general solutions for these four regions can be written as y1 ðtÞ ¼ edðtt0 Þ ½A1 cos bðt t0 Þ þ B1 sin bðt t0 Þ þ G1 sin ot þ H 1 cos ot _ for jkp y þ kd yjp1;
ð5aÞ
y2 ðtÞ ¼ emðtt1 Þ ½A2 cos qðt t1 Þ þ B2 sin qðt t1 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 _ for 1pðkp y þ kd yÞpu s;
ð5bÞ
y3 ðtÞ ¼ edðtt2 Þ ½A3 cos bðt t2 Þ þ B3 sin bðt t2 Þ þ G1 sin ot þ H 1 cos ot þ Y 2 _ for ðkp y þ kd yÞXu s;
ð5cÞ
y4 ðtÞ ¼ emðtt3 Þ ½A4 cos qðt t3 Þ þ B4 sin qðt t3 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 _ for 1pðkp y þ kd yÞpu s;
ð5dÞ
where, Ai ; Bi ; and ti1 (i ¼ 1; 2; 3; 4) are the unknown constants to be determined, pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ¼ 1 d2 ; q ¼ 1 þ k m2 ; Y 1 ¼ f d =ð1 þ kÞ sgn ðyÞ; Y 2 ¼ f s sgnðyÞ; G1 ¼
2dof ð1 o2 Þf ; H ¼ ; 1 ð1 o2 Þ2 þ ð2doÞ2 ð1 o2 Þ2 þ ð2doÞ2
G2 ¼
2mof ð1 þ k o2 Þf ; H ¼ : 2 ð1 þ k o2 Þ2 þ ð2moÞ2 ð1 þ k o2 Þ2 þ ð2moÞ2
ARTICLE IN PRESS 1814
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
The eight constants Ai and Bi ; and four crossing times ti1 ; can be determined by applying a set of initial conditions and periodicity, continuity as well as symmetry conditions, which can be written as; at point 1;
kp y1 ðt0 Þ þ kd y_ 1 ðt0 Þ ¼ 1;
(6a)
at point 2;
kp y2 ðt1 Þ þ kd y_ 2 ðt1 Þ ¼ 1;
(6b)
y1 ðt1 Þ ¼ y2 ðt1 Þ; at point 3;
(6e)
y_ 2 ðt2 Þ ¼ y_ 3 ðt2 Þ;
(6f,g)
kp y4 ðt3 Þ þ kd y_ 4 ðt3 Þ ¼ us ;
(6h)
y_ 3 ðt3 Þ ¼ y_ 4 ðt3 Þ;
(6i,j)
y3 ðt3 Þ ¼ y4 ðt3 Þ; at point 5;
(6c,d)
kp y3 ðt2 Þ þ kd y_ 3 ðt2 Þ ¼ us ;
y2 ðt2 Þ ¼ y3 ðt2 Þ; at point 4;
y_ 1 ðt1 Þ ¼ y_ 2 ðt1 Þ;
y4 ðt4 Þ ¼ y1 ðt0 Þ;
y_ 4 ðt4 Þ ¼ y_ 1 ðt0 Þ;
and t4 t0 ¼ p=o:
(6k,l) (6m)
For brevity, Eqs. 6(a), (b), (e) and (h) are subsequently referred to as the initial conditions, Eqs. 6(c), (f), (i) and (k) as the displacement matching conditions, and equations 6(d), (g), (j) and (l) as the velocity matching conditions, respectively. The constants, Ai and Bi ; can be determined by solving the equations that result from implementing four initial conditions and four displacement matching conditions. For the clarity of notation, these constants are expressed as Ai ¼ Ai0 þ Aic cos ot0 þ Ais sin ot0 ; Bi ¼ Bi0 þ Bic cos ot0 þ Bis sin ot0 ;
ð7Þ
where the coefficients Ai0 ; Aic ; Ais ; Bi0 ; Bic and Bis are given in Appendix A1. Implementing the four velocity matching conditions gives rise to a set of four transcendental equations for three time intervals ðti1 t0 Þ and time t0 as follows: Di1 cos ot0 þ Di2 sin ot0 ¼ Di0 ;
i ¼ 1; 2; 3; 4;
(8a2d)
where the coefficients ,Di1 ; Di2 and Di0 ; are listed in Appendix A2. Eq. (8) involves the system parameters and four unknown crossing times only. No analytic solutions can be found, thus a numerical mean is resorted to solve three time intervals ðti1 t0 Þ and time t0 : Then constants Ai and Bi can be obtained from Eq. (7) by back substitutions. After getting Ai ; Bi and crossing times ti1 ; the corresponding histories of yi can be determined by Eq. (5). For an asymmetric motion as shown in Fig. 3c or Fig. 3d, a determination of its solution can be performed by the same procedure. In these cases, eight general solutions are needed to determine for eight distinct segments. The unknowns are sixteen constants Ai ; Bi and eight crossing times ti1 (where i ¼ 1; 2; . . . 8:), which can be determined by imposing a set of initial, matching and periodicity conditions. A set of 24 conditions similar to the forms given by Eq. (6) can be written
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1815
at nine points, where the damping and stiffness discontinuities occur. An implementation of these conditions may lead to a set of eight transcendental equations for the eight crossing times, which are solved using a numerical procedure. The constants Ai ; Bi can be then evaluated by simple back substitutions, and hence the corresponding solutions yi can be determined. 3.3. Asymmetric motions with a doubly crossing dead zone and singly crossing saturation region The asymmetric motions like those shown in Figs. 3e and f, consist of six distinct segments, i.e., six intervals in which the damping and stiffness properties remain constant. It is thus needed to determine six general solutions similar to the forms given by Eq. (5) for the six segments. The unknowns are twelve constants Ai ; Bi and six crossing times ti1 ; where i ¼ 1; 2; 6: Similarly, these constants can be determined by implementing a set of initial and matching conditions. Totally eighteen independent conditions, which are not given in the present paper for brevity, can be written at seven points where the discontinuities occur. Imposing these conditions leads to an equal number of transcendental equations, which can be further reduced to a set of six transcendental equations in a form similar to Eq. (8). Numerical solutions of these six equations determine the six crossing times. Then the constants, Ai and Bi ; and the corresponding solutions yi can be determined by simple back substitutions. 3.4. Motions with a doubly crossing dead zone For a symmetric motion like that shown in Fig. 3g, its solution can be constructed based on two distinct segments. Two general solutions for the two segments are needed to determine. The corresponding constants Ai ; Bi and two crossing times ti1 (where i ¼ 1; 2) can be determined by implementing a set of six conditions. For asymmetric solutions as shown in Figs. 3h and i, the motions consist of four distinct segments. It is naturally required to determine four general solutions for the four segments. The unknowns are eight constants Ai ; Bi and four crossing times ti1 ; where i ¼ 1; 2; 3; 4: Totally a set of twelve independent conditions can be written at five points. Then a similar procedure developed in Section 3.2 can be performed to determine these solutions.
4. Stability of the period-one motions For a small amplitude motion like that Discussed in Section 3.1, its stability analysis can be performed using a classical method. Evidently, the motion is stable, provided that the damping coefficient in Eq. 3(a) is positive. For other types of periodic motions, their stability characteristics can be examined by considering the asymptotic behavior of the perturbations to the steady state solutions. In this section, a stability analysis is presented for the symmetric periodic solution with a doubly crossing dead zone and saturation region. A similar procedure can be performed to investigate the stability of other periodic solutions. The symmetry of the motion examined here ensures that its stability can be fully determined in one half of the solution.
ARTICLE IN PRESS 1816
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Let zi ðtÞ denote the corresponding perturbed solutions to yi ðtÞ (where i ¼ 1; 2; 3; 4), respectively, then the perturbed solutions can be written as
z1 ðtÞ ¼ edðtt0 Dt0 Þ P1 cos bðt t0 Dt0 Þ þ Q1 sin bðt t0 Dt0 Þ þ G1 sin ot þ H 1 cos ot; for jkp z þ kd z_jp1; ð9aÞ
z2 ðtÞ ¼ emðtt1 Dt1 Þ P2 cos qðt t1 Dt1 Þ þ Q2 sin qðt t1 Dt1 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 ; for 1pðkp z þ kd z_Þpus ;
ð9bÞ
z3 ðtÞ ¼ edðtt2 Dt2 Þ P3 cos bðt t2 Dt2 Þ þ Q3 sin bðt t2 Dt2 Þ þ G1 sin ot þ H 1 cos ot þ Y 2 ; for ðkp z þ kd z_ÞXus ;
ð9cÞ
z4 ðtÞ ¼ emðtt3 Dt3 Þ P4 cos qðt t3 Dt3 Þ þ Q4 sin qðt t3 Dt3 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 ; for 1pðkp z þ kd z_Þpus ;
ð9dÞ
where Pi and Qi ; (i ¼ 1; 2; 3; 4), are the constants to be determined. The initial conditions for the perturbed motions of the first segment can be written as z1 ðt0 þ Dt0 Þ ¼ z0 þ Dz0 ; kp z1 ðt0 þ Dt0 Þ þ kd z_1 ðt0 þ Dt0 Þ ¼ 1; where operator D indicates a small perturbation of the operand. At moment ðt1 þ Dt1 Þ; the input signal of the actuator reaches a boundary of dead zone and enters the region 1pðkp z þ kd z_Þpus thereafter. The perturbed response at the moment of departure from dead zone is assumed to be z1 ðt1 þ Dt1 Þ ¼ z1 þ Dz1 ; kp z1 ðt1 þ Dt1 Þ þ kd z_1 ðt1 þ Dt1 Þ ¼ 1:
(10)
Since the perturbations in the initial conditions are assumed to be small, it is expected that the ¯ 1 and coefficients P1 and Q1 will assume the values close to those of the unperturbed solutions, A ¯ 1 ; respectively. To the first order, they can be expressed as B ¯ 1 þ P11 Dt0 þ P12 Dz0 ; P1 ¼ A ¯ 1 þ Q11 Dt0 þ Q12 Dz0 ; Q1 ¼ B
ð11Þ
¯ 1 ¼ z0 G 1 sin ot0 H 1 cos ot0 ;P11 ¼ oH 1 sin ot0 oG 1 cos ot0 ; P12 ¼ 1; where, A ¯ 1 ¼ 1 þ kd b kp z0 dH 1 þ oG 1 cos ot0 þ oH 1 dG 1 sin ot0 ; B kd b kd b b b Q11 ¼ ½oðoG 1 þ dH 1 Þ sin ot0 þ oðoH 1 dG1 Þ cos ot0 =b; Q12 ¼ ðkd d kp Þ=kd b: Substituting Eq. (9a) into Eq. (10), and keeping only the first-order terms yields Dz1 ¼ M 11 Dt1 þ M 12 Dt0 þ M 13 Dz0 ; M 21 Dt1 þ M 22 Dt0 þ M 23 Dz0 ¼ 0; where the coefficients M ij are defined in Appendix A3. For brevity, Eq. (12) can be put in a matrix form as Dt1 Dt0 ¼ R1 ; Dz1 Dz0
ð12Þ
(13)
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1817
where R1 is a 2 2 matrix with the elements r11 ¼ M 22 =M 21 ; r12 ¼ M 23 =M 21 ; r21 ¼ M 12 M 11 M 22 =M 21 ; r22 ¼ M 13 M 11 M 23 =M 21 ; respectively. The asymptotic behavior of the perturbed solutions for other segments from ðt1 þ Dt1 Þ to ðt2 þ Dt2 Þ (for the second segment), from ðt2 þ Dt2 Þ to ðt3 þ Dt3 Þ (for the third segment), and from ðt3 þ Dt3 Þ to ðt4 þ Dt4 Þ (for the fourth segment), can be investigated in a similar procedure. After a lengthy algebraic manipulation, the following equations retaining the first order terms can be obtained Dti1 Dti ¼ Ri ; i ¼ 2; 3; 4; (14) Dzi Dzi1 where Ri are three 2 2 matrices, whose elements are not given here for brevity. The perturbation to the steady state solution during the first half period is then obtained by combining Eq. (13) and (14) to form an equation as Dt0 Dt4 ¼J ; (15) Dz4 Dz0 where J ¼ R1 R2 R3 R4 ; represents the transition matrix of the response from ðt0 þ Dt0 Þ to ðt4 þ Dt4 Þ: The eigenvalues of matrix J determines the stability of the solution. The symmetric periodone solution is asymptotically stable if both eigenvalues of J have modulus less than unity. When either of the two eigenvalues has modulus greater than one the solution is unstable. For some combinations of the system parameters, the modulus may take the value of unity, where a bifurcation occurs. To validate the present analysis, the analytical solutions of the symmetric period-one motion with a doubly crossing dead zone and saturation region were compared with those obtained from direct numerical integration of Eq. (3) for several cases. It was found that the solutions obtained by these two different means were identical. As an illustrative example, consider a system with parameters in Eq. (3) as follows; d ¼ 0:1; m ¼ 0:1; k ¼ 4:0; kp ¼ 4:0; us ¼ 4:0; o ¼ 2:0; and f ¼ 1:3: The crossing times obtained via the analysis are t0 ffi 0:405825; t1 ffi 0:535464; t2 ffi 0:740362; and t3 ffi 1:771968: The time intervals between two crossing times given by numerical integration of Eq. (3) were nearly identical to those predicted by the analysis. A very small difference only in the sixth decimal place was found between the values obtained by two approaches. Fig. 4 shows the phase portrait of the symmetric solutions obtained by the analysis and numerical integration of Eq. (3). The analytical solution is in excellent agreement with the result of numerical integration.
5. Other types of motions The analysis presented in Sections 3 and 4 can be extended to construct sub-harmonic periodic solutions. However, the procedure for finding such solutions is exceptionally lengthy and involves much more algebraic manipulation. In this section, results of numerical integration are given to show other behavior of the system. The widely used software package MATLAB/Simulink is employed to model system (3). The dead zone block and the saturation block in the library of the built-in Simulink nonlinear
ARTICLE IN PRESS 1818
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Fig. 4. Phase portrait of the symmetric motion with a doubly crossing dead zone and saturation region; solid curve indicates results of numerical integration and circles represent analytical solutions.
components are used to impose the starting and end values of dead zone and the upper and lower bounds of saturation on the input signal, respectively. The fourth-order Runge–Kutta procedure is chosen to be the integration method of the ODE solvers. The initial conditions for the integrator blocks are all set zero, unless otherwise stated. The simulation results are observed using scopes and other display blocks during the simulation. Extensive numerical simulations have been performed in order to explore the complex dynamical behavior. For the results given in this paper, the parameters are chosen as follows: d ¼ 0:025; ka ¼ 1:2; kp ¼ 3:0; kd ¼ 0:05; o ¼ 0:9; us ¼ 4:0; unless otherwise specified. In addition to period-one motions, a variety of sub-harmonic and chaotic motions may be observed. The minimum period of a sub-harmonic motion is an integral multiple of the forcing period. This integer indicates the order of the sub-harmonic motion. Sub-harmonic periodic motions of order two, three, four, five, seven and eight were found to exist for some combinations of the system parameters. Sub-harmonic motions may appear not only in the regions where the input of the actuator does not reach the boundaries of saturation regions, but also in the regions where the input signal exceeds the saturation value. Sub-harmonic motions seemed to disappear and period-one motions appear at higher damping levels. Chaotic motions, which result from either a period-one motion or a subharmonic motion losing its stability, are another widely observed complex dynamical behavior [18]. By taking the excitation amplitude as the bifurcation parameter, the bifurcation diagrams of the response are shown in Fig. 5. The system parameters and mapping conditions for sampling data are identical for two diagrams except initial conditions. The data for plotting Fig. 5a were obtained via a usual way to set initial conditions, while for Fig. 5b the initial conditions ðy0 ; y_ 0 Þ ¼ ð1:0; 0:6Þ remain unchanged at each step of the excitation amplitude. It can be easily seen in Fig. 5 that chaotic motions exist at many values of the excitation amplitude. Fig. 6 displays a typical chaotic motion at f ¼ 5:4; where (a) is the Poincare map at the excitation phase p4 and (b) is the phase portrait. Another interesting dynamical behavior, which widely exists in the system response, is a coexistence of different types of motions starting from different sets of initial conditions for a given set of system parameters. The particular dependence of motion on the initial conditions can
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1819
Fig. 5. Bifurcation diagrams of the response; (a) the usual way to set initial conditions, (b) the initial conditions at each step are ðy0 ; y_ 0 Þ ¼ ð1:0; 0:6Þ:
Fig. 6. The chaotic motion of the system at f ¼ 5:4 under initial conditions ðy0 ; y_ 0 Þ ¼ ð0:0; 1:0Þ; (a) Poincare map, (b) phase portrait.
be also observed in Fig. 5 for many values of the excitation amplitude, for example, when f ¼ 3:25; 4:15 4:7; 6:0: It has been numerically found that a period-one motion may coexist with one or two sub-harmonic motions, or that two asymmetric period-one motions coexist with a chaotic motion. For instance, two sub-harmonic motions of order two may coexist with a periodone motion. A sub-harmonic motion of order three or five was found to coexist with a period-one motion. Due to laborious calculation being involved, a determination of the domain of attraction
ARTICLE IN PRESS 1820
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Fig. 7. Phase portraits of the response for the co-existence of a period-one motion and a subharmoinc motion of order seven as well as of order five; (a) period-one motion starting from ðy0 ; y_ 0 Þ ¼ ð0:0; 1:0Þ; (b) subharmonic motion of order seven from (1.0,0.6), and (c) subharmonic motion of order five from ð0:6; 0:3Þ:
of various solutions is not considered in the present paper. Alternatively, some examples are given to illustrate the coexisting solutions. A sub-harmonic motion of order five and a subharmonic motion of order seven may coexist with a period-one motion when f ¼ 6:0; as shown in Fig. 7. Fig. 8 illustrates the coexistence of two asymmetric period-one motions and a chaotic motion for f ¼ 4:55: Each of the asymmetric motions is the other’s symmetric image about the origin. Starting from initial conditions ðy0 ; y_ 0 Þ ¼ ð0:0; 1:0Þ leads to the chaotic motion, while from ð0:6; 0:3Þ and ð1:0; 0:6Þ result in two asymmetric motions, respectively. The presence of coexistence of multiple stable motions and of chaotic motions causes a difficulty in predicting the system dynamics.
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1821
Fig. 8. Coexistence of a pair of asymmetric period-one motions and a chaotic motion; (a) phase portraits of two asymmetric period-one motions starting from ðy0 ; y_ 0 Þ ¼ ð0:6; 0:3Þ for the solid curve and ð1:0; 0:6Þ for the dashed curve, (b) Poincare map of the chaotic motion starting from ð0:0; 1:0Þ:
6. Conclusions The dynamics of a periodically forced single degree-of-freedom piecewise linear system with quintuple linear segments has been investigated. A detailed procedure has been performed to construct a solution for the symmetric period-one motion with a doubly entering dead zone and saturation region per cycle. The stability of the solution has been studied by examining the asymptotic behavior of the corresponding perturbed solution. It was found that the response of the system may exhibit symmetric and asymmetric period-one motions, sub-harmonic periodic motions, chaotic motions, and coexistence of three stable motions.
Appendix: List of the coefficients A1. The coefficients in Eq. (7) A10 ¼ ða81 h60 þ a82 h70 þ b80 Þ=ad; A1c ¼ ða81 h61 þ a82 h71 þ b81 Þ=ad; A1s ¼ ða81 h62 þ a82 h72 þ b82 Þ=ad; A20 ¼ h20 þ g2 A10 ; A2c ¼ h21 þ g2 A1c ; A2s ¼ h22 þ g2 A1s ; A30 ¼ h40 þ g4 A10 ; A3c ¼ h41 þ g4 A1c ; A3s ¼ h42 þ g4 A1s ; A40 ¼ h60 þ g6 A10 ; A4c ¼ h61 þ g6 A1c ; A4s ¼ h62 þ g6 A1s ; Bk0 ¼ hj0 þ gj A10 ; Bkc ¼ hj1 þ gj A1c ; Bks ¼ hj2 þ gj A1s ; k ¼ 1; 2; 3; 4; j ¼ 2k 1; where, e110 ¼ edðt1 t0 Þ ; e120 ¼ edðt2 t0 Þ ; e130 ¼ edðt3 t0 Þ ; e210 ¼ emðt1 t0 Þ ; e220 ¼ emðt2 t0 Þ ;
ARTICLE IN PRESS 1822
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825 p=o
e230 ¼ emðt3 t0 Þ ; e20 ¼ em ; s2 ¼ sin q
p p ; c2 ¼ cos q ; s10 ¼ sin oðt1 t0 Þ; o o
c10 ¼ cos oðt1 t0 Þ; s20 ¼ sin oðt2 t0 Þ; c20 ¼ cos oðt2 t0 Þ; s30 ¼ sin oðt3 t0 Þ; c30 ¼ cos oðt3 t0 Þ; s110 ¼ sin bðt1 t0 Þ; c110 ¼ cos bðt1 t0 Þ; s120 ¼ sin bðt2 t0 Þ; c120 ¼ cos bðt2 t0 Þ; s130 ¼ sin bðt3 t0 Þ; c130 ¼ cos bðt3 t0 Þ; s210 ¼ sin qðt1 t0 Þ; c210 ¼ cos qðt1 t0 Þ; s220 ¼ sin qðt2 t0 Þ; c220 ¼ cos qðt2 t0 Þ; s230 ¼ sin qðt3 t0 Þ; c230 ¼ cos qðt3 t0 Þ; a11 ¼ kp kd d; a12 ¼ kd b; b11 ¼ kp H 1 þ okd G1 ; b12 ¼ kp G 1 okd H 1 ; b10 ¼ 1; a21 ¼ kp kd m; a22 ¼ qkd ; b21 ¼ kp H 2 þ okd G 2 ; b22 ¼ kp G 2 okd H 2 ; b20 ¼ 1 kp Y 1 ; a31 ¼ e110 c110 ; a32 ¼ e110 s110 ; b31 ¼ ðG1 G2 Þs10 þ ðH 1 H 2 Þc10 ; b32 ¼ ðG1 G2 Þc10 ðH 1 H 2 Þs10 ; b30 ¼ Y 1 ; a41 ¼ kp kd d; a42 ¼ kd b; b41 ¼ ðkp G 1 okd H 1 Þs20 þ ðkp H 1 þ okd G 1 Þc20 ; b42 ¼ ðkp G1 okd H 1 Þc20 ðkp H 1 þ okd G 1 Þs20 ; b40 ¼ us kp Y 2 ; a51 ¼ e220 ðc220 c210 þ s220 s210 Þ=e210 ; a52 ¼ e220 ðs220 c210 c220 s210 Þ=e210 ; b51 ¼ ðG2 G1 Þs20 þ ðH 2 H 1 Þc20 ; b52 ¼ ðG2 G1 Þc20 ðH 2 H 1 Þs20 ; b50 ¼ Y 2 Y 1 ; a61 ¼ kp mkd ; a62 ¼ qkd ; b61 ¼ ðkp G 2 okd H 2 Þs30 þ ðkp H 2 þ okd G 2 Þc30 ; b62 ¼ ðkp G2 okd H 2 Þc30 ðkp H 2 þ okd G 2 Þs30 ; b60 ¼ us kp Y 1 ; a71 ¼ e130 ðc130 c120 þ s130 s120 Þ=e120 ; a72 ¼ e130 ðs130 c120 c130 s120 Þ=e120 ; b71 ¼ ðG1 G2 Þs30 þ ðH 1 H 2 Þc30 ; b72 ¼ ðG1 G2 Þc30 ðH 1 H 2 Þs30 ; b70 ¼ Y 1 Y 2 ; a81 ¼ e20 ðc2 c230 þ s2 s230 Þ=e230 ; a82 ¼ e20 ðs2 c230 c2 s230 Þ=e230 ; b81 ¼ H 1 H 2 ; b82 ¼ G 1 G 2 ; b80 ¼ Y 1 ; l 11 ¼ e110 ðc110 d þ s110 bÞ; l 12 ¼ e110 ðc110 b s110 dÞ; m11 ¼ oðG 1 G2 Þc10 oðH 1 H 2 Þs10 ; m12 ¼ oðG2 G 1 Þs10 þ oðH 2 H 1 Þc10 ; l 21 ¼ e220 ½mðc220 c210 þ s220 s210 Þ þ qðs220 c210 c220 s210 Þ =e210 ; l 22 ¼ e220 ½qðc220 c210 þ s220 s210 Þ mðs220 c210 c220 s210 Þ =e210 ; m21 ¼ oðG1 G 2 Þc20 þ oðH 1 H 2 Þs20 ; m22 ¼ oðG 1 G 2 Þs20 þ oðH 1 H 2 Þc20 ; l 31 ¼ e130 ½dðc130 c120 þ s130 s120 Þ þ bðs130 c120 c130 s120 Þ =e120 ; l 32 ¼ e130 ½bðc130 c120 þ s130 s120 Þ dðs130 c120 c130 s120 Þ =e120 ;
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
m31 ¼ oðG 1 G2 Þc30 þ oðH 2 H 1 Þs30 ; m32 ¼ oðG2 G 1 Þs30 þ oðH 2 H 1 Þc30 ; l 41 ¼ e20 ½mðc2 c230 þ s2 s230 Þ þ qðs2 c230 c2 s230 Þ =e230 ; l 42 ¼ e20 ½qðc2 c230 þ s2 s230 Þ mðs2 c230 c2 s230 Þ =e230 ; m41 ¼ oðG 1 G2 Þ; m42 ¼ oðH 2 H 1 Þ; g1 ¼ a11 =a12 ; h10 ¼ b10 =a12 ; h11 ¼ b11 =a12 ; h12 ¼ b12 =a12 ; g2 ¼ a31 þ a32 g1 ; h20 ¼ a32 h10 b30 ; h21 ¼ a32 h11 þ b31 ; h22 ¼ a32 h12 þ b32 ; g3 ¼ a21 g2 =a22 ; h30 ¼ ðb20 a21 h20 Þ=a22 ; h31 ¼ ðb21 þ a21 h21 Þ=a22 ; h32 ¼ ðb22 þ a21 h22 Þ=a22 ; g4 ¼ a51 g2 þ a52 g3 ; h40 ¼ a51 h20 þ a52 h30 b50 ; h41 ¼ a51 h21 þ a52 h31 þ b51 ; h42 ¼ a51 h22 þ a52 h32 þ b52 ; g5 ¼ a41 g4 =a42 ; h50 ¼ ðb40 a41 h40 Þ=a42 ; h51 ¼ ðb41 þ a41 h41 Þ=a42 ; h52 ¼ ðb42 þ a41 h42 Þ=a42 ; g6 ¼ a71 g4 þ a72 g5 ; h60 ¼ a71 h40 þ a72 h50 b70 ; h61 ¼ a71 h41 þ a72 h51 þ b71 ; h62 ¼ a71 h42 þ a72 h52 þ b72 ; g7 ¼ a61 g6 =a62 ; h70 ¼ ðb60 a61 h60 Þ=a62 ; h71 ¼ ðb61 þ a61 h61 Þ=a62 ; h72 ¼ ðb62 þ a61 h62 Þ=a62 ; ad ¼ 1 þ a81 g6 þ a82 g7 :
A2. The coefficients in Eq. (8) D11 ¼ l 11 A1c þ l 12 B1c þ m11 qB2c þ mA2c ; D12 ¼ l 11 A1s þ l 12 B1s þ m12 qB2s þ mA2s ; D10 ¼ qB20 mA20 l 11 A10 l 12 B10 ; D21 ¼ l 21 A2c þ l 22 B2c þ m21 B3c b þ A3c d; D22 ¼ l 21 A2s þ l 22 B2s þ m22 B3s b þ A3s d; D20 ¼ B30 b A30 d l 21 A20 l 22 B20 ; D31 ¼ l 31 A3c þ l 32 B3c þ m31 qB4c þ mA4c ; D32 ¼ l 31 A3s þ l 32 B3s þ m32 qB4s þ mA4s ; D30 ¼ qB40 mA40 l 31 A30 l 32 B30 ; D41 ¼ l 41 A4c þ l 42 B4c þ m41 A1c d þ B1c b; D42 ¼ l 41 A4s þ l 42 B4s þ m42 A1s d þ B1s b; D40 ¼ A10 d B10 b l 41 A40 l 42 B40 :
1823
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1824
A3. The coefficients in Eq. (12) M 11 ¼ e110 n11 e110 n10 d þ oG1 cos ot1 oH 1 sin ot1 ; M 12 ¼ e110 n12 þ e110 n10 d; M 13 ¼ e110 n13 ; M 21 ¼ kp M 11 þ kd N 11 ; M 22 ¼ kp M 12 þ kd N 12 ; M 23 ¼ kp M 13 þ kd N 13 ; where, ¯ 1 c110 þ B ¯ 1 s110 b; ¯ 1 s110 ; n11 ¼ B ¯ 1 c110 b A n10 ¼ A ¯ 1 s110 b B ¯ 1 c110 b þ P11 c110 þ Q11 s110 Þ; n13 ¼ P12 c110 þ Q12 s110 ; n12 ¼ A ¯ 1 dÞc110 ðA ¯ 1b þ B ¯ 1d B ¯ 1b þ B ¯ 1b A ¯ 1 dÞs110 ; n21 ¼ bðA ¯ 1 bÞs110 bðA ¯ 1 dÞc110 ; n20 ¼ ðB ¯ 1 dÞs110 þ ðQ11 b P11 dÞc110 þ bðA ¯ 1b þ B ¯ 1 dÞc110 ðP11 b þ Q11 dÞs110 ; n22 ¼ bðB¯ 1 b A n23 ¼ ðQ12 b P12 dÞc110 ðP12 b þ Q12 dÞs110 ; N 11 ¼ e110 n21 e110 n20 d o2 ðG 1 sin ot1 þ H 1 cos ot1 Þ; N 12 ¼ e110 n22 þ e110 n20 d; N 13 ¼ e110 n23 :
References [1] Corradini ML, Orlando G. Robust stabilization of nonlinear uncertain plants with backlash or dead zone in the actuator. IEEE Transactions on Control Systems Technology 2000;10(1):158–66. [2] Tao G, Kokotovic PV. Adaptive control of systems with actuator and sensor nonlinearities. New York: Wiley; 1996. [3] Shaw SW, Holmes PJ. A periodically forced piecewise linear oscillator. Journal of Sound and Vibration 1983;90:129–55. [4] Thompson JM, Bokain AB, Ghaffari R. Subharmonic resonances and chaotic motions of the bilinear oscillator. IMA Journal of Applied Mathematics 1983;31:207–34. [5] Karyeaclis MP, Caughey TK. Stability of a semi-active impact damper: part II, periodic solutions. ASME Journal of Applied Mechanics 1989;56:930–40. [6] Nordmark AB. Non-periodic motion caused by grazing incidence in an impact oscillator. Journal of Sound and Vibration 1991;145:279–97. [7] Foale S, Bishop SR. Dynamical complexities of forced impacting systems. Philosophical Transactions of Royal Society of London 1992;A338:547–56. [8] Chatterjee S, Mallik AK, Ghosh A. Period response of piecewise nonlinear oscillators under harmonic excitation. Journal of Sound and Vibration 1996;191:129–44. [9] Hinrichs N, Oestreich M, Popp K. Dynamics of oscillators with impact and friction. Chaos, Solitons and Fractals 1997;8(4):535–58. [10] Blazejczyk-Okolewska B, Kapitaniak T. Coexisting attractors of impact oscillator. Chaos, Solitons and Fractals 1998;9(8):1439–43. [11] Ouyang H, Mottershed JE, Cartmell MP, Brookfield DJ. Friction induced vibration of an elastic slider on a vibrating disc. International Journal of Mechanical Sciences 1999;41:325–36. [12] Lenci S, Rega G. Periodic solutions and bifurcations in an impact inverted pendulum under impulsive excitation. Chaos, Solitons and Fractals 2000;11:2453–72.
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1825
[13] Wiercigroch M. Mathematical models of mechanical systems with discontinuities. In: Wiercigroch M, de Kraker B, editors. Applied nonlinear dynamics and chaos of mechanical systems with discontinuities. Singapore: Word Scientific; 2000. p. 17–38. [14] Natsiavas S. Dynamics of piecewise linear oscillators. In: Wiercigroch M, de Kraker B, editors. Applied nonlinear dynamics and chaos of mechanical systems with discontinuities. Singapore: Word Scientific; 2000. p. 127–53. [15] Hogan SJ. Damping in rigid block dynamics contained between side-walls. Chaos, Solitons and Fractals 2000;11:495–506. [16] Wiercigroch M. Modelling of dynamical systems with motion dependent discontinuities. Chaos, Solitons and Fractals 2000;11:2429–42. [17] Karpenko EV, Wiercigroch M, Pavlovskaia EE, Cartmell MP. Piecewise approximate analytical solutions for a Jeffcott rotor with a snubber ring. International Journal of Mechanical Sciences 2002;44:475–88. [18] Thompson JMT, Stewart HB. Nonlinear dynamics and Chaos. New York: Wiley; 1987.
International Journal of Mechanical Sciences 46 (2004) 1807–1825 www.elsevier.com/locate/ijmecsci
Periodic and chaotic motions of a harmonically forced piecewise linear system J.C. Ji, A.Y.T. Leung Department of Building and Construction, City University of Hong Kong, Kowloon, Hong Kong Received 1 July 2003; received in revised form 13 October 2004; accepted 19 October 2004
Abstract The dynamics of a harmonically excited single degree-of-freedom linear system with a feedback control, in which the actuator is subjected to dead zone and saturation constraints, is investigated in detail. The controlled system is mathematically modeled by a set of three piecewise linear equations. It is found that the system may exhibit nine types of symmetric and asymmetric period-one motions, which are characterized by a different number of crossing dead zone and saturation region per cycle. A solution for the symmetric period-one motion with a doubly crossing dead zone and saturation region is analytically constructed and its stability characteristics is examined. Other types of dynamic response such as sub-harmonic periodic motions and chaotic motions, found through numerical simulations, are also presented. r 2004 Elsevier Ltd. All rights reserved. Keywords: Periodic motions; Chaotic motions; Piecewise linear system; Dead zone; Saturation
1. Introduction Various non-linearities of the structural and control elements are frequently encountered in modeling many of controlled systems. Typical nonlinearities in such systems are the non-smooth nature of actuator devices such as dead zone, saturation, relay and hysteresis. These Corresponding author. Department of Mechanical Engineering, The University of Adelaide, Adelaide, SA 5005,
Australia. Tel.: +61 8 8303 6941; fax: +61 8 8303 4367. E-mail addresses: [email protected] (J.C. Ji), [email protected] (A.Y.T. Leung). 0020-7403/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2004.10.005
ARTICLE IN PRESS 1808
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
non-linearities, which were often neglected in designing and modeling of the control system, may deteriorate the performance of a feedback control system and lead to a poor dynamical behavior of the controlled system [1]. Dead zone, which describes the component’s insensitivity to small signals, is a static input–output relationship which for a range of input values gives no output. Once the output appears, the slope between the input and output is constant. Saturation may be caused by the limitation of the physical devices. When an input signal reaches a certain value, further increases in the input signal produce only the same output, which is usually called the saturation value. Investigation of the effect of dead zone and saturation is of great importance for the control systems designers due to the fact that these nonlinearities exist in a large number of actuators [2]. On the other hand, a thorough understanding of the dynamical behavior can be of prime importance to design the system and to provide useful information for the fault diagnosis. The main objective of the present paper is to investigate the dynamic behavior of the controlled system with dead zone and saturation non-linearities. For simplicity, the non-linearities considered in the present paper are assumed to be approximated by piecewise linear characteristics. Due to the presence of dead zone and saturation constraints, equations of motion of the controlled system appear in a strongly nonlinear form with piecewise linear characteristic, which is switched by the displacement and velocity constraints. Though the dynamical behavior of piecewise linear systems with a bilinear or tri-linear characteristic has been an active topic of intensive research e.g. [3–17], the dynamics of a piecewise linear system with quintuple linear elements has not yet been reported in the literature to the authors’ knowledge. It is found that the system considered here may exhibit symmetric and asymmetric period-one motions. Period-one motions may lose their stability and give their places to sub-harmonic or chaotic motions. In addition, sub-harmonic motions may coexist with a period-one motion, and a chaotic motion may coexist with two asymmetric period-one motions. The present paper is organized as follows; in Section 2, the equations of motion are derived. An analytical method is developed to determine period-one motions in Section 3, and their stability characteristics are examined in Section 4. Other types of motions such as sub-harmonic periodic motions and chaotic motions, are discussed in Section 5. Conclusions are presented in Section 6.
2. Equations of motion Consider a simple linear feedback control system in which a periodic excitation (disturbance) is presented, as shown in Fig. 1. The equation governing the motion of the plant (i.e., a singledegree-of-freedom linear oscillator) can be written as Mx00 þ Cx0 þ Kx ¼ F cos OT U;
Fig. 1. Block diagram of a linear system with feedback control.
(1)
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1809
where M is the mass, C is the damping coefficient, K is the stiffness, F and O are the amplitude and frequency of the excitation, U denotes the output of the actuator, and a prime denotes differentiation with respect to the time T. In the presence of dead zone and saturation constraints, the output of the actuator, as shown in Fig. 2, is of the form 8 if V XV s ; K a ðV s V d Þ; > > > > > if V d pVpV s ; > < K a ðV V d Þ; if V d pVpVd ; U ¼ 0; (2) > > > K a ðV þ V d Þ; if V s pV p V d ; > > > : K a ðV s V d Þ; if V p V s ; where K a denotes the actuator gain, V d and V s denote the maximum threshold values of dead zone and saturation. For simplicity, the feedback control is assumed to employ a PD-controller with the form V ¼ K p x þ K d x0 ; where K p and K d are the proportional and derivative gains, respectively. Substituting Eq. (2) into Eq. (1) and then introducing nondimensional variables x ¼ yV d and t ¼ O1 T; yields the equations of motion in non-dimensional form y€ þ 2dy_ þ y ¼ f cos ot;
_ for jkp y þ kd yjp1;
y€ þ 2my_ þ ð1 þ kÞy f d sgnðyÞ ¼ f cos ot; y€ þ 2dy_ þ y þ f s sgnðyÞ ¼ f cos ot;
_ for 1pjkp y þ kd yjpu s;
_ for jkp y þ kd yjXu s;
(3a) (3b) (3c)
where, O21 ¼ K=M; 2d ¼ C=MO1 ; o ¼ O=O1 ; f ¼ F =MV d O21 ; f d ¼ K a =M O21 ; f s ¼ K a ðV s V d Þ=MO21 V d ; ka ¼ K a =M O21 ; kp ¼ K p ; kd ¼ K d O1 ; k ¼ ka kp ; us ¼ ðV s V d Þ=V d ; 2m ¼ 2d þ ka kd ; sgn( ) denotes the sign function, and a superscript dot indicates differentiation with respect to the time t.
Fig. 2. Output of the actuator subjected to dead zone and saturation constraints.
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1810
It is easy to note that system (3) is continuous and satisfies the Lipschitz condition. Therefore, the solution of Eq. (3) exists and is unique. However, except for a small amplitude response, the analytical solutions of Eq. (3) cannot be obtained in closed form, since, as will be seen, the exact times when the input of the actuator crosses the boundaries of dead zone and saturation region cannot be found explicitly. These time instances will be subsequently referred to as the crossing times. In the next section, a method is developed to determine various types of period-one motions.
3. Period-one motions The response of the overall system given by Eq. (3) may exhibit nine types of period-one motions, which are categorized by the characteristics of motions and the number of crossing dead zone and saturation region per cycle, as shown in Fig. 3. The period-one response will be subsequently referred to as small amplitude motion without crossing dead zone, symmetric and asymmetric motions with a doubly crossing dead zone and saturation region, asymmetric motions with a doubly crossing dead zone and singly crossing saturation region, symmetric and asymmetric motions with a doubly crossing dead zone, respectively. 3.1. Small amplitude motion For a small amplitude motion as shown in Fig. 3a, the actuator does not give any output. The response of the oscillator is linear and the steady state solution is of the form yðtÞ ¼ f 1 cosðot yÞ;
(4)
where
2do f1 ¼ ; y ¼ arctan : 1 o2 ½ð1 o2 Þ2 þ 4d2 o2 1=2 f
The critical value of the excitation amplitude at which the actuator gives an output, can be _ which is given by f c ¼ ð½ð1 o2 Þ2 þ obtained by applying the condition jkp y þ kd yjX1; 2 2 1=2 2 2 2 1=2 4d o Þ=ð½kp þ o kd Þ: When the excitation amplitude is larger than f c ; the feedback control will apply a control output to the plant. The motion of the system may cross the boundaries of dead zero and/or saturation regions and then leave and enter dead zone and/or saturation regions at least once over a period of period-one response. 3.2. Motions with a doubly crossing dead zone and saturation region A symmetric or an asymmetric periodic motion with a doubly crossing dead zone and saturation region, as shown in Figs. 3b, c, and d, may appear when the excitation amplitude is much larger than f c : The input signal of the actuator and the time trajectory of the response (as shown in Fig. 3b) are entirely composed of eight segments according to the time intervals: ½t0 ; t1 ; ½t1 ; t2 ; ½t2 ; t3 ; ½t3 ; t4 ; ½t4 ; t5 ; ½t5 ; t6 ; ½t6 ; t7 and ½t7 ; t8 ; where t8 ¼ t0 þ 2p o and ti represents the time instant that the input signal passes the boundaries of dead zero and saturation regions. An
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1811
Fig. 3. Time histories of the system response (left subfigure) and the corresponding input signal of the actuator (right subfigure) for nine types of period-one motions: (a) a small amplitude motion without crossing dead zone; (b) a symmetric motion with a doubly crossing dead zone and saturation region (where 1; 2; 9 denote the points at which the stiffness and damping discontinuities occur); (c) and (d) asymmetric motions with a doubly crossing dead zone and saturation region; (e) and (f) asymmetric motions with a doubly crossing dead zone and singly crossing saturation region; (g) asymmetric motion with a doubly crossing dead zone; (h) and (i) asymmetric motions with a doubly crossing dead zone.
ARTICLE IN PRESS 1812
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Fig. 3. (Continued)
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1813
Fig. 3. (Continued)
analytical construction of the solution for each interval is based on the fact that Eq. (3) is linear in _ _ _ 1pjkp y þ kd yjpu each of the regions jkp y þ kd yjp1; s ; and jk p y þ kd yjXu s ; respectively. In the present paper, emphasis is placed on the analytical determination of the symmetric period-one motion with a doubly crossing dead zone and saturation region. Other types of period-one solutions can be obtained in a similar procedure. Because of symmetry of the motion examined, only four segments of the periodic solution need to be considered. _ y2 ðtÞ from t1 in Let y1 ðtÞ denote the solution starting from t0 in the region jkp y þ kd yjp1; _ _ _ 1pðkp y þ kd yÞpu s ; y3 ðtÞ from t2 in ðk p y þ k d yÞXu s ; and y4 ðtÞ from t3 in 1pðkp y þ k d yÞpu s; respectively, then the general solutions for these four regions can be written as y1 ðtÞ ¼ edðtt0 Þ ½A1 cos bðt t0 Þ þ B1 sin bðt t0 Þ þ G1 sin ot þ H 1 cos ot _ for jkp y þ kd yjp1;
ð5aÞ
y2 ðtÞ ¼ emðtt1 Þ ½A2 cos qðt t1 Þ þ B2 sin qðt t1 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 _ for 1pðkp y þ kd yÞpu s;
ð5bÞ
y3 ðtÞ ¼ edðtt2 Þ ½A3 cos bðt t2 Þ þ B3 sin bðt t2 Þ þ G1 sin ot þ H 1 cos ot þ Y 2 _ for ðkp y þ kd yÞXu s;
ð5cÞ
y4 ðtÞ ¼ emðtt3 Þ ½A4 cos qðt t3 Þ þ B4 sin qðt t3 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 _ for 1pðkp y þ kd yÞpu s;
ð5dÞ
where, Ai ; Bi ; and ti1 (i ¼ 1; 2; 3; 4) are the unknown constants to be determined, pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ¼ 1 d2 ; q ¼ 1 þ k m2 ; Y 1 ¼ f d =ð1 þ kÞ sgn ðyÞ; Y 2 ¼ f s sgnðyÞ; G1 ¼
2dof ð1 o2 Þf ; H ¼ ; 1 ð1 o2 Þ2 þ ð2doÞ2 ð1 o2 Þ2 þ ð2doÞ2
G2 ¼
2mof ð1 þ k o2 Þf ; H ¼ : 2 ð1 þ k o2 Þ2 þ ð2moÞ2 ð1 þ k o2 Þ2 þ ð2moÞ2
ARTICLE IN PRESS 1814
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
The eight constants Ai and Bi ; and four crossing times ti1 ; can be determined by applying a set of initial conditions and periodicity, continuity as well as symmetry conditions, which can be written as; at point 1;
kp y1 ðt0 Þ þ kd y_ 1 ðt0 Þ ¼ 1;
(6a)
at point 2;
kp y2 ðt1 Þ þ kd y_ 2 ðt1 Þ ¼ 1;
(6b)
y1 ðt1 Þ ¼ y2 ðt1 Þ; at point 3;
(6e)
y_ 2 ðt2 Þ ¼ y_ 3 ðt2 Þ;
(6f,g)
kp y4 ðt3 Þ þ kd y_ 4 ðt3 Þ ¼ us ;
(6h)
y_ 3 ðt3 Þ ¼ y_ 4 ðt3 Þ;
(6i,j)
y3 ðt3 Þ ¼ y4 ðt3 Þ; at point 5;
(6c,d)
kp y3 ðt2 Þ þ kd y_ 3 ðt2 Þ ¼ us ;
y2 ðt2 Þ ¼ y3 ðt2 Þ; at point 4;
y_ 1 ðt1 Þ ¼ y_ 2 ðt1 Þ;
y4 ðt4 Þ ¼ y1 ðt0 Þ;
y_ 4 ðt4 Þ ¼ y_ 1 ðt0 Þ;
and t4 t0 ¼ p=o:
(6k,l) (6m)
For brevity, Eqs. 6(a), (b), (e) and (h) are subsequently referred to as the initial conditions, Eqs. 6(c), (f), (i) and (k) as the displacement matching conditions, and equations 6(d), (g), (j) and (l) as the velocity matching conditions, respectively. The constants, Ai and Bi ; can be determined by solving the equations that result from implementing four initial conditions and four displacement matching conditions. For the clarity of notation, these constants are expressed as Ai ¼ Ai0 þ Aic cos ot0 þ Ais sin ot0 ; Bi ¼ Bi0 þ Bic cos ot0 þ Bis sin ot0 ;
ð7Þ
where the coefficients Ai0 ; Aic ; Ais ; Bi0 ; Bic and Bis are given in Appendix A1. Implementing the four velocity matching conditions gives rise to a set of four transcendental equations for three time intervals ðti1 t0 Þ and time t0 as follows: Di1 cos ot0 þ Di2 sin ot0 ¼ Di0 ;
i ¼ 1; 2; 3; 4;
(8a2d)
where the coefficients ,Di1 ; Di2 and Di0 ; are listed in Appendix A2. Eq. (8) involves the system parameters and four unknown crossing times only. No analytic solutions can be found, thus a numerical mean is resorted to solve three time intervals ðti1 t0 Þ and time t0 : Then constants Ai and Bi can be obtained from Eq. (7) by back substitutions. After getting Ai ; Bi and crossing times ti1 ; the corresponding histories of yi can be determined by Eq. (5). For an asymmetric motion as shown in Fig. 3c or Fig. 3d, a determination of its solution can be performed by the same procedure. In these cases, eight general solutions are needed to determine for eight distinct segments. The unknowns are sixteen constants Ai ; Bi and eight crossing times ti1 (where i ¼ 1; 2; . . . 8:), which can be determined by imposing a set of initial, matching and periodicity conditions. A set of 24 conditions similar to the forms given by Eq. (6) can be written
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1815
at nine points, where the damping and stiffness discontinuities occur. An implementation of these conditions may lead to a set of eight transcendental equations for the eight crossing times, which are solved using a numerical procedure. The constants Ai ; Bi can be then evaluated by simple back substitutions, and hence the corresponding solutions yi can be determined. 3.3. Asymmetric motions with a doubly crossing dead zone and singly crossing saturation region The asymmetric motions like those shown in Figs. 3e and f, consist of six distinct segments, i.e., six intervals in which the damping and stiffness properties remain constant. It is thus needed to determine six general solutions similar to the forms given by Eq. (5) for the six segments. The unknowns are twelve constants Ai ; Bi and six crossing times ti1 ; where i ¼ 1; 2; 6: Similarly, these constants can be determined by implementing a set of initial and matching conditions. Totally eighteen independent conditions, which are not given in the present paper for brevity, can be written at seven points where the discontinuities occur. Imposing these conditions leads to an equal number of transcendental equations, which can be further reduced to a set of six transcendental equations in a form similar to Eq. (8). Numerical solutions of these six equations determine the six crossing times. Then the constants, Ai and Bi ; and the corresponding solutions yi can be determined by simple back substitutions. 3.4. Motions with a doubly crossing dead zone For a symmetric motion like that shown in Fig. 3g, its solution can be constructed based on two distinct segments. Two general solutions for the two segments are needed to determine. The corresponding constants Ai ; Bi and two crossing times ti1 (where i ¼ 1; 2) can be determined by implementing a set of six conditions. For asymmetric solutions as shown in Figs. 3h and i, the motions consist of four distinct segments. It is naturally required to determine four general solutions for the four segments. The unknowns are eight constants Ai ; Bi and four crossing times ti1 ; where i ¼ 1; 2; 3; 4: Totally a set of twelve independent conditions can be written at five points. Then a similar procedure developed in Section 3.2 can be performed to determine these solutions.
4. Stability of the period-one motions For a small amplitude motion like that Discussed in Section 3.1, its stability analysis can be performed using a classical method. Evidently, the motion is stable, provided that the damping coefficient in Eq. 3(a) is positive. For other types of periodic motions, their stability characteristics can be examined by considering the asymptotic behavior of the perturbations to the steady state solutions. In this section, a stability analysis is presented for the symmetric periodic solution with a doubly crossing dead zone and saturation region. A similar procedure can be performed to investigate the stability of other periodic solutions. The symmetry of the motion examined here ensures that its stability can be fully determined in one half of the solution.
ARTICLE IN PRESS 1816
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Let zi ðtÞ denote the corresponding perturbed solutions to yi ðtÞ (where i ¼ 1; 2; 3; 4), respectively, then the perturbed solutions can be written as
z1 ðtÞ ¼ edðtt0 Dt0 Þ P1 cos bðt t0 Dt0 Þ þ Q1 sin bðt t0 Dt0 Þ þ G1 sin ot þ H 1 cos ot; for jkp z þ kd z_jp1; ð9aÞ
z2 ðtÞ ¼ emðtt1 Dt1 Þ P2 cos qðt t1 Dt1 Þ þ Q2 sin qðt t1 Dt1 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 ; for 1pðkp z þ kd z_Þpus ;
ð9bÞ
z3 ðtÞ ¼ edðtt2 Dt2 Þ P3 cos bðt t2 Dt2 Þ þ Q3 sin bðt t2 Dt2 Þ þ G1 sin ot þ H 1 cos ot þ Y 2 ; for ðkp z þ kd z_ÞXus ;
ð9cÞ
z4 ðtÞ ¼ emðtt3 Dt3 Þ P4 cos qðt t3 Dt3 Þ þ Q4 sin qðt t3 Dt3 Þ þ G2 sin ot þ H 2 cos ot þ Y 1 ; for 1pðkp z þ kd z_Þpus ;
ð9dÞ
where Pi and Qi ; (i ¼ 1; 2; 3; 4), are the constants to be determined. The initial conditions for the perturbed motions of the first segment can be written as z1 ðt0 þ Dt0 Þ ¼ z0 þ Dz0 ; kp z1 ðt0 þ Dt0 Þ þ kd z_1 ðt0 þ Dt0 Þ ¼ 1; where operator D indicates a small perturbation of the operand. At moment ðt1 þ Dt1 Þ; the input signal of the actuator reaches a boundary of dead zone and enters the region 1pðkp z þ kd z_Þpus thereafter. The perturbed response at the moment of departure from dead zone is assumed to be z1 ðt1 þ Dt1 Þ ¼ z1 þ Dz1 ; kp z1 ðt1 þ Dt1 Þ þ kd z_1 ðt1 þ Dt1 Þ ¼ 1:
(10)
Since the perturbations in the initial conditions are assumed to be small, it is expected that the ¯ 1 and coefficients P1 and Q1 will assume the values close to those of the unperturbed solutions, A ¯ 1 ; respectively. To the first order, they can be expressed as B ¯ 1 þ P11 Dt0 þ P12 Dz0 ; P1 ¼ A ¯ 1 þ Q11 Dt0 þ Q12 Dz0 ; Q1 ¼ B
ð11Þ
¯ 1 ¼ z0 G 1 sin ot0 H 1 cos ot0 ;P11 ¼ oH 1 sin ot0 oG 1 cos ot0 ; P12 ¼ 1; where, A ¯ 1 ¼ 1 þ kd b kp z0 dH 1 þ oG 1 cos ot0 þ oH 1 dG 1 sin ot0 ; B kd b kd b b b Q11 ¼ ½oðoG 1 þ dH 1 Þ sin ot0 þ oðoH 1 dG1 Þ cos ot0 =b; Q12 ¼ ðkd d kp Þ=kd b: Substituting Eq. (9a) into Eq. (10), and keeping only the first-order terms yields Dz1 ¼ M 11 Dt1 þ M 12 Dt0 þ M 13 Dz0 ; M 21 Dt1 þ M 22 Dt0 þ M 23 Dz0 ¼ 0; where the coefficients M ij are defined in Appendix A3. For brevity, Eq. (12) can be put in a matrix form as Dt1 Dt0 ¼ R1 ; Dz1 Dz0
ð12Þ
(13)
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1817
where R1 is a 2 2 matrix with the elements r11 ¼ M 22 =M 21 ; r12 ¼ M 23 =M 21 ; r21 ¼ M 12 M 11 M 22 =M 21 ; r22 ¼ M 13 M 11 M 23 =M 21 ; respectively. The asymptotic behavior of the perturbed solutions for other segments from ðt1 þ Dt1 Þ to ðt2 þ Dt2 Þ (for the second segment), from ðt2 þ Dt2 Þ to ðt3 þ Dt3 Þ (for the third segment), and from ðt3 þ Dt3 Þ to ðt4 þ Dt4 Þ (for the fourth segment), can be investigated in a similar procedure. After a lengthy algebraic manipulation, the following equations retaining the first order terms can be obtained Dti1 Dti ¼ Ri ; i ¼ 2; 3; 4; (14) Dzi Dzi1 where Ri are three 2 2 matrices, whose elements are not given here for brevity. The perturbation to the steady state solution during the first half period is then obtained by combining Eq. (13) and (14) to form an equation as Dt0 Dt4 ¼J ; (15) Dz4 Dz0 where J ¼ R1 R2 R3 R4 ; represents the transition matrix of the response from ðt0 þ Dt0 Þ to ðt4 þ Dt4 Þ: The eigenvalues of matrix J determines the stability of the solution. The symmetric periodone solution is asymptotically stable if both eigenvalues of J have modulus less than unity. When either of the two eigenvalues has modulus greater than one the solution is unstable. For some combinations of the system parameters, the modulus may take the value of unity, where a bifurcation occurs. To validate the present analysis, the analytical solutions of the symmetric period-one motion with a doubly crossing dead zone and saturation region were compared with those obtained from direct numerical integration of Eq. (3) for several cases. It was found that the solutions obtained by these two different means were identical. As an illustrative example, consider a system with parameters in Eq. (3) as follows; d ¼ 0:1; m ¼ 0:1; k ¼ 4:0; kp ¼ 4:0; us ¼ 4:0; o ¼ 2:0; and f ¼ 1:3: The crossing times obtained via the analysis are t0 ffi 0:405825; t1 ffi 0:535464; t2 ffi 0:740362; and t3 ffi 1:771968: The time intervals between two crossing times given by numerical integration of Eq. (3) were nearly identical to those predicted by the analysis. A very small difference only in the sixth decimal place was found between the values obtained by two approaches. Fig. 4 shows the phase portrait of the symmetric solutions obtained by the analysis and numerical integration of Eq. (3). The analytical solution is in excellent agreement with the result of numerical integration.
5. Other types of motions The analysis presented in Sections 3 and 4 can be extended to construct sub-harmonic periodic solutions. However, the procedure for finding such solutions is exceptionally lengthy and involves much more algebraic manipulation. In this section, results of numerical integration are given to show other behavior of the system. The widely used software package MATLAB/Simulink is employed to model system (3). The dead zone block and the saturation block in the library of the built-in Simulink nonlinear
ARTICLE IN PRESS 1818
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Fig. 4. Phase portrait of the symmetric motion with a doubly crossing dead zone and saturation region; solid curve indicates results of numerical integration and circles represent analytical solutions.
components are used to impose the starting and end values of dead zone and the upper and lower bounds of saturation on the input signal, respectively. The fourth-order Runge–Kutta procedure is chosen to be the integration method of the ODE solvers. The initial conditions for the integrator blocks are all set zero, unless otherwise stated. The simulation results are observed using scopes and other display blocks during the simulation. Extensive numerical simulations have been performed in order to explore the complex dynamical behavior. For the results given in this paper, the parameters are chosen as follows: d ¼ 0:025; ka ¼ 1:2; kp ¼ 3:0; kd ¼ 0:05; o ¼ 0:9; us ¼ 4:0; unless otherwise specified. In addition to period-one motions, a variety of sub-harmonic and chaotic motions may be observed. The minimum period of a sub-harmonic motion is an integral multiple of the forcing period. This integer indicates the order of the sub-harmonic motion. Sub-harmonic periodic motions of order two, three, four, five, seven and eight were found to exist for some combinations of the system parameters. Sub-harmonic motions may appear not only in the regions where the input of the actuator does not reach the boundaries of saturation regions, but also in the regions where the input signal exceeds the saturation value. Sub-harmonic motions seemed to disappear and period-one motions appear at higher damping levels. Chaotic motions, which result from either a period-one motion or a subharmonic motion losing its stability, are another widely observed complex dynamical behavior [18]. By taking the excitation amplitude as the bifurcation parameter, the bifurcation diagrams of the response are shown in Fig. 5. The system parameters and mapping conditions for sampling data are identical for two diagrams except initial conditions. The data for plotting Fig. 5a were obtained via a usual way to set initial conditions, while for Fig. 5b the initial conditions ðy0 ; y_ 0 Þ ¼ ð1:0; 0:6Þ remain unchanged at each step of the excitation amplitude. It can be easily seen in Fig. 5 that chaotic motions exist at many values of the excitation amplitude. Fig. 6 displays a typical chaotic motion at f ¼ 5:4; where (a) is the Poincare map at the excitation phase p4 and (b) is the phase portrait. Another interesting dynamical behavior, which widely exists in the system response, is a coexistence of different types of motions starting from different sets of initial conditions for a given set of system parameters. The particular dependence of motion on the initial conditions can
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1819
Fig. 5. Bifurcation diagrams of the response; (a) the usual way to set initial conditions, (b) the initial conditions at each step are ðy0 ; y_ 0 Þ ¼ ð1:0; 0:6Þ:
Fig. 6. The chaotic motion of the system at f ¼ 5:4 under initial conditions ðy0 ; y_ 0 Þ ¼ ð0:0; 1:0Þ; (a) Poincare map, (b) phase portrait.
be also observed in Fig. 5 for many values of the excitation amplitude, for example, when f ¼ 3:25; 4:15 4:7; 6:0: It has been numerically found that a period-one motion may coexist with one or two sub-harmonic motions, or that two asymmetric period-one motions coexist with a chaotic motion. For instance, two sub-harmonic motions of order two may coexist with a periodone motion. A sub-harmonic motion of order three or five was found to coexist with a period-one motion. Due to laborious calculation being involved, a determination of the domain of attraction
ARTICLE IN PRESS 1820
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
Fig. 7. Phase portraits of the response for the co-existence of a period-one motion and a subharmoinc motion of order seven as well as of order five; (a) period-one motion starting from ðy0 ; y_ 0 Þ ¼ ð0:0; 1:0Þ; (b) subharmonic motion of order seven from (1.0,0.6), and (c) subharmonic motion of order five from ð0:6; 0:3Þ:
of various solutions is not considered in the present paper. Alternatively, some examples are given to illustrate the coexisting solutions. A sub-harmonic motion of order five and a subharmonic motion of order seven may coexist with a period-one motion when f ¼ 6:0; as shown in Fig. 7. Fig. 8 illustrates the coexistence of two asymmetric period-one motions and a chaotic motion for f ¼ 4:55: Each of the asymmetric motions is the other’s symmetric image about the origin. Starting from initial conditions ðy0 ; y_ 0 Þ ¼ ð0:0; 1:0Þ leads to the chaotic motion, while from ð0:6; 0:3Þ and ð1:0; 0:6Þ result in two asymmetric motions, respectively. The presence of coexistence of multiple stable motions and of chaotic motions causes a difficulty in predicting the system dynamics.
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1821
Fig. 8. Coexistence of a pair of asymmetric period-one motions and a chaotic motion; (a) phase portraits of two asymmetric period-one motions starting from ðy0 ; y_ 0 Þ ¼ ð0:6; 0:3Þ for the solid curve and ð1:0; 0:6Þ for the dashed curve, (b) Poincare map of the chaotic motion starting from ð0:0; 1:0Þ:
6. Conclusions The dynamics of a periodically forced single degree-of-freedom piecewise linear system with quintuple linear segments has been investigated. A detailed procedure has been performed to construct a solution for the symmetric period-one motion with a doubly entering dead zone and saturation region per cycle. The stability of the solution has been studied by examining the asymptotic behavior of the corresponding perturbed solution. It was found that the response of the system may exhibit symmetric and asymmetric period-one motions, sub-harmonic periodic motions, chaotic motions, and coexistence of three stable motions.
Appendix: List of the coefficients A1. The coefficients in Eq. (7) A10 ¼ ða81 h60 þ a82 h70 þ b80 Þ=ad; A1c ¼ ða81 h61 þ a82 h71 þ b81 Þ=ad; A1s ¼ ða81 h62 þ a82 h72 þ b82 Þ=ad; A20 ¼ h20 þ g2 A10 ; A2c ¼ h21 þ g2 A1c ; A2s ¼ h22 þ g2 A1s ; A30 ¼ h40 þ g4 A10 ; A3c ¼ h41 þ g4 A1c ; A3s ¼ h42 þ g4 A1s ; A40 ¼ h60 þ g6 A10 ; A4c ¼ h61 þ g6 A1c ; A4s ¼ h62 þ g6 A1s ; Bk0 ¼ hj0 þ gj A10 ; Bkc ¼ hj1 þ gj A1c ; Bks ¼ hj2 þ gj A1s ; k ¼ 1; 2; 3; 4; j ¼ 2k 1; where, e110 ¼ edðt1 t0 Þ ; e120 ¼ edðt2 t0 Þ ; e130 ¼ edðt3 t0 Þ ; e210 ¼ emðt1 t0 Þ ; e220 ¼ emðt2 t0 Þ ;
ARTICLE IN PRESS 1822
J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825 p=o
e230 ¼ emðt3 t0 Þ ; e20 ¼ em ; s2 ¼ sin q
p p ; c2 ¼ cos q ; s10 ¼ sin oðt1 t0 Þ; o o
c10 ¼ cos oðt1 t0 Þ; s20 ¼ sin oðt2 t0 Þ; c20 ¼ cos oðt2 t0 Þ; s30 ¼ sin oðt3 t0 Þ; c30 ¼ cos oðt3 t0 Þ; s110 ¼ sin bðt1 t0 Þ; c110 ¼ cos bðt1 t0 Þ; s120 ¼ sin bðt2 t0 Þ; c120 ¼ cos bðt2 t0 Þ; s130 ¼ sin bðt3 t0 Þ; c130 ¼ cos bðt3 t0 Þ; s210 ¼ sin qðt1 t0 Þ; c210 ¼ cos qðt1 t0 Þ; s220 ¼ sin qðt2 t0 Þ; c220 ¼ cos qðt2 t0 Þ; s230 ¼ sin qðt3 t0 Þ; c230 ¼ cos qðt3 t0 Þ; a11 ¼ kp kd d; a12 ¼ kd b; b11 ¼ kp H 1 þ okd G1 ; b12 ¼ kp G 1 okd H 1 ; b10 ¼ 1; a21 ¼ kp kd m; a22 ¼ qkd ; b21 ¼ kp H 2 þ okd G 2 ; b22 ¼ kp G 2 okd H 2 ; b20 ¼ 1 kp Y 1 ; a31 ¼ e110 c110 ; a32 ¼ e110 s110 ; b31 ¼ ðG1 G2 Þs10 þ ðH 1 H 2 Þc10 ; b32 ¼ ðG1 G2 Þc10 ðH 1 H 2 Þs10 ; b30 ¼ Y 1 ; a41 ¼ kp kd d; a42 ¼ kd b; b41 ¼ ðkp G 1 okd H 1 Þs20 þ ðkp H 1 þ okd G 1 Þc20 ; b42 ¼ ðkp G1 okd H 1 Þc20 ðkp H 1 þ okd G 1 Þs20 ; b40 ¼ us kp Y 2 ; a51 ¼ e220 ðc220 c210 þ s220 s210 Þ=e210 ; a52 ¼ e220 ðs220 c210 c220 s210 Þ=e210 ; b51 ¼ ðG2 G1 Þs20 þ ðH 2 H 1 Þc20 ; b52 ¼ ðG2 G1 Þc20 ðH 2 H 1 Þs20 ; b50 ¼ Y 2 Y 1 ; a61 ¼ kp mkd ; a62 ¼ qkd ; b61 ¼ ðkp G 2 okd H 2 Þs30 þ ðkp H 2 þ okd G 2 Þc30 ; b62 ¼ ðkp G2 okd H 2 Þc30 ðkp H 2 þ okd G 2 Þs30 ; b60 ¼ us kp Y 1 ; a71 ¼ e130 ðc130 c120 þ s130 s120 Þ=e120 ; a72 ¼ e130 ðs130 c120 c130 s120 Þ=e120 ; b71 ¼ ðG1 G2 Þs30 þ ðH 1 H 2 Þc30 ; b72 ¼ ðG1 G2 Þc30 ðH 1 H 2 Þs30 ; b70 ¼ Y 1 Y 2 ; a81 ¼ e20 ðc2 c230 þ s2 s230 Þ=e230 ; a82 ¼ e20 ðs2 c230 c2 s230 Þ=e230 ; b81 ¼ H 1 H 2 ; b82 ¼ G 1 G 2 ; b80 ¼ Y 1 ; l 11 ¼ e110 ðc110 d þ s110 bÞ; l 12 ¼ e110 ðc110 b s110 dÞ; m11 ¼ oðG 1 G2 Þc10 oðH 1 H 2 Þs10 ; m12 ¼ oðG2 G 1 Þs10 þ oðH 2 H 1 Þc10 ; l 21 ¼ e220 ½mðc220 c210 þ s220 s210 Þ þ qðs220 c210 c220 s210 Þ =e210 ; l 22 ¼ e220 ½qðc220 c210 þ s220 s210 Þ mðs220 c210 c220 s210 Þ =e210 ; m21 ¼ oðG1 G 2 Þc20 þ oðH 1 H 2 Þs20 ; m22 ¼ oðG 1 G 2 Þs20 þ oðH 1 H 2 Þc20 ; l 31 ¼ e130 ½dðc130 c120 þ s130 s120 Þ þ bðs130 c120 c130 s120 Þ =e120 ; l 32 ¼ e130 ½bðc130 c120 þ s130 s120 Þ dðs130 c120 c130 s120 Þ =e120 ;
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
m31 ¼ oðG 1 G2 Þc30 þ oðH 2 H 1 Þs30 ; m32 ¼ oðG2 G 1 Þs30 þ oðH 2 H 1 Þc30 ; l 41 ¼ e20 ½mðc2 c230 þ s2 s230 Þ þ qðs2 c230 c2 s230 Þ =e230 ; l 42 ¼ e20 ½qðc2 c230 þ s2 s230 Þ mðs2 c230 c2 s230 Þ =e230 ; m41 ¼ oðG 1 G2 Þ; m42 ¼ oðH 2 H 1 Þ; g1 ¼ a11 =a12 ; h10 ¼ b10 =a12 ; h11 ¼ b11 =a12 ; h12 ¼ b12 =a12 ; g2 ¼ a31 þ a32 g1 ; h20 ¼ a32 h10 b30 ; h21 ¼ a32 h11 þ b31 ; h22 ¼ a32 h12 þ b32 ; g3 ¼ a21 g2 =a22 ; h30 ¼ ðb20 a21 h20 Þ=a22 ; h31 ¼ ðb21 þ a21 h21 Þ=a22 ; h32 ¼ ðb22 þ a21 h22 Þ=a22 ; g4 ¼ a51 g2 þ a52 g3 ; h40 ¼ a51 h20 þ a52 h30 b50 ; h41 ¼ a51 h21 þ a52 h31 þ b51 ; h42 ¼ a51 h22 þ a52 h32 þ b52 ; g5 ¼ a41 g4 =a42 ; h50 ¼ ðb40 a41 h40 Þ=a42 ; h51 ¼ ðb41 þ a41 h41 Þ=a42 ; h52 ¼ ðb42 þ a41 h42 Þ=a42 ; g6 ¼ a71 g4 þ a72 g5 ; h60 ¼ a71 h40 þ a72 h50 b70 ; h61 ¼ a71 h41 þ a72 h51 þ b71 ; h62 ¼ a71 h42 þ a72 h52 þ b72 ; g7 ¼ a61 g6 =a62 ; h70 ¼ ðb60 a61 h60 Þ=a62 ; h71 ¼ ðb61 þ a61 h61 Þ=a62 ; h72 ¼ ðb62 þ a61 h62 Þ=a62 ; ad ¼ 1 þ a81 g6 þ a82 g7 :
A2. The coefficients in Eq. (8) D11 ¼ l 11 A1c þ l 12 B1c þ m11 qB2c þ mA2c ; D12 ¼ l 11 A1s þ l 12 B1s þ m12 qB2s þ mA2s ; D10 ¼ qB20 mA20 l 11 A10 l 12 B10 ; D21 ¼ l 21 A2c þ l 22 B2c þ m21 B3c b þ A3c d; D22 ¼ l 21 A2s þ l 22 B2s þ m22 B3s b þ A3s d; D20 ¼ B30 b A30 d l 21 A20 l 22 B20 ; D31 ¼ l 31 A3c þ l 32 B3c þ m31 qB4c þ mA4c ; D32 ¼ l 31 A3s þ l 32 B3s þ m32 qB4s þ mA4s ; D30 ¼ qB40 mA40 l 31 A30 l 32 B30 ; D41 ¼ l 41 A4c þ l 42 B4c þ m41 A1c d þ B1c b; D42 ¼ l 41 A4s þ l 42 B4s þ m42 A1s d þ B1s b; D40 ¼ A10 d B10 b l 41 A40 l 42 B40 :
1823
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1824
A3. The coefficients in Eq. (12) M 11 ¼ e110 n11 e110 n10 d þ oG1 cos ot1 oH 1 sin ot1 ; M 12 ¼ e110 n12 þ e110 n10 d; M 13 ¼ e110 n13 ; M 21 ¼ kp M 11 þ kd N 11 ; M 22 ¼ kp M 12 þ kd N 12 ; M 23 ¼ kp M 13 þ kd N 13 ; where, ¯ 1 c110 þ B ¯ 1 s110 b; ¯ 1 s110 ; n11 ¼ B ¯ 1 c110 b A n10 ¼ A ¯ 1 s110 b B ¯ 1 c110 b þ P11 c110 þ Q11 s110 Þ; n13 ¼ P12 c110 þ Q12 s110 ; n12 ¼ A ¯ 1 dÞc110 ðA ¯ 1b þ B ¯ 1d B ¯ 1b þ B ¯ 1b A ¯ 1 dÞs110 ; n21 ¼ bðA ¯ 1 bÞs110 bðA ¯ 1 dÞc110 ; n20 ¼ ðB ¯ 1 dÞs110 þ ðQ11 b P11 dÞc110 þ bðA ¯ 1b þ B ¯ 1 dÞc110 ðP11 b þ Q11 dÞs110 ; n22 ¼ bðB¯ 1 b A n23 ¼ ðQ12 b P12 dÞc110 ðP12 b þ Q12 dÞs110 ; N 11 ¼ e110 n21 e110 n20 d o2 ðG 1 sin ot1 þ H 1 cos ot1 Þ; N 12 ¼ e110 n22 þ e110 n20 d; N 13 ¼ e110 n23 :
References [1] Corradini ML, Orlando G. Robust stabilization of nonlinear uncertain plants with backlash or dead zone in the actuator. IEEE Transactions on Control Systems Technology 2000;10(1):158–66. [2] Tao G, Kokotovic PV. Adaptive control of systems with actuator and sensor nonlinearities. New York: Wiley; 1996. [3] Shaw SW, Holmes PJ. A periodically forced piecewise linear oscillator. Journal of Sound and Vibration 1983;90:129–55. [4] Thompson JM, Bokain AB, Ghaffari R. Subharmonic resonances and chaotic motions of the bilinear oscillator. IMA Journal of Applied Mathematics 1983;31:207–34. [5] Karyeaclis MP, Caughey TK. Stability of a semi-active impact damper: part II, periodic solutions. ASME Journal of Applied Mechanics 1989;56:930–40. [6] Nordmark AB. Non-periodic motion caused by grazing incidence in an impact oscillator. Journal of Sound and Vibration 1991;145:279–97. [7] Foale S, Bishop SR. Dynamical complexities of forced impacting systems. Philosophical Transactions of Royal Society of London 1992;A338:547–56. [8] Chatterjee S, Mallik AK, Ghosh A. Period response of piecewise nonlinear oscillators under harmonic excitation. Journal of Sound and Vibration 1996;191:129–44. [9] Hinrichs N, Oestreich M, Popp K. Dynamics of oscillators with impact and friction. Chaos, Solitons and Fractals 1997;8(4):535–58. [10] Blazejczyk-Okolewska B, Kapitaniak T. Coexisting attractors of impact oscillator. Chaos, Solitons and Fractals 1998;9(8):1439–43. [11] Ouyang H, Mottershed JE, Cartmell MP, Brookfield DJ. Friction induced vibration of an elastic slider on a vibrating disc. International Journal of Mechanical Sciences 1999;41:325–36. [12] Lenci S, Rega G. Periodic solutions and bifurcations in an impact inverted pendulum under impulsive excitation. Chaos, Solitons and Fractals 2000;11:2453–72.
ARTICLE IN PRESS J.C. Ji, A.Y.T. Leung / International Journal of Mechanical Sciences 46 (2004) 1807–1825
1825
[13] Wiercigroch M. Mathematical models of mechanical systems with discontinuities. In: Wiercigroch M, de Kraker B, editors. Applied nonlinear dynamics and chaos of mechanical systems with discontinuities. Singapore: Word Scientific; 2000. p. 17–38. [14] Natsiavas S. Dynamics of piecewise linear oscillators. In: Wiercigroch M, de Kraker B, editors. Applied nonlinear dynamics and chaos of mechanical systems with discontinuities. Singapore: Word Scientific; 2000. p. 127–53. [15] Hogan SJ. Damping in rigid block dynamics contained between side-walls. Chaos, Solitons and Fractals 2000;11:495–506. [16] Wiercigroch M. Modelling of dynamical systems with motion dependent discontinuities. Chaos, Solitons and Fractals 2000;11:2429–42. [17] Karpenko EV, Wiercigroch M, Pavlovskaia EE, Cartmell MP. Piecewise approximate analytical solutions for a Jeffcott rotor with a snubber ring. International Journal of Mechanical Sciences 2002;44:475–88. [18] Thompson JMT, Stewart HB. Nonlinear dynamics and Chaos. New York: Wiley; 1987.