Free vibration of structures with cubic non-linearity-remarks on amplitude equation and Rayleigh quotient

Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709 www.elsevier.com/locate/cma

Free vibration of structures with cubic non-linearity-remarks on amplitude equation and Rayleigh quotient R. Lewandowski

*

Institute of Structural Engineering, Poznan University of Technology, 5 Piotrowo Street, 60-965 Poznan, Poland Received 24 January 2002; received in revised form 23 July 2002

Abstract In the paper, several aspects concerning the formulation of and solution to amplitude equations for free vibration of systems with cubic non-linearity are discussed in detail. The use of several possible constraint equations to achieve the non-trivial solution to the amplitude equation is also presented. The well-known concept of the Rayleigh quotient is extended to the non-linear case of interest. Some general remarks on the formulation of the non-linear amplitude equation and the non-linear Rayleigh quotient are also presented. Moreover, the introduced Rayleigh quotient is used to develop a new and suitably efficient method to solve amplitude equations. The accuracy of solutions obtained by means of some widely used methods of analysis of free vibration problems is also presented. It is found that the Galerkin method gives the most accurate results.
2003 Elsevier Science B.V. All rights reserved.

1. Introduction An analysis of the problem of free vibration of non-linear oscillating systems leads to determination of frequencies and modes of vibration. Since these quantities are the essential characteristics of each vibrating system, a lot of attention is given to the problem of free vibration. In the case of linear vibrating systems, the exact solution to the motion equation is known and the problem of determination of frequencies and modes of vibration is reduced to finding a solution to a generalized eigenvalue problem. The generalized eigenvalue problem could be solved using many well documented and numerically efficient methods, like the subspace iteration method, the Lanczos method, the vector iteration methods and other ones. In these methods the Rayleigh quotient is often used. A detailed description of the above mentioned methods could be found in many monographs (see, for example [1–3]). The problems of free vibration of non-linear systems are much more complex. Beside a few exceptions, the exact solutions to motion equations are not known and only approximate solutions are available. Very

*

Tel.: +48-61-6652472. E-mail address: [email protected] (R. Lewandowski).

0045-7825/03/$ - see front matter
2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0045-7825(03)00189-0

1682

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

often in the case of non-linear continuous systems, the approximate solutions to motion equations are written as a sum of products of time and space functions. To give an example, for one-dimensional structures like beams these solutions can be written in the following form: wðx; tÞ ¼ vi ðxÞqi ðtÞ;

ð1Þ

where vi ðxÞ and qi ðtÞ are the functions of space x and time t, respectively. Moreover, i ¼ 1; 2; . . . ; n and the summation convention for repeated indices is used in (1) and throughout this work. Usually, the approximate solution in the time domain is assumed in the form of the truncated Fourier series and then qi ðtÞ ¼ cos ixt, where x denotes the non-linear frequency of vibration. Additionally, in many instances only the fundamental harmonic of the Fourier series is taken into account, i.e. n ¼ 1 in (1). The space functions vi ðxÞ are usually treated as unknown quantities or it is assumed they are a linear combination of trial functions. In the latter case the following can be written: vi ðxÞ ¼ aij uj ðxÞ;

ð2Þ

where j ¼ 1; 2; . . . ; m, and symbols uj ðxÞ and aij denote the trial function and the unknown factor, respectively. Quite frequently, the modes of vibration of linearised system are chosen as the trial functions. A similar description is used in the finite element method and then relation (2) refers to the typical finite element. Moreover, in many cases it is assumed that m ¼ 1. Using well-known methods, like the harmonic balance method, the Galerkin method or the Ritz method, it is possible to reduce the considered problem to finding a solution for a set of non-linear, ordinary differential equations (if the unknown quantities are vi ðxÞ and x) or for a system of non-linear algebraic equations (if aij and x are unknowns). Later in this paper, the above-mentioned system of non-linear algebraic equations will be referred to as the amplitude equations. From the mathematical point of view, these equations can be considered as a non-linear eigenvalue problem. In this paper the analysis of the second approach, reducing the problem to the amplitude equations, is of particular interest. Our considerations are restricted to non-linear dynamic systems with a cubic characteristic of restoring forces. In Section 2, the existing methods of derivation of amplitude equations are briefly summarized. The problem of normalisation of the vector of amplitudes is discussed in Section 3. In Section 4, several possible ways to extend the Rayleigh quotient to the non-linear systems with cubic characteristic are presented. Some remarks concerning the amplitude equations and the possible extensions of the Rayleigh quotient are formulated in Section 5. In Section 6, the application of the vector iteration methods and the Rayleigh quotient for solving the non-linear eigenvalue problems are presented. The results of exemplary calculations for beams are briefly described and discussed in Section 7. Concluding remarks are formulated in Section 8.

2. Derivation of the matrix amplitude equation Strings, beams, membranes and plates are examples of structures for which the geometrically non-linear theory must be used to correctly describe the responses of such structures under sufficiently large external forces. Moreover, the internal forces of these structures can be written as cubic functions of displacements. In the case of free vibration, the motion equations of undamped structures with cubic characteristic, treated as the discrete systems, can be written in the following matrix form:
 M€ wðtÞ þ K0 þ 12K2 ðwðtÞÞ wðtÞ ¼ 0; ð3Þ where M, K0 and K2 ðwðtÞÞ denote the global mass matrix, the global, linear stiffness matrix and the global, non-linear stiffness matrix, respectively. Moreover, wðtÞ is the global vector of nodal parameters and dots

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1683

means differentiation with respect to a time variable. It is useful to note that, on a finite element level, the non-linear stiffness matrix Ke2 ðwe ðtÞÞ can be written in the following form (see Refs. [4,5]): Z BT1 ðwe ÞEB1 ðwe Þ dV ; ð4Þ Ke2 ðwe Þ ¼ V

where E denotes the matrix of elasticity and V is the volume of the finite element. The matrix B1 ðwe Þ is the linear and homogenous function of nodal parameters. It means that the non-linear stiffness matrix Ke2 ðwe ðtÞÞ is a quadratic and homogenous function of nodal parameters. It is assumed that the displacements version of the finite element method is used to discretize the continuous systems. The large displacements and small rotations were assumed during the process of derivation of the matrix equation of motion (3). Derivation of Eq. (3) is given in detail in Refs. [4,5]. Moreover, the kinetic energy K and the strain energy U of the structure can be written in the following form (see [4,6]): K ¼ 12w_ T ðtÞMw_ ðtÞ;

ð5Þ

U ¼ 12wT ðtÞK0 wðtÞ þ 18wT ðtÞK2 ðwðtÞÞwðtÞ:

ð6Þ

In this paper, the beam structures are considered as an example of structures with cubic characteristic. The non-linear motion equation for beams with immovable ends can be written in the form (see [7]) M€ wðtÞ þ KwðtÞ þ

EA BwðtÞwT ðtÞBwðtÞ ¼ 0; 2l

ð7Þ

where now K0 ¼ K;

K2 ðwðtÞÞ ¼

EA BwðtÞwT ðtÞB l

ð8Þ

and E, A, l denote the YoungÕs constant, the area of beam cross-section and the beam length, respectively. On an element level the transverse displacements wðx; tÞ are written in the form wðx; tÞ ¼ NT ðxÞwe ðtÞ;

ð9Þ

where NðxÞ ¼ colðN1 ðxÞ; N2 ðxÞ; N3 ðxÞ; N4 ðxÞÞ is the vector of shape functions and we ðtÞ ¼ colðw1 ; u1 ; w2 ; u2 Þ is the vector of nodal parameters. The Hermite polynomials and two-node element are used in this paper. The definitions of M, K and B matrices on an element level are Z le Z le Z le Me ¼ mNðxÞNT ðxÞ dx; Ke ¼ EJ N;xx ðxÞNT;xx ðxÞ dx; Be ¼ N;x ðxÞNT;x ðxÞ dx; ð10Þ 0

0

0

where m is the mass of beam per unit length, J is the moment of inertia of the cross-section and ðÞ;x denotes differentiation with respect to the space variable. The axial force N ðtÞ and the strain energy U of the beam is given by (see [7]) Z EA l 2 EA T w ðtÞBwðtÞ; N ðwðx; tÞÞ ¼ w;x dx ¼ ð11Þ 2l 0 2l 2 1 EA  T w ðtÞBwðtÞ : U ¼ wT ðtÞKwðtÞ þ 2 8l

ð12Þ

The simplest method for derivation of the amplitude equation is proposed by Mei in [8]. Subsequently, this method was used in many papers devoted to non-linear dynamic problems of beams and plates [9–15]. The method assumes that the solution of motion equations fulfils the following conditions:

1684

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

€ max ðt1 Þ ¼ x2 wmax ðt1 Þ; w_ ðt1 Þ ¼ 0 ð13Þ w at the point of the reversal of motion (i.e. the maximum amplitude point at t ¼ t1 ). It is easy to see that the above conditions are fulfilled if only the fundamental harmonic is taken into account in the approximate solution of motion equations i.e. wðtÞ ¼ a cos xt:

ð14Þ

Introducing relation (14) into the equation of motion (3) or (7) and using conditions (13), the matrix amplitude equation is obtained in the form 
ð15aÞ K0  x2 M þ 12K2 ða; aÞ a ¼ 0; 

 EA T 2 Baa B  x M a ¼ 0 Kþ 2l

ð15bÞ

for a general case and for beams, respectively. The notation like K2 ða; aÞ underlines that this matrix is the quadratic function of amplitudes of nodal parameters. Conditions (13) can also be understood as collocation conditions. The collocation method is applied to solve the problem of forced vibration in [16] and general considerations concerning the collocation method can be found in [17]. Among very popular methods of derivation of amplitude equations, there are the harmonic balance method, the Galerkin method and the Ritz method. All of the above mentioned methods make it possible to analyze strongly non-linear systems. Moreover, as it is demonstrated in [4,6], all of these methods give us identical amplitude equations. The identical amplitude equation is also obtained if the so-called incremental harmonic balance method suggested by Lau and Cheung [18,19] is used. Also the method proposed in Refs. [20,21] results in the same amplitude equation. Sometimes, the perturbation methods are used to investigate the free or forced vibrations of non-linear systems with many degrees of freedom. A computational formulation of the perturbation method is given in [22]. However, application of these methods is restricted to weakly non-linear systems and for this reason they are out of the scope of this paper. Below, the matrix amplitude equation will be briefly derived with the help of the Galerkin method. An approximate solution to the motion equations both (3) and (7) is assumed in the following form: wðtÞ ¼ ai cos zi xt;

ð16Þ

where i ¼ 1; 2; . . . ; n and symbols ai and zi denote the unknown vector of amplitudes and the appropriately chosen natural number, respectively. The assumed solution is approximate and after introducing (16) into the motion equation (3) or (7) the residual vector rðtÞ is obtained. The Galerkin conditions have the following form: Z 2 T rðtÞ cos zl xt dt ¼ 0; ð17Þ T 0 where l ¼ 1; 2; . . . ; n and T ¼ 2p=x is the fundamental period of free vibration. From conditions (17) the following matrix amplitude equations can be derived ðK0  zi x2 MÞdil ai þ 12aijkl K2 ðai ; aj Þak ¼ 0;

ð18aÞ

EA aijkl Bai aTj Bak ¼ 0 2l

ð18bÞ

ðK  z2i x2 MÞdil ai þ

for the general system and beams, respectively. The coefficients dil and aijkl appearing in Eq. (18) are defined in the Appendix A.

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1685

After introducing notations Hkl ¼ 12aijkl K2 ðai ; aj Þ;

ð19aÞ

EA aijkl Bai aTj B; 2l

ð19bÞ

Hkl ¼

both Eqs. (19a) and (19b) can be rewritten in the identical form: ðKdlk  z2k x2 dlk M þ Hlk ðaÞÞak ¼ 0:

ð20Þ

Moreover, introducing subsequently the following notations a ¼ colða1 ; a2 ; . . . ; ak ; . . . ; an Þ;

ð21Þ

f M ¼ diag½z21 M; z22 M; . . . ; z2k M; . . . ; z2n M;

ð22Þ

e ¼ diag½K; K; . . . ; K; . . . ; K; K 2 H11 ðaÞ H12 ðaÞ    6 H21 ðaÞ H22 ðaÞ    6 6 .. .. 6 . . e H ðaÞ ¼ 6 6 Hl1 ðaÞ Hl2 ðaÞ    6 6 . .. 4 .. . Hn1 ðaÞ Hn2 ðaÞ   

ð23Þ H1k ðaÞ    H2k ðaÞ    .. . Hlk ðaÞ .. .



Hnk ðaÞ   

3 H1n ðaÞ H2n ðaÞ 7 7 .. 7 . 7 7; Hln ðaÞ 7 7 .. 7 . 5 Hnn ðaÞ

ð24Þ

the set of Eq. (20) can be rewritten in the form of a single matrix amplitude equation: e þH e ðaÞ  x2 f ðK M Þa ¼ 0:

ð25Þ

In many instances, only the fundamental harmonic is taken into account in the assumed solution to the e ¼ K, f M ¼ M, a ¼ a1 , a1111 ¼ 3=4, motion equation. In this case, z1 ¼ 1, n ¼ 1, K e ðaÞ ¼ H11 ðaÞ ¼ 3K2 ða; aÞ; H 8

ð26aÞ

e ðaÞ ¼ H11 ¼ 3EA BaaT B H 8l

ð26bÞ

and the amplitude equation takes the following form ðK0  x2 MÞa þ 38K2 ða; aÞa ¼ 0; ðK  x2 MÞa þ

3EA BaaT Ba ¼ 0 8l

ð27aÞ ð27bÞ

for the general case and beams, respectively. It is easy to observe that Eqs. (15) and (27) differ by the factor 3/4 in the non-linear term. If the dimension of wðtÞ vector is m, the assumed solution to the equation of motion contains n harmonics, then the dimension of a vector is mn and is equal to the number of the equations resulting from the Galerkin conditions. However, there are mn þ 1 unknowns in the amplitude equations because apart from the vector of amplitudes, the frequency of vibration x is also unknown. The additional condition must be added to Eq. (25). These conditions are discussed in the section that follows.

1686

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

3. Conditions of normalization of amplitude vector The amplitude equation (25) is treated as the matrix equation with a parameter, and frequency x is chosen as the main parameter. It is obvious that a trivial solution to the amplitude equation exists for all x 2 R. A non-trivial solution exist if eþH e ðaÞ  x2 f detð K M Þ ¼ 0: ð28Þ In the linear case, natural frequencies of vibration can be determined from condition (28) but in the nonlinear case the same is not possible because of the a vector appearing in (28). However, it is possible to use condition (28) as the additional equation to the amplitude equation (25). To date, this possibility has not been investigated. In literature, one can also find other proposals for the above-mentioned condition, which can be understood simply as a norm of the amplitude vector. According to the simplest and most popular condition, the freely chosen element of a vector (say r and usually the largest one) has a prescribed value, which can be written in the form ar ¼ a; ð29Þ where a denotes an assumed maximum norm of the amplitude vector. Condition (29) was used in [8–15]. In papers by Lewandowski [4,7], the additional equation is taken in the form a T a ¼ a2 ;

ð30Þ

where a is the assumed Euclidean norm of the amplitudes vector. In papers by Wellford et al. [20] and by Ghabrial et al. [21], the amplitude vector is normalised in such a way that the average kinetic energy Ka of the system of interest has a given value R. If the average kinetic energy is defined as Z T 2 KðtÞ dt; ð31Þ Ka ¼ T x2 0 then the above-mentioned condition can be written in the form 1 2 T z a Mai 2 i i

¼ R:

ð32Þ

In papers [20,21] this condition is used to determine vibration amplitudes when only one harmonic is taken into account in the solution for the motion equation. A different condition, based on the principle of energy conservation, has recently been proposed and used in a number of papers by Benamar and his co-workers [23–29]. This condition states that the maximum value of the strain energy Umax ðt1 Þ and the maximum value of kinetic energy Kmax ðt2 Þ are equal and the following can be written: Umax ðt1 Þ ¼ Kmax ðt2 Þ:

ð33Þ

After inserting the assumed form of solution (16) into relations (5) and (6) or (5) and (12), condition (33) can be written in the form x2 ¼

dij aTi K0 aj þ 14 lijkl aTi K2 ðaj ; ak Þal zi zj dij aTi Maj

ð34aÞ

in the general case and x2 ¼

dij aTi Kaj þ EA l aT Baj aTk Bal 4l ijkl i zi zj dij aTi Maj

for beams, where lijkl ¼ 1 for all i, j, k, l.

ð34bÞ

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1687

In the papers mentioned above, relation (34) is used to remove the natural frequency of vibration from the amplitudes equation (25). However, after such modification, the amplitude equation is still homogenous. For this reason, condition (29) is included in the set of problem equations and simultaneously the amplitude equation corresponding to amplitude r is removed from the matrix amplitude equation. Please note that condition (34) looks like the Rayleigh quotient extended to the non-linear case. This problem will be discussed in Section 4.

4. Non-linear Rayleigh quotient The classical Rayleigh quotient can be derived in a few different ways, which will be briefly discussed below. The amplitude equation for the linear system has the form ðK  x2 MÞa ¼ 0:

ð35Þ T

By pre-multiplying Eq. (35) by a we obtain the well-known Rayleigh quotient x2 ¼

aT Ka : aT Ma

ð36Þ

From the mathematical point of view, the Rayleigh quotient can be understood as the solution for the overdetermined set of equations with respect to x (here the matrix amplitude equation) if the approximation of the amplitude vector is known. The identical Rayleigh quotient is obtained from the principle of energy conservation, written in the form of relation (33) because for linear systems Kmax ¼ 12x2 aT Ma;

Umax ¼ 12aT Ka:

ð37Þ

The considered systems are autonomous and the principle of energy conservation must be satisfied in every time instance. Therefore, the integral of action Z 2 T J¼ ½U ðtÞ  KðtÞ dt ð38Þ T 0 can be understood as a functional which expresses the principle of energy conservation in a weak sense. Introducing the solution for motion equation of linear systems wðtÞ ¼ a cos xt

ð39Þ

into the integral of action (38), we obtain J ða; x2 Þ ¼ 12x2 aT Ma  12aT Ka: 2

ð40Þ 2

If J ða; x Þ is considered as a functional of a where x is a parameter, then from the stationary condition dJ ¼ 0 the amplitude equation (35) is obtained. Moreover, it is assumed that the value of functional J ða; x2 Þ in the point of stationarity is equal to zero (i.e. J ð~a; x2 Þ ¼ 0), which means that the condition of conservation of average energy is exactly fulfilled. Once again, from the condition J ð~a; x2 Þ ¼ 0, where ~a is the eigenvector, we obtain the Rayleigh quotient. In conclusion, we have three ways of derivation of the Rayleigh quotient for linear dynamic systems. Up to now, the Rayleigh quotient for non-linear systems has very rarely been considered. In Ref. [30] a special kind of the non-linear eigenvalue problem is considered and the corresponding Rayleigh quotient is introduced. Below, the possibilities of extending this quotient for non-linear dynamic systems with cubic characteristics will be discussed.

1688

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

By pre-multiplying the amplitude equation (18) by aTl the following quotient x2 ¼

dil aTl K0 ai þ 12 aijkl aTl K2 ðai ; aj Þak ; zi dil aTl Mai

ð41aÞ

x2 ¼

dil aTl Kai þ EA a aT Bai aTj Bak 2l ijkl l z2i dil aTl Mai

ð41bÞ

is obtained. If only one harmonic is included in the assumed solution for the motion equation, then x2 ¼

aT K0 a þ 38 aT K2 ða; aÞa ; aT Ma

x2 ¼

aT Ka þ 3EA aT BaaT Ba 8l aT Ma

ð42Þ

for the general case and beams, respectively. Following, in the same way, with the amplitude equation (15), resulting from the collocation condition, we obtain x2 ¼

aT K0 a þ 12 aT K2 ða; aÞa ; aT Ma

x2 ¼

aT Ka þ EA aT BaaT Ba 2l aT Ma

ð43Þ

for the general case and for beams, respectively. Quotients (41), (42) and (43) will be called the non-linear Rayleigh quotients of first kind. Possibilities resulting from the condition of energy conservation will now be analyzed in detail. Please note that Eq. (34) resulting from the condition of maximal energies (33) can be understood as the second possible definition of the non-linear Rayleigh quotient. It is easy to see that this quotient differs from the one given by (41) because, in general, lijkl 6¼ aijkl . Therefore, quotients resulting from energy conditions will be called the non-linear Rayleigh quotients of second kind. If one harmonic solution for the motion equation is of interest then n ¼ 1, a1 ¼ a, z1 ¼ 1 and from (34) we obtain the following quotients x2 ¼

aT K0 a þ 14 aT K2 ða; aÞa ; aT Ma

x2 ¼

aT Ka þ EA aT BaaT Ba 4l aT Ma

ð44Þ

which differ both from quotient (42) and (43). Furthermore, the non-linear Rayleigh quotient can be defined on the basis of the principle of energy conservation in an average sense. This principle is written in the form Ua ¼ Ka ;

ð45Þ

where the symbols Ka and Ua denote the average kinetic and the average strain energy, respectively. These energies are defined as Z Z 2 T 2 T KðtÞ dt; Ua ¼ U ðtÞ dt: ð46Þ Ka ¼ T 0 T 0 Please note that the definition of the average kinetic energy given above differs by the factor 1=x2 from the definition given by relation (31). It is obvious that the considered system is autonomous and that the principle of energy conservation can be satisfied in every time instance t. However, we must take into account that the exact solution to the motion equation is not known. Neither the principle of energy conservation nor the equation of motion are fulfilled by the approximate solution in every time instance t. It seems to be reasonable that for the approximate solutions, the principle of energy conservation must be satisfied at least in an average sense.

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1689

After introducing the solution to the motion equation in the form of (16) into relations (46), (5), (6) or (12) we obtain Ka ¼ 12zi zj dij x2 aTi Maj ; ð47Þ Ua ¼ 12dij aTi K0 aj þ 18aijkl aTi K2 ðaj ; ak Þal ; 1 EA aijkl aTi Baj aTk Bal : Ua ¼ dij aTi Kaj þ 2 8l From condition (45), the following non-linear Rayleigh quotients dij aTi K0 aj þ 14 aijkl aTi K2 ðaj ; ak Þal x2 ¼ ; dij aTi Maj x2 ¼

dij aTi Kaj þ EA a aT Baj aTk Bal 4l ijkl i zi zj dij aTi Maj

ð48aÞ ð48bÞ

ð49aÞ

ð49bÞ

are obtained. If the solution to the motion equation contains only one harmonic, then relation (49b) can be rewritten as x2 ¼

aT Ka þ 3EA aT BaaT Ba 16l : aT Ma

ð50Þ

Once again, we obtain another formula for the non-linear Rayleigh quotient. In the linear case, the Rayleigh quotient can be derived from the condition J ða; x2 Þ ¼ 0. Now, we will analyse the consequences of this condition in a non-linear case. Our consideration is restricted to the one harmonic approximate solution for the motion equation, i.e. when wðtÞ ¼ a cos xt. Two cases are considered. In the first case, the functional J ða; x2 Þ results from the principle of energy conservation in an average sense i.e. J ða; x2 Þ  Ja ða; x2 Þ ¼ Ua ðaÞ  Ka ða; x2 Þ;

ð51Þ

while in the second case the functional is defined as J ða; x2 Þ  Jm ða; x2 Þ ¼ Umax ðaÞ  Kmax ða; x2 Þ:

ð52Þ

Inserting Eq. (14) into (5) and (12) and using relations (51) and (52) we can write 1 EA T 1 ja BaaT Ba  x2 aT Ma; J ða; x2 Þ ¼ aT Ka þ 2 4l 2

ð53Þ

where j ¼ 3=4 in the first case and j ¼ 1 in the second. Now a is the trial vector and x2 is treated as a parameter. The first variation of functional J ða; x2 Þ is dJ ða; x2 Þ ¼ daT

oJ : oa

ð54Þ

At a stationary point of J ða; x2 Þ the variation dJ ða; x2 Þ must be eliminated, which yields the following stationary condition Ka þ

EA T ja BaaT Ba  x2 Ma ¼ 0 2l

ð55Þ

which agrees with the previously derived amplitudes equation (27b) (first case) or (15b) (second case). If we additionally introduce the condition that at the point of stationarity J ða; x2 Þ ¼ 0 then the parameter x2 can be calculated from the following relation

1690

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

x2 ¼

aT Ka þ EA aT BaaT Ba 2l aT Ma

ð56Þ

which for j ¼ 3=4 agrees with (42b) and for j ¼ 1 is identical with (43b).

5. General remarks on formulation of amplitude equation and non-linear Rayleigh quotient The following conclusions result from the above considerations: • The amplitudes equation (15b) resulting from the collocation condition is simultaneously the condition of stationarity of the functional Jm ða; x2 Þ while the amplitude equation (27b) derived with the help of the Galerkin method is also the stationary condition of the functional Ja ða; x2 Þ resulting from the principle of conservation energy in an average sense. • We have two possible ways of extending the linear Rayleigh quotient concept in a non-linear case. The first way, leading to the non-linear Rayleigh quotient of first kind, can be understood as the condition of fulfilling the amplitude equation with respect to x2 if some approximation of the a vector is available. The proof of this statement is described in Appendix B. In the second way, the frequency of vibrations is calculated from the principle of energy conservation (in the averaging sense or in the sense of conservation of maximal energies). • From the mathematical point of view, the considered problem is fully defined by the amplitude equations and the condition of existence of non-trivial solutions (28). Condition (28) is hard to use in computing procedures, that is why different conditions, such as (29) and (30) or (32), are used in practice. • Definitions of the non-linear Rayleigh quotient of first and second kind, given by relations (41a) and (34a), differ from each other in the non-linear term. Both have clear mathematical or physical interpretations. The Rayleigh quotient of second kind is an additional physical condition which makes the considered problem overdetermined because now we have n þ 2 equations (i.e. n amplitude equations, the condition of existence of a non-trivial solution and the condition of energy conservation) with only n þ 1 unknowns. The Rayleigh quotient of first kind can be interpreted as an additional condition from which the best approximation of frequency of vibration can be determined for a given approximation of non-linear eigenvector a. This condition does not introduce any additional equation into problem formulation and is only an auxiliary equation which can be used in the procedure for solving the matrix amplitude equation. • The principle of energy conservation is an essential that should be fulfilled for every autonomous dynamic system. Differences in the proposed non-linear Rayleigh quotients suggest that the approximate solution for the motion equation (16) do not fulfil the principle of energy conservation. The value of energy error will be shown later for typical beam structures. • In this paper, two formulations of amplitude equation are discussed. The first formulation uses the collocation condition while the second formulation utilises some averaging procedure. Both formulations must be internally consistent. This will be achieved if the functional, the amplitude equation and the non-linear Rayleigh quotient are used in the forms shown in Table 1. Other formulations, including the formulation proposed in papers [23–29] are internally inconsistent. Obviously, this conclusion applies also to approximate solutions with many harmonics. • Both internally consistent formulations are approximate, which means that all solutions are burdened with errors. The rational choice of formulation requires an analysis of the accuracy of the results obtained. It is very difficult in the general case, therefore, only numerical results of accuracy analysis for beam structures will be presented and briefly discussed in this paper.

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1691

Table 1 Summary of the Galerkin and collocation formulation Type of formulation

Collocation formulation

Galerkin formulation

Type of functional Type of amplitude equation

J ða; x2 Þ ¼ Umax ðaÞ  Kmax ða; x2 Þ Resulting from the condition € max ðt1 Þ ¼ x2 wmax ðt1 Þ w Formula (43) Resulting from condition Umax ðaÞ ¼ Kmax ða; x2 Þ

J ða; x2 Þ ¼ Ua ðaÞ  Ka ða; x2 Þ Resulting from the conditions R 2 T rðtÞ cos zl xt dt ¼ 0 T 0 Formula (41) Resulting from condition Ua ðaÞ ¼ Ka ða; x2 Þ

Rayleigh quotient of first kind Rayleigh quotient of second kind

6. Vectors iteration method In this section the method to solve amplitude equations using the non-linear Rayleigh quotient is described. The method will be referred to as the Ôvectors iteration methodÕ because of similarities between the algorithm of the proposed method and the well-known method to solve linear eigenvalue problems. As it was shown above, analysis of non-linear free vibrations problem requires solution of the matrix amplitude equation which can be rewritten here in the following form: Fða; aÞ ¼ kMa;

ð57Þ

where k ¼ x2 and Fða; aÞ is the vector whose elements are non-linear functions of the elements of a and the parameter a (a is the amplitude of vibration at the point p). It is assumed that the following conditions are fulfilled: Fð0; aÞ ¼ 0;

kak1 ¼ max aj ¼ ap ¼ 1; j

Gða; aÞ ¼

oF oa

ð58Þ

and the Gða; aÞ matrix is positive definiteness. Moreover, Fða; lÞ and Gða; lÞ are continuous functions of a and a. The solution to the amplitude equation presented in this section has been obtained by means of the vectors iteration method. It is a well-known method to solve linear eigenproblems [1], and is extended here for the non-linear case. The methods for determining the fundamental and arbitrarily chosen eigenvalues, as well as the results of sample calculation with information concerning the convergence of the proposed method are discussed successively. All the procedures given below concern the case a ¼ const. but it is possible to find the backbone curve kðaÞ by repeating the calculations for a set of a values. 6.1. Determining the fundamental eigenvalue and eigenvector The algorithm of the method has an iterative character. It is assumed that certain approximation for kðaÞ and a is known and denoted as ki and ai (i means the number of iteration). The good first approximation is a non-linear eigenpair obtained for the previous a value or linear ones obtained by solving the linear eigenproblem as follows: ðGð0; aÞ  kMÞa ¼ 0:

ð59Þ

Now, for the given ki and ai we can solve the following non-linear algebraic system, for instance, by means of the Newton–Raphson method: Fðciþ1 ; aÞ ¼ ki Mai ; where ciþ1 is now the unknown vector.

ð60Þ

1692

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

As the next approximation of the eigenvector a, we take over the vector proportional to ciþ1 and whose norm is kaiþ1 k1 ¼ maxj aj;iþ1 ¼ ap;iþ1 ¼ 1. A new approximation of k is obtained from kai k1 kiþ1 ¼ ki ðfirst wayÞ; ð61aÞ kciþ1 k1 or from the relation kiþ1 ¼

aTiþ1 Fðaiþ1 ; aÞ aTiþ1 MaTiþ1

ðsecond wayÞ:

ð61bÞ

In the second way, the above-defined non-linear Rayleigh quotient is used. The iterative process is finished if the following inequalities are satisfied: jkiþ1  ki j 6 e1 ; kiþ1

kaiþ1  ai k2 6 e2 ; kaiþ1 k2

ð62Þ

where k  k2 denotes the Euclidean norm and e1 and e2 are the assumed accuracies of calculation. As it is easy to observe, the above methods are very similar to the well-known inverse iteration and the Rayleigh quotient method (compare [1]), respectively. 6.2. Determining the largest eigenvalue and eigenvector It is possible to write Eq. (57) in the form Ma ¼ hFða; aÞ;

ð63Þ

where h ¼ 1=k. For certain values hi and ai , the following equation Mciþ1 ¼ hi Fðai ; aÞ

ð64Þ

is solved. The new approximation for hiþ1 is obtained from hiþ1 ¼ hi

kai k1 kciþ1 k1

ðfirst wayÞ

ð65aÞ

or hiþ1 ¼

aTiþ1 MaTiþ1 aTiþ1 Fðaiþ1 ; aÞ

ðsecond wayÞ;

ð65bÞ

where we take over ciþ1 (after its normalisation) as aiþ1 . The iteration process is finished if jhiþ1  hi j 6 e1 ; hiþ1

kaiþ1  ai k2 6 e2 : kaiþ1 k2

ð66Þ

The above algorithm required only one factorisation of matrix M and it is very similar to the forward iteration method (see [1]). 6.3. Algorithm for determining the arbitrarily chosen eigenpair First, the shift u for Fða; aÞ is performed by introducing Rða; aÞ ¼ Fða; aÞ  uMa

ð67Þ

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1693

and then the following eigenproblem is considered: Rða; aÞ ¼ gMa;

ð68Þ

where g ¼ k  u. If we take a suitable u, then g (which is k decreased by u) is the fundamental eigenvalue of non-linear eigenproblem (68), which now can be solved using the method given in Section 6.1. The searched eigenvalue k is obtained from k¼gþu

ð69Þ

and the corresponding eigenvector is the same as the one of the eigenproblem (68). However, from a wide experience with the vector iteration method it results that the method is sensitive to round-off errors like the well-known version used to solve the linear eigenproblem (compare Ref. [1]). Depending on the magnitude of these errors, the actual value of shift and the starting amplitude vector, the iterative process can converge to the desired non-linear eigenpair or to the fundamental frequency and mode of vibration. This is an important disadvantage of the method and a natural question must be how to improve the convergence rate in the vector iterations. Please notice, that the application of the Gram– Schmidt orthogonalization process and the concept of the matrix deflation are not obvious because of the problem non-linearity. In the linear eigenvalue problem, the vector iteration converges to the desired eigenpair if the shift is reasonably close to the eigenvalue of interest. In a non-linear case, the frequency of vibration depends on amplitude of vibration and during the process of determination of the backbone curve the distance from the shift and the frequency of interest can growth due to change of frequency. If this distance is large enough the iteration process fails. For the above reasons and on a base of wide numerical experiences it is proposed to calculate the shift ui (i.e. the shift in the ith iteration) from ui ¼ cki1 :

ð70Þ

A value of c parameter must equals approximately 1 to ensure that the shift is close to current approximation of searched eigenvalue. It is especially important if the strongly non-linear system is analysed or if the large increment of amplitude of vibration a is used to determine the backbone curve. c ¼ 0:98 is chosen in our own calculation. 6.4. Discussion of convergence properties of iteration processes Taking into account that ci ai ¼ ¼ fi ðci Þci : kci k1

ð71Þ

Eq. (64) can be written as ciþ1 ¼ hi ðci ÞM1 Fðfðci Þci ; aÞ ¼ zðci Þ:

ð72Þ

Likewise, by substituting (71) into (60) the following is obtained Fðciþ1 ; aÞ ¼ fi ðci Þki ðci Þci

ð73Þ

ciþ1 ¼ ci þ G1 ðci Þðfi ðci Þki ðci Þci  Fðci ; aÞÞ ¼ zðci Þ

ð74Þ

or

after linearization of the left-hand side of (73).

1694

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

Eqs. (72) and (73) have a form of the successive iteration scheme, which is convergent if (see, for example [3]) kzðc1i Þ  zðci Þk ¼ kci  ciþ1 k 6 bkci1  ci k;

ð75Þ

where 0 6 b 6 1 and the error of the solution is estimated by kci  ak 6

bi kc2  c1 k: 1b

ð76Þ

6.5. Numerical examples 6.5.1. Mass-spring system The free vibration problem of the dynamic system pictured in Fig. 1 is considered. The motion equation of the system can be written in the form € r  k1 ðwr1  2wr þ wrþ1 Þ þ k3 ðwr  wr1 Þ3  k3 ðwrþ1  wr Þ3 ¼ 0; Mw

ð77Þ

where r ¼ 1; 2; . . . ; n and w0 ¼ wnþ1 ¼ 0. Assuming that the system oscillates periodically with one mode, the following can be written: wr ¼ aar cos xt:

ð78Þ

By substituting (78) into (77) and applying the Galerkin method, the following amplitude equation written in the form of a difference equation is obtained: 3a2 k3 3 3 ½ðar  ar1 Þ  ðarþ1  ar Þ  ¼ 0; ð79Þ 4 where k ¼ x2 and again r ¼ 1; 2; . . . ; n. After defining the non-zero entries of a three-diagonal matrix HðaÞ as kMar  k1 ðar1  2ar þ arþ1 Þ þ

2

2

hrr ¼ 2k3 a2 ½ðar  ar1 Þ þ ðarþ1  ar Þ ; hr;rþ1 ¼ 2k3 a2 ðarþ1  ar Þ

2

hr;r1 ¼ 2k3 a2 ðar  ar1 Þ ;

2

ð80Þ

and introducing the following notation a ¼ colða1 ; a2 ; . . . ; ar ; . . . ; an Þ;

M ¼ diag½M; M; . . . ; M;

K0 ¼

;

the amplitude equation (79) can be written in the form of Eq. (27a).

Fig. 1. The mass-spring system.

ð81Þ

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1695

Table 2 The non-linear frequencies of the mass-spring system a (m)

1st eigenvalue p ¼ 6

2nd eigenvalue p ¼ 3

3rd eigenvalue p ¼ 2

10th eigenvalue p ¼ 5

linear 0.03 0.06 0.09 0.12 0.15

8.10141 8.11834 8.16889 8.25234 8.36760 8.51339

31.7493 32.0084 32.7706 33.9991 35.6512 37.6895

69.0279 70.2435 73.7346 79.1849 86.3428 95.0997

391.899 426.368 518.077 674.722 897.622 1186.340

p means that kak1 ¼ ap ¼ 1.

Table 3 Statistics of the vector iteration procedure––the mass-spring system Mode

a (m)

Da (m)

Number of iterations for a

Number of solutions to linear algebraic equations

1

0.03 0.15

0.03 0.03

4 (1), 3 (2) 11 (1), 3 (2)

5 (1), 4 (2) 18 (1), 6 (2)

2

0.03 0.09 0.09 0.15

0.03 0.03 0.09 0.03

3 3 4 4

5 6 8 9

10

0.03 0.15 0.15

0.03 0.03 0.15

35 (1,2) 9 (1,2) 14 (2)

(2) (2) (2) (2)

(2) (2) (2) (2)

35 (1,2) 9 (1,2) 14 (2)

The number in brackets means the way of calculation of kiþ1 . Moreover, Da is the applied increment of a.

The backbone curves of the system of interest are calculated and the following data are adopted: n ¼ 10, M ¼ 10:0 kg, k1 ¼ 1000:0 N/m, k3 ¼ 50000:0 N/m3 , e1 ¼ e2 ¼ 0:0001. The backbone curves for x1 , x2 , x3 and x10 are determined. As the first approximation for k ¼ x2 and ar , the appropriate solution of the linear eigenproblem was chosen. The results of calculation and information concerning the convergence of the method are given in Tables 2 and 3. 6.5.2. Multispan beam The three span beam shown in Fig. 2 is considered. The beam is divided into twenty identical, two node finite elements. The one harmonic solution to the motion equation given by Eq. (14) is assumed in this case. The first and second backbone curves are calculated using the above described vector iteration method. Calculation results are shown in Fig. 3. In the figure the non-dimensional, first and second frequencies of vibration x=xl versus the non-dimensional amplitude a=i in the middle of the first span are presented. Here i denotes the radius of inertia and xl is the linear frequency of vibration (first or second, respectively). The solid curve is the fundamental backbone curve while the second backbone curve is shown as the dashed line. Information concerning the convergence of the method is given in Table 4. In both considered cases and for different increments of non-dimensional amplitudes Da=i the iterative process converges very fast. 6.5.3. Fully clamped square plate The main aim of this example is to illustrate the power of the proposed vector iteration method when applied to the solution of the two-dimensional in space problem. However, some observation concerning different formulations is also made.

1696

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

Fig. 2. The three span beam.

Fig. 3. The first (––) and second (- - -) non-dimensional backbone curves of the three span beam.

Table 4 Statistics of the vector iteration procedure––the three span beam Mode

Increment of non-dimensional amplitude Da=i

Total number of increments

Average number of changes of frequency per one amplitude increment

Average number of Newton– Raphson iterations per one change of frequency

1

0.05 0.10 0.15 0.20

101 51 34 25

1.00 1.16 1.38 1.48

1.05 1.41 1.43 1.57

2

0.05 0.10 0.15 0.20

101 51 34 25

1.00 1.04 1.32 1.40

1.14 1.62 1.71 1.89

The fundamental backbone curve for the fully clamped square isotropic plate is determined by means of the vector iteration method. The formulation presented in Ref. [31] is used to describe the plate behaviour. In this formulation the contribution of the in-plane displacements to the non-linear strain energy is neglected because the terms involving the in-plane displacements are considered insignificant compared with the terms involving only out-of-plane displacements. Moreover, the axial and rotary inertia forces are neglected. Taking into account the above assumptions, the strain energy U and the kinetic energy K of the plate can be written as Z m K¼ ð82Þ w_ 2 dS; 2 S

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

U¼

D 2

Z

2

½ðw;xx þ w;yy Þ þ 2ð1  mÞw2;xy  w;xx w;yy dS þ S

3D 2h2

Z

2

ðw2;x þ w2;y Þ dS;

1697

ð83Þ

S

where ðÞ;x ¼ d=dx, D is the bending stiffness of the plate, h is the plate thickness, m is the Poisson ratio and wðx; yÞ denotes the plate transverse displacement. It is assumed that the periodic free vibrations of the plate can be described by wðx; yÞ ¼ ai wi ðx; yÞ cos xt;

ð84Þ

where ai ði ¼ 1; 2; . . . ; nÞ are the unknown amplitudes of global shape functions wi ðx; yÞ. Substituting relations (82)–(84) into (38) leads to the following functional: J ¼ 12kij ai aj þ 163 hijkl ai aj ak al  12x2 mij ai aj ; where mij ¼ m

Z

ð85Þ

wi ðx; yÞwj ðx; yÞ dS;

ð86Þ

½ðwi;xx þ wi;yy Þðwj;xx þ wj;yy Þ þ 2ð1  mÞwi;xy wj;xy  wi;xx wj;yy  dS;

ð87Þ

S

kij ¼ D

Z S

hijkl

6D ¼ 2 h

Z

ðwi;x wj;x þ wi;y wj;y Þðwk;x wl;x þ wk;y wl;y Þ dS:

ð88Þ

S

The non-linear Rayleigh quotient of first kind is given by x2 ¼

kij ai aj þ 38 hijkl ai aj ak al : mij ai aj

ð89Þ

The first variation of J ðai ; x2 Þ leads us to the following amplitude equation kij aj þ 34hijkl aj ak al  x2 mij aj ¼ 0:

ð90Þ

Each global shape function wi ðx; yÞ is chosen as a product of eigenfunctions ur ðxÞ and us ðxÞ of the clamped–clamped beam i.e. wi ðx; yÞ ¼ ur ðxÞus ðyÞ:

ð91Þ

These functions satisfy all plate boundary conditions. To obtain the fundamental non-linear mode shape twenty-five global shape functions obtained as products of the first five clamped–clamped beam modes are used. The discussion of convergence of the above spectral expansion, given by relations (84) and (91), is presented in Ref. [31] for both the linear and non-linear case. The results of our calculation combined with findings taken from literature are presented in Table 5. In the table the fundamental non-dimensional non-linear frequency x=xlin versus the non-dimensional amplitudes of vibrations in the middle of the plate w=h is presented. It is obvious that there are some differences between the compared results. In Ref. [32] the influence of inplane displacements is taken into account while it is neglected in our formulation. Therefore, this reason the axial strains in the formulation used here are smaller in comparison with the formulation presented by Lau et al., our model is stiffer and in consequence the non-linear frequency is higher compared with results taken from Ref. [32]. These differences are not negligible, which means that the effects of the in-plane displacement cannot be ignored. The differences between the results taken from Ref. [31] and presented here have a different source. Now, the formulation is identical, the identical shape functions are used and the identical amplitude equations are obtained. The difference is only that the additional equation, resulting from the

1698

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

Table 5 Comparison of results for the square clamped–clamped plate x=xlin

w=h 0.2 0.4 0.6 0.8 1.0 1.5

w=h

Results from Ref. [31]

Results from Ref. [32]

1.0070 1.0276 1.0607 1.1044 1.1573 1.3499

1.0073 1.0291 1.0648 1.1138 1.1762

x=xlin Present results

0.1951 0.3892 0.5815 0.7718 0.9602 1.4250

1.0100 1.0390 1.0843 1.1425 1.2106 1.4082

principle of energy conservation, is added to the amplitude equations in Ref. [31]. The consequences of this change will be discussed in the next section. The rapid convergence of the vector iteration procedure is also observed in this case. In the incrementaliterative procedure for the non-dimensional amplitude increment Da=h ¼ 0:1 only three iterations for the non-linear frequency and five iterations in the Newton–Raphson method are sufficient to obtain the convergent solution for a=h ¼ 2:0. Here a denotes the amplitude of first global shape function in the middle of the plate. These results have shown that the proposed method is sufficiently effective. No modification of the standard finite element procedure except one, generating the non-linear matrix HðaÞ, is necessary. The method is self-starting because it may be used to solve the linear eigenproblem too. In the case of strongly non-linear solutions, the procedure is very fast convergent, especially when calculating the fundamental frequency and mode of vibration. 7. Discussion of results of exemplary calculations In this paper, conclusions concerning the accuracy of approximate solutions to equations of motion are drawn from the values of the energy error introduced below. This definition follows from the integral of action Z 2 T L¼ ½KðtÞ  U ðtÞ þ W ðtÞ dt; ð92Þ T 0 where W ðtÞ is the potential of external loads. If a system is autonomous, then W ðtÞ ¼ 0 and L ¼ Ka  Ua . From the principle of energy conservation it follows that Ka ¼ Ua and then L ¼ 0. After introducing the approximate solution into the equation of motion the vector of residuals is obtained, which is denoted by rðtÞ. The rðtÞ vector is treated as the vector of fictitious forces that must load the system if its oscillations are exactly described by the above-mentioned approximate solution to the motion equation. The potential of the fictitious forces is given by W ðtÞ ¼ wT ðtÞrðtÞ ð93Þ and the average potential Wa is Z 2 T T w ðtÞrðtÞ dt: ð94Þ Wa ¼  T 0 Taking into account the above observations, the following non-dimensional energy error DEa ¼ is proposed.

La Ka  Ua þ Wa ¼ Ea Ka þ Ua

ð95Þ

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1699

If the approximate solution of motion equation has the form wðtÞ ¼ al cos zl xt then Wa ¼ aTl sl ;

ð96Þ

where l ¼ 1; 2; . . . ; n, and Z 2 T sl ¼ rðtÞ cos zl xt dt: T 0

ð97Þ

If the Galerkin method is used, then in all cases Wa ¼ 0. However, Wa 6¼ 0 if the method proposed in papers [23–29] is used because in this case one of the equations from the matrix amplitude equation is not fulfilled. If it is the equation number j associated with the harmonic l, then Wa ¼ alj slj ;

ð98Þ

where the symbols alj and slj denote the jth entries of al and sl vectors, respectively. If the amplitude equation is derived with the help of the collocation method, then sl 6¼ 0. For beam structures and for n ¼ 1, we can write a1 ¼ a, s1 ¼ s and Wa ¼ aT s where s ¼ ðK  x2 MÞa þ

3EA BaaT Ba: 8l

ð99Þ

7.1. Simply supported beam with immovable ends The motion equation of beams treated as a continuous system has the form: m€ wðx; tÞ þ EJw;xxxx ðx; tÞ  N ðwðx; tÞÞw;xx ðx; tÞ ¼ 0;

ð100Þ

where symbols m, EJ and w denote the mass per unit length of beam, the beam rigidity and the transverse displacements of beam, respectively. Moreover, Z EA l 2 N ðwðx; tÞÞ ¼ w dx; ð101Þ 2l 0 ;x where l is the beam length and EA is the axial rigidity of the beam. If the solution to the motion equation is assumed in the form px wðx; tÞ ¼ a sin cos xt; l then the amplitude equation is " #  2 x j  a 2 1 þ a ¼ 0; xl 4 i

ð102Þ

ð103Þ

where xl is the fundamental, linear frequency of vibration, and i denotes the radius of inertia of the crosssection. Furthermore, j ¼ 3=4 if the amplitude equation results from the Galerkin condition and j ¼ 1 if the amplitude equation follows from the collocation condition. It is easy to observe that from condition (28) and from the Rayleigh quotient of first kind, an identical backbone curve equation is obtained in the form  2 x j  a 2 ¼1þ : ð104Þ xl 4 i Moreover, the frequency of vibration determined with the help of the collocation method differs from the frequency of vibration resulting from the Galerkin method.

1700

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

The energies defined in the previous section can be written in the following form: Ka ¼ Kmax ¼ 14ma2 lx2 ;

ð105Þ

1  p 4 1  p 4 EJa2 l þ EAa4 l; 4 l 32 l 1  p 4 3  p 4 EJa2 l þ EAa4 l: Ua ¼ 4 l 128 l

Umax ¼

The average potential of fictitious external forces is given by 1 1  p 4 3  p 4 Wa ¼ ma2 lx2  EJa2 l  EAa4 l: 2 2 l 32 l

ð106Þ ð107Þ

ð108Þ

In the method proposed by Benamar et al. [24] the frequency of vibration is determined from the condition Kmax ¼ Umax . Using this condition, the following is obtained:  2 x 1  a 2 ¼1þ : ð109Þ xl 8 i The exact backbone curve is described by (see [33])  2      x p 2 1  a 2 p ¼ 1þ F 2 k2; ; xl 2 4 i 2 where 1  a 2 k ¼ 8 i 2



1  a 2 1þ 4 i

ð110Þ

 ð111Þ

and F ðÞ is the elliptic integral of first kind. In Fig. 4 the ratios of frequency obtained by the approximate method and the exact frequency 2 ðxapprox =xexact Þ are shown versus the non-dimensional amplitude ratio (a=i). In addition, Fig. 5 presents the backbone curves describes by Eqs. (104), (109) and (110). In both figures results obtained from the Galerkin

Fig. 4. The ratios of the approximate frequency and the exact frequency versus the non-dimensional amplitude a=i for the simply supported beam.

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1701

Fig. 5. Fundamental backbone curves for the simply supported beam.

method are shown as a solid line, the results of the collocation method are presented as the solid line with small crosses whereas the results of the Benamar method are pictured as the dashed line. The solid line with triangles in Fig. 5 is the exact backbone curve. It is easy to observe that the frequency of vibration determined with the help of the Galerkin method is burdened with the smallest errors. Errors of the two other methods are considerably larger. Using the collocation method, the frequency of vibration is determined with excess while the Benamar method gives an underestimated frequency of vibration. The energy error of the approximate method will now be considered. Introducing Eqs. (105), (107) and (108) into (95), the non-dimensional average energy error of the Galerkin method can be written in the following form:   3  a 2 9  a 2 DEa ¼  2þ : ð112Þ 32 i 32 i Likewise, the following is obtained:   9  a 2 11  a 2 DEa ¼ 2þ ð113Þ 32 i 32 i for the collocation method and   3  a 2 7  a 2 DEa ¼  2þ 32 i 32 i

ð114Þ

for the Benamar method. From relations (112) and (113) it follows that the average energy error of the solution obtained by the Galerkin method and the Benamar method is negative while the energy error of the collocation method is positive. In Fig. 6 the modulus of average energy error is shown versus the non-dimensional amplitude ratio (a=i). Once again the energy error of the Galerkin method is the smallest. 7.2. Fixed–fixed beam––comparison of results obtained by Galerkin and Benamar methods It is assumed, as in [24], that the solution to the motion equation has the form wðx; tÞ ¼ ða1 v1 ðxÞ þ a2 v2 ðxÞÞ cos xt;

ð115Þ

where v1 ðxÞ and v2 ðxÞ denote the first and second symmetrical, linear mode of vibration, respectively.

1702

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

Fig. 6. The energy error of approximate solutions versus the non-dimensional amplitude a=i.

Using the Galerkin method with respect to space and time variables, the following amplitude equations are obtained: 

  3EA  c11 a21 þ 2c12 a1 a2 þ c22 a22 ðc11 a1 þ c12 a2 Þ ¼ 0; k11  x2 m11 a1 þ 8    3EA  2 c11 a21 þ 2c12 a1 a2 þ c22 a22 ðc12 a1 þ c22 a2 Þ ¼ 0; k22  x m22 a2 þ 8

ð116Þ

where mii ¼

m l

Z

l

v2i ðxÞ dx; 0

kii ¼

EJ l

Z

l

v2i;xx ðxÞ dx; 0

cij ¼

1 l

Z

l

vi;x vj;x ðxÞ dx:

In this case, the maximal and average energies can be determined from the following relations:   Ka ¼ Kmax ¼ 12x2 m11 a21 þ m22 a22 ; Umax ¼ Ua ¼

ð117Þ

0

 EA  2 1 k11 a21 þ k22 a22 þ c11 a21 þ 2c12 a1 a2 þ c22 a22 ; 2 8

 3EA  2 1 k11 a21 þ k22 a22 þ c a2 þ 2c12 a1 a2 þ c22 a22 : 2 32 11 1

ð118Þ ð119Þ ð120Þ

Moreover, the condition of energy conservation Umax ¼ Kmax introduced in paper [24] can be rewritten in the form  1  EA  2 1 2 x m11 a21 þ m22 a22  k11 a21 þ k22 a22  c11 a21 þ 2c12 a1 a2 þ c22 a22 ¼ 0: 2 2 8 Because in the Benamar method Eq. (116) is not satisfied, then    3EA  a1 c11 a21 þ 2c12 a1 a2 þ c22 a22 ðc11 a1 þ c12 a2 Þ: Wa ¼  k11  x2 m11 a21  8

ð121Þ

ð122Þ

The fundamental backbone curves obtained by the Galerkin method and the Benamar method are shown in Fig. 7 as the solid and dashed lines, respectively. In Fig. 7 the non-dimensional amplitude means the ratio

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1703

Fig. 7. Fundamental backbone curves of the fixed–fixed beam.

Table 6 The non-dimensional frequency of vibration––comparison of results obtained by the Galerkin method and the Benamar method The Galerkin method

The Benamar method

x=xl

a1 =i

a2 =i

x=xl

a1 =i

a2 =i

1.00 1.05 1.10 1.15 1.20 1.25 1.30

0.00000 1.52259 2.17846 2.70677 3.17089 3.59580 3.99409

0.0000000 )0.0037843 )0.0108308 )0.0201679 )0.0314380 )0.0444276 )0.0589805

1.00000 1.00373 1.01483 1.03297 1.05768 1.08837 1.12441 1.16516 1.21000

0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00

0.0000000 )0.0001389 )0.0010995 )0.0036494 )0.0084564 )0.0160570 )0.0268427 )0.0410608 )0.0588266

a1 =i. Moreover, results of calculation are presented in Table 6. Differences between the results obtained with the help of both methods are significant. In Fig. 8 the average energy error of the compared methods is also shown. As previously, all errors are negative and the modulus of energy error of the Galerkin method is smaller than the energy error of the Benamar method. The backbone curve for the fixed–fixed beam is also determined using the finite element method. The beam is divided into ten identical, two-node elements. The Hermite polynomials of third degree are adopted as shape functions. One harmonic solution is assumed. The amplitude equation is given by Eq. (27b) and the additional condition needed in the Benamar formulation has the form of relation (34). Calculation results are shown in Fig. 7 as small circles and small crosses for the Galerkin and the Benamar methods, respectively. Now the non-dimensional amplitude of vibration means a=i where a denotes the amplitude in the middle of the beam. Good agreement between the semi-analytical solution and the finite element results could be observed. The non-linear algebraic equations arising in the finite element method are solved using the incremental-iterative procedure for the gradually increasing amplitude in the middle of the beam. The Newton–Raphson method is used to solve a system of non-linear equations. Usually, only two or three iterations are enough to fulfil convergence conditions.

1704

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

Fig. 8. Energy errors of the approximate solution for the fixed–fixed beam.

7.3. Two mass system The solutions for a two mass system similar to that shown in Fig. 1 are also calculated. In this case, the amplitude equations resulting from (77) have the following form: i 3k3 h 3 3 a1  ða2  a1 Þ ¼ 0; 4 i 3k3 h 3 a2 þ ða2  a1 Þ3 ¼ 0: x2 Ma2 þ k1 ð2a2  a1 Þ þ 4

x2 Ma1 þ k1 ð2a1  a2 Þ þ

ð123Þ ð124Þ

Because the maximal kinetic and potential energies of the system of interest can be written as Kmax ¼ 12x2 Mða21 þ a22 Þ;

ð125Þ

i k h i 1 h 3 a41 þ ða2  a1 Þ4 þ a42 ; Umax ¼ k1 a21 þ ða2  a1 Þ2 þ a22 þ 2 4

ð126Þ

the condition Kmax ¼ Umax takes the form h i k h i   3 2 4 ð127Þ k1 a21 þ ða2  a1 Þ þ a22 þ a41 þ ða2  a1 Þ þ a42  x2 M a21 þ a22 ¼ 0: 2 The linear frequencies of vibration are x21;lin ¼ k1 =M and x22;lin ¼ 3k1 =M and the corresponding first and second modes of vibration are described by a1 ¼ a2 ¼ a and a1 ¼ a2 ¼ a, respectively. The non-linear modes of vibration determined by the Galerkin method are also symmetrical and antisymmetrical. The backbone curves corresponding to the symmetrical and antisymmetrical modes are given by x21 ¼

k1 3k3 2 þ a; M 4M

ð128Þ

x22 ¼

3k1 27k3 2 þ a: M 4M

ð129Þ

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1705

The Benamar method requires finding the solutions to equations (124) and (127) for a set of gradually increased values of a1 . The incremental-iterative procedure (with respect to a1 ) and the Newton–Raphson method are used to solve these equations. Calculation results in the form of the fundamental backbone curve are presented in Fig. 9. As previously, differences between the results of the Benamar method, shown as the dashed line, and those obtained with the help of the Galerkin method (the solid line in Fig. 9) are significant. Moreover, the non-linear mode of vibration obtained by the Benamar method does not preserve the property of symmetry. This is illustrated in Fig. 10 where the amplitude a2 versus the amplitude a1 is shown. This fact needs some comments. In the case of the Benamar method, the symmetry condition a1 ¼ a2 ¼ a leads us to equations in the form x2 Ma þ k1 a þ

3k3 3 a ¼ 0; 4

Fig. 9. Fundamental backbone curves for two mass system.

Fig. 10. The second mass amplitude versus the first mass amplitude for a two mass system.

ð130Þ

1706

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

Fig. 11. Fundamental backbone curves of the three span beam.

2k1 a2 þ k3 a4  2x2 Ma2 ¼ 0:

ð131Þ

Now, it is clear that the above set of equations has a trivial solution a ¼ 0 only for all x. This means that non-trivial symmetrical solutions do not exist, which cannot be accepted from the physical point of view. Similarly, it can be proved that a purely antisymmetrical solution does not exist either. Unfortunately, the second backbone curve cannot be determined with the help of the above-mentioned incremental-iterative procedure because the iteration process of the Newton–Raphson method does not converge. 7.4. Multispan beam – comparison of results obtained using the finite element method and the vector iteration method The three span beam shown in Fig. 2 is considered. The beam is divided into twenty identical, two node finite elements. The one harmonic solution to the motion equation given by Eq. (14) is assumed in this case. Calculation results are shown in Fig. 11. The non-dimensional, fundamental frequency of vibration x=xl versus the non-dimensional amplitude of vibration a=i in the middle of the first span of the beam is presented. The results obtained by means of the Galerkin and the collocation methods are shown as the solid and the dashed line, respectively. Moreover, the vector iteration method is used to determine the mentioned backbone curve. If the Rayleigh quotation of first kind is used, then the method gives results also shown in Fig. 11. The small circles present the results of solution of the amplitude equations derived with the help of the Galerkin method, while the small crosses show the results of solution of the amplitude equations resulting from the collocation method. It is obvious that both approaches to the solution of amplitude equation give identical results. However, the iteration process of the vector iteration method does not converge if the Rayleigh quotient of second kind is used. These results are in agreement with the remarks presented in Section 5.

8. Concluding remarks The following concluding remarks can be formulated on the basis of the considerations presented in this paper.

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1707

• Several aspects concerning the formulation of and solution to the amplitude equation arising in problems of non-linear free vibration of a system with cubic non-linearity are discussed in this paper. • The concept of the Rayleigh quotient is extended to the non-linear case. The Rayleigh quotient of first kind is proposed as a correct definition in the non-linear case. Moreover, the quotient is used to develop a new and sufficiently efficient method for solving the non-linear eigenvalue problem. The method is an extension of the well-known vector iteration method to the non-linear case. • Calculation results suggest that the Galerkin method gives the most accurate results, when compared with the collocation method and the method proposed in [24]. Moreover, some results obtained by the latter method cannot be accepted as a physically correct solution.

Acknowledgement The author acknowledges the financial support by Poznan University of Technology (Grant No. BW 11307/2003) related to this work. Appendix A The coefficients dil and aijkl are defined as Z 1 T cosðzi xtÞ cosðzl xtÞ dt; dil ¼ T 0 aijkl

1 ¼ T

Z

ðA:1Þ

T

cosðzi xtÞ cosðzj xtÞ cosðzk xtÞ cosðzl xtÞ dt:

ðA:2Þ

0

It can be proved that (see [4]) dli ¼ Iðzi þ zl Þ þ Iðzi  zl Þ;

alijk ¼ 14ðI1 þ I2 þ I3 þ I4 þ I5 þ I6 þ I7 þ I8 Þ;

ðA:3Þ

where I1 ¼ Iðzi þ zj þ zk þ zl Þ; I3 ¼ Iðzi  zj þ zk þ zl Þ;

I2 ¼ Iðzi þ zj  zk  zl Þ; I4 ¼ Iðzi  zj  zk  zl Þ;

I5 ¼ Iðzi þ zj þ zk  zl Þ; I7 ¼ Iðzi  zj þ zk  zl Þ;

I6 ¼ Iðzi þ zj  zk þ zl Þ; I8 ¼ Iðzi  zj  zk þ zl Þ:

ðA:4Þ

Moreover, the function IðzÞ ¼ 0 if z 6¼ 0 and IðzÞ ¼ 1 if z ¼ 0. Now, it is easy to prove that dli can be understood as the Kronecker symbol. Appendix B In this Appendix it is proved that the non-linear Rayleigh quotient of first kind is an approximate solution for amplitude equations with respect to x2 in the sense of the smallest squares method if some approximation of the amplitude vector a is known. Let us consider the general amplitude equation (25). The symmetrical and non-singular matrix f M can be written as a product of two triangular matrices i.e. bT C b: f M¼C

ðB:1Þ

1708

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

Furthermore, after introducing the following transformations b T a; ^¼C x bx ^ a¼C

ðB:2Þ ðB:3Þ

b T , the equation can be rewritten in the form and pre-multiplying the amplitude equation (25) by C b ð^ b  x2 c ^rðx2 Þ ¼ ð K M Þ^ xþG xÞ ¼ 0;

ðB:4Þ

where b T K b 1 ; b ¼C eC K

b T f b 1 ; c M¼C MC

2

b ð^ b T H b 1 x b 1 x e ðC ^Þ C ^: G xÞ ¼ C

ðB:5Þ

2

^ are the solution for the amplitude equation. In the opposite The vector ^rðx Þ is eliminated only if x and x case, we have ^rðx2 Þ 6¼ 0. Let us introduce the error functional in the form Rðx2 Þ ¼ 12^rT ðx2 Þ^rðx2 Þ:

ðB:6Þ 2

The stationary condition dRðx Þ ¼ 0 gives us the following equation: b ð^ b  x2 c ^T bð K x M Þ^ xþG xÞc ¼ 0:

ðB:7Þ

Using the transformation (B.2), Eq. (B.7) can be rewritten in the form x2 ¼

eþH e ðaÞÞa aT ð K T f a Ma

ðB:8Þ

which is simply the non-linear Rayleigh quotient of first kind.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

K.J. Bathe, Finite Element Procedures in Engineering Analysis, Prentice Hall, New York, 1982. L. Meirovitch, Dynamics and Control of Structures, Wiley, New York, 1990. J. Stoer, R. Bulirsch, Wstez p do metod numerycznych, Pa nstwowe Wydawnictwo Naukowe, Warsaw, 1980. R. Lewandowski, Computational formulation for periodic vibration of geometrically nonlinear structures––Part I: Theoretical background, Int. J. Solids Struct. 34 (1997) 1925–1947. O. Zienkiewicz, The Finite Element Method, McGraw-Hill, New York, 1976. R. Lewandowski, Periodic Vibrations of Geometrically Nonlinear Structures, Wydawnictwo Politechniki Pozna nskiej, Pozna n, 1993 (in Polish). R. Lewandowski, Non-linear, steady state analysis of multispan beams by the finite element method, Comput. Struct. 39 (1991) 89–93. C. Mei, Non-linear vibration of beams by matrix displacement method, AIAA J. 10 (1972) 355–357. K. Decha-Umphai, C. Mei, Finite element for non-linear forced vibrations of circular plates, Int. J. Numer. Meth. Engrg 23 (1986) 1715–1726. C. Mei, Finite element displacement method for large amplitude free flexural vibrations of beams and plates, Comput. Struct. 3 (1972) 163–174. C. Mei, K. Decha-Umphai, A finite element method for non-linear forced vibrations of beams, J. Sound Vib. 102 (1985) 369–380. G. Prathap, G.R. Bhashyam, Galerkin finite element method for nonlinear beam vibrations, J. Sound Vib. 72 (1980) 191–203. G.V. Rao, K.K. Raju, I.S. Raju, Finite element formulation for the large amplitude vibrations of beams and orthotropic circular plates, Comput. Struct. 6 (1976) 169–172. B.S. Sarma, T.K. Varadan, Lagrange-type formulation for finite element analysis of non-linear beam vibrations, J. Sound Vib. 86 (1983) 61–70. B.S. Sarma, T.K. Varadan, Ritz finite element approach to nonlinear vibrations of beams, Int. J. Numer. Meth. Engrg. 20 (1984) 353–367. A.N. Jean, H.D. Nelson, Periodic response investigation of large order non-linear rotordynamic systems using collocation, J. Sound Vib. 143 (1990) 473–489.

R. Lewandowski / Comput. Methods Appl. Mech. Engrg. 192 (2003) 1681–1709

1709

[17] A.M. Samoilenko, N.I. Ronto, Numerical-Analytic Methods of Investigating Periodic Solutions, Mir, Moscow, 1979. [18] S.L. Lau, Y.K. Cheung, Amplitude incremental variational principle for nonlinear vibration of elastic systems, ASME J. Appl. Mech. 48 (1981) 959–964. [19] A.A. Ferri, On the equivalence of the incremental harmonic balance-Newton Raphson method, ASME J. Appl. Mech. 53 (1986) 455–457. [20] L. Wellford, G.M. Dib, W. Mindle, Free and steady state vibration of non-linear structures using a finite element-non-linear eigenvalue technique, Earthquake Engrg. Struct. Dyn. 8 (1980) 97–115. [21] M.A.E. Ghabrial, L.C. Wellford, An averaged LagrangeÕ-finite element technique for the solution of nonlinear vibration problems, Comput. Struct. 16 (1983) 207–214. [22] J. Padovan, Non-linear vibrations of general structures, J. Sound Vib. 72 (1980) 427–441. [23] R. Benamar, Non-linear dynamic behaviour of fully clamped beams and rectangular isotropic and laminated plates, Ph.D. Thesis, Institute of Sound and Vibration Research, 1990. [24] R. Benamar, M.M.K. Bennouna, R.G. White, The effects of large vibration amplitudes on the mode shapes and natural frequencies of thin elastic structures. Part I: simply supported and clamped–clamped beams, J. Sound Vib. 149 (1991) 179–195. [25] M. El Kadiri, R. Benamar, R.G. White, Improvement of the semi-analytical method, for determining of the geometrically nonlinear response of thin straight structures. Part I: application to clamped–clamped and simply supported-clamped beams, J. Sound Vib. 249 (2002) 263–305. [26] B. Harras, R. Benamar, R.G. White, Experimental and theoretical investigation of the linear and non-linear dynamic behaviour of a glare 3 hybrid composite panel, J. Sound Vib. 252 (2002) 281–315. [27] B. Harras, R. Benamar, R.G. White, Geometrically non-linear free vibration analysis of fully clamped composite laminated plates, Part I: general theory, experimental and theoretical investigations of the fundamental non-linear mode shape, J. Sound Vib. 251 (2002) 579–619. [28] R. Benamar, M.M.K. Bennouna, R.G. White, The effects of large vibration amplitudes on the fundamental mode shape of thin elastic structures. Part II: fully clamped rectangular isotropic plates, J. Sound Vib. 164 (1993) 295–316. [29] M. El Kadiri, R. Benamar, R.G. White, The non-linear free vibration of fully clamped rectangular plates: second non-linear mode for various plate aspect ratios, J. Sound Vib. 228 (1999) 333–358. [30] P. Lindqvist, Note on a nonlinear eigenvalue problem, Rocky Mt. J. Math. 23 (1993) 281–288. [31] R. Benamar, M.M.K. Bennouna, R.G. White, The effects of large vibration amplitudes on the mode shapes and natural frequencies of thin elastic structures, Part II: fully clamped rectangular isotropic plates, J. Sound Vib. 164 (1993) 295–316. [32] S.L. Lau, Y.K. Cheung, S.Y. Wu, Nonlinear vibrations of thin elastic plates, Part I: generalized incremental HamiltonÕs principle and element formulation, J. Appl. Mech. 51 (1984) 837–844. [33] S. Woinowsky-Krieger, The effect of axial force on the vibration of hinged bars, J. Appl. Mech. 17 (1950) 35–36.