Application of the incremental harmonic balance method to cubic non-linearity systems

Journal of Sound and Vibration (1990) 140(2), 273-286

APPLICATION BALANCE

OF THE

METHOD Y.

TO K.

INCREMENTAL

CUBIC

CHEUNG

HARMONIC

NON-LINEARITY

AND

S. H.

SYSTEMS

CHEN

Department of Civil and Structural Engineering, University qf Hong Kong, Hong Kong AND s.

L.

LAU

Department of Civil and Structural Engineering, Hong Kong Polytechnic, Hong Kong (Received 10 November 1988, and in revised form

26 September 1989)

In this paper, the formulation of the incremental harmonic balance (IHB) method IS derived for a general system of differential equations with cubic non-linearity, which governs a wide range of engineering problems such as large-amplitude vibration of beams or plates. An incremental arc-length method combined with a cubic extrapolation technique is adopted to trace the response curve automatically. The stability of its periodic solutions can also be evaluated from the IHB formulation by multi-variable Floquet theory. Hsu’s method is adopted for computing the transition matrix at the end of one period. Two numerical examples are presented which demonstrate the effectiveness of the IHB method and Hsu’s method in treating the non-linear vibration of a strongly non-linear multiple degree of freedom system. 1.

INTRODUCTION

Usually, in non-linear analysis, the elastic system is first reduced to a discrete system by using the linear, free vibration modes as the basic functions to express the deflection as an expansion and a set of non-linear ordinary differential equations are obtained. Such equations, which usually contain quadratic and cubic non-linearities, have been investigated in the past by many researchers. Nayfeh and Mook have given an extensive bibliography up to 1979 in their book [ 11. Recently, a large amount of research has been devoted to the multiple degree of freedom (DOF) systems under various conditions such as parametric excitation [2,3], subharmonic or combination resonance [4-61. However, in most of the literature perturbation methods, especially the multiple scales method, have been used. It is well known that the common weakness of classical perturbation methods restricts them to solving problems with weak nonlinearity. Bajkowski and Szemplinska-Stupnicka [7] investigated the internal resonance effects in a two-degree-offreedom system by using both the analytical method and an analog computer simulation method. It shows that the amplitude-frequency response curves obtained by the averaging method and the simulation method are very different when the amplitude or the frequency shift becomes large. This example illustrates again that when dealing with strongly non-linear systems researchers must be very cautious about employing various perturbation methods. A practical weakness of perturbation methods is that carrying out the expansion to higher order is very cumbersome, especially for multiple DOF systems. In practice one seldom goes beyond the third order unless the algebraic manipulations are performed by a computer. Accordingly, it is desirable to establish a realistic and systematic computer method for analyzing strongly non-linear periodic vibrations of multiple DOF 273 0022-460X/90/140273+

14 $03.00/O

@

1990 Academic

Press Limited

274

Y.

K. CHEUNC

ET

AL,

systems. For such cases the incremental harmonic balance (IHB) method may be one of the most suitable approaches. Since the IHB method was presented in 1981 [8], it has been developed further and was successfully applied to the analysis of periodic and almost periodic non-linear structural vibrations and related problems [9-191. It is interesting to note that the IHB method has also been applied to the stationary responses of a multi-dimensional non-linear system with random input [20] and to the construction of chaotic regions [21]. Many advantages of the IHB method have been described in detail in the literature mentioned above. The most important one is that it is capable of dealing with strongly non-linear systems to any desired accuracy. The other one, which is worth pointing out, is that the IHB method is in fact a general computer method for finding the periodic solution of non-linear differential equations. The IHB method is a combination of the incremental method (Newton-Raphson procedure) with the harmonic balance method (Ritz and Galerkin averaging method). It is exactly equivalent to a Galerkin procedure followed by a Newton-Raphson method (HBNR). Ferri [22] has given a detailed comparison and discussion of the two methods. The HBNR method differs from the IHB method only in the order. Both methods produce the same set of linear equations. However, in many cases the linearized equations obtained by the IHB method would be easier to generate than those of the HBNR method. Moreover, in order to find a satisfactory approximation, not only the residue corresponding to the prescribed harmonics should be estimated, but also the residue corresponding to higher harmonics must be controlled. By using the IHB method, it is easier to select and adjust the proper harmonic terms automatically depending on the unbalanced harmonics. Nevertheless, we believe that the two basically equivalent methods (IHB and HBNR) do possess their own advantages and both have great potential for further development. Ling and Wu [23] presented a Fast Galerkin method, in which the Fast Fourier Transform (FFT) technique is introduced into the Galerkin method to reduce the amount of computation work and to obtain a high order approximation. Recently, Lau [24] has further developed the IHB method by introducing a time transformation so that it is able to obtain very satisfactory solutions with many fewer harmonic terms for those problems in which the vibration inherently contains a large number of harmonic terms. In this paper, efforts are devoted to the application of the IHB method to the cubic non-linearity systems including the determination of stability of the periodic solutions. For those systems containing quadratic and cubic non-linearities, the IHB method described can be easily applied by adding a quadratic term to the original differential equations. An incremental arc-length method combined with a cubic extrapolation technique is presented to trace the response curve automatically. After the steady state periodic solutions (both stable and unstable) are obtained, the determination of the stability of the solutions can be taken as a following exercise. Taking increments as disturbance of the solutions for the ordinary differential equations, one obtains a set of linear variational equations with periodic coefficients which can be used to determine the stability characteristics. Based on the Floquet theorem, the stability of the system is determined by a transition matrix, which can be evaluated by using Hsu’s method [25-271. Both the procedures of solving equations and determining the stability of solutions can be performed point by point along the frequency-response curves, so that response curves with stable

and unstable

characteristics

can be obtained.

2. INCREMENTAL HARMONIC BALANCE METHOD For a multi-degree-of-freedom system with cubic non-linearities, equations, in general, have the form

the non-linear

INCREMENTAL

co2 i A2&jj+ w f &jj+ 111 /=I

HARMONIC

i

BALANCE

&q, + i j=C

J=I

i

i

k=l

I=1

275

METHOD

Kj;;lqjqkqr =J; cos (2m - l)T, i=1,2

,...)

n,

(I)

where k”,‘k:=

r=Wt,

k(3) = /q(3) d,k

(2,3)

tkf,r

the 4, are the unknowns of the system, the dots denote derivatives with respect to dimensionless time T, and iii,, Cij, K,, Kj$‘, ,A and w are coefficients of the mass, damping, linear stiffness, cubic stiffness, excitation amplitude and exciting frequency respectively. Equation (1) can be rewritten in matrix form as w2~ii+w~‘il+(K+K’~))q=Fcos(2m-1)7,

(4)

where q=[q1,q2,...,qnlT, F=[f,,h,.. . ,fnlT, A, c and g are mass, damping, linear stiffness matrices, with elements are denoted by fiU, C, and Kij respectively, and g(i) is the cubic non-linear stiffness matrix, its element K!;,‘, being taken in the form

The first step of the IHB method is a Newton-Raphson denote a state of vibration; the neighboring state can corresponding increments to them as follows: 9, = qjo+ Aqj,

j = 1,2, . . . , ny

Substituting expansions (S)-(7) order, one obtains the following

i-1,2,.

f;=&+Adf;,

procedure. Let qiO, Jo and wg be expressed by adding the

.., n,

w = w,)+ Aw.

in to equation (1) and neglecting small terms of higher linearized incremental equation in matrix form:

w,:~A~+w,~A(l+(~+3~‘3’)Aq=~-(2w,~~,+~~o)Ao+cos(2m-1)~AF,

(8)

R=F,cos(2m-1)7-((w~~iii_ii~,+w”C~,+Kq”+K’~~q,),

in which qo, Aq, Fo, AF which goes to zero when The second step of the is odd and the excitation qj(b= z

3

(9)

and RF’ are given in the Appendix. R is a corrective vector the solution is reached. IHB method is the Galerkin procedure. Because equation (1) force is periodic, one can assume, for steady state response, U,kCOS(2k-l)T+

?

k=l

Aqj=

(5-7)

b,,sin(2k--l)r=C,A,,

(10)

Abjksin(2k-l)r=C,AAi,

(11)

h=l

AUjkCOS(Zk-l)T+?

k-l

k-l

where c, =

[COS 7, COS 3T, . . . , COs (2N,.-

l)~,

sin T, sin

3r,.

_. , sin (2N,, -

l)T],

Ai=[~,,,~j2,...,~iN,,~j~,~j2,...,b,~~lT,

AA, = [Aaj,)

Aa,>, . . . , AajN,, Ab, I, Abj,,

. . . , Abj,\]“.

Hence, the vectors of unknowns and their increments coefficients vector A and its increment AA as follows: qo =

S-4

Aq = SAA,

can be expressed

by the Fourier

(12,13)

276

Y.

K.

CHEUNG

ET

AL.

where S, A and AA are given in the Appendix. Substituting equations equation (8) and using the Galerkin procedure gives ?rr G(Aq)T[w;$iA~+w,CA~+(R+3jZ”‘)Aq] d7 I0

(12) and (13) into

2rr

s(Aq)*[R-(2w&I&,+~&,)Aw+cos(2m-l)~AF]dr.

=

I0 One can easily obtain

a set of linear

equations

in terms of AA, Ao and AF,

K,,AA=R-R,,Ao+R,AF,

(14)

in which K,, = c&M + w,C + K+ 3Kc3’ R=F-(~;M+w&+K+K(~))A,

R,,. = (2w,M+C)A.

(15) (1617)

M, C, K, Kc3), F and R, are given in detail in the Appendix. It is worth mentioning that in equation (14) the number of incremental unknowns is greater than the number of equations available due to the existence of AF and Aw. However, since one is primarily interested in the frequency-response curves of the system for a constant force level, F is fixed as a parameter vector, which implies AF = 0. Hence equation (14) is reduced to K,,AA

= R - R,,Aw.

(18)

The solution process starts from a suggested solution (in general, from a corresponding known linear solution), and then the non-linear amplitude-frequency response is solved point by point by incrementing frequency w or incrementing one component of the amplitudes A. The Newton-Raphson iteration can be applied within an incremental step. In the incremental process, an increment which is prescribed a priori is called a control or active increment. If Aw is specified as a control increment, then w remains constant are solved from the through the iterative process: i.e., Aw = 0, while other increments equation K,,AA

= R.

(19)

The process is repeated until the magnitude of the corrected vector R is acceptably small-in which case a solution is obtained. This process is called an iteration. The value of w is then augmented an increment Aw artificially, and a new iteration is repated with the new value of w until a new solution is obtained. The above process is called an augmentation. The whole solution process is an alternative application of augmentation and iteration. The above incremental process in which Aw is taken as active increment is called w-incrementation. Similarly, it is equally possible to have amplitude incrementation. In then aJk this case, one component of AA, say Aajk, is specified as the control increment; remains constant: i.e., Aaj, = 0 through the iteration and one has to solve equation (18) to obtain other increments of AA and Aw. After the amplitude of R has reached the desired accuracy, the iteration is terminated and a new augmentation can be started by adding an increment on the ajk. This process is called aj,-incrementation. In practice, the active increment is chosen as the one that varies faster and therefore the w-incrementation or the aj,-incrementation can be adopted along the response curves. If one is interested in the forcing amplitude, say A, response curves of the system for a constant frequency level, then Aw = 0 and AA = 0, j # i, and hence equation (14) is reduced to K,,.AA = R+R,AF.

(20)

INCREMENTAL

HARMONIC

BALANCE

METHOD

277

In this case, ajk-incrementation or &incrementation can be adopted along the forcing amplitude response curves. It is worth emphasizing that the method is successful in predicting the response curves in general by incrementation in amplitude or frequency. However, it may fail to converge at the sharp peaks (e.g., points B and C in Fjgure 5 of section 4.2) if the alternative incrementation is simply used. To overcome this shortcoming and to reduce the number of iterations required for the solution to converge, an incremental arc-length method combined with a cubic extrapolation technique is adopted. Let x,,, x, , x2 and x3 be known points; the next point xq will be predicted by increasing a prescribed arc-length As to extrapolate (see Figure 1). If the arc length t is used as parameter, then, with ti corresponding to xi, to = 0, t, = S, , f2 = t, + sz, f3 = t2 + sj, f4 = t, + As, one has

x4=; i=”

(

3 t4-

n

J=(J

lX t,-tj

j’l

t.

)

”

in which

Nf is the number of increments: i.e., the number of total degrees of freedom in equation (19).

Figure

1. Cubic

extrapolation

to predict

the next solution

point.

x0, x, ,

x2 and x3 are known points.

The choice of the magnitude of the increment length As is determined by the principal curvature of the hyperspace curve. The larger the curvature, the smaller will be the increment As. In practice, the following very simple approach is adopted [28]: As, = Asj_,Nd/Zj_, . Here AS,_, was adopted in the (j- 1)th augmentation step; 4-i is the number of iterations at the (j - 1)th augmentation step and Nd is the desired number of iterations which is selected based on experience. This procedure will automatically lead to small increment lengths in the areas of large curvature and long increment lengths in the areas of small curvature. Once the position of x is determined, the components of the increment vector AA are predicted. By comparing the magnitudes of these components, one can choose the largest one as the control increment and fix its value in the iteration process, while the others are taken as initial values of iteration. The control increment

278

Y.

K.

ET AL

C‘HEUNG

usually remains the same before the corresponding response curve reaches the peak. However, as the peak is approached, there will be another component coming up which will be chosen as the new control increment. As a result, automatic generation of the response curve is entirely possible. In our numerical work, all the response curves (e.g., shown in Figures 3-6 of section 4) have been automatically generated from the beginning to the end without any difficulty. 3. STABILITY OF STEADY STATE SOLUTIONS Once the steady state solution is determined, its stability is usually investigated superposing a small perturbation Aq on q,: that is, q=qo+Aq.

by (21)

Substituting equation (21) into equation (4), linearizing the resulting equation in terms of Aq and noting that q. satisfies equation (4), one finally obtains w2kA~+oCAtj+(ii+3ii’3’)Aq=0.

(22)

Equation (22) is called the perturbed equation: i.e., perturbed from the known solution equation (8) with ii = 0, do = 0 and A F = 0. Therefore the stability of the steady state solutions corresponds to the stability of the solutions of equation (22) which is a linear ordinary differential equation with periodic coefficients in ifc3). The stability characteristics can be studied by multi-variable Floquet theory. Equation (22) can be rewritten as qo. It can be obtained from the incremental

X = Q( r)X, where

1’

I

X = 14, Ad’,

(23)

- (l/W)WC

Qzl = -3

lt+‘(ii+

3iic3’), (24-26)

0 is a null matrix, and I is a unit matrix. Since each component of q. is a periodic function of 7 with a period T= 2r/w, each element of Q2, is also a periodic function with the same period T. For equation (23), there exists a fundamental set of solutions yk = FYI,,

where N = 2n. This fundamental matrix solution

y=

Yll

Y2k,.

k=l,2

. . , YNklT,

,...,

N,

set can be expressed in a matrix called a fundamental

Y.21

[I:YNI Clearly, Y satisfies the matrix equation

:I

Yl2

. . .

YIN

y22 .

*..

Y2N

YN~

“*

. . .

.

(28)

YNN

Y = Q( T)Y. Since Q(T+ T) = Q(T), Y( T+ T) is also a fundamental be expressed by

(29) matrix solution; therefore it can

Y(r+ T) =PY(T), where P is a non-singular

(27)

constant matrix called the transition matrix.

(30)

INCREMENTAL

HARMONIC

BALANCE

279

METHOD

According to the Floquet theory, the stability criteria for the system is related to the eigenvalues of the matrix P. If all the moduli of the eigenvalues of P are less than unity, the motion is bounded and therefore the solution is stable. Otherwise the motion is unbounded and the solution is unstable. One can numerically calculate a fundamental set of solutions of equation (29) by using the initial condition Y(0) = I.

(31)

According to equation (30), one obtains the transition matrix P = Y(T).

(32)

From the preceding discussion, it is clear that an efficient numerical computation of P is the key to efficient numerical treatment of stability analysis. Hsu has developed an efficient method for approximating the transition matrix (or the fundamental matrix) during one period [25-271. This method consists of evaluating the transition matrix by dividing a period into a number of equal parts and considering the equations over each interval to be a set of equations with constant coefficients. Friedmann summarized this method and gave a formulation clearly and concisely [29]. Here we just give the final formula. Suppose each period T is divided into Nk intervals denoted by rk, and the kth interval is denoted by Ak = Q - TV_, . In the kth interval, the periodic coefficient matrix Q(T) is replaced by a constant matrix Qk defined by Q&[’ I

Q(i')di-.

(33)

?-I Finally, the transition matrix is given in the form

1. 4.

(34)

NUMERICAL EXAMPLES

To illustrate the power of the IHB method, two numerical examples are presented in this section. The first one is a two-degree-of-freedom system with non-linear cubic-type characteristics subjected to a harmonic load. This example is selected from reference [7] to show the accuracy of the IHB method by comparing the results with those obtained by other numerical methods and analog computer simulation method. The internal resonance effect of the ratio of two natural frequencies in the same system is further investigated as the second example. The results for this example reveal some interesting behavior of internal resonances which have not been noted by other researchers before. 4.1.

NON-LINEAR CUBIC

VIBRATION

OF

TWO-DEGREE-OF-FREEDOM

SYSTEM

WITH

NONLINEARITIES

Consider a two-degree-of-freedom system consisting of two point masses and two springs with a linear damper, under a harmonic excitation shown in Figure 2. One of the springs is linear with the stiffness coefficient k ,. and the other has a cubic non-linear type. Its restoring force is defined by

fiz=k&q1- qd+ dq,

- qd3.

(35) Bajkowski and his co-workers investigated the internal resonance effects of this system [7]. They employed the analog computer simulation method and compared the results with those produced by analytical methods: the averaging method and the Ritz method. It is shown in their paper that the averaging method fails to predict the behavior of the

280

ET AL.

Y. K. CHEUNG

Figure

2. Mechanical

model

of two-degree-of-freedom

system.

system, but the Ritz method gives results close to those obtained by the simulation method. To demonstrate the effectiveness of IHB method, this system is investigated as the example. The equilibrium differential equations of motion of the system can easily be written in a non-dimensional form as ii, + k2q, + Y(41&+

q2) + PYl(41-

q2 - a-

A41

42) + CL7491 - 92j3 = P cos - 42) - A491

(37)

where t=i&JSG,

Y = m1/m2,

P=cLIk,z, l2=~~,

(36)

fit,

- q213 = 0,

idlQ&,

I=

k2 = k,oylk,,,

p=Isrlb,

Q=dq/dt.

_ _

q1 and q2 are displacements of point-masses, i is time, and m,, m2, k,,, k,2, @, 1, 0 and p are the masses of the system, coefficient of linear stiffness, coefficient of non-linear stiffness, coefficient of damping, excitation frequency and excitation amplitude respectively. Equations (36) and (37) are special forms of equations (4) in the case of n = 2, and the matrices $I, c, ii and P are of the following forms:

lz;:‘,,= #c;;:, =Eg:, = /.Ly 9 r;:‘,:,

=

E(3) _ 73) 2212 - k2,22

-

-

CL,

j73) _ EC31 1112 1121p3’ 1111 -PY,

p;‘,,

= _ cLy

Q3’ 2222 -

Pu,

pi\,

9 G%2

= p31 _ k(3) 2121 2112 -

= -

PY,

CL,

E(3) __~ 2111-

.

In Figures 3 and 5 are shown the response curves for k2 = l-5, y = 1.582, p = 2, 1= 2, p = 0.01 and R = 30 which are described in reference [7]. In this case, the second natural frequency of the system is three times the first natural frequency and therefore internal resonance will occur. In the solution process, the number of harmonic terms is taken as NC = N, = 2: i.e., q,=a,,cos~+a,,cos3~+bl,sin~+bl2sin3~=A,,cos(~+~,,)+A,,cos(3~+~,,), (38) q2=u~2cos~+a22cos3~+b2,sin~+b22sin3~=A2,cos(~+~2,)+A22~~s(3~+~22). (39) where A,=-,

&=tan-‘(-bJa,),

i=l,2,

j=1,2.

There exist two types of non-trivial solutions: (a) fundamental resonance only, i.e., 41 = A12 cos (37+ 4121, 42 = A,, cos (37+ +22); (b) both fundamental resonance and subharmonic resonance occur simultaneously: i.e., q1 and q2 take the forms of equations (38) and (39). In Figure 3 is shown the Ar2-0 fundamental response curve as Al, = A2, = 0 in equations (38) and (39): i.e., the first type of solutions. Figure 4(a) shows the A,,4 response curve and Figure 4(b) shows the A,,-0 response curve: i.e., Figure 4 shows the second type of solutions. The results obtained by the Ritz method and the simulation method [7] are

INCREMENTAL

HARMONIC

BALANCE

METHOD

281

30, 25 1

Frequency, Figure 3. Frequency response p = 2, L = 2 and /1 = 0.01. -,

D

curve of equations (36) and (37) when A,, = A,, = 0, with k’ = 1.5, y = l-582, 0, numerical [7]; 0, simulation [7]. Stable; - - -, unstable; - - -, backbone;

25 -

20: q g

15-

1 ‘a E a

lo-

5-

01

$ 2 = F a

4

0 196

200

2-04 Frequency, 0

2.08

Figure 4. Frequency-response curves of equations (36) and (37) with k2 = 1.5, y = 1.582, p = 2, L = 2 and cc = 0.01. -, Stable; - - -, unstable; 0, numerical [7]; 0, simulation [7].

plotted in the figures with dots and little circles respectively for comparison. It is readily seen that the IHB method gives results, both stable and unstable, almost identical to those obtained by Ritz method and close to those obtained by the simulation method. The physical interpretation of the shape of the subharmonic resonance curve can be found in reference [7]. It should be noted that it is important to choose the appropriate numbers of harmonic terms N, and N, in equation (lo), the number of integral intervals Nk and the number of Nj-the terms of the series of the matrix exponential in equation (34). Although there is no theoretical limitation to these numbers, there is obviously a CPU time, and hence

Y. K. CHEUNG

282

ETAL.

a computer cost limitation. Practically, in our examples N, = 2, N, = 2, NI = 200 and N, =4 are sufficient to obtain very good accuracy at a reasonable cost. 4.2. THE

INTERNAL

RESONANCE

Consider the system are given by

VERSUS

THE

of the first example

RATIO

again.

2

wlo,o:o=;(l+kL+y)~

OF

NATURAL

The natural

FREQUENC‘IE:S

frequencies

t(l+k’+y)-k2.

w,(, and w?,)

(40)

Obviously, wlo and w20 depend on the parameters k’ and y of the original system. According to earlier results, k2 and y are selected in such a way that w20 = 3w,,. Here, we further investigate the internal resonance influenced by the ratio of two natural frequencies. The free vibration of the undamped system is of primary importance for investigation of the resonance phenomena, so only free vibration is considered in this example. Four cases corresponding to different values of k2 and y as shown in Table I are investigated. Since the response curves are symmetrical with respect to the frequency axis, only the positive amplitudes are plotted in the figures. Figure 5(a)-(d) correspond to case IV in Table 1. They show the typical amplitudefrequency response curves A ,l-w, A,z-w, A,,+ and A,,-o respectively. A,, and Alz, defined in equation (38), are the amplitudes of the first and third harmonic terms of the of the first variable ql, while A,, and A,, , defined in equation (39), are the amplitudes first and third harmonic terms of the second variable q2, respectively. As in Figures 3 and 4, the solid lines represent stable solutions and the dashed lines represent unstable solutions. One can find from the figures that the response curves have the looping characteristics which are similar to those of the backbone curves of clamped-hinged beams [30-311 obtained by using the perturbation method. Moreover, these looping resonances are all unstable. Comparing Figures 5(a) and (c) to Figures 5(b) and (d), one can also see the internal resonance on the amplitudes of the responding modes. At the beginning, all the response curves start from the point A corresponding to the first natural frequency w,~, and then they climb up along the curve with increasing frequency. At the end of the stable response curves, both modes would jump down to their respective steady state values. After this, the first mode becomes larger and larger rapidly with increasing frequency, while the second mode increases very slowly until the saturation phenomenon occurs. The full features of free vibration for different parameters k2 and y are shown in Figures 6(a) and (b). The four A,,-@ and A12- w response curves are for the four cases I-IV in Table 1. The other two families of curves A,,-w and A,,-w are similar in shape and are not shown. From Figure 6, some to those of A,,-w and A22- w respectively interesting phenomena can be listed as follows. disappears. (1) In case I, w~~/w,~= 3.0, the looping resonance

TABLE

1

Four cases corresponding to dzxerent parameters k2

Y

I II

5.875 6.272

1.205 1.3

0.8988 0.8988

III IV

6.700 6.893

1.4 1448

0.8988 0.8988

Case

WI0

R=

w20/~1,,

3.0 3.1 3.2 3.25

INCREMENTAL

HARMONIC

BALANCE

283

METHOD

al

a.0 i

0

I It3

1

I

0

ul

I 0

a

284

Y.

K. C’HEUNG

ET

Al.

11

10 Frequency,

Figure

6. Free vibration

response

curves with different

12

w

R: 1, R =3.0;

II, R = 3.1; III,

R =3.2;

IV, R = 3.25.

(2) The looping resonance becomes larger with larger OJ~~/C.O,~ ratios, and the looping response curves are farther away from the starting point. However, these looping resonances are unstable. (3) No matter what the value of the ratio w~,,/w,~ is, the stable response curves invariably start from the same point corresponding to the first natural frequency CO,,,,and then proceed along the same track to their corresponding critical points, and after that the jump phenomenon occurs. (4) After crossing the looping resonance, all the response curves are under the response curve I.

CONCLUDING

REMARKS

1. The incremental harmonic balance method has been shown to be a straightforward, efficient and reliable method for treating the vibration of strongly non-linear systems. 2. In the incrementation process, the use of the arc-length method combined with the extrapolation technique can automatically trace complicated response curves with sharp peaks, and can considerably reduce the number of iterations required for the solution to converge. 3. The stability of steady state solutions can be analyzed by using the multi-variable Floquet theory. Moreover, it can be performed in association with the solution process at the same time in each augmentation. In other words, the computer program can be easily integrated to cover both the solution process and stability analysis in each augmentation step.

REFERENCES 1. A. H. NAYFEH and D. T. MOOK 1979 Nonlinear Oscikztions. New York: Wiley-Interscience. 2. A. H. NAYFEH 1983 Journal of Sound and Vibration 88,547-X7. The response of two-degree-offreedom systems with quadratic non-linearities to a parametric excitation. 3. A. H. NAYFEH and L. D. ZAVODNEY 1986 Journal of Sound and Vibration 107, 329-350. The response of two-degree-of-freedom systems with quadratic non-linearities to a combination parametric excitation.

INCREMENTAL

HARMONIC

BALANCE

METHOD

285

4. D. T. MOOK, R. H. PLAUT and N. HAQUANG 1985 Journal of Sound and Vibration 102, 473-492. The influence of an internal resonance on non-linear structural vibrations under subharmonic resonance conditions. 5. D. T. MoOK, N. HAQUANG and R. H. PLAUT 1986 Journalof Sound and Vibration 104, 229-241. The influence of an internal resonance on non-linear structural vibrations under combination resonance conditions. 6. R. H. PLAUT, N. HAQUANG and D. T. MOOK 1986 Journal of Sound and Vibration 107, 309-319. The influence of an internal resonance on non-linear structural vibrations under two-frequency excitation. 7. J. BAJKOWSKI and W. SZEMPLINSKA-STUPNICKA 1986 Journal of Sound and Vibration 104, 259-275. Internal resonances effects-simulation versus analytical methods results. 8. S. L. LAU and Y. K. CHEUNG 1981 American Society of Mechanical Engineers, Journal of Applied Mechanics 48, 959-964. Amplitude incremental variational principle for nonlinear vibration of elastic systems. 9. Y. K. CHEUNG and S. L. LAU 1982 Earthquake Engineering and Structural Dynamics 10, 239-253. Incremental time-space finite strip method for nonlinear structural vibrations. 10. S. L. LAU, Y. K. CHEUNG and S. Y. WU 1984 American Society of Mechanical Engrneers, Journal of Applied Mechanics 51, 837-844. Nonlinear vibration of thin elastic plates, part 1: generalized incremental Hamilton’s principle and element formulation. 11. S. L. LAU, Y. K. CHEUNG and S. Y. Wu 1984 American Society of Mechanical Engrneers, Journal of Applied Mechanics 51, 845-851. Nonlinear vibration of thin elastic plates, part 2: internal resonance by amplitude-incremental finite element. 12. V. P. Iu, Y. K. CHEUNG and S. L. LAUJ 1985 Journal ofSound and Vibration 100,359-372. Nonlinear vibration analysis of multilayer beams by incremental finite elements, part 1: theory and numerical formulation. 13. V. P. Iu, Y. K. CHEUNG and S. L. LAU 1985 Journal of Sound and Vibration 100,373-382. Nonlinear vibration analysis of multilayer beams by incremental finite elements, part II: damping and forced vibrations. 14. V. P. Iu and Y. K. CHEUNG 1986 Engineering Computation 3, 36-42. Non-linear vibration analysis of multilayer sandwich plates by incremental finite elements, part I: theoretical development. 15. V. P. IU and Y. K. CHEUNG 1986 Engineering Computation 3, 43-52. Non-linear vibration analysis of multilayer sandwich plates by incremental finite elements, part II: solution techniques and examples. 16. S. L. LAU, Y. K. CHEUNG and S. Y. WV 1983 American Society of Mechanical Engineers, Journal ofApplied Mechanics 50,871-876. Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems 17. S. L. LAU, Y. K. CHEUNG and S. Y. WU 1982 American Society of Mechanical Engineers, method for Journal of Applied Mechanics 49, 849-853. A variable parameter incrementation dynamic instability of linear and nonlinear elastic systems. 18. C. PIERRE and E. H. DOWELL 1985 American Society of Mechanical Engineers, Journal of Applied Mechanics 52, 693-697. A study of dynamic instability of plates by an extended incremental harmonic balance method. 19. C. PIERRE, A. A. FERRI and E. H. DOWELL 1985 American Society of Mechanical Engineers, Journal of Applied Mechanics 52, 958-964. Multi-harmonic analysis of dry friction damped systems using an incremental harmonic balance method. 20. R. BOUc and M. DEFILIPPI 1987 International Journal of Engineering Science 25, 723-733. A Galerkin multi-harmonic procedure for nonlinear multi-dimensional random vibration. 121.A. Y. T. LEUNG andT. C. FUNG 1989 JournalofSoundand Vibration 131,445-455. Construction of chaotic regions. 122. A. A. FERRI 1986 American Society of Mechanical Engineers, Journal of Applied Mechanics 53, 455-457. On the equivalence of the incremental harmonic balance method and the harmonic balance - Newton-Raphson method. 23. F. H. LING and X. X. WU 1987 International Journal of Non-linear Mechanics 22, 89-98. Fast Galerkin method and its application to determine periodic solutions of non-linear oscillators. 24. S. L. LAU 1989 Incremental harmonic balance method with time transformation for nonlinear vibrations (submitted) 25. C. S. HSU 1972 American Society of Mechanical Engineers, Journal of Applied Mechanics 39, 551-558. Impulsive parametric excitation: theory.

Y.

286

K.C‘HEUNG

ET

AL.

26. C. S. Hsu and W. H. CHENG 1973 American Society of Mechanical Engineers, Journal qfApplied of the theory of impulsive parametric excitation and new Mechanics 40, 78-86. Applications treatments of general parametric excitation problems. 27. C. S. Hsu 1974 Journal of Mathematical Analysis Applications 45, 234-251. On approximating a general linear periodic system. solution 28. M. A. CRISFIELD 1981 Computers and Structures 13,55-62. A fast incremental/iterative procedure that handles “snap-through”. C. E. HAMMOND and T. H. WOO 1977 International Journal of Numerical 29. P. FRIEDMANN, Methods in Engineering 11, 1117-l 136. Efficient numerical treatment of periodic systems with application to stability problems. 30. S. L. LAU, Y. K. CHEUNG and S. H. CHEN 1989 American Society of Mechanical Engineers, procedure of multiple Journal of Applied Mechanics 56, 667-675. An alternative perturbation scales for nonlinear dynamic systems. 31. S. H. CHEN, Y. K. CHEUNG and S. L. LAU 1989 Journal ofSound and Vibration 128, 13-24. On the internal resonance of multi-degree-of-freedom systems with cubic nonlinearity. APPENDIX Qo= Fo=

r910,420, . . . , qnolT, [fio,_&or .

Aq=[Aq,,Aq2,...,Aq,lT, ~F=W,,~h,...,4LlT,

. . .LolT,

Pi lip=

f k=l

i

H$o,q:qkOq,O,

Ol CT

s= I

A=MI,&,...,JLI~,

9

*.

I=1

0

-1 C*

AA=[AA,,AA,,..

.,AA,,lr, 2?r

M=

STtis

d r,

c=

S’c$

dq

K=

STk3

d 7,

I0

STgC3)S dq J0

l-zr

f2m

r2n

K(3)=

F=

ST& cos (2m - 1)~ dr, J0

R,=

STcos(2m-l)~d~. J0