Stability Chart

Available online at www.sciencedirect.com

International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

Stability chart of parametric vibrating systems using energy-rate method G. Nakhaie Jazar∗ Department of Mechanical Engineering and Applied Mechanics, North Dakota State University, Fargo, ND 58105-5285, USA Received 5 March 2003; accepted 6 August 2003

Abstract Based on the integral of energy and numerical integration, we introduce, develop, and apply a general algorithm to predict parameters of a parametric equation to produce a periodic response. Using the new method, called energy-rate, we are able to 3nd not only stability chart of a parametric equation which indicates the boundaries of stable and unstable regions, but also periodic responses that are embedded in stable or unstable regions. There are three main important advantages in energy-rate method. It can be applied not only to linear but also to non-linear parametric equations; most of the perturbation methods cannot. It can be applied to large values of parameters; most of the perturbation methods cannot. Depending on the accuracy of numerical integration method, it can also 3nd the value of parameters for a periodic response more accurate than classical methods, no matter if the periodic response is on the boundary of stability and instability or it is a periodic response within the stable or unstable region. In order to introduce the energy-rate method and indicate its advantages we apply the method to the standard Mathieu’s equation, x7 + ax − 2bx cos(2t) = 0 and show how to 3nd its stability chart for the large values of b in a–b plane. The results are compared with McLachlan’s report (Theory and Application of Mathieu Function, Clarendon, Oxford, 1947). ? 2003 Elsevier Ltd. All rights reserved. Keywords: Mathieu stability chart; Energy-rate method; Parametric vibrations

1. Introduction In this paper, an applied algorithm based on integral energy and numerical integration will be introduced, developed, and presented to determine the parameters of parametric equations corresponding to periodic solutions.

∗

Tel.: +1-701-231-8303; fax: +1-701-231-8913. E-mail address: [email protected] (G.N. Jazar).

0020-7462/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2003.08.009

The algorithm is based on a method called energy-rate method that can be applied to the following general 3rst-order diEerential equation: x7 + f(x) + g(x; x; ˙ t) = 0:

(1)

There are some important advantages in energy-rate method comparing to the perturbation methods and Floquet theory. First, it can be applied to non-linear parametric equations as well as linear equations. Floquet theory, similar to the most of perturbation methods cannot be applied on non-linear equations.

1320

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

Second, it can be applied not even to the small values of parameters, but also to large values. Most of the approximation and perturbation methods are stick to the very small values of parameters. Third, depending on the accuracy of numerical integration method, it can 3nd the value of parameters for a periodic response faster and more accurate than other classical methods. The accuracy of most classical methods depend on smallness of the value of parameters, and the number of terms in perturbed solution. Their accuracy cannot sometimes be increased by increasing the degree of polynomial or the number of terms in perturbed solution. Fourth, the algorithm of detection periodic solution using energy-rate method can be adjusted to detect periodic response within a stable or unstable region, as well as periodic solutions on the boundary of stability and instability regions. I study the well-known Mathieu’s equation as an example to show the applicability of the energy-rate method in identifying the stability chart, and some of the advantages of method. The results are interesting because they uncover some hidden facts about the Mathieu’s equation. The Mathieu’s equation x7 + ax − 2bx cos(2t) = 0;

x˙ =

dx dt

(2)

is the simplest and the most widely known parametric diEerential equation, in which a and b are constant parameters. One of the most interesting characteristics of this linear equation is that depending on the coeIcients, it may have bounded or unbounded solutions. Fig. 1 depicts the stability chart of Mathieu’s equation (2) based on McLachlan’s [1] results. Parametric equations govern problems of the greatest diversity in astronomy and theoretical physics and stability of the oscillatory processes in non-linear systems. They have accordingly been the subject of a vast number of investigations since the beginning of the last century [2,3]. Parametric equation arises in applied mathematics in three main groups of problems. The 3rst group is the transverse vibrations of a taut elastic member. The second group of problems arises typically in the modulation of radio carrier wave [4]. The third group of problems come from the standard procedure of investigating the stability of a periodic motion in a

Fig. 1. Mathieu stability chart based on the numerical values, generated by McLachlan [1].

non-linear system. A linearization in neighbourhood of a periodic motion results in a linear diEerential equation with periodic coeIcients, i.e. Floquet theory. The paradigm example is given by Mathieu’s equation [5]. The Mathieu’s equation is commonly known as the equation of the elliptic cylinder functions. The contribution to the theory of diEerential equations, related to the third group, which has been stimulated by problems in celestial mechanics, has come chieLy from Hill and Poincare [6]. Examples of vibrating elastic bodies are: the lateral vibration of stretched strings and thin rods, which are perhaps the most amenable to theoretical and experimental treatment [7]. The linear and various modes of vibration of bars and beams, which was originated back by Daniel Bernoulli, was investigated completely by Rayleigh [8] and Love [9]. It is well known that classical linearized analysis of the vibrating strings and rods could lead to some results related to the Mathieu’s equation, which are reasonably accurate if only the tension and displacements are assumed to be small [10]. 2. Historical background and literature review The 3rst recorded demonstration of parametric behaviour belongs to Faraday [11] in 1831, in which he produced wave motion in Luids by vibrating a plate in contact with the Luid [12]. Shortly afterwards, many scientists discovered parametric phenomena. Among

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

them are, Melde (1860), Mathieu (1860), Matthiessen (1870), Floquet (1883), Hill (1886), Rayleigh (1887), Stephenson (1906), Meissner (1918), Strutt (1927), van der Pol (1928), and Barrow (1923), in the late 19th and the beginning of the 20th century [13]. It is known that the 3rst detailed theory relevant to the study of periodically time varying systems was given by Mathieu [14], Hill [15], and Lord Rayleigh [16] in 1860 –1890. One of the most important earlier works on the behaviour of periodically time-varying systems was that by Hill in 1886, in which he laid down the very mathematical foundations of the stability theory of parametric systems [15]. The 3rst analysis of parametric resonance of a structural con3guration (a pinned perfect column) was presented by Beliaev [17]. Many studies started appearing in the literature in the late 1940s and 1950s. Reviews on the subject, including historical sketches, may be found in the works of Beilin and Dzhanelidze [18]. For more recent contributions on the subject, the reader is referred to the work done by Yao [19], Hsu [42], Esmailzadeh and Nakhaie Jazar [20,21], Luo [22], and Zounes and Rand [23]. Solution of the Mathieu’s equation depends on two independent parameters, a and b. It is known from Floquet theory that its solution is, in general, of the following structure: x(t) = C1 e t ’(t) + C2 e− t ’(−t)

(3)

in which the function ’(t) is periodic. The characteristic exponent, , is a time invariant constant that depends in an intricate way upon the parameters a and b. If it is real, the solution becomes exponentially in3nite, i.e., a so called unstable solution. If the exponent is purely imaginary the solutions remain bounded along the real axis. The intermediate case in which

= 0 is of especial importance because it includes a solution known as a Mathieu function which is periodic [24]. Broadly speaking, the determination of the Mathieu functions is the matter of prime importance in the applications of the equation. The most comprehensive treatment of classical methods for analyzing the Mathieu’s equation has been given by McLachlan [1]. It is known that the Mathieu’s equation (2), could have periodic solution depending on parameters a and b, and independent of the initial conditions, due to the linearity of the equation. The relationship between a and b for periodic

1321

solution generates a graph which is called a stability chart [4]. The relation could be developed analytically by Fourier coeIcients or perturbation methods, but all of these methods fail to predict the stability chart when b is large. Many authors investigated the stability chart of the time-varying systems including Mathieu’s equation using a variety of concepts such as Lyapunov exponents [25], Poincare Mapping [26], Lyapunov–Floquet transformation, and perturbative Hamiltonian normal forms [27]. None of these analysis could give accurate prediction of the stability boundaries when b is large and/or some non-linearity is presented. Due to the nature of the solutions given by the classical approaches, extensive calculation is generally required to ensure suIcient convergence to give accurate answers. Determination of stability, classically, rests upon the use of the Hill determinant procedure [28]. For values of a and b signi3cantly larger than unity, quite large determinant need to be calculated which is also a time consuming task. Further, these may diverge initially prior to acceptable convergence [13]. Taylor and Narendra [29] and later on Gunderson et al. [30] found the stability boundaries of the damped Mathieu’s equation using Laplace transformation and Lyapunov function. Their stability theory was only correct for small multiplier of both damping and periodic terms. Within the limits of numerical accuracy, most of these procedures have the appeal that they give exact results since they do not rely upon algebraic approximations to the solutions of the equation. Their shortcomings, however, lie with their numerical error, especially when solutions near stability boundaries are required. Accuracy is also governed by the number of iteration points chosen for iterative techniques per period of the periodic coeIcient. Numerical solution of the Mathieu’s equation has been extensively used in the past to model systems which contain a sinusoidal time-varying parameter [31]. Rand and Hastings [32] used numerical integration to determine regions of stability for a quasi-periodic Mathieu’s equation. They assumed that the equation is stable for a given initial condition r(0), if r(t) ¡ 106 r(0) for 1000 time
units and for all t between 0 and 1000, where r(t) = x2 (t) + x˙2 (t). Using a numerical integration, they could develop and present some useful stability charts.

1322

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

3. The investigation method Stability analysis for even a single, second order, linear diEerential equation with periodic coeIcients, such as Mathieu’s equation, is rather cumbersome, but diEerent methods are available. Among them are the method of in3nite determinants [33], the perturbation methods [34], the Galerkin method [35], and the classic Floquet method [36]. The method of continued fractions is the most strong method that McLachlan used to 3nd the stability chart of Mathieu’s equation. This method starts by knowing that the periodic solution of (2) which admit the period 2 falls into four classes. Each class depends on its Fourier series which involve, respectively, cosines or sines of even or odd multiples of t. They are de3ned by Ince [37] as ce2n =

∞ 

(2n) C2r

r=0 ∞ 

se2n+1 =

r=0

cos(2rx)

have period ;

(4)

(2n+1) S2r+1 sin((2r + 1)x)

have period 2;

(5)

[1,37,38],  aC0 + bC2 = 0; r=0    (a − 4)C2 + b(2C0 + C4 ) = 0; r=1 (8)    r ¿ 2 (a − 4r 2 )C2r + b(2C2r−2 + C2r+2 ) = 0;    (a − 1 + b)C1 + bC3 = 0; r=0  (a − (2r + 1)2 )C2r+1 r¿1    +b(C2r−1 + C2r+3 ) = 0;

(9)

   (a − 1 − b)S1 + bS3 = 0; r=0  (a − (2r + 1)2 )S2r+1 r¿1    +b(S2r−1 + S2r+3 ) = 0;

(10)

   (a − 4)S2 + bS4 = 0; r=0  (a − (2r + 2)2 )S2r+2 r¿1    +b(S2r + S2r+4 ) = 0:

(11)

Using the recurrence relations (8)–(11) we can 3nd a continued fraction equation between a and b such as 2

2

ce2n+1 =

∞  r=0

(2n+1) C2r+1 cos((2r + 1)x)

∞  r=0

2

b b − b4 1 1 16 36 1 − 4a 4 1 − 16a 16 1 − 36a 2

have period 2; se2n+2 =

a=2 (6)

(2n+1) S2r+1 sin((2r + 2)x)

have period ;

b 1 1 64 × ··· 36 1 − 64a 4(r − 1)2 1 −

a 4r 2

b2 4r 2

−

1 4r 2 Zr+1

(12)

2

(7)

which are corresponding to characteristic numbers denoted by a2n ; b2n+1 ; a2n+1 , and b2n+2 . In every case, satisfaction of the diEerential equation necessitates recurrence formulae connecting three successive coeIcients. The continued fraction forms the basis of the technique developed by Ince [37] and Goldstein [38] for the computation of the characteristics numbers. If each series of (4)–(7) is substituted in turn in the Mathieu’s equation (2), and the coeIcients of cos(2rt); cos((2r +1)t); sin((2r +1)t); sin((2r +2)t) equated to zero for r = 0; 1; 2; : : :, the following recurrence relations are obtained, respectively

− 4rb 2 Zr = ; 1 − 4ra2 − 4r12 Zr+1

lim Zr = 0

r→∞

for cos(2rt), to draw the stability chart of Fig. 1. A similar method was used by McLachlan to report Table 1 for stability chart of the Mathieu’s equation (2). In Table 1, the numerical values related to aci and asi are the stability boundaries for ith cosine elliptic (cei ) and ith sine elliptic (sei ) functions, respectively. Fig. 1 is a graphical representation of Table 1. It can be seen in Table 1 and Fig. 1 that continued fractions method is strong enough to detect the stability boundaries for large values of parameters. The method of continued fractions was used by later investigators, and it is shown that the continued fraction

Table 1 McLachlan’s report for Mathieu stability chart ac0

as1

ac1

as2

ac2

as2

ac3

as4

ac4

as5

ac5

as6

0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 24, 28, 32, 36,

0, −0.4551386, −1.5139569, −2.8343919, −4.2805188, −5.8000460, −7.3688308, −8.9737425, −10.6067292, −12.2624142, −13.93698, −17.3320660, −20.7760553, −24.2586795, −27.7728422, −31.3133901, −38.4589732, −45.6733696, −52.9422230, −60.2555679,

1, −0.1102488, −1.3906765, −2.7853797, −4.2591829, −5.7900806, −7.3639110, −8.9712024, −10.6053681, −12.2616617, −13.9365525, −17.3319184, −20.7760004, −24.2586578, −27.7728332, −31.3133862, −38.4589724, −45.6733694, −52.9422229, −60.2555679,

1, 1.8591081, 2.3791999, 2.5190391, 2.3180082, 1.8581875, 1.2142782, 0.4383491, −0.4359436, −1.3867016, −2.3991424, −4.5701329, −6.8934005, −9.33523671, −11.8732425, −14.4913014, −19.9225956, −25.5617471, −31.3651544, −37.3026391,

4, 3.9170248, 3.6722327, 3.2769220, 2.7468810, 2.0994604, 1.3513812, 0.5175454, −0.3893618, −1.3588101, −2.3821582, −4.5635399, −6.8907007, −9.3341097, −11.8727265, −14.4910633, −19.9225403, −25.5617329, −31.3651505, −37.3026380,

4, 4.3713010, 5.1726651, 6.0451969, 6.8290748, 7.4491097, 7.8700645, 8.0866231, 8.1152388, 7.9828432, 7.7173698, 6.8787369, 5.7363123, 4.3712326, 2.8330567, 1.15422829, −2.5397657, −6.5880630, −10.9143534, −15.4667703,

9, 9.0477393, 9.14062277, 9.2231328, 9.2614461, 9.2363277, 9.1379058, 8.9623855, 8.7099144, 8.3831192, 7.9860691, 7.0005668, 5.7926295, 0.3978962, 2.8459917, 1.1607057, −2.5380779, −6.5875850, −10.9142090, −15.4667243,

9, 9.0783688, 9.3703225, 9.91155063, 10.6710271, 11.5488320, 12.4656007, 13.3584213, 14.1818804, 14.9036797, 15.5027844, 16.3015349, 16.5985405, 16.4868843, 16.0619754, 15.3958109, 13.5228427, 11.1110798, 8.2914962, 5.1456363,

16, 16.0329701, 16.1276880, 16.2727012, 16.4520353, 16.6482199, 16.8446016, 17.0266608, 17.1825278, 17.3030110, 17.3813807, 17.3952497, 17.2071153, 16.8186837, 16.2420804, 15.4939776, 13.5527965, 11.1206227, 8.2946721, 5.1467375,

16, 16.0338323, 16.1412038, 16.3387207, 16.6468189, 17.0965817, 17.6887830, 18.4166087, 19.2527051, 20.1609264, 21.1046337, 22.9721275, 24.6505951, 26.0086783, 26.9877664, 27.5945782, 27.8854408, 27.2833082, 26.0624482, 24.3785094,

25, 25.0208408, 25.0833490, 25.1870798, 25.3305449, 25.5108160, 25.7234107, 25.9624472, 26.2209995, 26.4915472, 26.7664264, 27.3000124, 27.7697667, 28.136359, 28.3738582, 28.4682213, 28.2153594, 27.4057488, 26.1083526, 24.3960665,

25, 25.0208543, 25.0837778, 25.1902855, 25.3437576, 25.5499717, 25.8172720, 26.1561202, 26.5777533, 27.0918661, 27.7037687, 29.2080550, 31.0000508, 32.9308951, 34.8530587, 36.6449897, 39.5125519, 41.2349503, 41.9535112, 41.9266646,

36, 36.0142899, 36.0572070, 36.1288712 36.22944114 36.3588668 36.5170667 36.7035027 36.9172131 37.1566950 37.4198588 38.0060087 38.6484719 39.3150108 39.9723511 40.5896641 41.6057099 42.2248415 42.3939428 42.1183561

1323

method can be applied to all linear parametric diEerential equations, but it is not applicable to non-linear parametric equations [39]. Due to this problem, we need a more reliable and accurate method to be applicable not only for linear parametric equation with large coeIcients, but also to non-linear equations. In the next section, we use the principle of the energy integral combined with numerical integration to develop an applied method to determine the transition curves.

4. Energy-rate method and stability chart algorithm

(13)

(14)

(15)

Consider that we are looking for a boundary of stability for the system

x7 + f(x) + g(x; x; ˙ t) = 0;

where

g(0; 0; t) = 0

(f(x)x) ˙ dt (16)

and f(x) is a single variable function, and g(x; x; ˙ t) is a non-linear time-dependent function. The functions f and g are depending on a set of parameters. Writing Eq. (13) in the following form:

1 2 x˙ + 2

x7 + f(x) = −g(x; x; ˙ t)

shows that the system is a model of a unit mass particle attached to a conservative spring, acted upon by a non-conservative force −g(x; x; ˙ t). The free motion of the system is governed by x7 +f(x)=0. In general, the applied force generates or absorbs energy depending on the value of parameters, position and velocity of the particle, and also time. If f(x) is a restoring force and g(x; x; ˙ t) is a periodic function of time, then we should expect a tendency to oscillate [40]. De3ning a kinetic, potential, and mechanical energy function  for the system by V (x) = f(x) d x; T (x) ˙ = 21 x˙2 , and E =T (x)+V ˙ (x), respectively, we can write an integral of energy as follows: 



d d E˙ = (E) = dt dt

= −xg(x; ˙ x; ˙ t):

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

b

1324

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

This expression represents the rate of energy generated or absorbed by the term −g(x; x; ˙ t). Suppose that for some set of parameters and a non-zero response x(t); E˙ = −xg(x; ˙ x; ˙ t) is negative then, E continuously decreases along the path of x(t). The eEect of g(x; x; ˙ t) resembles damping or resistance; energy is continuously withdrawn from the system, and this produces a general decrease in amplitude until the rate of initial amplitude depended energy runs out and a new solution make E˙ = 0. On the other hand, if for a set of parameters and a non-zero response x(t); E˙ = −xg(x; ˙ x; ˙ t) is positive then, the amplitude increases so long as the path of x(t) runs away. Evaluating a suitable numerical integration can show if a set of parameters belongs to a ˙ stable or unstable region by evaluating E. It is known that time derivative of mechanical energy E must be zero over one period for conservative and autonomous systems in a steady state periodic cycle [41]. We may integrate the Mathieu’s equation (2), to get the following equation: 1 d 2 ˙ (17) (x˙ + ax2 ) = 2bxx˙ cos(2t) = E; 2 dt where E is the mechanical energy of the system. In order to 3nd a set of a and b to indicate a steady-state periodic response, we may choose a pair of parameters a and b and integrate Eq. (17) numerically. We may evaluate the averaged energy over a period 1 T 1 T ˙ Eav = E dt = (2bxx˙ cos(2t)) dt; T 0 T 0 T = 2

(18)

to compare with zero. If Eav is greater than zero, (a; b) belongs to a region that energy being inserted to the system and then Eq. (2) is unstable. However, if it is less than zero, then (a; b) belongs to a region that energy being extracted from the system and Eq. (2) is stable. On the common boundary of these two regions, Eav =0. Assume b is 3xed and we are looking for a to be on a boundary that its left-hand side is stable and its right-hand side is unstable then, if the chosen parameters show that Eav is less than zero, increasing a increases Eav . On the other hand, if Eav is greater than zero, decreasing a decreases Eav . Using this strategy we may 3nd the appropriate value of a such that Eav equates zero. In this way, we will 3nd a point on the

boundary that has a stable region on its left-hand side and unstable region on its right-hand side. This group of points constitutes a branch of -periodic boundary corresponding to ce2n or se2n+2 . By changing the strategy and looking for a to be on a boundary that its right-hand side is stable and its left-hand side is unstable then, we can generate the 2-periodic branches corresponding to ce2n+1 or se2n+1 . Now we may change b by some increment and repeat the procedure. This procedure could be arranged in an algorithm called “stability chart algorithm”.

4.1. Stability chart algorithm 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.

set a, equal to one of its special values; set b, equal to some arbitrary small value; solve the diEerential equation numerically; evaluate Eav ; decrease (increase) a, if Eav ¿ 0 (Eav ¡ 0) by some small increment; the increment of a must be decreased if Eav (ai ) · Eav (ai−1 ) ¡ 0; save a and b when Eav 1; while b ¡ b3nal , increase b and go to step 3; set a, equal to another special value and go to step 2; reverse the decision in step 5 and go to step 1.

It is known that the boundary between stable and unstable regions in the stability chart starts from some special values of a. So, it is better to start the algorithm from one of these characteristic values.

5. Applying the stability chart algorithm In order to 3nd more accurate values than those presented in Table 1, we have used the elements of McLachlan’s table as an initial value. Then, the Mathieu’s equation was solved numerically for 0 ¡ t ¡ T; T =2, in order to include both -periodic and 2-periodic solutions. Two solutions of identical period, even  or 2, bound the region of instability, while two solutions of diEerent periods bound the region of stability.

Table 2 Stability chart for Mathieu’s equation generated by the algorithm ac0

as1

ac1

as2

ac2

as2

ac3

as4

ac4

as5

ac5

as6

0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 24, 28, 32, 36,

0, −0.455138604, −1.5139568851, −2.834391889959, −4.280518818368, −5.8000460207, −7.36883083214, −8.97374250574, −10.60672923566, −12.2624142182, −13.93697995675, −17.33206603507, −20.77605531238, −24.25867947476, −27.77284216366, −31.31339007051, −38.4589731694, −45.673369663942, −52.94222296412, −60.25556789136,

1, −0.110249001, −1.3906767, −2.7853799, −4.2591831, −5.7900808, −7.3639112, −8.9712026, −10.6053683, −12.2616619, −13.9365527, −17.3319186, −20.7760006, −24.258658, −27.7728334, −31.3133864, −38.4589726, −45.673369545, −52.9422229395, −60.2555681,

1, 1.8591084, 2.3791999, 2.519039087, 2.31800817, 1.8581875414, 1.2142781642, 0.43834908996, −0.43594360159, −1.386701566696, −2.399142400351, −4.570132851088, −6.8934005,3343 −9.33523701, −11.87324251702, −14.49130142564, −19.92259564555, −25.561747709501, −31.36515444857, −37.30263912327,

4, 3.9170245, 3.6722324, 3.2769217, 2.7468807, 2.0994601, 1.3513809, 0.5175451, −0.3893621, −1.3588104, −2.3821585, −4.5635402, −6.890701, −9.33411, −11.8727268, −14.4910636, −19.9225406, −25.5617332, −31.3651508, −37.3026383,

4, 4.3713010, 5.1726648, 6.0451972, 6.8290745, 7.4491094, 7.8700648, 8.086623145, 8.11523883, 7.982843163, 7.7173698494, 6.87873685481, 5.736312349375, 4.371232605296, 2.833056732297, 1.1542285900001, −2.539765704347, −6.588062973934, −10.91435338774, −15.46677033703,

9, 9.04773900, 9.14062247, 9.2231325, 9.2614458, 9.2363274, 9.1379055, 8.9623852, 8.7099141, 8.3831189, 7.9860688, 7.0005665, 5.7926292, 0.3978959, 2.8459914, 1.1607054, −2.5380782, −6.5875853, −10.9142093, −15.4667246,

9, 9.0783688, 9.3703225, 9.91155033, 10.6710271, 11.5488317, 12.465601, 13.358421, 14.1818807, 14.90368, 15.5027847, 16.301534946, 16.598540469, 16.4868842565, 16.06197536045, 15.39581091191, 13.52284271453, 11.1110798375, 8.291496150361, 5.145636255447,

16, 16.0329698, 16.1276877, 16.2727009, 16.452035, 16.6482196, 16.8446013, 17.0266605, 17.1825275, 17.3030107, 17.3813804, 17.3952494, 17.207115, 16.8186834, 16.2420801, 15.4939773, 13.5527962, 11.1206224, 8.2946718, 5.1467372,

16, 16.0338323, 16.1412038, 16.3387207, 16.64681817, 17.0965817, 17.688783, 18.416609, 19.2527054, 20.1609267, 21.1046334, 22.9721278, 24.6505954, 26.008678, 26.98776644, 27.594578154, 27.8854407976, 27.28330817057, 26.062448443, 24.37850942577,

25, 25.0208405, 25.0833487, 25.1870795, 25.3305446, 25.5108157, 25.7234104, 25.9624469, 26.2209992, 26.4915469, 26.7664261, 27.3000121, 27.7697664, 28.1363587, 28.3738579, 28.468221, 28.2153591, 27.4057485, 26.1083523, 24.3960662,

25, 25.0208543, 25.0837778, 25.1902855, 25.3437576, 25.5499717, 25.8172720, 26.1561202, 26.5777533, 27.0918661, 27.7037687, 29.2080550, 31.0000505, 32.9308948, 34.8530584, 36.6449894, 39.5125519, 41.234950267, 41.9535111621, 41.92666456905,

36 36.0142899, 36.0572070, 36.1288709 36.22944084 36.3588665 36.5170664 36.7035024 36.9172128 37.1566947 37.4198585 38.0060084 38.6484716 39.3150105 39.9723508 40.5896638 41.6057096 42.2248412 42.3939425 42.1183558

0

Fig. 2. Plot of TEav as a function of a for b = 1.

E˙ dt =

n−1

h ˙ ˙ [E(kh) + E((k + 1)h)]; 2

k=0

1325

The averaged energy Eav is evaluated by the simple trapezoidal integration method T T (2bxx˙ cos(2t)) dt

0

=

T (19) h = ; n = 1000; n using 1000 segments over the period T . Following the algorithm presented in the previous section, the value of a must be decreased (or increased) if Eav ¿ 0 (or Eav ¡ 0). If Eav at step i was greater than zero and Eav at step i − 1 was less than zero, then the value of increment which must be added to a (or subtracted from a) would be decreased at step i + 1 in order to prevent cyclic repetition. We have applied the algorithm and found the results tabulated in Table 2, which could be compared with Table 1. Distribution of regions of instability will be more clear if we let b → 0 then, using relations (8)–(11), we 3nd that the solutions with a 2-period lie in pairs near the characteristic values a=(2k +1)2 ; k =0; 1; 2; 3; : : :, and the solutions with a -period lie in pairs near the characteristic values a = (2k)2 ; k = 0; 1; 2; 3; : : : . Both cases can be combined in one formula to indicate characteristic values, a = k 2 ; k = 0; 1; 2; 3; : : : . Indeed, 2-periodic transition curves start from a = 1; 32 ; 52 ; : : :, and -periodic transition curves start from a = 0; 22 ; 42 ; : : :. In order to investigate the applicability of the method, we have plotted TEav versus a in Fig. 2, for

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

b

1326

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

the line b = 1 shown in Fig. 1. It can be seen that the curve TEav cuts the a axis exactly at the same points where line b = 1 cuts the transition curves. Therefore, Fig. 1 could be assumed as intersection of a three-dimensional surface T q(a; b) = TEav = (2bxx˙ cos(2t)) dt 0

and the plane TEav = 0. Due to physical consideration, the 3rst instability region, started at a = 1, is the most dangerous and has therefore the greatest practical importance. Bolotin [33] calls this region the “principal region of dynamic instability”. In the next section, we compare Tables 1 and 2, and show the advantages of the stability chart algorithm. 6. Comparison, application, and modi%ed algorithm How important the determined transient curve of a parametric equation is, depends on the physical system, real period of free oscillation, closeness of design parameters to the transient curves, mechanism and rate of dissipation system. If mistakenly design parameters belong to an unstable region instead of periodic or stable region, then parametric resonance might occur. Table 3 shows the value of TEav for the 3rst ten digits of -periodic branch ac0 of Tables 1 and 2. The last three rows indicate that how important is the accuracy of 3nding the transition lines. On the 10th row, decreasing 0.00000004325 unit could reduce TEav from

Fig. 3. Time response of the diEerence of amplitude.

104.6426231398 to −0:07985525604 signifying that point (−13:93697995675; 10) is much more closer to the transient curve than point (−13:93698; 10). We compare Tables 1 and 2 by study the Mathieu’s equation (2) for two pairs of a and b equal to (−0:4551386; 1) and (−0:455138604; 1); denoting by x1 and x2 , respectively. The pairs of a and b belong to two identical elements of Tables 1 and 2, and Table 3 says that TEav ¡ 0 for both of them. Although we cannot recognize any diEerence between the plots of time responses of the Mathieu’s equation for x1 and x2 , the function y =|x1 −x2 | for 0 ¡ t ¡ 10 shown in Fig. 3 indicates that the response of Mathieu’s equation for x1 and x2 are not identical. It can be seen that the diEerence of the solutions is growing up, which shows that the solution of Mathieu’s equation for point x1 deviates from periodic response faster than its solution for x2 .

Table 3 The values of TEav for 3rst 10 digits of ac0 for a and b from Tables 1 and 2 b

a, Table 1

a, Table 2

TEav for Table 1

TEav for Table 2

1 2 3 4 5 6 7 8 9 10

−0.4551386 −1.5139569 −2.8343919 −4.2805188 −5.800046 −7.3688308 −8.9737425 −10.6067292 −12.2624142 −13.93698

−0.455138604 −1.5139568851 −2.83439189959 −4.280518818368 −5.8000460207 −7.36883083214 −8.97374250574 −10.60672923566 −12.2624142182 −13.93697995675

−0.41036782785e − 6 0.000033303989524 0.000269337200715 −0.004127213067066 −0.03087386648124 −0.2604694332158 −0.2235618874348 −5.978866028905 −12.03902889408 104.6426231398

−0.1213741795e − 10 0.147415342e − 7 0.58511601e − 7 −0.382606859e − 6 −0.000483913357634 0.000306168820171 0.000289910400872 0.002771638776375 −0.0061065122282 −0.07985525604

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

Fig. 4. Time response for a point on the boundary of stability.

Fig. 5. TEav as a function of time for a point on the boundary of stability.

The time response shown in Fig. 4 depicts that the period of oscillation is . It is also indicated in the plots of time history of TEav shown in Fig. 5. The plot of TEav in Fig. 5 has an interesting behaviour. The time integral of TEav is zero for 0 ¡ t ¡ , but the value of time integral of TEav for 0 ¡ t ¡ =2 is minus of the time integral of TEav for =2 ¡ t ¡ . If time integral of Eav for 0 ¡ t ¡ T is minus of the time integral of Eav for T ¡ t ¡ 2T we call that Eav has skew-symmetric property. Whenever Eav has skew-symmetric property then, the system is 2T -periodic instead of T -periodic. Hence, searching for responses having period T = n, and =n; n = 1; 2; 3; : : :, also includes responses having period T = 2n and 2=n; n = 1; 2; 3; : : : . Now another advantage of energy-rate method will be revealed that are not easily available by traditional

1327

Fig. 6. Plot of TEav as a function of a for b = 5 integrated for diEerent period of time.

perturbations and approximation methods. We analyze super and sub-harmonic oscillation, and open new doors of research to discover periodic responses of the parametric equations. Although mathematical theories can provide some analytic basis, we may use a cross section of the three-dimensional energy plots, similar to Fig. 2, to 3nd the possibility of super and sub-harmonic oscillations. Integration of the Mathieu’s equation for a constant b, and plotting Eav versus a indicates where Eav is negative (stability) and positive (instability). If the plot of Eav versus a, integrated for T = n or =n; n = 1; 2; 3; : : :, crosses the axis Eav = 0, then there is a transition or periodicity curve separating stable and unstable regions. Fig. 6 depicts a plot of TEav for the Mathieu’s equation as a function of a. The value of b is set to 5. The line b = 5 is plotted in Fig. 1 to give a view of the domain of stability and instability. The period of integration of TEav in Fig. 6 is diEerent for each curve. It is seen that the curve for T = =3 and =2 has just one zero for positive a, close to a = 9:1 and 7.4, respectively. Time response of the Mathieu’s equation at those points would be periodic, although both of them belong to instability regions of the Mathieu’s stability chart. Fig. 7 shows time response of the Mathieu’s equation for (a = 7:45; b = 5). The curve for T = in Fig. 6 intersects the horizontal axis exactly at points where the line b = 5 intersects the transition curves in Fig. 2, but the curve for T =2 not only intersects at the -periodic and 2-periodic transition curves of the Mathieu’s equation, but also

1328

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

Fig. 7. Time history of -periodic response of Mathieu’s equation for (a = 7:45; b = 5), a periodic response in instability region.

Fig. 9. Time history of 3- and 6-periodic response of Mathieu’s equation for (a = 12:57; b = 5) and (a = 14:54; b = 5).

tion? Figs. 6, 8 and 9 indicate that it is possible. Now we may have for instance, a 2-periodic curve, where both left- and right-hand sides of that curve are stable. It is a periodic line “splitting” a stable region. We may also have a 2-periodic curve, where both left- and right-hand sides are unstable; a periodic line “splitting” an unstable region. The same situation could be seen for other sub and super harmonics. We may now modify the stability chart algorithm to be able to detect splitting lines. The following algorithm is prepared to detect splitting lines of the Mathieu’s equation. Fig. 8. Time history of 4-periodic response of Mathieu’s equation for (a = 8:13; b = 5) and (a = 13:47; b = 5).

touches the horizontal axis at some other points. For instance, the 3rst touching point is around a = 8:13 which is a point in stable region of Mathieu’s stability chart. Fig. 8 depicts two 4-periodic time responses using the 3rst two touching points (a = 8:13; b = 5) and (a = 13:47; b = 5) of the curve T = 2 in Fig. 6. These two periodic solutions were found when we were looking for 2-periodic responses due to skew symmetric characteristic of Eav . A 3-periodic and a 6-periodic responses are shown in Fig. 9 for two touching points (a = 12:57; b = 5) and (a = 14:54; b = 5) of the curve T = 3 in Fig. 6. These two super harmonic responses belong to the stability region of the Mathieu’s equation. Is it possible to have other periodicity lines in the stability or instability is regions of Mathieu’s equa-

6.1. Splitting chart algorithm 0. set period of integration T equal to a multiple of ; 1. set a, equal to a1 , and a2 , two arbitrary values in a stability region; 2. set b, equal to some arbitrary small value; 3. solve the diEerential equation numerically for both pairs of a, and b; 4. evaluate (Eav )1 and (Eav )2 ; 5. set a1 = a2 ; 6. increase a2 , by some small increment if |(Eav )1 | ¿ |(Eav )2 | and (a2 − a1 ) ¿ 0 or |(Eav )1 | ¡ |(Eav )2 | and (a2 − a1 ) ¡ 0 else decrease a2 , if |(Eav )1 | ¿ |(Eav )2 | and (a2 − a1 ) ¡ 0 or |(Eav )1 | ¡ |(Eav )2 | and (a2 − a1 ) ¿ 0;

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

b

12

12

10

10

8

8 b

6

4

2

2 0

b

2

4

8

10

12

0

10

10

8

8 b

6

4

2

2 4

6 a

8

10

12

0

12

10

10

8

8 b

6

4

2

2 4

6 a

8

10

12

2

4

6 a

8

10

12

2

4

6 a

8

10

12

6

4

2

6 a

(d)

12

0

4

6

4

2

2

(b)

12

0

(e)

6 a

12

(c)

b

6

4

(a)

8

10

12

0 (f)

1329

Fig. 10. (a) Mathieu’s -periodicity chart and splitting curves, based on the “splitting chart algorithm”; (b) Mathieu’s 2-periodicity chart and splitting curves, based on the “splitting chart algorithm”; (c) Mathieu’s 3-periodicity chart and splitting curves, based on the “splitting chart algorithm”; (d) Mathieu’s 4-periodicity chart and splitting curves, based on the “splitting chart algorithm”; (e) Mathieu’s 5-periodicity chart and splitting curves, based on the “splitting chart algorithm”; and (f) Mathieu’s periodicity chart and splitting curves, for ; 2; 3; 4, and 5-periodic solutions.

1330

7. 8. 9. 10.

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331

the increment of a must be decreased if a2 repeats twice; save a and b when (a2 − a1 )0; while b ¡ b3nal , increase b and go to step 1; set a, equal to a1 , and a2 , two arbitrary values in another stability region and go to step 2.

It should be mentioned that the algorithms can be extended to cover the other parametric diEerential equations whose stability chart is two dimensional or could be reduced to two dimensional. Applying the algorithms on parametric equation with more than two parameters produces a set of stability charts by varying only two parameters and 3xing the others. We have applied the “splitting chart algorithm” to 3nd the splitting curves of Mathieu’s equation, and Fig. 10(a) – (e) illustrates the results for ; 2; 3; 4, and 5-periodic solutions on the second and third stable regions. Fig. 10(f) depicts all splitting curves related to ; 2; 3; 4, and 5-periodic solutions. The algorithm can easily be applied to a wider range, and also to the other periodic splitting curves. 7. Closure We have developed a stability chart algorithm based on energy-rate method, to 3nd the stability boundaries of the parametric diEerential equations. The algorithm provides more informative results than perturbation and approximation methods. The algorithm has been applied to the Mathieu’s equation, and its stability chart has been found, comparable to previously reported results. Although the procedure is sometimes time consuming, it generates a more eIcient and more accurate stability chart. The idea may be used to detect other periodicity curves not necessarily dividing stability and instability regions. It may be applied to detect splitting curves, a periodicity curve embedded in a stable or unstable region. The algorithm can be applied to any linear parametric diEerential equation as well as non-linear time-varying equations. Acknowledgements The author acknowledges M. Rastgaar Aagaah as an assistant who developed the computer program to apply the energy-rate Algorithm.

References [1] N.W. McLachlan, Theory and Application of Mathieu Functions, Clarendon, Oxford, 1947, U.P., Reprinted by Dover, New York, 1964. [2] E.T. Whittaker, G.N. Watson, A Course of Modern Analysis, 4th Edition, Cambridge University Press, Cambridge, 1927. [3] B. Van der Pol, M.J.O. Strutt, On the stability of the solution of Mathieu’s equation, Philos. Mag. 5 (1928) 18. [4] W.G. Bickley, The tabulation of Mathieu equations, Math. Tables Other Aids to Comput. 1 (11) (1945) 409–419. [5] C. Hayashi, Nonlinear Oscillations in Physical Systems, Princeton University Press, Princeton, 1964. [6] F.R. Moulton, W.D. MacMillan, On the solutions of certain types of linear diEerential equations with periodic coeIcients, Amer. J. Math. 33 (10) (1911) 65–96. [7] G.F. Carrier, On the Nonlinear vibration problem of the elastic string, Quart. Appl. Math. 3 (2) (1945) 157–165. [8] J.W.S. Rayleigh, The Theory of Sound, Vol. 1, Dover Publications, New York, 1945 (Macmillan Company, 1894). [9] A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, Dover, New York, 1927. [10] L. Meirovitch, Analytical Methods in Vibrations, Macmillan, New York, 1967. [11] M. Faraday, On a particular class of acoustical 3gures; and on certain forms assumed by a group of particles upon vibrating elastic surfaces, Philos. Trans. R. Soc. London 121 (1831) 299–318. [12] E.R.M. Iwanowski, On the parametric response of structures, Appl. Mech. Rev. 18 (9) (1965) 699–702. [13] J.A. Richards, Analysis Periodically Time Varying Systems, Springer, New York, 1983. [14] E. Mathieu, Memoire sur le mouvement vibratoire d’une membrane de forme elliptique, J. Math. Pure Appl. 13 (1868) 137. [15] G.W. Hill, On the part of the moon’s motion which is a function of the mean motions of the sun and the moon, Acta Math. 8 (1886) 1–36. [16] J.W. Strutt, (Lord Rayleigh), On the crispations of Luid resting upon a vibrating support, Philos. Mat. 16 (1883) 50 –53. [17] N.M. Beliaev, Stability of prismatic rods subject to variable longitudinal forces, Collection of Papers: Engineering Constructions and Structural Mechanics, Leningrad, 1924, pp. 149 –167. [18] E.A. Beilin, G.H. Dzhanelidze, Survey of work on the dynamic stability of elastic systems, Prikl. Math. I Mekh. 16 (5) (1952) 635–648. [19] J.C. Yao, Dynamic stability of cylindrical shells under static and periodic axial and radial loads, AIAA J. 1 (6) (1963) 1391–1396. [20] E. Esmailzadeh, G. Nakhaie Jazar, Periodic solution of a Mathieu–DuIng type equation, Internat. J. Nonlinear Mech. 32 (5) (1997) 905–912. [21] E. Esmailzadeh, G. Nakhaie Jazar, Periodic behaviour of a cantilever with end mass subjected to harmonic base excitation, Internat. J. Nonlinear Mech. 33 (4) (1998) 567–577.

G.N. Jazar / International Journal of Non-Linear Mechanics 39 (2004) 1319 – 1331 [22] A.C.J. Luo, Chaotic motion in the generic separatrix band of a Mathieu–DuIng oscillator with a twin-well potential, J. Sound Vib. 248 (3) (2001) 521–532. [23] R.S. Zounes, R.H. Rand, Global behaviour of a nonlinear quasiperiodic Mathieu equation, Proceedings of DETC01, ASME 2001 Design Engineering Technical Conference, Pittsburgh, PA, September 2001, pp. 1–11. [24] R.E. Langer, The solution of the Mathieu equation with a complex variable and at least one parameter large, Trans. Amer. Math. Soc. 36 (3) (1934) 637–695. [25] F. Colonius, W. Kliemann, Stability of time varying systems, ASME Design Engineering Technical Conferences DE. 84-1, Vol. 3, Part A, 1995, pp. 365 –373. [26] R.S. Guttalu, H. Flashner, An analytical study of stability of periodic systems by Poincare mappings, ASME Des. Eng. Tech. Conf. DE 84-1 Part A 3 (1995) 387–398. [27] E.A. Butcher, S.C. Sinha, On the analysis of time-periodic nonlinear hamiltonian dynamical systems, ASME Des. Eng. Tech. Conf. DE 84-1 Part A 3 (1995) 375–386. [28] D.W. Jordan, P. Smith, Nonlinear Ordinary DiEerential Equations, Oxford University Press, Oxford, 1999. [29] J. Taylor, K. Narendra, Stability regions for the damped Mathieu equation, SIAM J. Appl. Math. 17 (1969) 343–352. [30] H. Gunderson, H. Rigas, F.S. VanVleck, A technique for determining stability regions for the damped Mathieu equation, SIAM J. Appl. Math. 26 (2) (1974) 345–349.

1331

[31] P.H. Dawson, N.R. Whetten, Ion storage in three dimensional, rotationally symmetric, quadrupole 3eld. I. Theoretical treatment, J. Vac. Sci. Technol. 5 (1968) 1–6. [32] R. Rand, R. Hastings, A quasiperiodic Mathieu equation, ASME Des. Eng. Tech. Conf. DE 84-1 Part A 3 (1995) 747–758. [33] V.V. Bolotin, The Dynamic Stability of Elastic Systems, Holden-Day, San-Francisco, USA, 1964. [34] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillation, Wiley, New York, USA, 1979. [35] P. Pederson, On stability diagrams for damped Hill equations, Quart. Appl. Math. 42 (1985) 477–495. [36] L. Cesari, Asymptotic Behaviour and Stability Problems in Ordinary DiEerential Equations, 2nd Edition, Academic Press, New York, USA, 1964. [37] E.L. Ince, Tables of the elliptic-cylinder functions, R. Soc. Edinburgh Proc. 52 (1932) 355–423. [38] S. Goldstein, Mathieu functions, Trans. Cambridge Philos. Soc. 23 (1927) 303–336. [39] G. Nakhaie Jazar, Analysis of nonlinear parametric vibrating systems, Ph.D. Thesis, Mechanical Engineering Department, Sharif University of Technology, 1997. [40] E. Esmailzadeh, B. Mehri, G. Nakhaie Jazar, Periodic solution of a second order, autonomous, nonlinear system, Nonlinear Dyn. 10 (1996) 307–316. [41] W. Szemplinska-Stupnicka, The Behaviour of Nonlinear Vibrating Systems, Kluwer Publishers, Dordrecht, 1990. [42] C.S. Hsu, Impulsive parametric excitation: Theory, Journal of Applied Mechanics 39 (1972) 551–559.