Nonlinear oscillations and chaotic dynamics of a supported FGM rectangular plate system under mixed excitations

Yasser S. Hamed1

1Department of Mathematics and Statistics, Faculty of Science, Taif University, El-Taif, El-Haweiah, P.O. Box 888, Zip Code 21974, Kingdom of Saudi Arabia

1Department of Engineering Mathematics, Faculty of Electronic Engineering, Menouf, 32952, Egypt

Journal of Vibroengineering, Vol. 16, Issue 7, 2014, p. 3218-3235.
Received 3 December 2013; received in revised form 6 June 2014; accepted 20 June 2014; published 15 November 2014

Copyright © 2014 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

In this paper, we investigated the nonlinear oscillations, vibrations and chaotic dynamics of a simply supported rectangular plate made of functionally graded materials (FGMs) through a temperature field subjected to mixed excitations. The rectangular FGM plate system described by a coupled of nonlinear differential equations (two degree of freedom) including the quadratic and cubic nonlinear terms. The mathematical solutions of the governing equations of motion for the FGM plate derived using the perturbation method based on the power series expansion up to and including the second order approximation. The numerical simulations investigated using Runge-Kutta of fourth order using MATLAB and MAPLE programs. All different resonance cases reported and studied numerically. We applied both frequency response equations and phase-plane technique and also Lyapunov’s first method near the worst resonance cases to analyze the stability of the steady state solution of vibrating system. The effects of the different parameters of the rectangular plate system studied numerically. Results compared to previously published work. In the future work, we can deal with the same system subjected multi-external and tuned and parametric excitation forces. Also, the system can be studied at another worst different resonance cases, active and passive controller.

Keywords: FGM rectangular plate, vibrations, chaotic dynamics, stability, saturation.

1. Introduction

The functionally graded materials (FGMs) have received a considerable attention recently as one class of inhomogeneous composite materials. Functionally graded materials (FGMs) are inhomogeneous composites within which materials properties vary continuously. This is achieved by gradually changing the composition of the constituent materials, usually in one direction, to obtain smooth variation in material properties [1]. Reddy [2] developed both theoretical and finite element formulations for thick FGM plates according to the higher-order shear deformation plates theory, and studied the non-linear dynamic response of FGM plates subjected to a sudden applied uniform pressure. Yang et al. [3] studied the large amplitude vibration of pre-stressed FGM laminated plates that were composed of an FGM layer and two surface-mounted piezoelectric actuator layers. Huang and Shen [4] dealt with the nonlinear vibration and dynamic response of FGM plates in heat environment. Heat conduction and temperature-dependent material properties were both considered. More recently, Yang and his coworkers studied the sensitivity of nonlinear vibration and dynamic response of FGM plates to initial geometrical imperfections of arbitrary shape [5, 6]. Ye et al. [7] dealt with nonlinear dynamic behavior of a parametrically excited, simply supported, symmetric cross-ply laminated rectangular thin plate. The geometric nonlinearity and nonlinear damping were included in the governing equations of motion. FGM were initially designed as thermal barriers for aerospace structures and fusion reactors. Usually, they are made from a mixture of metals and ceramics. Due to their unique advantage of being able to withstand sever high-temperature, FGM have found a wide range of applications in many industries, especially in space vehicles and aircrafts, where they are often subjected to high levels of thermal and dynamic loading, such as large temperature gradients and acoustic pressure, which may result in a complicated nonlinear vibration behavior of the FGM plate structures [8]. Hao et al. [9] established the nonlinear governing equations of motion for a simply supported FGM rectangular plate subjected to the transversal and in-plane excitations in thermal environment and investigated the nonlinear oscillations, bifurcations and chaos of the functionally graded material plate. Zhang et al. [10] studied the chaotic vibrations of an orthotropic FGM rectangular plate based on third-order shear deformation theory. Zhang et al. [10] investigated the nonlinear vibrations and chaotic dynamics of a simply supported orthotropic FGM rectangular plate subjected to the in-lane and transverse excitations together with thermal loading. The present study is focused on the case of 1:2:4 internal resonance, primary parametric resonance and sub-harmonic resonance of order 1/2. Yang et al. [11] analyzed the nonlinear dynamic responses of a functionally graded plate with a through-width surface crack. Li et al. [12] used the extended Melnikov method to investigate the multi-pulse chaotic dynamics of a functionally graded material rectangular plate with one-to-one internal resonance. Malekzadeh [13] studied free vibration of functionally graded (FG) thin-to-moderately thick annular plates subjected to thermal environment and supported on two-parameter elastic foundation. Akbarzadeh et al. [14] investigated the dynamic response of a simply supported functionally graded rectangular plate subjected to a lateral thermomechanical by using the hybrid Fourier-Laplace transform method. Zhang et al. [15] investigated the nonlinear periodic and chaotic motions of a simply supported FGM rectangular plate under both in plane and transverse excitations in thermal environment. The present analysis is focused on the case of 1:2 internal resonance and primary parametric resonance, 1/2 sub-harmonic resonance. Eissa and Sayed [16-18], studied the effects of different active controllers on simple and spring pendulum at the primary resonance via negative velocity feedback or it’s square or cubic. Hamed et al. [19-21] studied USM model subject to multi-external or both multi-external and multi-parametric and both multi-external and tuned excitation forces. The model consists of multi-degree-of-freedom system consisting of the tool holder and absorbers (tools) simulating ultrasonic machining process. The advantages of using multi-tools are to machine different materials and different shapes at the same time. This leads to time saving and higher machining efficiency. Sayed and Hamed [22] studied the response of a two-degree-of-freedom system with quadratic coupling under parametric and harmonic excitations. The method of multiple scale perturbation technique is applied to solve the non-linear differential equations and obtain approximate solutions up to and including the second-order approximations. Hamed et al [23] presented the behavior of the nonlinear string beam coupled system subjected to external, parametric and tuned excitations for case 1:1 internal resonance. The stability of the system studied using frequency response equations and phase-plane method. It is found from numerical simulations that there are obvious jumping phenomena in the frequency response curves. Sayed and et al. [24] investigated the non-linear dynamics of a two-degree-of freedom vibration system including quadratic and cubic non-linearities subjected to external and parametric excitation forces. There exist multi-valued solutions which increase or decrease by the variation of some parameters. The numerical simulations show the system exhibits periodic motions and chaotic motions. Amer and Sayed [25] studied the response of one-degree-of freedom, non-linear system under multi-parametric and external excitation forces simulating the vibration of the cantilever beam. The solution of this system up to and including the second order approximation is determined applying the multiple time scale perturbation. The steady state solution and its stability are determined. Sayed and Kamel [26, 27] investigated the effect of different controllers on the vibrating system and the saturation control of a linear absorber to reduce vibrations due to rotor blade flapping motion. The stability of the obtained numerical solution is investigated using both phase plane methods and frequency response equations. Sayed and Mousa [28] investigated the influence of the quadratic and cubic terms on non-linear dynamic characteristics of the angle-ply composite laminated rectangular plate with parametric and external excitations. They focused on the case of 1:2 internal resonance and sub-harmomic resonance. Sayed and Mousa [29] investigated the perturbation method and stability of the composite laminated piezoelectric rectangular plate under simultaneous transverse and in-plane excitations. They focused on the case of 1:1:3 internal resonance and primary resonance. The stability of the system and the effects of different parameters on system behavior have been studied using frequency response curves. Stability is performed of figures by solid and dotted lines. The analytical results given by the method of multiple time scale is verified by comparison with results from numerical integration of the modal equations. It is quite clear that some of the simultaneous resonance cases are undesirable in the design of such system.

2. Mathematical analysis

A simply supported FGM rectangular plate subjected to a through-thickness temperature field together with tuned, parametric and external excitations is considered in this study. The plate is of length a, width b, and thickness h as shown in Fig. 1. Assume that (uvw) represent the displacements of an arbitrary point of the FGM rectangular plate in the x, y, and z directions respectively. Along the plate edges x= 0, x=a, the plate is subjected to an in-plane excitation p=p0-p1cosγ2t where p0, p1 are the external excitation amplitude in y direction and a transverse excitation F(x,y)(cosγ1t+sinγ3t*cosγ4t), where F(x,y) are the parametric and tuned excitation amplitudes in z direction. Here ϕx and ϕy represent the mid-plane rotations of a transverse normal about the y and x axes, γ1, γ2, γ3 and γ4 are the excitation frequencies of the transverse and in-plane excitations, respectively. According to the Reddy’s third-order shear deformation plate theory (TSDT) [30], the displacement field of the plate is assumed as:

(1a)
u x , y , t = u 0 x , y , t + z ϕ x x , y , t - c 1 z 3 ϕ x + w 0 x ,
(1b)
v ( x , y , t ) = v 0 ( x , y , t ) + z ϕ y ( x , y , t ) - c 1 z 3 ϕ y + w 0 y ,
(1c)
w x , y , t = w 0 x , y , t ,

where (u0,v0,w0) represent the displacements of a point in middle plane of the FGM rectangular plate in the x, y, and z directions respectively. According to the von Karman strain-displacements relationship and the Hamilton’s principle, the nonlinear governing equations of motion for the FGM rectangular plate are derived as:

(2a)
N x x , x + N x y , y = I 0 u ¨ 0 + I 1 - c 1 I 3 ϕ ¨ x - c 1 I 3 w ¨ 0 x ,
(2b)
N y y , y + N x y , y = I 0 v ¨ 0 + ( I 1 - c 1 I 3 ) ϕ ¨ y - c 1 I 3 w ¨ 0 y ,
(2c)
N y y , y w 0 y + N y y 2 w 0 y 2 + N x y , y w 0 x + N x y , x w 0 y + 2 N x y 2 w 0 y x + N x x , x w 0 x + N x x 2 w 0 x 2
            + c 1 p x x , x x + 2 p x y , x y + p y y , y y + Q x , x - c 2 R x , x + Q y , y - c 2 R y , y + F - γ w ˙ 0
            = I 0 w ¨ 0 + c 1 I 3 u ¨ 0 x + v ¨ 0 y + c 1 I 4 - c 1 2 I 6 ϕ ¨ x x + ϕ ¨ y y - c 1 2 I 6 2 w ¨ 0 x 2 + 2 w ¨ 0 y 2 ,
(2d)
M x x , x + M x y , y - c 1 p x x , x + p x y , y - Q x - c 2 R x = I 1 - c 1 I 3 u ¨ 0
            + I 2 - 2 c 1 I 4 + c 1 2 I 6 ϕ ¨ x - c 1 I 4 - c 1 2 I 6 w ¨ 0 x ,
(2e)
M y y , y + M x y , y - c 1 p y y , y + p x y , x - Q y - c 2 R y = I 1 - c 1 I 3 v ¨ 0
            + I 2 - 2 c 1 I 4 + c 1 2 I 6 ϕ ¨ y - ( c 1 I 4 - c 1 2 I 6 ) w ¨ 0 y .

The boundary conditions for the simply supported plate under current consideration require at:

(3a)
x = 0 ,       x = a ,       w = ϕ y = M x x = p x x = N x y = 0 ,
(3b)
y = 0 ,       y = b ,       w = ϕ x = M y y = p y y = N x y = 0 ,
(3c)
N y y | y = 0 , b = 0 ,       0 b N x x | x = 0 , a d y = - 0 b ( p 0 - p 1 c o s γ 2 t ) d y .

Fig. 1. The model of a rectangular FGM plate and the coordinate system

 The model of a rectangular FGM plate and the coordinate system

Applying the Galerkin procedure yield the dimensionless governing differential equation of transverse motion of the FGM rectangular plate with two degrees of freedom as follows:

The governing differential equation of transverse motion of the FGM rectangular plate with two degrees of freedom as follows:

(4a)
x ¨ + ω 1 2 x = ε - β 1 x 2 - β 2 y 2 - β 3 x 3 - β 4 x y 2 - β 5 x y
            + ε 2 ρ 1 c o s γ 1 t - ρ 11 x c o s γ 2 t + ρ 21 s i n γ 3 t * c o s γ 4 t - μ 1 x ˙ ,
(4b)
y ¨ + ω 2 2 y = ε - η 1 y 2 - η 2 x 2 - η 3 y 3 - η 4 x 2 y - η 5 x y
            + ε 2 ρ 2 c o s γ 1 t - ρ 12 y c o s γ 2 t + ρ 22 s i n γ 3 t * c o s γ 4 t - μ 2 y ˙ .

The initial conditions for the simply supported plate are:

(5)
x t = 0 = 0 ,       y t = 0 = 0 ,       x ˙ t = 0 = 0 ,       y ˙ ( t = 0 ) = 0 ,

where x and y are the amplitudes of the FGM rectangular plate for the first and second modes, ω1 and ω2 the natural frequencies of the rectangular plate modes, μ1 and μ2 are the modal damping coefficients, and γ1, γ2, γ3 and γ4 are the external and parametric and tuned excitation frequencies respectively. ρ1, ρ2, ρ11, ρ12, ρ21 and ρ22 are related to the magnitude of the excitation forces of the FGM rectangular plate modes, ε is a small perturbation parameter and 0<ε<<1. We introducing three time scales:

(6)
T 0 = t ,       T 1 = ε t ,       T 2 = ε 2 t .

We seek an asymptotic solution of Eq. (4) for x and y in the form:

(7a)
x ( t , ε ) = x 10 ( T 0 , T 1 , T 2 ) + ε x 11 ( T 0 , T 1 , T 2 ) + ε 2 x 12 ( T 0 , T 1 , T 2 ) + O ( ε 3 ) ,
(7b)
y t , ε = y 10 T 0 , T 1 , T 2 + ε y 11 T 0 , T 1 , T 2 + ε 2 y 12 T 0 , T 1 , T 2 + O ε 3 .

Derivatives with respect to t are then transformed into:

(8a)
d d t = D 0 + ε D 1 + ε 2 D 2 ,
(8b)
d 2 d t 2 = D 0 2 + 2 ε D 0 D 1 + ε 2 D 1 2 + 2 D 0 D 2 .

The derivatives x˙, y˙, x¨ and y¨ are in the form:

(9a)
x ˙ = D 0 x 10 + ε ( D 0 x 11 + D 1 x 10 ) + ε 2 ( D 0 x 12 + D 1 x 11 + D 2 x 10 ) ,
(9b)
x ¨ = D 0 2 x 10 + ε ( D 0 2 x 11 + 2 D 0 D 1 x 10 ) + ε 2 ( D 0 2 x 12 + 2 D 0 D 1 x 11 + D 1 2 x 10 + 2 D 0 D 2 x 10 ) ,
(9c)
y ˙ = D 0 y 10 + ε ( D 0 y 11 + D 1 y 10 ) + ε 2 ( D 0 y 12 + D 1 y 11 + D 2 y 10 ) ,
(9d)
y ¨ = D 0 2 y 10 + ε D 0 2 y 11 + 2 D 0 D 1 y 10 + ε 2 D 0 2 y 12 + 2 D 0 D 1 y 11 + D 1 2 y 10 + 2 D 0 D 2 y 10 ,

where Dn=/Tn, n= 0, 1, 2. Terms of O(ε3) and higher orders are neglected. Substituting Eqs. (7)-(9) into Eq. (4) we get:

(10a)
D 0 2 x 10 + ε D 0 2 x 11 + 2 D 0 D 1 x 10 + ε 2 D 0 2 x 12 + 2 D 0 D 1 x 11 + D 1 2 x 10 + 2 D 0 D 2 x 10
            + ω 1 2 x 10 + ε x 11 + ε 2 x 12
            = ε - β 1 x 10 + ε x 11 + ε 2 x 12 2 - β 2 y 10 + ε y 11 + ε 2 y 12 2 - β 3 x 10 + ε x 11 + ε 2 x 12 3 - β 4 x 10 + ε x 11 + ε 2 x 12 y 10 + ε y 11 + ε 2 y 12 2 - β 5 x 10 + ε x 11 + ε 2 x 12 y 10 + ε y 11 + ε 2 y 12
            + ε 2 ρ 1 c o s γ 1 t - ρ 11 x 10 + ε x 11 + ε 2 x 12 c o s γ 2 t + ρ 21 s i n γ 3 t c o s γ 4 t - μ 1 D 0 x 10 + ε D 0 x 11 + D 1 x 10 + ε 2 D 0 x 12 + D 1 x 11 + D 2 x 10 ,
(10b)
D 0 2 y 10 + ε D 0 2 y 11 + 2 D 0 D 1 y 10 + ε 2 D 0 2 y 12 + 2 D 0 D 1 y 11 + D 1 2 y 10 + 2 D 0 D 2 y 10
            + ω 2 2 y 10 + ε y 11 + ε 2 y 12
              = ε - η 1 y 10 + ε y 11 + ε 2 y 12 2 - η 2 x 10 + ε x 11 + ε 2 x 12 2 - η 3 y 10 + ε y 11 + ε 2 y 12 3 - η 4 x 10 + ε x 11 + ε 2 x 12 2 y 10 + ε y 11 + ε 2 y 12 - η 5 x 10 + ε x 11 + ε 2 x 12 y 10 + ε y 11 + ε 2 y 12
              + ε 2 ρ 2 c o s γ 1 t - ρ 12 y 10 + ε y 11 + ε 2 y 12 c o s γ 2 t + ρ 22 s i n γ 3 t c o s γ 4 t - μ 2 D 0 y 10 + ε D 0 y 11 + D 1 y 10 + ε 2 D 0 y 12 + D 1 y 11 + D 2 y 10 .

We expand the two sides of Eqs. (10) and equating the coefficients of similar powers of ε, one obtains the following set of ordinary differential equations:

Order ε0:

(11a)
( D 0 2 + ω 1 2 ) x 10 = 0 ,
(11b)
( D 0 2 + ω 2 2 ) y 10 = 0 .

Order ε1:

(12a)
D 0 2 + ω 1 2 x 11 = - 2 D 0 D 1 x 10 - β 1 x 10 2 - β 2 y 10 2 - β 3 x 10 3 - β 4 x 10 y 10 2 - β 5 x 10 y 10 ,
(12b)
D 0 2 + ω 2 2 y 11 = - 2 D 0 D 1 y 10 - η 1 y 10 2 - η 2 x 10 2 - η 3 y 10 3 - η 4 y 10 x 10 2 - η 5 y 10 x 10 .

Order ε2:

(13a)
D 0 2 + ω 1 2 x 12 = - 2 D 0 D 2 x 10 - D 1 2 x 10 - 2 D 0 D 1 x 11 - 2 β 1 x 10 x 11 - 2 β 2 y 10 y 11
            - 3 β 3 x 10 2 x 11 - μ 1 D 0 x 10 - β 4 x 11 y 10 2 + 2 x 10 y 10 y 11 - β 5 x 10 y 11 + x 11 y 10
            - ρ 11 x 10 c o s γ 2 t + ρ 1 c o s γ 1 t + ρ 21 s i n γ 3 t c o s γ 4 t ,
(13b)
D 0 2 + ω 2 2 y 12 = - 2 D 0 D 2 y 10 - D 1 2 y 10 - 2 D 0 D 1 y 11 - 2 η 1 y 10 y 11 - 2 η 2 x 10 x 11
            - 3 η 3 y 10 2 y 11 - μ 2 D 0 y 10 - η 4 y 11 x 10 2 + 2 y 10 x 10 x 11 - η 5 y 10 x 11 + y 11 x 10
            - ρ 12 y 10 c o s γ 2 t + ρ 2 c o s γ 1 t + ρ 22 s i n γ 3 t c o s γ 4 t .

The general solutions of Eq. (11), can be written in the form:

(14a)
x 10 = A 1 e x p i ω 1 T 0 + A - 1 e x p - i ω 1 T 0 ,
(14b)
y 10 = A 2 e x p i ω 2 T 0 + A - 2 e x p - i ω 2 T 0 ,

where A1 and A2 are a complex function in T1, T2. Substituting Eq. (14) into Eq. (12) we get:

(15a)
D 0 2 + ω 1 2 x 11 = - 2 D 0 D 1 A 1 - 2 β 4 A 1 A 2 A - 2 - 3 β 3 A 1 2 A - 1 exp i ω 1 T 0
            + - 2 D 0 D 1 A - 1 - 2 β 4 A - 1 A - 2 A 2 - 3 β 3 A - 1 2 A 1 exp - i ω 1 T 0 - β 1 A 1 2 exp 2 i ω 1 T 0
            - β 1 A - 1 2 exp - 2 i ω 1 T 0 - β 2 A 2 2 exp 2 i ω 2 T 0 - β 2 A - 2 2 exp - 2 i ω 2 T 0
            - β 3 A 1 3 e x p 3 i ω 1 T 0 - β 3 A - 1 3 e x p - 3 i ω 1 T 0 - β 5 A 1 A 2 e x p i ω 1 + ω 2 T 0
            - β 5 A - 1 A - 2 exp - i ω 1 + ω 2 T 0 - β 5 A 1 A - 2 exp i ω 1 - ω 2 T 0
            - β 5 A - 1 A 2 exp - i ω 1 - ω 2 T 0 - β 4 A 1 A 2 2 exp i ω 1 + 2 ω 2 T 0
            - β 4 A - 1 A - 2 2 exp - i ω 1 + 2 ω 2 T 0 - β 4 A 1 A - 2 2 exp i ω 1 - 2 ω 2 T 0
            - β 4 A - 1 A 2 2 e x p ( - i ( ω 1 - 2 ω 2 ) T 0 ) - 2 β 1 A 1 A - 1 - 2 β 2 A 2 A - 2 ,
(15b)
D 0 2 + ω 2 2 y 11 = - 2 D 0 D 1 A 2 - 2 η 4 A 1 A - 1 A 2 - 3 η 3 A 2 2 A - 2 e x p i ω 2 T 0
            + - 2 D 0 D 1 A - 2 - 2 η 4 A - 1 A 1 A - 2 - 3 η 3 A 2 A - 2 2 e x p - i ω 2 T 0 - η 2 A 1 2 e x p 2 i ω 1 T 0
            - η 1 A - 1 2 e x p - 2 i ω 1 T 0 - η 1 A 2 2 e x p 2 i ω 2 T 0 - η 1 A - 2 2 e x p - 2 i ω 2 T 0
            - η 3 A 2 3 e x p 3 i ω 2 T 0 - η 3 A - 2 3 e x p - 3 i ω 2 T 0 - η 5 A 1 A 2 e x p i ω 2 + ω 1 T 0
            - η 5 A - 1 A - 2 e x p - i ω 2 + ω 1 T 0 - η 5 A 2 A - 1 e x p i ω 2 - ω 1 T 0
            - η 5 A - 2 A 1 e x p - i ω 2 - ω 1 T 0 - η 4 A 2 A 1 2 e x p i ω 2 + 2 ω 1 T 0
            - η 4 A - 2 A - 1 2 e x p - i ω 2 + 2 ω 1 T 0 - η 4 A 2 A - 1 2 e x p i ω 2 - 2 ω 1 T 0
            - η 4 A - 2 A 1 2 e x p - i ω 2 - 2 ω 1 T 0 - 2 η 1 A 2 A - 2 - 2 η 2 A 1 A - 1 .

Eliminating the secular terms e(±iω1T0) and e(±iω2T0) from Eqs. (15), then the first-order approximations are given by:

(16a)
x 11 = K 1 e x p 2 i ω 1 T 0 + K 2 e x p 2 i ω 2 T 0 + K 3 e x p 3 i ω 1 T 0 + K 4 e x p i ω 1 + ω 2 T 0
            + K 5 e x p i ω 1 - ω 2 T 0 + K 6 e x p i ω 1 + 2 ω 2 T 0 + K 7 e x p i ω 1 - 2 ω 2 T 0
            + K 8 + c c ,
(16b)
y 11 = G 1 e x p 2 i ω 1 T 0 + G 2 e x p 2 i ω 2 T 0 + G 3 e x p 3 i ω 2 T 0 + G 4 e x p i ω 2 + ω 1 T 0
            + G 5 e x p i ω 2 - ω 1 T 0 + G 6 e x p i ω 2 + 2 ω 1 T 0 + G 7 e x p i ω 2 - 2 ω 1 T 0
            + G 8 + c c ,

where Ki, (i= 1,..., 8) and Gi, (i= 1,..., 8) are a complex function in T1, T2 which are given in the appendix and cc stands for the complex conjugate of the preceding terms.

With the same steps used to get the first-order approximation, substituting Eqs. (14), (16) into Eq. (13) and eliminating the secular terms, then the second-order approximations are given by:

(17a)
x 12 T 0 , T 1 , T 2 =    E 1 e x p ω 2 T 0 + E 2 e x p 2 i ω 1 T 0 + E 3 e x p 2 i ω 2 T 0
            + E 4 e x p 3 i ω 1 T 0 + E 5 e x p 3 i ω 2 T 0 + E 6 e x p 4 i ω 1 T 0 + E 7 e x p 4 i ω 2 T 0
            + E 8 e x p 5 i ω 1 T 0 + E 9 e x p i ω 1 + ω 2 T 0 + E 10 e x p i ω 1 - ω 2 T 0
            + E 11 e x p i ω 1 + 2 ω 2 T 0 + E 12 e x p i ω 1 - 2 ω 2 T 0 + E 13 e x p i ω 1 + 3 ω 2 T 0
            + E 14 e x p i ω 1 - 3 ω 2 T 0 + E 15 e x p i ω 1 + 4 ω 2 T 0 + E 16 e x p i ω 1 - 4 ω 2 T 0
            + E 17 e x p i 2 ω 1 + ω 2 T 0 + E 18 e x p i 2 ω 1 - ω 2 T 0 + E 19 e x p i 2 ω 1 + 2 ω 2 T 0
            + E 20 e x p i 2 ω 1 - 2 ω 2 T 0 + E 21 e x p i 3 ω 1 + ω 2 T 0 + E 22 e x p i 3 ω 1 - ω 2 T 0
            + E 23 e x p i 3 ω 1 + 2 ω 2 T 0 + E 24 e x p i 3 ω 1 - 2 ω 2 T 0 + E 25 e x p i γ 2 + ω 1 T 0
            + E 26 e x p i γ 2 - ω 1 T 0 + E 27 e x p i γ 1 T 0 + E 28 e x p i γ 3 + γ 4 T 0
            + E 29 e x p i γ 3 - γ 4 T 0 + E 30 + c c ,
(17b)
y 12 T 0 , T 1 , T 2 =    H 1 e x p ω 1 T 0 + H 2 e x p 2 i ω 1 T 0 + H 3 e x p 2 i ω 2 T 0
            + H 4 e x p 3 i ω 1 T 0 + H 5 e x p 3 i ω 2 T 0 + H 6 e x p 4 i ω 1 T 0 + H 7 e x p 4 i ω 2 T 0
            + H 8 e x p 5 i ω 2 T 0 + H 9 e x p i ω 1 + ω 2 T 0 + H 10 e x p i ω 1 - ω 2 T 0
            + H 11 e x p i ω 1 + 2 ω 2 T 0 + H 12 e x p i ω 1 - 2 ω 2 T 0 + H 13 e x p i ω 1 + 3 ω 2 T 0
            + H 14 e x p i ω 1 - 3 ω 2 T 0 + H 15 e x p i 4 ω 1 + ω 2 T 0 + H 16 e x p i 4 ω 1 - ω 2 T 0
            + H 17 e x p i 2 ω 1 + ω 2 T 0 + H 18 e x p i 2 ω 1 - ω 2 T 0 + H 19 e x p i 2 ω 1 + 2 ω 2 T 0
            + H 20 e x p i 2 ω 1 - 2 ω 2 T 0 + H 21 e x p i 3 ω 1 + ω 2 T 0 + H 22 e x p i 3 ω 1 - ω 2 T 0
            + H 23 e x p i 2 ω 1 + 3 ω 2 T 0 + H 24 e x p i 2 ω 1 - 3 ω 2 T 0 + H 25 e x p i γ 2 + ω 2 T 0
            + H 26 e x p i γ 2 - ω 2 T 0 + H 27 e x p i γ 1 T 0 + H 28 e x p i γ 3 + γ 4 T 0
            + H 29 e x p i γ 3 - γ 4 T 0 + H 30 + c c ,

whereEj, (j= 1,..., 30) and Hj, (j= 1,..., 30) are complex functions in T1, T2 and cc represents the complex conjugates.

From the above-derived solutions, many resonance cases can be deduced. The reported resonance cases are classified into:

a) Primary resonance: γ1ω1, γ1ω2.

b) Sub-Harmonic resonance: γ22ω1, γ22ω2.

c) Internal or secondary resonance: ω2ω1,ω22ω1,ω23ω1,ω24ω1,ω21/2ω1, ω21/3ω1, ω21/4ω1 and ω22/3ω1, ω23/2ω1.

d) Combined resonance: ±γ3±γ4±ω1, ±γ3±γ4±ω2.

e) Simultaneous or incident resonance: Any combination of the above resonance cases is considered as simultaneous or incident resonance.

3. Stability of motion

Here, to investigate the stability the simultaneous primary γ1ω1, combined γ3-γ4ω1 and internal resonance ω21/2ω1 are considered. We introduce detuning parameters σ1, σ2 and σ3 such that:

(18)
γ 1 = ω 1 + ε σ 1 ,       γ 3 - γ 4 = ω 1 + ε σ 2 ,       2 ω 2 = ω 1 + ε σ 3 .

This case represents the system worst case. Substituting Eq. (18) into Eqs. (12) and Eq. (13) and eliminating the secular terms, leads to the solvability conditions for the first and second order approximation, we get:

(19a)
2 i ω 1 d A 1 d t = ε - 2 β 4 ω 1 A 1 A 2 A - 2 - 3 β 3 A 1 2 A - 1 - β 2 A 2 2 exp i σ 3 T 1
            + ε 2 - D 1 2 A 1 - i ω 1 μ 1 A 1 + Γ 1 A 1 3 A - 1 2 + Γ 2 A 1 2 A - 1 A 2 A - 2 + Γ 3 A 1 2 A - 1 + Γ 4 A 1 A 2 2 A - 2 2
            + Γ 5 A 1 A 2 A - 2 + Γ 6 A 1 A - 1 A 2 2 + Γ 7 A 2 3 A - 2 exp i σ 3 T 1 + Γ 8 A 1 2 A - 2 2 exp - i σ 3 T 1
            + ( Γ 9 A - 1 A 2 4 ) e x p ( 2 i σ 3 T 1 ) + ρ 1 2 e x p ( i σ 1 T 1 ) - i ρ 21 4 e x p ( i σ 2 T 1 ) ,
(19b)
2 i ω 2 d A 2 d t = ε - 2 η 4 A 2 A 1 A - 1 - 3 η 3 A 2 2 A - 2 - η 5 A - 2 A 1 exp - i σ 3 T 1
            + ε 2 - D 1 2 A 2 - i ω 2 μ 2 A 2 + Γ 10 A 2 3 A - 2 2 + Γ 11 A 2 2 A - 2 A 1 A - 1 + Γ 11 A 2 2 A - 2 + Γ 12 A 2 A 1 2 A - 1 2
            + Γ 13 A 2 A 1 A - 1 + Γ 14 A - 1 A 2 3 e x p i σ 3 T 1 + ( Γ 15 A 1 2 A - 2 3 ) e x p ( - 2 i σ 3 T 1 ) + ( Γ 16 A 1 A 2 A - 2 2
            + Γ 17 A 1 2 A - 1 A - 2 ) e x p ( - i σ 3 T 1 ) ,

where Γm, (m= 1,..., 17) are constants.

To analyze the solutions of Eqs. (19), we express A1(T1,T2) and A2(T1,T2) in the polar form:

(20)
A 1 T 1 , T 2 = a 1 2 e i φ 1 ,       A 2 T 1 , T 2 = a 2 2 e i φ 2 ,

where a1, a2 and φ1, φ2 are the steady state amplitudes and phases of the motion. Substituting Eqs. (20) into Eq. (19) and equating the real and imaginary parts we obtain the following equations describing the modulation of the amplitudes and phases of the first and second modes of FGM rectangular plate response:

(21a)
a ˙ 1 = - μ 1 2 a 1 + ρ 1 2 ω 1 s i n θ 1 - ρ 21 4 ω 1 c o s θ 2 - Γ 23 ω 1 a 1 2 a 2 2 s i n θ 3 + Γ 24 ω 1 a 1 a 2 4 s i n 2 θ 3
            + Γ 25 ω 1 a 2 2 s i n θ 3 + Γ 26 ω 1 a 2 4 s i n θ 3 + Γ 27 ω 1 a 1 2 a 2 2 s i n θ 3 ,
(21b)
a 1 φ ˙ 1 = - Γ 18 ω 1 a 1 3 - Γ 19 ω 1 a 1 5 - Γ 20 ω 1 a 1 a 2 2 - Γ 21 ω 1 a 1 a 2 4 - Γ 22 ω 1 a 1 3 a 2 2 - ρ 1 2 ω 1 c o s θ 1 - ρ 21 4 ω 1 s i n θ 2
            - Γ 23 ω 1 a 1 2 a 2 2 c o s θ 3 - Γ 24 ω 1 a 1 a 2 4 c o s 2 θ 3 - Γ 25 ω 1 a 2 2 c o s θ 3 - Γ 26 ω 1 a 2 4 c o s θ 3 - Γ 27 ω 1 a 1 2 a 2 2 c o s θ 3 ,
(22a)
a ˙ 2 = - μ 2 2 a 2 - Γ 33 ω 2 a 1 a 2 s i n θ 3 - Γ 34 ω 2 a 1 a 2 3 s i n θ 3 - Γ 35 ω 2 a 1 3 a 2 s i n θ 3 - Γ 36 ω 2 a 1 a 2 3 s i n θ 3
            + Γ 37 ω 2 a 1 a 2 3 s i n θ 3 - Γ 38 ω 2 a 1 2 a 2 3 s i n 2 θ 3 ,
(22b)
a 2 φ ˙ 2 = - Γ 28 ω 2 a 2 3 - Γ 29 ω 2 a 2 5 - Γ 30 ω 2 a 1 2 a 2 - Γ 31 ω 2 a 1 4 a 2 - Γ 32 ω 2 a 1 2 a 2 3 - Γ 33 ω 2 a 1 a 2 c o s θ 3
            - Γ 34 ω 2 a 1 a 2 3 c o s θ 3 - Γ 35 ω 2 a 1 3 a 2 c o s θ 3 - Γ 36 ω 2 a 1 a 2 3 c o s θ 3 - Γ 37 ω 2 a 1 a 2 3 c o s θ 3
            - Γ 38 ω 2 a 1 2 a 2 3 c o s 2 θ 3 ,

where Γr(r=18,..., 38) are constants, θ1=σ1T1-φ1,θ2=σ2T1-φ1 and θ3=σ3T1-φ1+2φ2. Steady-state solutions of the system correspond to the fixed points of Eqs. (21)-(22), which in turn correspond to:

(23)
φ ˙ 1 = σ 1 ,       φ ˙ 2 = ( σ 1 - σ 3 ) 2 .

Hence, the fixed points of Eqs. (21)-(22) are given by:

(24a)
- μ 1 2 a 1 + ρ 1 2 ω 1 s i n θ 1 - ρ 21 4 ω 1 c o s θ 2 - Γ 23 ω 1 a 1 2 a 2 2 s i n θ 3 + Γ 24 ω 1 a 1 a 2 4 s i n 2 θ 3 + Γ 25 ω 1 a 2 2 s i n θ 3
            + Γ 26 ω 1 a 2 4 s i n θ 3 + Γ 27 ω 1 a 1 2 a 2 2 s i n θ 3 = 0 ,
(24b)
a 1 σ 1 + Γ 18 ω 1 a 1 3 + Γ 19 ω 1 a 1 5 + Γ 20 ω 1 a 1 a 2 2 + Γ 21 ω 1 a 1 a 2 4 + Γ 22 ω 1 a 1 3 a 2 2 + ρ 1 2 ω 1 c o s θ 1 + ρ 21 4 ω 1 s i n θ 2
            + Γ 23 ω 1 a 1 2 a 2 2 c o s θ 3 + Γ 24 ω 1 a 1 a 2 4 c o s 2 θ 3 + Γ 25 ω 1 a 2 2 c o s θ 3 + Γ 26 ω 1 a 2 4 c o s θ 3
            + Γ 27 ω 1 a 1 2 a 2 2 c o s θ 3 = 0 ,
(25a)
- ζ 2 2 a 2 - Γ 33 ω 2 a 1 a 2 s i n θ 3 - Γ 34 ω 2 a 1 a 2 3 s i n θ 3 - Γ 35 ω 2 a 1 3 a 2 s i n θ 3 - Γ 36 ω 2 a 1 a 2 3 s i n θ 3
            + Γ 37 ω 2 a 1 a 2 3 s i n θ 3 - Γ 38 ω 2 a 1 2 a 2 3 s i n 2 θ 3 = 0 ,
(25b)
a 2 ( σ 1 - σ 3 ) 2 + Γ 28 ω 2 a 2 3 + Γ 29 ω 2 a 2 5 + Γ 30 ω 2 a 1 2 a 2 + Γ 31 ω 2 a 1 4 a 2 + Γ 32 ω 2 a 1 2 a 2 3 + Γ 33 ω 2 a 1 a 2 c o s θ 3
            + Γ 34 ω 2 a 1 a 2 3 c o s θ 3 + Γ 35 ω 2 a 1 3 a 2 c o s θ 3 + Γ 36 ω 2 a 1 a 2 3 c o s θ 3 + Γ 37 ω 2 a 1 a 2 3 c o s θ 3
            + Γ 38 ω 2 a 1 2 a 2 3 c o s 2 θ 3 = 0 .

Solving the resulting algebraic equations for the fixed points for practical case a10, a20 we obtained:

(26a)
( a 1 2 ) σ 1 2 + 2 Γ 22 a 1 4 a 2 2 ω 1 + 2 Γ 21 a 1 2 a 2 4 ω 1 + 2 Γ 20 a 1 2 a 2 2 ω 1 + 2 Γ 19 a 1 6 ω 1 + 2 Γ 18 a 1 4 ω 1 σ 1
            + μ 1 2 a 1 2 4 + Γ 18 2 a 1 6 ω 1 2 + Γ 19 2 a 1 10 ω 1 2 + Γ 20 2 a 1 2 a 2 4 ω 1 2 + Γ 21 2 a 1 2 a 2 8 ω 1 2 + Γ 22 2 a 1 6 a 2 4 ω 1 2 - Γ 23 2 a 1 4 a 2 4 ω 1 2 - Γ 24 2 a 1 2 a 2 8 ω 1 2
            - Γ 25 2 a 2 4 ω 1 2 - Γ 26 2 a 2 8 ω 1 2 - Γ 27 2 a 1 4 a 2 8 ω 1 2 + 2 Γ 18 Γ 19 a 1 8 ω 1 2 + 2 Γ 18 Γ 20 a 1 4 a 2 2 ω 1 2 + 2 Γ 18 Γ 21 a 1 4 a 2 4 ω 1 2
            + 2 Γ 18 Γ 22 a 1 6 a 2 2 ω 1 2 + 2 Γ 19 Γ 20 a 1 6 a 2 2 ω 1 2 + 2 Γ 19 Γ 21 a 1 6 a 2 4 ω 1 2 + 2 Γ 19 Γ 22 a 1 8 a 2 2 ω 1 2 + 2 Γ 20 Γ 21 a 1 2 a 2 6 ω 1 2 + 2 Γ 20 Γ 22 a 1 4 a 2 4 ω 1 2 + 2 Γ 21 Γ 22 a 1 4 a 2 6 ω 1 2 - 2 Γ 23 Γ 24 a 1 3 a 2 6 ω 1 2 - 2 Γ 23 Γ 25 a 1 2 a 2 4 ω 1 2 - 2 Γ 23 Γ 26 a 1 2 a 2 6 ω 1 2
            - 2 Γ 23 Γ 27 a 1 4 a 2 4 ω 1 2 - 2 Γ 24 Γ 25 a 1 a 2 6 ω 1 2 - 2 Γ 24 Γ 26 a 1 a 2 8 ω 1 2 - 2 Γ 24 Γ 27 a 1 3 a 2 6 ω 1 2 - 2 Γ 25 Γ 26 a 2 6 ω 1 2
            - 2 Γ 25 Γ 27 a 1 2 a 2 4 ω 1 2 - 2 Γ 26 Γ 27 a 1 2 a 2 6 ω 1 2 - ρ 1 2 4 ω 1 2 - ρ 21 2 16 ω 1 2 - ρ 1 Γ 23 a 1 2 a 2 2 ω 1 2 - ρ 1 Γ 24 a 1 a 2 4 ω 1 2
            - ρ 1 Γ 25 a 2 2 ω 1 2 - ρ 1 Γ 26 a 2 4 ω 1 2 - ρ 1 Γ 27 a 1 2 a 2 2 ω 1 2 = 0 ,
(26b)
a 2 2 4 σ 3 2 + - a 2 2 σ 1 2 - Γ 28 a 2 4 ω 2 - Γ 29 a 2 6 ω 2 - Γ 30 a 1 2 a 2 2 ω 2 - Γ 31 a 1 4 a 2 2 ω 2 - Γ 32 a 1 2 a 2 4 ω 2 σ 3
            + μ 2 2 a 2 2 4 + σ 1 2 a 2 2 4 + Γ 28 2 a 2 6 ω 2 2 + Γ 29 2 a 2 10 ω 2 2 + Γ 30 2 a 1 4 a 2 2 ω 2 2 + Γ 31 2 a 1 8 a 2 2 ω 2 2 + Γ 32 2 a 1 4 a 2 6 ω 2 2 - Γ 33 2 a 1 2 a 2 2 ω 2 2
            - Γ 34 2 a 1 2 a 2 6 ω 2 2 - Γ 35 2 a 1 6 a 2 2 ω 2 2 - Γ 36 2 a 1 2 a 2 6 ω 2 2 - Γ 37 2 a 1 2 a 2 6 ω 2 2 - Γ 38 2 a 1 4 a 2 6 ω 2 2 + Γ 28 σ 1 a 2 4 ω 2 + Γ 29 σ 1 a 2 6 ω 2
            + Γ 30 σ 1 a 1 2 a 2 2 ω 2 + Γ 31 σ 1 a 1 4 a 2 2 ω 2 + Γ 32 σ 1 a 1 2 a 2 4 ω 2 + 2 Γ 28 Γ 29 a 2 8 ω 2 2 + 2 Γ 28 Γ 30 a 1 2 a 2 8 ω 2 2 + 2 Γ 28 Γ 31 a 1 4 a 2 4 ω 2 2
            + 2 Γ 28 Γ 32 a 1 2 a 2 6 ω 2 2 + 2 Γ 29 Γ 30 a 1 2 a 2 6 ω 2 2 + 2 Γ 29 Γ 31 a 1 4 a 2 6 ω 2 2 + 2 Γ 29 Γ 32 a 1 2 a 2 8 ω 2 2 + 2 Γ 30 Γ 31 a 1 6 a 2 2 ω 2 2
            + 2 Γ 30 Γ 32 a 1 4 a 2 4 ω 2 2 + 2 Γ 31 Γ 32 a 1 6 a 2 4 ω 2 2 - 2 Γ 33 Γ 34 a 1 2 a 2 4 ω 2 2 - 2 Γ 33 Γ 35 a 1 4 a 2 2 ω 2 2 - 2 Γ 33 Γ 36 a 1 2 a 2 4 ω 2 2
            - 2 Γ 33 Γ 37 a 1 2 a 2 4 ω 2 2 - 2 Γ 33 Γ 38 a 1 3 a 2 4 ω 2 2 - 2 Γ 34 Γ 35 a 1 4 a 2 4 ω 2 2 - 2 Γ 34 Γ 36 a 1 2 a 2 6 ω 2 2 - 2 Γ 34 Γ 37 a 1 2 a 2 6 ω 2 2
            - 2 Γ 34 Γ 38 a 1 3 a 2 6 ω 2 2 - 2 Γ 35 Γ 36 a 1 4 a 2 4 ω 2 2 - 2 Γ 35 Γ 37 a 1 4 a 2 4 ω 2 2 - 2 Γ 35 Γ 38 a 1 3 a 2 6 ω 2 2 - 2 Γ 36 Γ 37 a 1 2 a 2 6 ω 2 2
            - 2 Γ 36 Γ 38 a 1 3 a 2 6 ω 2 2 - 2 Γ 37 Γ 38 a 1 3 a 2 6 ω 2 2 = 0 .

In the frequency response curves, the stable (unstable) steady-state solutions have been represented by solid (dashed) lines.

3.1. Stability of non-linear solution

To determine the stability of the fixed points, one lets:

(27)
a 1 =   a 10 +   a 11 ,       a 2 =   a 20 +   a 21 ,       θ S = θ S 0 + θ S 1       a n d       θ 3 = θ 30 + θ 31 ,       s   =   1 ,   2 ,

where a10, a20, θS0 and θ30 are the solutions of Eqs. (24)-(25) and where a11, a21, θS1 and θ31 are perturbations which are assumed to be small compared to where a10, a20, θS0 and θ30. Substituting Eq. (27) into Eqs. (21)-(22), using Eqs. (24)-(25) and keeping only the linear terms in where a11, a21, θS1 and θ31 we obtain:

(28a)
a 11 ' = - μ 1 2 - 2 Γ 23 ω 1 a 10 a 20 2 s i n θ 30 + Γ 24 ω 1 a 20 4 s i n 2 θ 30 + 2 Γ 27 ω 1 a 10 a 20 2 s i n θ 30 a 11
            + ( ρ 1 2 ω 1 c o s θ S 0 + ρ 21 4 ω 1 s i n θ S 0 ) θ S 1 + - 2 Γ 23 ω 1 a 10 2 a 20 s i n θ 30 + 4 Γ 24 ω 1 a 10 a 20 3 s i n 2 θ 30
            + 2 Γ 25 ω 1 a 20 s i n θ 30 + 4 Γ 26 ω 1 a 20 3 s i n θ 30 + 2 Γ 27 ω 1 a 10 2 a 20 s i n θ 30 a 21
            + - 2 Γ 23 ω 1 a 10 2 a 20 2 c o s θ 30 + 2 Γ 24 ω 1 a 10 a 20 4 c o s 2 θ 30 + Γ 25 ω 1 a 20 2 c o s θ 30 + Γ 26 ω 1 a 20 4 c o s θ 30
            + Γ 27 ω 1 a 10 2 a 20 2 c o s θ 30 θ 31 ,
(28b)
θ S 1 ' = σ a 10 + 3 Γ 18 ω 1 a 10 + 5 Γ 19 ω 1 a 10 3 + Γ 20 a 10 ω 1 a 20 2 + Γ 21 a 10 ω 1 a 20 4 + 3 Γ 22 ω 1 a 10 a 20 2
            + 2 Γ 23 ω 1 a 20 2 c o s θ 30 + Γ 24 a 10 ω 1 a 20 4 c o s 2 θ 30 + 2 Γ 27 ω 1 a 20 2 c o s θ 30 a 11 + ρ 21 4 ω 1 a 10 c o s θ S 0
            - ρ 1 2 ω 1 a 10 s i n θ S 0 θ S 1 + 2 Γ 20 ω 1 a 20 + 4 Γ 21 ω 1 a 20 3 + 2 Γ 22 ω 1 a 10 2 a 20 + 2 Γ 23 ω 1 a 10 a 20 c o s θ 30
            + 4 Γ 24 ω 1 a 20 3 c o s 2 θ 30 + 2 Γ 25 a 10 ω 1 a 20 c o s θ 30 + 4 Γ 26 a 10 ω 1 a 20 3 c o s θ 30 + 2 Γ 27 ω 1 a 10 a 20 c o s θ 30 a 21
            - Γ 23 ω 1 a 10 a 20 2 s i n θ 30 + Γ 24 ω 1 a 20 4 s i n 2 θ 30 + Γ 25 a 10 ω 1 a 20 2 s i n θ 30 + Γ 26 a 10 ω 1 a 20 4 s i n θ 30
            + Γ 27 ω 1 a 10 a 20 2 s i n θ 30 θ 31 ,
(29a)
a 21 ' = - Γ 33 ω 2 a 20 s i n θ 30 - Γ 34 ω 2 a 20 3 s i n θ 30 - 2 Γ 35 ω 2 a 20 a 10 2 s i n θ 30 - Γ 36 ω 2 a 20 3 s i n θ 30
            + Γ 37 ω 2 a 20 3 s i n θ 30 - 2 Γ 38 ω 2 a 10 a 20 3 s i n 2 θ 30 a 11 - μ 2 2 + Γ 33 ω 2 a 10 s i n θ 30
            + 3 Γ 34 ω 2 a 10 a 20 2 s i n θ 30 + Γ 35 ω 2 a 10 3 s i n θ 30 + 3 Γ 36 ω 2 a 10 a 20 2 s i n θ 30 + 3 Γ 37 ω 2 a 10 a 20 2 s i n θ 30
            + 3 Γ 38 ω 2 a 10 2 a 20 2 s i n 2 θ 30 a 21 - Γ 33 ω 2 a 10 a 20 c o s θ 30 + Γ 34 ω 2 a 10 a 20 3 c o s θ 30
            + Γ 35 ω 2 a 10 3 a 20 c o s θ 30 + Γ 36 ω 2 a 10 a 20 3 c o s θ 30 + Γ 37 ω 2 a 10 a 20 3 c o s θ 30
            + Γ 38 ω 2 a 10 2 a 20 3 c o s 2 θ 30 θ 31 ,
(29b)
θ 31 ' = - 4 Γ 30 ω 2 a 10 + 8 Γ 31 ω 1 a 10 3 + 4 Γ 32 ω 2 a 10 a 20 2 + 2 Γ 33 ω 2 c o s θ 30 + 2 Γ 34 ω 2 a 20 2 c o s θ 30
            + 6 Γ 35 ω 2 a 10 2 c o s θ 30 + 2 Γ 36 ω 2 a 20 2 c o s θ 30 + 2 Γ 37 ω 2 a 20 2 c o s θ 30 + 4 Γ 38 ω 2 a 10 a 20 2 c o s θ 30 a 11
            - σ - σ 3 a 20 + 6 Γ 28 ω 2 a 20 + 10 Γ 29 ω 2 a 20 3 + 2 Γ 30 a 20 ω 2 a 10 2 + 6 Γ 32 ω 2 a 10 2 a 20 + 2 Γ 33 a 20 ω 2 a 10 c o s θ 30
            + 6 Γ 34 ω 2 a 10 a 20 c o s θ 30 + 2 Γ 35 a 20 ω 2 a 10 3 c o s θ 30 + 6 Γ 36 ω 2 a 10 a 20 c o s θ 30 + 6 Γ 37 ω 2 a 10 a 20 c o s θ 30
            + 6 Γ 38 ω 2 a 10 2 a 20 c o s 2 θ 30 a 21 + 2 Γ 33 ω 2 a 10 s i n θ 30 + 2 Γ 34 ω 2 a 10 a 20 2 s i n 2 θ 30
            + 2 Γ 35 ω 2 a 10 3 s i n θ 30 + 2 Γ 36 ω 2 a 10 a 20 2 s i n θ 30 + 2 Γ 37 ω 2 a 10 a 20 2 s i n θ 30
            + 4 Γ 38 ω 2 a 10 2 a 20 2 s i n 2 θ 30 θ 31 .

To study the stability of the fixed points corresponding to the practical case, we let a110, θS10, a210 and θ310 in Eqs. (28)-(29), and obtain the eigenvalues from the Jacobian matrix of the right hand sides. The zeros of the characteristic equation are given by:

(30)
λ 4 + L 1 λ 3 + L 2 λ 2 + L 3 λ + L 4 = 0 ,

where, L1, L2, L3 and L4 are functions of the parameters (a1, a2,ω1,ω2, μ1, μ2, β1, β2, β3, β4, β5, η1, η2, η3, η4, η5, ρ1, ρ2, ρ11, ρ12, ρ21, ρ22, θs, θ3, σ, σ3). According to the Routh-Hurwitz criterion the necessary and sufficient conditions for all the roots of Eq. (30) to posses negative real parts, such that:

(31)
L 1 > 0 ,       L 1 L 2 - L 3 > 0 ,       L 3 L 1 L 2 - L 3 - L 1 2 L 4 > 0 ,       L 4 > 0 .

4. Numerical results

Fig. 2. Non-resonance system behavior (basic case)

 Non-resonance system behavior (basic case)
 Non-resonance system behavior (basic case)
 Non-resonance system behavior (basic case)
 Non-resonance system behavior (basic case)

To determine the numerical solution and response of the given system of Eq. (4), the Runge-Kutta of fourth order method was applied. Fig. 2 illustrates the response and the phase-plane for the non-resonant system (basic case) where γ1γ2γ3γ4ω1ω2 at some practical values of the equation parameters μ1=0.02, β1=0.03, β2=0.3, β3=0.2, β4=0.05, β5=0.3, μ2=0.02, η1=0.3, η2=0.03, η3=0.5, η4=0.03, η5=0.3, ρ1=8.6, ρ2=8.6, ρ11=0.02, ρ12=0.02, ρ21=0.005, ρ22=0.005, η1=0.3, η2=0.03, η3=0.5, η4=0.03, η5=0.3, ρ1=8.6, ρ2=8.6, ρ11=0.02, ρ12=0.02, ρ21=0.005, ρ22=0.005. It is observed from this figure that the response of the first and second modes of the FGM rectangular plate start with increasing amplitude with some chaotic and tuned oscillation respectively, the oscillation of the two modes becomes stable and the steady state amplitudes x and y are about 0.35 and 0.4 respectively and the phase plane shows limit cycle, denoting that the system is free from chaos. Table 1 shows the results of the worst resonance conditions.

Fig. 3 shows that the time response and phase-plane of the simultaneous primary, combined and internal resonance (γ1ω1, γ3-γ4ω1 and ω21/2ω1) which is one of the worst resonance cases. From this figure we have that the amplitude of the first mode of the FGM rectangular plate is increased to about 1950 % of that values shown in Fig. 2, while the amplitude of the second mode is increased to about 1150 % and becomes stable and the phase plane shows multi-limit cycle.

Fig. 3. Simultaneous primary combined and internal resonance case: μ1=0.02, β1=0.03, β2=0.3, β3=0.2, β4=0.05, β5=0.3, μ2=0.02, η1=0.3, η2=0.03, η3=0.5, η4=0.03, η5=0.3, ρ1=8.6, ρ2=8.6, ρ11=0.02, ρ12=0.02, ρ21=0.005, ρ22=0.005 (γ1ω1, γ3-γ4ω1 and ω21/2ω1)

 Simultaneous primary combined and internal resonance case:  μ1=0.02, β1=0.03, β2=0.3, β3=0.2, β4=0.05, β5=0.3, μ2=0.02, η1=0.3, η2=0.03,  η3=0.5, η4=0.03, η5=0.3, ρ1=8.6, ρ2=8.6, ρ11=0.02, ρ12=0.02, ρ21=0.005,  ρ22=0.005 ( γ1≅ω1,  γ3-γ4≅ω1 and ω2≅1/2ω1)
 Simultaneous primary combined and internal resonance case:  μ1=0.02, β1=0.03, β2=0.3, β3=0.2, β4=0.05, β5=0.3, μ2=0.02, η1=0.3, η2=0.03,  η3=0.5, η4=0.03, η5=0.3, ρ1=8.6, ρ2=8.6, ρ11=0.02, ρ12=0.02, ρ21=0.005,  ρ22=0.005 ( γ1≅ω1,  γ3-γ4≅ω1 and ω2≅1/2ω1)
 Simultaneous primary combined and internal resonance case:  μ1=0.02, β1=0.03, β2=0.3, β3=0.2, β4=0.05, β5=0.3, μ2=0.02, η1=0.3, η2=0.03,  η3=0.5, η4=0.03, η5=0.3, ρ1=8.6, ρ2=8.6, ρ11=0.02, ρ12=0.02, ρ21=0.005,  ρ22=0.005 ( γ1≅ω1,  γ3-γ4≅ω1 and ω2≅1/2ω1)
 Simultaneous primary combined and internal resonance case:  μ1=0.02, β1=0.03, β2=0.3, β3=0.2, β4=0.05, β5=0.3, μ2=0.02, η1=0.3, η2=0.03,  η3=0.5, η4=0.03, η5=0.3, ρ1=8.6, ρ2=8.6, ρ11=0.02, ρ12=0.02, ρ21=0.005,  ρ22=0.005 ( γ1≅ω1,  γ3-γ4≅ω1 and ω2≅1/2ω1)

Table 1. Summary of the worst resonance cases

Resonance cases
x %
y %
Remarks
Non-resonance case
100 %
100 %
Limit cycle
γ 1 ω 1 ,   γ 1 ω 2 ,   ω 2 ω 1
1450 %
1000 %
Multi Limit cycle
γ 1 ω 1 ,   γ 2 2 ω 2 ,   ω 2 ω 1
1450 %
900 %
Multi Limit cycle
γ 1 ω 1 ,  γ 3 - γ 4 ω 2 ,  ω 2 ω 1
1450 %
1000 %
Multi Limit cycle
γ 1 ω 2 ,  γ 2 2 ω 1 ,  ω 2 ω 1
1450 %
900 %
Multi Limit cycle
γ 1 ω 2 ,  γ 3 - γ 4 ω 1 ,  ω 2 ω 1
1450 %
1000 %
Multi Limit cycle
γ 1 ω 1 ,  γ 1 ω 2 ,   γ 3 - γ 4 ω 1
1450 %
1000 %
Multi Limit cycle
γ 1 ω 1 ,  γ 2 2 ω 2 ,  γ 3 - γ 4 ω 1
1400 %
900 %
Multi Limit cycle
γ 1 ω 2 ,  γ 2 2 ω 1 ,  γ 3 - γ 4 ω 1
1450 %
800 %
Multi Limit cycle
γ 1 ω 1 ,  γ 3 - γ 4 ω 1 ,   ω 2 1 / 2 ω 1
1950 %
1150 %
Multi Limit cycle
γ 1 ω 1 ,  γ 3 - γ 4 ω 1 ,   ω 2 1 / 3 ω 1
1900 %
1100 %
Multi Limit cycle
γ 1 ω 1 ,  γ 3 - γ 4 ω 1 ,   ω 2 1 / 4 ω 1
1900 %
1100 %
Multi Limit cycle

4.1. Response curves and effects of different parameters

In this section, the steady state response of the given system at various parameters near the simultaneous primary, combined and internal resonance case is investigated and studied. The frequency response equations given by Eq. (26) are solved numerically at the same values of the parameters shown in Fig. 3. In all figures, the solid lines stand for the stable solution and the dashed lines for the unstable solution.

Fig. 4. The steady state amplitudes of the first mode of the FGM rectangular plate

The steady state amplitudes of the first mode of the FGM rectangular plate

a) Effects of the detuning parameter σ1

The steady state amplitudes of the first mode of the FGM rectangular plate

b) Effects of the damping coefficient μ1

The steady state amplitudes of the first mode of the FGM rectangular plate

c) Effect of the natural frequencies ω1, ω2

The steady state amplitudes of the first mode of the FGM rectangular plate

d) Effects of the non-linear parameter β1

The steady state amplitudes of the first mode of the FGM rectangular plate

e) Effects of the non-linear parameter β2

The steady state amplitudes of the first mode of the FGM rectangular plate

f) Effects of the non-linear parameter β3

The steady state amplitudes of the first mode of the FGM rectangular plate

g) Effects of the non-linear parameter β4

The steady state amplitudes of the first mode of the FGM rectangular plate

h) Effects of the non-linear parameter β5

The steady state amplitudes of the first mode of the FGM rectangular plate

i) Effects of the excitation amplitude ρ1

The steady state amplitudes of the first mode of the FGM rectangular plate

j) Effects of the excitation amplitude ρ21

Fig. 4(a) shows the steady state amplitudes of the first mode of the FGM rectangular plate against the detuning parameter σ1 at the practical case, where a10, a20. Fig. 4(b) shows that the steady state amplitude of the first mode of the FGM rectangular plate is a monotonic decreasing function in the linear damping coefficients μ1. Fig. 4(c) show that the steady state amplitude of the first mode is a monotonic decreasing function in the natural frequencies ω1 and ω2, in this Figure, the response curve is bent to the right and has harding spring type and there exists jump phenomena. Figs. 4(d), (e) shows that the steady state amplitude of the first mode is a monotonic decreasing function in nonlinear parameter β1,β2 and the response curves are bent to the left leading to the occurrence of the jump phenomena and multi-valued amplitudes. The steady state amplitude of the first mode is a monotonic increasing function in the nonlinear parameter β3, β4 and the excitation amplitudes ρ1 as shown in Figs. 4(f), (g), (i). For increasing value of the nonlinear parameter β5, the amplitude of the first mode is decreased and the curve is shifted to the left as shown in Fig. 4(h). The effect of increasing or decreasing the excitation amplitudes ρ21 on the steady state amplitude of the FGM is trivial due to saturation occurrence as shown in Fig. 4(j).

Fig. 5. The steady state amplitudes of the second mode of the FGM rectangular plate

 The steady state amplitudes of the second mode of the FGM rectangular plate

a) Effects of the detuning parameter σ3

 The steady state amplitudes of the second mode of the FGM rectangular plate

b) Effects of the damping coefficient μ2

 The steady state amplitudes of the second mode of the FGM rectangular plate

c) Effect of the natural frequencies ω1, ω2

 The steady state amplitudes of the second mode of the FGM rectangular plate

d) Effects of the non-linear parameter η1

 The steady state amplitudes of the second mode of the FGM rectangular plate

e) Effects of the non-linear parameter η2

 The steady state amplitudes of the second mode of the FGM rectangular plate

f) Effects of the non-linear parameter η3

 The steady state amplitudes of the second mode of the FGM rectangular plate

g) Effects of the non-linear parameter η4

 The steady state amplitudes of the second mode of the FGM rectangular plate

h) Effects of the non-linear parameter η5

Fig. 5(a) shows the steady state amplitudes of the second mode of the FGM rectangular plate against the detuning parameter σ3 at the practical case, where a10, a20. In this Figure, the response amplitude consists of a continuous curve which is bent to the right and there exists jump phenomena. Fig. 5(b) shows that the effect of increasing or decreasing the linear damping coefficients μ2 on the steady state amplitude of the second mode of FGM is trivial due to saturation occurrence. Figs. 5(c), (d), (f), (g) show that the steady state amplitude of the second mode is a monotonic decreasing function in nonlinear parameter η1,η3, η4 and the natural frequencies ω1, ω2 and the response curves are bent to the right leading to the occurrence of the jump phenomena and multi-valued amplitudes. For increasing value of the nonlinear parameter η2, the amplitude of the second mode is shifted to the right as shown in Fig. 5(e). The steady state amplitude of the first mode is a monotonic increasing function in the nonlinear parameter η5, as shown in Fig. 5(h).

4.2. Comparison between the numerical and analytical simulation

Fig. 6. Comparison between numerical solution (using RKM) and analytical solution (using perturbation method) of the system at ρ1=8.6, ρ2=8.6

 Comparison between numerical solution (using RKM) and analytical solution  (using perturbation method) of the system at ρ1=8.6, ρ2=8.6

a)

 Comparison between numerical solution (using RKM) and analytical solution  (using perturbation method) of the system at ρ1=8.6, ρ2=8.6

b)

Fig. 7. Comparison between numerical solution (using RKM) and analytical solution (using perturbation method) of the system at ρ1=1.6, ρ2=1.6

 Comparison between numerical solution (using RKM) and analytical solution  (using perturbation method) of the system at ρ1=1.6, ρ2=1.6

a)

 Comparison between numerical solution (using RKM) and analytical solution  (using perturbation method) of the system at ρ1=1.6, ρ2=1.6

b)

In this section, the system of nonlinear differential equations given by Eqs. (21) and (22) are solved numerically at the simultaneous primary, combined and internal resonance case where γ1ω1, γ3-γ4ω1 and ω21/2ω1. Figs. 6 and 7 show the comparison between numerical integration for the system Eqs. 4(a), (b) and the amplitude-phase modulating Eqs. (21) and (22) at different values of excitation force amplitudes ρ1 and ρ2. We found that all predictions from analytical solutions are in good agreement with the numerical simulation.

5. Conclusions

Multiple time scale perturbation method is useful to determine approximate solutions for the coupled nonlinear differential equations describing the FGM rectangular plate system up to and including the second order approximation. Both the frequency response equations and the phase-plane technique are applied to study the stability of the system. The effect of the different parameters of the system is studied numerically. From the above study the following may be concluded:

1) The amplitude of the first mode FGM rectangular plate is increased to about 0.35, while the amplitude of the second mode is increased to about 0.4 at the non-resonant case (γ1γ2γ3γ4ω1ω2).

2) The amplitude of the first mode FGM rectangular plate is increased to about 1950 % of that values shown in Fig. 2, while the amplitude of the second mode is increased to about 1150 % at the worst resonance case (the simultaneous primary, combined and internal resonance γ1ω1, γ3-γ4ω1 and ω21/2ω1).

3) The amplitude of the first and second mode FGM rectangular plate are increased to about 1900 % and 1100 % respectively at the worst resonance case (γ1ω1, γ3-γ4ω1 and ω21/3ω1) and (γ1ω1, γ3-γ4ω1 and ω21/4ω1).

4) The steady state amplitudes of the FGM rectangular plate are a monotonic decreasing function the linear damping coefficients μ1, the nonlinear parameters (β1, β2, β5, η1, η3, η4) and the natural frequencies ω1 and ω2.

5) The steady state amplitudes of FGM rectangular plate are a monotonic increasing function in the nonlinear parameters (β3, β4, η5) and the excitation amplitude ρ1.

6) The effects of increasing or decreasing the linear damping coefficients μ2 and the excitation amplitude ρ21 on the steady state amplitude of the FGM are trivial due to saturation occurrence.

7) The analytical solutions are in good agreement with the numerical simulation at different values of excitation amplitude forces ρ1 and ρ2.

Acknowledgements

The author would like to thank the anonymous reviewers for their valuable comments and suggestions for improving the quality of this paper.

References

  1. Reddy J. N., Cheng Z. Q. Three-dimensional thermomechanical deformations of functionally graded rectangular plates. European Journal of Mechanics, Vol. 20, 2001, p. 841-855. [Search CrossRef]
  2. Reddy J. N. Analysis of functionally graded plates. International Journal of Numerical Methods Engineering, Vol. 47, 2000, p. 66 3-684. [Search CrossRef]
  3. Yang J., Kitipomchai S., Liew K. M. Large amplitude vibration of thermo-electro mechanically stressed FGM laminated plates. Computer Methods in Applied Mechanics and Engineering, Vol. 192, 2003, p. 3861-3885. [Search CrossRef]
  4. Huang X. L., Shen H. S. Nonlinear vibration and dynamic response of functionally graded plates in thermal environments. International Journal of Solids and Structures, Vol. 41, 2004, p. 2403-2407. [Search CrossRef]
  5. Kitipornchai S., Yang J., Liew M. K. Semi-analytical solution for nonlinear vibration of laminated FGM plates with geometric imperfections. International Journal of Solids and Structures, Vol. 41, 2004, p. 2235-2257. [Search CrossRef]
  6. Yang J., Huang X. L. Nonlinear transient response of functionally graded plates with general imperfections in thermal environments. Computer Methods in Applied Mechanics and Engineering, Vol. 196, 2007, p. 2619-2630. [Search CrossRef]
  7. Ye M., Lu J., Zhang W., Ding Q. Local and global nonlinear dynamics of a parametrically excited rectangular symmetric cross-ply laminated composite plate. Chaos Solitons and Fractals, Vol. 26, 2005, p. 195-213. [Search CrossRef]
  8. Haddadpour H., Navazi H. M., Shadmehri F. Nonlinear oscillations of a fluttering functionally graded plate. Composite Structure, Vol. 79, 2007, p. 242-250. [Search CrossRef]
  9. HaoY. X., Chen L. H., Zhang W., Lei J. G. Nonlinear oscillations, bifurcations and chaos of functionally graded materials plate. Journal of Sound and Vibration, Vol. 312, 2008, p. 862-892. [Search CrossRef]
  10. Zhang W., Yang J., Hao Y. X. Chaotic vibrations of an orthotropic FGM rectangular plate based on third-order shear deformation theory. Nonlinear Dynamics, Vol. 59, 2010, p. 619-660. [Search CrossRef]
  11. Yang J., Hao Y. X., Zhang W., Kitipornchai S. Nonlinear dynamic response of a functionally graded plate with a through-width surface crack. Nonlinear Dynamics, Vol. 59, 2010, p. 207-219. [Search CrossRef]
  12. Li S. B., Zhang W., Hao Y. X. Multi-pulse chaotic dynamics of a functionally graded material rectangular plate with one-to-one internal resonance. International Journal of Nonlinear Science and Numerical Simulation, Vol. 11, 2010, p. 351-362. [Search CrossRef]
  13. Malekzadeh P., Golbahar Haghighi M. R., Atashi M. M. Free vibration analysis of elastically supported functionally graded annular plates subjected to thermal environment. Meccanica, Vol. 46, Issue 5, 2011, p. 893-913. [Search CrossRef]
  14. Akbarzadeh A. H., Abbasi M., Hosseini S. K., Eslami M. R. Dynamic analysis of functionally graded plates using the hybrid Fourier-Laplace transform under thermomechanical loading. Meccanica, Vol. 46, Issue 6, 2011, p. 1373-1392. [Search CrossRef]
  15. Zhang W., Hao Y., Guo X., Chen L. Complicated nonlinear responses of a simply supported FGM rectangular plate under combined parametric and external excitations. Meccanica, Vol. 47, 2012, p. 985-1014. [Search CrossRef]
  16. Eissa M., Sayed M. A comparison between passive and active control of non-linear simple pendulum Part-I. Mathematical and Computational Applications, Vol. 11, 2006, p. 137-149. [Search CrossRef]
  17. Eissa M., Sayed M. A comparison between passive and active control of non-linear simple pendulum Part-II. Mathematical and Computational Applications, Vol. 11, 2006, p. 151-162. [Search CrossRef]
  18. Eissa M., Sayed M. Vibration reduction of a three DOF non-linear spring pendulum. Communication in Nonlinear Science and Numerical Simulation, Vol. 13, 2008, p. 465-488. [Search CrossRef]
  19. Hamed Y. S., El-Ganaini W. A., Kamel M. M. Vibration suppression in ultrasonic machining described by non-linear differential equations. Journal of Mechanical Science and Technology, Vol. 23, Issue 8, 2009, p. 2038-2050. [Search CrossRef]
  20. Hamed Y. S., EL-Ganaini W. A., Kamel M. M. Vibration suppression in multi-tool ultrasonic machining to multi-external and parametric excitations. Acta Mechanica Sinica, Vol. 25, 2009, p. 403-415. [Search CrossRef]
  21. Hamed Y. S., El-Ganaini W. A., Kamel M. M. Vibration reduction in ultrasonic machine to external and tuned excitation forces. Applied Mathematical Modeling, Vol. 33, 2009, p. 2853-2863. [Search CrossRef]
  22. Sayed M., Hamed Y. S. Stability and response of a nonlinear coupled pitch-roll ship model under parametric and harmonic excitations. Nonlinear Dynamics, Vol. 64, 2011, p. 207-220. [Search CrossRef]
  23. Hamed Y. S., Sayed M., Cao D.-X., Zhang W. Nonlinear study of the dynamic behavior of a string-beam coupled system under combined excitations. Acta Mechanica Sinica, Vol. 27, Issue 6, 2011, p. 1034-1051. [Search CrossRef]
  24. Sayed M., Hamed Y. S., Amer Y. A. Vibration reduction and stability of non-linear system subjected to external and parametric excitation forces under a non-linear absorber. International Journal of Contemporary Mathematical Sciences, Vol. 6, Issue 22, 2011, p. 1051-1070. [Search CrossRef]
  25. Amer Y. A., Sayed M. Stability at principal resonance of multi-parametrically and externally excited mechanical system. Advances in Theoretical and Applied Mechanics, Vol. 4, 2011, p. 1-14. [Search CrossRef]
  26. Sayed M., Kamel M. Stability study and control of helicopter blade flapping vibrations. Applied Mathematical Modelling, Vol. 35, 2011, p. 2820-2837. [Search CrossRef]
  27. Sayed M., Kamel M. 1:2 and 1:3 internal resonance active absorber for non-linear vibrating system. Applied Mathematical Modelling, Vol. 36, 2012, p. 310-332. [Search CrossRef]
  28. Sayed M., Mousa A. A. Second-order approximation of angle-ply composite laminated thin plate under combined excitations. Communication in Nonlinear Science and Numerical Simulation, Vol. 17, 2012, p. 5201-5216. [Search CrossRef]
  29. Sayed M., Mousa A. A. Vibration, stability, and resonance of angle-ply composite laminated rectangular thin plate under multi excitations. Mathematical Problems in Engineering, Vol. 2013, 2013, p. 1-26. [Search CrossRef]
  30. Reddy N. A refined nonlinear theory of plates with transverse shear deformation. International Journal of Solids and Structures, Vol. 20, Issue 9-10, 1984, p. 881-896. [Search CrossRef]