# Vibration, Stability, and Resonance of Angle-Ply Composite Laminated Rectangular Thin Plate under Multiexcitations

Authors: M. Sayed; A. A. Mousa
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

An analytical investigation of the nonlinear vibration of a symmetric cross-ply composite laminated piezoelectric rectangular plate under parametric and external excitations is presented. The method of multiple time scale perturbation is applied to solve the nonlinear differential equations describing the system up to and including the second-order approximation. All possible resonance cases are extracted at this approximation order. The case of 1 : 1 : 3 primary and internal resonance, whereΩ3≅ω1, ω2≅ω1 and ω3≅3ω1, is considered. The stability of the system is investigated using both phase-plane method and frequency response curves. The influences of the cubic terms on nonlinear dynamic characteristics of the composite laminated piezoelectric rectangular plate are studied. The analytical results given by the method of multiple time scale is verified by comparison with results from numerical integration of the modal equations. Reliability of the obtained results is verified by comparison between the finite difference method (FDM) and Runge-Kutta method (RKM). It is quite clear that some of the simultaneous resonance cases are undesirable in the design of such system. Such cases should be avoided as working conditions for the system. Variation of the parameters μ1,μ2, α7,β8, ω1,ω2, f1, and f2 leads to multivalued amplitudes and hence to jump phenomena. Some recommendations regarding the different parameters of the system are reported. Comparison with the available published work is reported.


---

## Body

## 1. Introduction

Composite laminated plates that are widely used in several engineering fields such as machinery, shipbuilding, aircraft, automobiles, robot arm, watercraft-hydropower, and wings of helicopters are made of the angle-ply composite laminated plates. Several researchers have focused their attention on studying the nonlinear dynamics, bifurcations, and chaos of the composite laminated plates.Internal resonance has been found in many engineering problems in which the natural frequencies of the system are commensurable. Ye et al. [1] investigated the local and global nonlinear dynamics of a parametrically excited symmetric cross-ply composite laminated rectangular thin plate under parametric excitation. The study is focused on the case of 1 : 1 internal resonance and primary parametric resonance. Zhang [2] dealt with the global bifurcations and chaotic dynamics of a parametrically excited, simply supported rectangular thin plate. The method of multiple scales is used to obtain the averaged equations. The case of 1 : 1 internal resonance and primary parametric resonance is considered. Guo et al. [3] studied the nonlinear dynamics of a four-edge simply supported angle-ply composite laminated rectangular thin plate excited by both the in-plane and transverse loads. The asymptotic perturbation method is used to derive the four averaged equations under 1 : 1 internal resonance. Zhang et al. [4] investigated the local and global bifurcations of a parametrically and externally excited simply supported rectangular thin plate under simultaneous transversal and in-plane excitations. The studies are focused on the case of 1 : 1 internal resonance and primary parametric resonance. Tien et al. [5] applied the averaging method and Melnikov technique to study local, global bifurcations and chaos of a two-degrees-of-freedom shallow arch subjected to simple harmonic excitation for case of 1 : 2 internal resonances. Sayed and Mousa [6] investigated the influence of the quadratic and cubic terms on nonlinear dynamic characteristics of the angle-ply composite laminated rectangular plate with parametric and external excitations. Two cases of the subharmonic resonances cases in the presence of 1 : 2 internal resonances are considered. The method of multiple time scale perturbation is applied to solve the nonlinear differential equations describing the system up to and including the second-order approximation. Zhang et al. [7] gave further studies on the nonlinear oscillations and chaotic dynamics of a parametrically excited simply supported symmetric cross-ply laminated composite rectangular thin plate with the geometric nonlinearity and nonlinear damping. Zhang et al. [8] dealt with the nonlinear vibrations and chaotic dynamics of a simply supported orthotropic functionally graded material (FGM) rectangular plate subjected to the in-plane and transverse excitations together with thermal loading in the presence of 1 : 2 : 4 internal resonance, primary parametric resonance, and subharmonic resonance of order 1/2. Zhang et al. [9] investigated the bifurcations and chaotic dynamics of a simply supported symmetric cross-ply composite laminated piezoelectric rectangular plate subject to the transverse, in-plane excitations and the excitation loaded by piezoelectric layers. Zhang and Li [10] analyzed the resonant chaotic motions of a simply supported rectangular thin plate with parametrically and externally excitations using exponential dichotomies and an averaging procedure. Zhang et al. [11] analyzed the chaotic dynamics of a six-dimensional nonlinear system which represents the averaged equation of a composite laminated piezoelectric rectangular plate subjected to the transverse, in-plane excitations and the excitation loaded by piezoelectric layers. The case of 1 : 2 : 4 internal resonances is considered. Zhang and Hao [12] studied the global bifurcations and multipulse chaotic dynamics of the composite laminated piezoelectric rectangular plate by using the improved extended Melnikov method. The multipulse chaotic motions of the system are found by using numerical simulation, which further verifies the result of theoretical analysis. Guo and Zhang [13] studied the nonlinear oscillations and chaotic dynamics for a simply supported symmetric cross-ply composite laminated rectangular thin plate with parametric and forcing excitations.  The case of 1 : 2 : 3 internal resonance is considered. The method of multiple scales is employed to obtain the six-dimensional averaged equation. The numerical method is used to investigate the periodic and chaotic motions of the composite laminated rectangular thin plate. Eissa and Sayed [14–16] and Sayed [17] investigated the effects of different active controllers on simple and spring pendulum at the primary resonance via negative velocity feedback or its square or cubic. Sayed and Kamel [18, 19] investigated the effects of different controllers on the vibrating system and the saturation control to reduce vibrations due to rotor blade flapping motion. The stability of the system is investigated using both phase-plane method and frequency response curves. Sayed and Hamed [20] 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 nonlinear differential equations and obtain approximate solutions up to and including the second-order approximations. Amer et al. [21] investigated the dynamical system of a twin-tail aircraft, which is described by two coupled second-order nonlinear differential equations having both quadratic and cubic nonlinearities under different controllers. Best active control of the system has been achieved via negative acceleration feedback. The stability of the system is investigated applying both frequency response equations and phase-plane method. Hamed et al. [22] studied the nonlinear dynamic behavior of a string-beam coupled system subjected to external, parametric, and tuned excitations that are presented. The case of 1 : 1 internal resonance between the modes of the beam and string, and the primary and combined resonance for the beam is considered. The method of multiple scales is applied to obtain approximate solutions up to and including the second-order approximations. All resonance cases are extracted and investigated. Stability of the system is studied using frequency response equations and the phase-plane method. Awrejcewicz et al. [23–25] studied the chaotic dynamics of continuous mechanical systems such as flexible plates and shallow shells. The considered problems are solved by the Bubnov-Galerkin, Ritz method with higher approximations, and finite difference method. Convergence and validation of those methods are studied. Awrejcewicz et al. [26] investigated the chaotic vibrations of flexible nonlinear Euler-Bernoulli beams subjected to harmonic load and with various boundary conditions. Reliability of the obtained results is verified by the finite difference method and finite element method with the Bubnov-Galerkin approximation for various boundary conditions and various dynamic regimes.In this paper, the perturbation method and stability of the composite laminated piezoelectric rectangular plate under simultaneous transverse and in-plane excitations are investigated. The method of multiple scales are applied to obtain the second-order uniform asymptotic solutions. All possible resonance cases are extracted at this approximation order. The study is 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. Such cases should be avoided as working conditions for the system. Some recommendations regarding the different parameters of the system are reported. Comparison with the available published work is reported.


## 2. Mathematical Analysis

Consider a simply supported four edges composite laminated piezoelectric rectangular plate of lengthsa,b and thickness h, as shown in Figure 1. The composite laminated piezoelectric rectangular plate is considered as regular symmetric cross-ply laminates with n layers. A Cartesian coordinate system is located in the middle surface of the plate. Assume that u,  v, and w represent the displacements of an arbitrary point of the composite laminated piezoelectric rectangular plate in the x, y, and z directions, respectively. The in-plane excitations are loaded along the ydirection at x=0 in the form q0+qxcosΩ1t, and the excitations are loaded along the x direction at y=0 in the form q1+qycosΩ2t. The transverse excitation subjected to the composite laminated piezoelectric rectangular plate is represented by qcosΩ3t. The dynamic electrical loading is expressed as EzcosΩ4t. Based on Reddy's third-order shear deformation plate theory [27], the displacement field at an arbitrary point in the composite laminated plate is expressed as [13](1a)u(x,y,t)=u0(x,y,t)+zφx(x,y,t)-z343h2(φx+∂w0∂x),(1b)v(x,y,t)=v0(x,y,t)+zφy(x,y,t)-z343h2(φy+∂w0∂y),(1c)w(x,y,t)=w0(x,y,t),Figure 1
The model of the composite laminated piezoelectric rectangular plate.whereu0, v0, and w0 are the original displacement on the midplane of the plate in the x,y, and z directions, respectively. Let φx and φy represent the midplane rotations of transverse normal about the x and y axes, respectively. From the van Karman-type plate theory and Hamilton's principle, the nonlinear governing equations of motion of the composite laminated piezoelectric rectangular plate are given as follows [12]:(2a)∂Nxx∂x+∂Nxy∂y=I0u¨0+J1φ¨x-43h2I3∂w¨0∂x,(2b)∂Nxy∂x+∂Nyy∂y=I0v¨0+J1φ¨y-43h2I3∂w¨0∂y,(2c)∂Q-x∂x+∂Q-y∂y+∂∂x(Nxx∂w0∂x+Nxy∂w0∂y)+∂∂y(Nxy∂w0∂x+Nyy∂w0∂y)+43h2(∂2Pxx∂x2+2∂2Pxy∂x∂y+∂2Pyy∂y2)+q=I0w¨0-169h4I6(∂2w¨0∂x2+∂2w¨0∂y2)+43h2[I3(∂u¨0∂x+∂v¨0∂y)+J4(∂φ¨x∂x+∂φ¨y∂y)],(2d)∂M-xx∂x+∂M-xy∂y-Q-x=J1u¨0+k2φ¨x-43h2J4∂w¨0∂x,(2e)∂M-xy∂x+∂M-yy∂y-Q-y=J1v¨0+k2φ¨y-43h2J4∂w¨0∂y.The boundary conditions of the simply supported rectangular composite laminated plate are expressed as follows:(3a)atx=0,x=a:∂u∂x=0,w=u=φx=φy=Mxx=Nxy=Q-x=0,(3b)aty=0,y=b:∂v∂y=0,w=v=Nxy=Myy=Q-y=0,(3c)∫0bNxx|x=0,ady=∫0b(q0+qxcosΩ1t)dy,(3d)∫0aNyy|y=0,bdx=∫0a(q1+qycosΩ2t)dx.Applying the Galerkin procedure, we obtain that the dimensionless differential equations of motion for the simply supported symmetric cross-ply rectangular thin plate are shown as follows [11, 28]:(4a)u¨1+μ1u˙1+ω12u1+(f11cosΩ1t+f12cosΩ2t+f14cosΩ4t)u1+α1u12u2+α2u12u3+α3u22u1+α4u22u3+α5u32u1+α6u32u2+α7u13+α8u23+α9u33+α10u1u2u3=f1cosΩ3t,(4b)u¨2+μ2u˙2+ω22u2+(f21cosΩ1t+f22cosΩ2t+f24cosΩ4t)u2+β1u12u2+β2u12u3+β3u22u1+β4u22u3+β5u32u1+β6u32u2+β7u13+β8u23+β9u33+β10u1u2u3=f2cosΩ3t,(4c)u¨3+μ3u˙3+ω32u3+(f31cosΩ1t+f32cosΩ2t+f34cosΩ4t)u3+γ1u12u2+γ2u12u3+γ3u22u1+γ4u22u3+γ5u32u1+γ6u32u2+γ7u13+γ8u23+γ9u33+γ10u1u2u3=f3cosΩ3t,whereu1,u2, and u3 are the vibration amplitudes of the composite laminated piezoelectric rectangular plate for the first-order, second-order, and the third-order modes, respectively, μ1, μ2, and μ3 are the linear viscous damping coefficients, ω1, ω2, and ω3 are the natural frequencies of the rectangular plate, and Ω1, Ω2, Ω3, and Ω4 are the excitations frequencies. fn1, fn2, fn3, and fn(n=1,2,3) are the amplitudes of parametric and external excitation forces corresponding to the three nonlinear modes, and αi, βi, and γi(i=1,2,…,10) are the nonlinear coefficients. The linear viscous damping and exciting forces are assumed to be
(5)μn=εμ^n,fn1=εf^n1,fn2=εf^n2,fn3=εf^n3,fn=ε2f^n,n=1,2,3,
where ε is a small perturbation parameter and  0<ε≪1.The external excitation forcesf^n are of the order 2, and the linear viscous damping μ^n, parametric exciting forces f^n1, f^n2, and f^n3 are of the order 1.To consider the influence of the cubic terms on nonlinear dynamic characteristics of the composite laminated piezoelectric rectangular plate, we need to obtain the second-order approximate solution of (4a), (4b), and (4c). Method of multiple scales [29–31] is applied to obtain a second-order approximation for the system. For the second-order approximation, we introduce three time scales defined by
(6)T0=t,T1=εt,T2=ε2t.
In terms of these scales, the time derivatives become(7a)ddt=D0+εD1+ε2D2,(7b)d2dt2=D02+2εD0D1+ε2(D12+2D0D2),whereDn=∂/∂Tn, n=0,1,2. We seek a uniform approximation to the solution of (4a), (4b), and (4c) in the form:
(8)un(t,ε)=εun1(T0,T1,T2)+ε2un2(T0,T1,T2)+ε3un3(T0,T1,T2)+O(ε4),n=1,2,3.
Terms of O(ε4) and higher orders are neglected. Substituting (5) and (7a)–(8) into (4a), (4b), and (4c) and equating the coefficients of like powers of ε, we obtain the following:Order  ε(9a)(D02+ω12)u11=0,(9b)(D02+ω22)u21=0,(9c)(D02+ω32)u31=0.Order  ε2(10a)(D02+ω12)u12=-2D0D1u11-μ^1D0u11-(f^11cosΩ1T0+f^12cosΩ2T0+f^14cosΩ4T0)u11+f^1cosΩ3T0,(10b)(D02+ω22)u22=-2D0D1u21-μ^2D0u21-(f^21cosΩ1T0+f^22cosΩ2T0+f^24cosΩ4T0)u21+f^2cosΩ3T0,(10c)(D02+ω32)u32=-2D0D1u31-μ^3D0u31-(f^31cosΩ1T0+f^32cosΩ2T0+f^34cosΩ4T0)u31+f^3cosΩ3T0.Order  ε3(11a)(D02+ω12)u13=-D12u11-2D0D2u11-2D0D1u12-μ^1(D0u12+D1u11)-(f^11cosΩ1T0+f^12cosΩ2T0+f^14cosΩ4T0)u12-α1u112u21-α2u112u31-α3u212u11-α4u212u31-α5u312u11-α6u312u21-α7u113-α8u213-α9u313-α10u11u21u31,(11b)(D02+ω22)u23=-D12u21-2D0D2u21-2D0D1u22-μ^2(D0u22+D1u21)-(f^21cosΩ1T0+f^22cosΩ2T0+f^24cosΩ4T0)u22-β1u112u21-β2u112u31-β3u212u11-β4u212u31-β5u312u11-β6u312u21-β7u113-β8u213-β9u313-β10u11u21u31,(11c)(D02+ω32)u33=-D12u31-2D0D2u31-2D0D1u32-μ^3(D0u32+D1u31)-(f^31cosΩ1T0+f^32cosΩ2T0+f^34cosΩ4T0)u32-γ1u112u21-γ2u112u31-γ3u212u11-γ4u212u31-γ5u312u11-γ6u312u21-γ7u113-γ8u213-γ9u313-γ10u11u21u31.The general solutions of (9a), (9b), and (9c) can be written in the form(12a)u11=A1(T1,T2)exp(iω1T0)+cc,(12b)u21=A2(T1,T2)exp(iω2T0)+cc,(12c)u31=A3(T1,T2)exp(iω3T0)+cc,whereA1, A2, and A3 are a complex function in T1, T2 which can be determined from eliminating the secular terms at the next approximation and cc stands for the complex conjugate of the preceding terms. Substituting (12a), (12b), and (12c) into (10a), (10b), and (10c) and eliminating the secular terms, then the first-order approximations are given by (13a)u12=E1exp(iω1T0)+E2exp(i(Ω1+ω1)T0)+E3exp(i(Ω1-ω1)T0)+E4exp(i(Ω2+ω1)T0)+E5exp(i(Ω2-ω1)T0)+E6exp(i(Ω4+ω1)T0)+E7exp(i(Ω4-ω1)T0)+cc,(13b)u22=E8exp(iω2T0)+E9exp(i(Ω1+ω2)T0)+E10exp(i(Ω1-ω2)T0)+E11exp(i(Ω2+ω2)T0)+E12exp(i(Ω2-ω2)T0)+E13exp(i(Ω4+ω2)T0)+E14exp(i(Ω4-ω2)T0)+cc,(13c)u32=E15exp(iω3T0)+E16exp(i(Ω1+ω3)T0)+E17exp(i(Ω1-ω3)T0)+E18exp(i(Ω2+ω3)T0)+E19exp(i(Ω2-ω3)T0)+E20exp(i(Ω4+ω3)T0)+E21exp(i(Ω4-ω3)T0)+E22exp(iΩ3T0)+cc,whereEi(i=1,2,…,22) are the complex functions in T1 and T2. From (12a)–(13c) into (11a), (11b), and (11c) and eliminating the secular terms, the second-order approximation is given by(14a)u13(T0,T1,T2)=H1exp(iω2T0)+H2exp(iω3T0)+H3exp(3iω1T0)+H4exp(3iω2T0)+H5exp(3iω3T0)+H6exp(i(ω2±2ω1)T0)+H7exp(i(ω3±2ω1)T0)+H8exp(i(2ω2±ω1)T0)+H9exp(i(2ω3±ω1)T0)+H10exp(i(ω2±2ω3)T0)+H11exp(i(ω3±2ω2)T0)+H12exp(i(ω3±ω2±ω1)T0)+H13exp(iΩ3T0)+H14exp(i(Ω3±Ω1)T0)+H15exp(i(Ω3±Ω2)T0)+H16exp(i(Ω3±Ω4)T0)+H17exp(i(Ω1±ω1)T0)+H18exp(i(Ω2±ω1)T0)+H19exp(i(Ω4±ω1)T0)+H20exp(i(2Ω1±ω1)T0)+H21exp(i(2Ω2±ω1)T0)+H22exp(i(2Ω4±ω1)T0)+H23exp(i(Ω2±Ω1±ω1)T0)+H24exp(i(Ω4±Ω1±ω1)T0)+H25exp(i(Ω4±Ω2±ω1)T0)+cc,(14b)u23(T0,T1,T2)=H26exp(iω1T0)+H27exp(iω3T0)+H28exp(3iω1T0)+H29exp(3iω2T0)+H30exp(3iω3T0)+H31exp(i(ω2±2ω1)T0)+H32exp(i(2ω2±ω1)T0)+H33exp(i(ω3±2ω1)T0)+H34exp(i(2ω3±ω1)T0)+H35exp(i(ω2±2ω3)T0)+H36exp(i(2ω2±ω3)T0)+H37exp(i(ω3±ω2±ω1)T0)+H38exp(iΩ3T0)+H39exp(i(Ω3±Ω1)T0)+H40exp(i(Ω3±Ω2)T0)+H41exp(i(Ω4±Ω3)T0)+H42exp(i(Ω1±ω2)T0)+H43exp(i(Ω2±ω2)T0)+H44exp(i(Ω4±ω2)T0)+H45exp(i(2Ω1±ω2)T0)+H46exp(i(2Ω2±ω2)T0)+H47exp(i(2Ω4±ω2)T0)+H48exp(i(Ω2±Ω1±ω2)T0)+H49exp(i(Ω4±Ω1±ω2)T0)+H50exp(i(Ω4±Ω2±ω2)T0)+cc,(14c)u33(T0,T1,T2)=H51exp(3iω3T0)+H52exp(iω1T0)+H53exp(3iω1T0)+H54exp(iω2T0)+H55exp(3iω2T0)+H56exp(i(ω2±2ω1)T0)+H57exp(i(2ω2±ω1)T0)+H58exp(i(ω3±2ω1)T0)+H59exp(i(2ω3±ω1)T0)+H60exp(i(ω2±2ω3)T0)+H61exp(i(2ω2±ω3)T0)+H62exp(i(ω3±ω2±ω1)T0)+H63exp(iΩ3T0)+H64exp(i(Ω3±Ω1)T0)+H65exp(i(Ω3±Ω2)T0)+H66exp(i(Ω3±Ω4)T0)+H67exp(i(Ω1±ω3)T0)+H68exp(i(Ω2±ω3)T0)+H69exp(i(Ω4±ω3)T0)+H70exp(i(2Ω1±ω3)T0)+H71exp(i(2Ω2±ω3)T0)+H72exp(i(2Ω4±ω3)T0)+H73exp(i(Ω2±Ω1±ω3)T0)+H74exp(i(Ω4±Ω1±ω3)T0)+H75exp(i(Ω4±Ω2±ω3)T0)+cc,whereHi(i=1,2,…,75) are the complex functions in T1 and T2. From the above derived solutions, the reported resonance cases are the following.(i)
Primary resonance:   Ω1≅ωn, Ω2≅ωn, Ω3≅ωn, Ω4≅ωn, and n=1,2,3.(ii)
Subharmonic resonance:   Ω1≅2ωn, Ω2≅2ωn, Ω4≅2ωn, and n=1,2,3.(iii)
Internal or secondary resonance: ω1≅ω2, ω2≅ω3, ω3≅ω1, ω1≅3ωs, ω2≅3ωr, ω3≅3ωm,  s=2,3, r=1,3, and m=1,2.(iv)
Combined resonance: ω3±ω2≅2ω1, ω3±ω1≅2ω2, ω1±ω2≅2ω3, ω2±2ω3≅ω1, ω3±2ω2≅ω1, ω3±2ω1≅ω2, ω1±2ω3≅ω2, ω2±2ω1≅ω3, ω1±2ω2≅ω3, Ω3±Ωt≅ωn, Ω4±Ωm≅2ωn, Ω2±Ω1≅2ωn, t=1,2,4,  m=1,2, and   n=1,2,3.(v)
Simultaneous or incident resonance.Any combination of the previous resonance cases is considered as simultaneous resonance.


## 3. Stability Analysis

The behavior of such a system can be very complex, especially when the natural frequencies and the forcing frequency satisfy certain internal and external resonance conditions. The study is focused on the case of 1 : 1 : 3 primary resonance and internal resonance, whereΩ3≅ω1, ω2≅ω1, and ω3≅3ω1. To describe how close the frequencies are to the resonance conditions, we introduce detuning parameters as follows:
(15)Ω3=ω1+σ1=ω1+εσ^1,ω2=ω1+σ2=ω1+εσ^2,ω3=3ω1+σ3=3ω1+εσ^3,
where σ1 and σ2, σ3 are called the external and internal detuning parameters, respectively. Eliminating the secular terms leads to solvability conditions for the first- and second-order expansions as follows:(16a)2iω1D1A1=-iμ^1ω1A1+f^12exp(iσ^1T1),(16b)2iω2D1A2=-iμ^2ω2A2+f^22exp(i(σ^1-σ^2)T1),(16c)2iω3D1A3=-iμ^3ω3A3,(17a)
2
i
ω
1
D
2
A
1
=
-
D
1
2
A
1
-
μ
^
1
D
1
A
1
-
{
f
^
11
2
2
(
Ω
1
2
-
4
ω
1
2
)
+
f
^
12
2
2
(
Ω
2
2
-
4
ω
1
2
)
+
f
^
14
2
2
(
Ω
4
2
-
4
ω
1
2
)
}
A
1
-
{
2
α
1
A
1
A
-
1
+
2
α
6
A
3
A
-
3
+
3
α
8
A
2
A
-
2
}
A
2
exp
(
i
σ
^
2
T
1
)
-
α
1
A
1
2
A
-
2
exp
(
-
i
σ
^
2
T
1
)
-
α
2
A
-
1
2
A
3
exp
(
i
σ
^
3
T
1
)
-
{
2
α
3
A
2
A
-
2
+
2
α
5
A
3
A
-
3
+
3
α
7
A
1
A
-
1
}
A
1
-
α
4
A
-
2
2
A
3
exp
(
i
(
σ
^
3
-
2
σ
^
2
)
T
1
)
-
α
10
A
-
1
A
-
2
A
3
exp
(
i
(
σ
^
3
-
σ
^
2
)
T
1
)
,
(17b)
2
i
ω
2
D
2
A
2
=
-
D
1
2
A
2
-
μ
^
2
D
1
A
2
-
{
f
^
21
2
2
(
Ω
1
2
-
4
ω
2
2
)
+
f
^
22
2
2
(
Ω
2
2
-
4
ω
2
2
)
+
f
^
24
2
2
(
Ω
4
2
-
4
ω
2
2
)
}
A
2
-
{
2
β
3
A
2
A
-
2
+
2
β
5
A
3
A
-
3
+
3
β
7
A
1
A
-
1
}
×
A
1
exp
(
-
i
σ
^
2
T
1
)
-
β
3
A
2
2
A
-
1
exp
(
i
σ
^
2
T
1
)
-
β
1
A
1
2
A
-
2
exp
(
-
2
i
σ
^
2
T
1
)
-
{
2
β
6
A
3
A
-
3
+
2
β
1
A
-
1
A
1
+
3
β
8
A
2
A
-
2
}
A
2
-
β
2
A
-
1
2
A
3
exp
(
i
(
σ
^
3
-
σ
^
2
)
T
1
)
-
β
10
A
-
1
A
-
2
A
3
exp
(
i
(
σ
^
3
-
2
σ
^
2
)
T
1
)
-
β
4
A
-
2
2
A
3
exp
(
i
(
σ
^
3
-
3
σ
^
2
)
T
1
)
,
(17c)
2
i
ω
3
D
2
A
3
=
-
D
1
2
A
3
-
μ
^
3
D
1
A
3
-
{
f
^
31
2
2
(
Ω
1
2
-
4
ω
3
2
)
+
f
^
32
2
2
(
Ω
2
2
-
4
ω
3
2
)
+
f
^
34
2
2
(
Ω
4
2
-
4
ω
3
2
)
}
A
3
-
γ
10
A
1
A
-
2
A
3
exp
(
-
i
σ
^
2
T
1
)
-
γ
10
A
-
1
A
2
A
3
exp
(
i
σ
^
2
T
1
)
-
γ
7
A
1
3
exp
(
-
i
σ
^
3
T
1
)
-
{
2
γ
2
A
1
A
-
1
+
2
γ
4
A
2
A
-
2
+
3
γ
9
A
3
A
-
3
}
A
3
-
γ
1
A
1
2
A
2
exp
(
i
(
σ
^
2
-
σ
^
3
)
T
1
)
-
γ
3
A
2
2
A
1
exp
(
i
(
2
σ
^
2
-
σ
^
3
)
T
1
)
-
γ
8
A
2
3
exp
(
i
(
3
σ
^
2
-
σ
^
3
)
T
1
)
.From (7a), multiplying both sides by 2iωn, we get
(18)2iωndAndt=ε2iωnD1An+ε22iωnD2An,2iωndAndt=ε2iωnD1Ann=1,2,3.To analyze the solutions of (16a)–(17c), we express An in the polar form as follows:
(19)An(T1,T2)=a^n2exp(iφn),an=εa^n,(n=1,2,3),
where an and φn are the steady-state amplitudes and phases of the motion, respectively. Substituting (16a)–(17c) and (19) into (18) and equating the real and imaginary parts, we obtain the following equations describing the modulation of the amplitudes and phases of the response:
(20)a˙1=-μ12a1+{f12ω1-σ1f14ω12}sinθ1+μ1f18ω12cosθ1-α18ω1a12a2sinθ2-α28ω1a12a3sinθ3-α48ω1a22a3sin(θ3-2θ2)-α64ω1a32a2sinθ2-3α88ω1a23sinθ2-α108ω1a1a2a3sin(θ3-θ2),a1φ˙1={-μ128ω1+Γ12ω1}a1-{f12ω1-σ1f14ω12}cosθ1+μ1f18ω12sinθ1+3α18ω1a12a2cosθ2+α28ω1a12a3cosθ3+α34ω1a1a22+α48ω1a22a3cos(θ3-2θ2)+α54ω1a1a32+α64ω1a32a2cosθ2+3α78ω1a13+3α88ω1a23cosθ2+α108ω1a1a2a3cos(θ3-θ2),a˙2=-μ22a2+{f22ω2-(σ1-σ2)f24ω22}sin(θ1-θ2)+μ2f28ω22cos(θ1-θ2)+β18ω2a12a2sin2θ2-β28ω2a12a3sin(θ3-θ2)+β38ω2a22a1sinθ2-β48ω2a22a3sin(θ3-3θ2)+β54ω2a32a1sinθ2+3β78ω2a13sinθ2-β108ω2a1a2a3sin(θ3-2θ2),a2φ˙2={-μ228ω2+Γ22ω2}a2+{(σ1-σ2)f24ω22-f22ω2}cos(θ1-θ2)+μ2f28ω22sin(θ1-θ2)+β18ω2a12a2cos2θ2+β14ω2a12a2+β28ω2a12a3cos(θ3-θ2)+3β38ω2a22a1cosθ2+β48ω2a22a3cos(θ3-3θ2)+β54ω2a32a1cosθ2+β64ω2a32a2+3β78ω2a13cosθ2+3β88ω2a23+β108ω2a1a2a3cos(θ3-2θ2),a˙3=-μ32a3-γ18ω3a12a2sin(θ2-θ3)-γ38ω3a22a1sin(2θ2-θ3)+γ78ω3a13sinθ3-γ88ω3a23sin(3θ2-θ3),a3φ˙3={-μ328ω3+Γ32ω3}a3+γ18ω3a12a2cos(θ2-θ3)+γ24ω3a12a3+γ38ω3a22a1cos(2θ2-θ3)+γ44ω3a22a3+γ78ω3a13cosθ3+γ88ω3a23cos(3θ2-θ3)+3γ98ω3a33+γ104ω3a1a2a3cosθ2,
where
(21)Γn={fn122(Ω12-4ωn2)+fn222(Ω22-4ωn2)+fn422(Ω42-4ωn2)},{fn122(Ω12-4ωn2)+fn422(Ω42-4ωn2)}n=1,2,3,θ1=σ^1T1-φ1,θ2=σ^2T1+φ2-φ1,θ3=σ^3T1+φ3-3φ1.
Steady-state solutions of the system correspond to the fixed points of (20), which in turn correspond to
(22)φ˙1=σ1,φ˙2=σ1-σ2,φ˙3=3σ1-σ3.Hence, the fixed points of (20) are given by
(23)-μ12a1+{f12ω1-σ1f14ω12}sinθ1+μ1f18ω12cosθ1-α18ω1a12a2sinθ2-α28ω1a12a3sinθ3-α48ω1a22a3sin(θ3-2θ2)-α64ω1a32a2sinθ2-3α88ω1a23sinθ2-α108ω1a1a2a3sin(θ3-θ2)=0,a1σ1+{μ128ω1-Γ12ω1}a1+{f12ω1-σ1f14ω12}cosθ1-μ1f18ω12sinθ1-3α18ω1a12a2cosθ2-α28ω1a12a3cosθ3-α34ω1a1a22-α48ω1a22a3cos(θ3-2θ2)-α54ω1a1a32-α64ω1a32a2cosθ2-3α78ω1a13-3α88ω1a23cosθ2-α108ω1a1a2a3cos(θ3-θ2)=0,-μ22a2+{f22ω2-(σ1-σ2)f24ω22}sin(θ1-θ2)+μ2f28ω22cos(θ1-θ2)+β18ω2a12a2sin2θ2-β28ω2a12a3sin(θ3-θ2)+β38ω2a22a1sinθ2-β48ω2a22a3sin(θ3-3θ2)+β54ω2a32a1sinθ2+3β78ω2a13sinθ2-β108ω2a1a2a3sin(θ3-2θ2)=0,a2(σ2-σ1)-{μ228ω2-Γ22ω2}a2+{(σ1-σ2)f24ω22-f22ω2}cos(θ1-θ2)+μ2f28ω22sin(θ1-θ2)+β18ω2a12a2cos2θ2+β14ω2a12a2+β28ω2a12a3cos(θ3-θ2)+3β38ω2a22a1cosθ2+β48ω2a22a3cos(θ3-3θ2)+β54ω2a32a1cosθ2+β64ω2a32a2+3β78ω2a13cosθ2+3β88ω2a23+β108ω2a1a2a3cos(θ3-2θ2)=0,-μ32a3-γ18ω3a12a2sin(θ2-θ3)-γ38ω3a22a1sin(2θ2-θ3)+γ78ω3a13sinθ3-γ88ω3a23sin(3θ2-θ3)=0,a3(σ3-3σ1)-{μ328ω3-Γ32ω3}a3+γ18ω3a12a2cos(θ2-θ3)+γ24ω3a12a3+γ38ω3a22a1cos(2θ2-θ3)+γ44ω3a22a3+γ78ω3a13cosθ3+γ88ω3a23cos(3θ2-θ3)+3γ98ω3a33+γ104ω3a1a2a3cosθ2=0.There are six possibilities besides the trivial solution as follows:(1)
a1≠0, a2=0, a3=0  (single mode),(2)
a2≠0, a1=0, a3=0 (single mode),(3)
a1≠0, a2≠0, a3=0 (two modes),(4)
a1≠0, a3≠0, a2=0  (two modes),(5)
a2≠0, a3≠0, a1=0   (two modes),(6)
a1≠0, a2≠0, a3≠0 (three modes).Case 1.
In this case, wherea2=0, a3=0, the frequency response equation is given by
(24)9α7264ω12a16+[R3+3α7σ14ω1]a14+[R2+σ12+μ12σ14ω1-Γ1σ1ω1]a12-μ12f1264ω14-R12=0.Case 2.
In this case, wherea1=0, a3=0, the frequency response equation is given by
(25)9β8264ω22a26+[Q3+3β8(σ2-σ1)4ω2]a24+[Q2+(σ2-σ1)2-μ22(σ2-σ1)4ω2+Γ2(σ2-σ1)ω2]a22-μ22f2264ω24-Q12=0.Case 3.
In this case, wherea3=0, the frequency response equations are given by
(26)9α7264ω12a16+[R3+3α7σ14ω1]a14+[R2+σ12+μ12σ14ω1-Γ1σ1ω1]a12-μ12f1264ω14-R12-9α1264ω12a14a22-9α1α832ω12a12a24-9α8264ω12a26+3R1α84ω1a23+3R1α14ω1a12a2=0,9β8264ω22a26+[Q3+3β8(σ2-σ1)4ω2]a24+[Q2+(σ2-σ1)2-μ22(σ2-σ1)4ω2+Γ2(σ2-σ1)ω2]a22-μ22f2264ω24-Q12+3Q1β34ω2a22a1+3Q1β74ω2a13+Q1β14ω2a12a2-9β3264ω22a24a12-[β1264ω22+9β3β732ω22]a14a22-9β7264ω22a16-3β1β732ω22a15a2-3β1β332ω22a13a23=0.Case 4.
In this case, wherea2=0, the frequency response equations are given by
(27)9α7264ω12a16+[R3+3α7σ14ω1]a14+[R2+σ12+μ12σ14ω1-Γ1σ1ω1]a12-μ12f1264ω14-R12-α2264ω12a14a32+R1α24ω1a12a3=0,9γ9264ω32a36+[K2+3γ9(σ3-3σ1)4ω3]a34+[K1+(σ3-3σ1)2-μ32(σ3-3σ1)4ω3+Γ3(σ3-3σ1)ω3]a32-γ7264ω32a16=0.Case 5.
In this case, wherea1=0, the frequency response equations are given by
(28)9β8264ω22a26+[Q3+3β8(σ2-σ1)4ω2]a24+[Q2+(σ2-σ1)2-μ22(σ2-σ1)4ω2+Γ2(σ2-σ1)ω2]a22-μ22f2264ω24-Q12+Q1β44ω2a22a3-β4264ω22a24a32=0,9γ9264ω32a36+[K2+3γ9(σ3-3σ1)4ω3]a34+[K1+(σ3-3σ1)2-μ32(σ3-3σ1)4ω3+Γ3(σ3-3σ1)ω3]a32-γ8264ω32a26=0.Case 6.
In this case, wherea1≠0, a2≠0, and a3≠0, this is the practical case, the frequency response equations are given by
(29)9α7264ω12a16+[R3+3α7σ14ω1]a14+[R2+σ12+μ12σ14ω1-Γ1σ1ω1]a12-μ12f1264ω14-R12+R1α24ω1a12a3+R1α62ω1a2a32+R1α44ω1a22a3+3R1α84ω1a23+3R1α14ω1a12a2+R1α104ω1a1a2a3-α6216ω12a22a34-α4α616ω12a23a33-α2264ω12a14a32-9α1264ω12a14a22-9α1α832ω12a12a24-9α8264ω12a26-3α4α832ω12a25a3-3α8α1032ω12a1a24a3-α2α1032ω12a13a2a32-α4α1032ω12a1a23a32-α6α1016ω12a1a22a33-α2α616ω12a12a2a33-3α1α1032ω12a13a22a3-3α1α232ω12a14a2a3-R4a12a22a32-R5a24a32-R6a23a12a3=0,9β8264ω22a26+[Q3+3β8(σ2-σ1)4ω2]a24+[Q2+(σ2-σ1)2-μ22(σ2-σ1)4ω2+Γ2(σ2-σ1)ω2]a22-μ22f2264ω24-Q12+3Q1β34ω2a22a1+3Q1β74ω2a13+Q1β14ω2a12a2+Q1β44ω2a22a3+Q1β52ω2a1a32+Q1β24ω2a12a3+Q1β104ω2a1a2a3-3β1β732ω22a15a2-9β3264ω22a24a12-9β7264ω22a16-β4264ω22a24a32-3β1β332ω22a13a23-β5216ω22a12a34-β2β516ω22a13a33-3β3β432ω22a1a24a3-3β2β732ω22a15a3-β4β516ω22a1a22a33-β5β1016ω22a12a2a23-β4β1032ω22a1a23a32-Q4a12a22a32-Q5a14a22-Q6a14a32-Q7a13a22a3-Q8a12a23a3-Q9a13a2a32-Q10a14a2a3=0,9γ9264ω32a36+[K2+3γ9(σ3-3σ1)4ω3]a34+[K1+(σ3-3σ1)2-μ32(σ3-3σ1)4ω3+Γ3(σ3-3σ1)ω3]a32-γ8264ω32a26-γ7264ω32a16-γ1γ732ω32a15a2-γ3γ832ω32a1a25-K3a14a22-K4a12a24-K5a13a22=0,
where
(30)R1=[f12ω1-σ1f14ω12],R2=[μ124+μ1464ω12+Γ124ω12-Γ1μ128ω12],R3=[3α7μ1232ω12-3α7Γ18ω12],R4=[α10264ω12+3α1α616ω12+α2α432ω12],R5=[3α6α816ω12+α4264ω12],R6=[3α1α432ω12+3α2α832ω12],Q1=[f22ω2-(σ1-σ2)f24ω22],Q2=[μ224+μ2464ω22+Γ224ω22-Γ2μ228ω22],Q3=[3β8Γ28ω22-3β8μ2232ω22],Q4=[β10264ω22+3β3β516ω22+β2β432ω22],Q5=[β1264ω22+9β3β732ω22],Q6=[β2264ω22+3β5β716ω22],Q7=[3β2β332ω22+3β4β732ω22+β1β1032ω22],Q8=[3β3β1032ω22+β1β432ω22],Q9=[β2β1032ω22+β1β516ω22],Q10=[3β7β1032ω22+β1β232ω22],K1=[μ324+μ3464ω32+Γ324ω32-Γ3μ328ω32],K2=[3γ9Γ38ω32-3γ9μ3232ω32],K3=[γ1264ω32+γ3γ732ω32],K4=[γ3264ω32+γ1γ832ω32],K5=[γ1γ332ω32+γ7γ832ω32].
In the frequency response curves, the stable (unstable) steady-state solutions have been represented by solid (dotted) lines.


## 4. Numerical Results and Discussion

To study the behavior of the system, the Runge-Kutta fourth-order method (RKM) was applied to (4a), (4b), and (4c) governing the oscillating system. A good criterion of both stability and dynamic chaos is the phase plane trajectories. Figure 2 illustrates the response and the phase-plane for the nonresonant system at some practical values of the equation parameters: μ1=0.05, μ2=0.05, μ3=0.05, α1=β1=γ1=1.5, α2=β2=γ2=1.9, α3=β3=γ3=0.4, α4=β4=γ4=0.6, α5=β5=γ5=0.6, α6=β6=γ6=0.2, α7=β7=γ7=0.01, α8=β8=γ8=0.01, α9=β9=γ9=0.4, α10=β10=γ10=1.8, Ω1=6.4, Ω2=6.1, Ω3=5.7, Ω4=2, ω1=4.7, ω2=2.7, ω3=13.1, f1=f2=4, f3=20.5, f11=f21=f31=0.1, f12=f22=f32=0.2, and f14=f24=f34=0.5. Figure 3 shows the steady-state amplitudes and phase plane of the system at simultaneous resonance case Ω3≅ω1, ω2≅ω1, and ω3≅3ω1. It is clear from Figure 3 that the steady-state amplitude of the first, second, and third modes is increased to about 370%, 890%, and 260%, respectively, of its value shown in Figure 2. Also, it can be seen that the time response of the system is tuned with multilimit cycle.Figure 2
Time response and phase-plane diagrams of the system for nonresonance case.Figure 3
Time response and phase-plane diagrams of the system at  simultaneous resonance case,Ω3≅ω1, ω2≅ω1, and ω3≅3ω1.It is quite clear that such case is undesirable in the design of such system because it represents one of the worst behaviors of the system. Such case should be avoided as working condition for the system. It is advised for such system not to haveω2=ω1 or ω3=3ω1. Figure 4 shows the comparison between numerical integration for the system equation (4a), (4b), and (4c) solid lines and the amplitude-phase modulating equation (20) dashed lines. We found that all predictions from analytical solutions dashed lines are in good agreement with the numerical simulation solid lines.Figure 4
Comparison between numerical solution (using RKM) and analytical solution (using perturbation method) of the system at resonance case,Ω3≅ω1, ω2≅ω1, and ω3≅3ω1.


### 4.1. FDM with Approximation O(c2)

The infinite equations (2a)–(2e) will be solved via a finite difference method (FDM). We briefly describe the procedure here. More details are available in [23, 26]. The infinite dimensional equations (2a)–(2e) can be reduced to the finite dimensional one via the finite difference method with second-order approximation O(c2). Namely, at each mesh node the following system of ordinary differential equations is obtained:
(31)L1,c(Nxx,Nxy)=I0(u¨0)i,j+J1(ϕ¨x)i,j-23h2cI3((w¨0)i+1,j-(w¨0)i-1,j),L2,c(Nxy,Nyy)=I0(v¨0)i,j+J1(ϕ¨y)i,j-23h2cI3((w¨0)i,j+1-(w¨0)i,j-1),L3,c(Q-x,Q-y,Nxx,w0,Nxy)=I0(w¨0)i,j-169h4c2I6((w¨0)i+1,j-2(w¨0)i,j+(w¨0)i+1,j+(w¨0)i,j+1-2(w¨0)i,j+(w¨0)i,j+1)+43h2c2I3((u¨0)i+1,j-(u¨0)i-1,j+(v¨0)i,j+1-(v¨0)i,j-1)+J42c((ϕ¨x)i+1,j-(ϕ¨x)i-1,j+(ϕ¨y)i,j+1-(ϕ¨y)i,j-1),L4,c(M-xx,M-xy)=J1(u¨0)i,j+k2(ϕ¨x)i,j-2J43h2c((w¨0)i+1,j-(w¨0)i-1,j),L5,c(M-xy,M-yy)=J1(v¨0)i,j+k2(ϕ¨y)i,j-2J43h2c((w¨0)i,j+1-(w¨0)i,j-1),(i=1,2,…,n),(j=1,2,…,n),
where n denotes the partition number regarding a spatial coordinate, c is the computational step regarding spatial coordinate, and L1,c(·), L2,c(·), L3,c(·), L4,c(·), and L5,c(·) are the difference operators as follows:(32)L1,c(Nxx,Nxy)=12c((Nxx)i+1,j-(Nxx)i-1,j+(Nxy)i,j+1-(Nxy)i,j-1),L2,c(Nxy,Nyy)=12c((Nxy)i+1,j-(Nxy)i-1,j+(Nyy)i,j+1-(Nyy)i,j-1),L3,c(Q-x,Q-y,Nxx,w0,Nxy)=12c((Q-x)i+1,j-(Q-x)i-1,j+(Q-y)i,j+1-(Q-y)i,j-1)+Nxx((w0)i-1,j-2(w0)i,j+(w0)i+1,jc2)+((w0)i+1,j-(w0)i-1,j2c)((Nxx)i+1,j-(Nxx)i-1,j2c)+Nxy(((w0)j+1,i+1-2(w0)j+1,i-1)((w0)j-1,i+1-2(w0)j-1,i-1)4c2)+((w0)i,j+1-(w0)i,j-12c)((Nxy)i+1,j-(Nxy)i-1,j2c)+Nxy(((w0)i+1,j+1-(w0)i+1,j-1)-((w0)i-1,j+1-(w0)i-1,j-1)8c3)+((w0)i+1,j-(w0)i-1,j2c)((Nyy)i,j+1-(Nyy)i,j-12c)+Nyy((w0)i,j-1-2(w0)i,j+(w0)i,j+1c2)+((w0)i,j+1-(w0)i,j-12c)((Nyy)i,j+1-(Nyy)i,j-12c)+43h2((Pxx)i-1,j-2(Pxx)i,j+(Pxx)i+1,jc2)+43h2(2(((Pxy)i+1,j+1-(Pxy)i-1,j+1)-((Pxy)i+1,j-1-(Pxy)i-1,j-1)8c3)+(Pyy)i,j-1-2(Pyy)i,j+(Pyy)i,j+1c2)+q,L4,c(M-xx,M-xy)=12c((M-xx)i+1,j-(M-xx)i-1,j+(M-xy)i,j+1-(M-xy)i,j-1)-Q-x,L5,c(M-xy,M-yy)=12c((M-xy)i+1,j-(M-xy)i-1,j+(M-yy)i,j+1-(M-yy)i,j-1)-Q-y.The obtained system of (31) with the supplemented boundary conditions equation and the initial conditions equation is solved by the fourth-order Runge-Kutta method. Figure 5 shows a comparison between the time responses of the system equations (4a), (4b), and (4c) using the Runge-Kutta of fourth-order method and the time response of the problems (31), using the finite difference method at the same values of the parameters shown in Figure 2.Figure 5
Time response of composite laminated rectangular plate using RKM and FDM.


### 4.2. Frequency Response Curves

When the amplitude achieves a constant nontrivial value, a steady-state vibration exists. Using the frequency response equations we can assess the influence of the damping coefficients, the nonlinear parameters, and the excitation amplitude on the steady-state amplitudes. The frequency response equations (24)–(29) are nonlinear equations in a1, a2, and a3 which are solved numerically. The numerical results are shown in Figures 6–25, and in all figures the region of stability of the nonlinear solutions is determined by applying the Routh-Hurwitz criterion. The solid lines stand for the stable solution, and the dotted lines stand for the unstable solution. From the geometry of the figures, we observe that each curve is continuous and has stable and unstable solutions.Figure 6
Comparison between analytical prediction using multiple time scale and numerical integration of the first mode.Figure 7
Effects of the linear dampingμ1.Figure 8
Effects of the natural frequencyω1.Figure 9
Effects of the external excitationf1.Figure 10
Effects of the nonlinear parameterα7.Figure 11
Effects of the parametric excitationf14.Figure 12
Comparison between analytical prediction using multiple time scale and numerical integration of the second mode.Figure 13
Effects of the nonlinear parameterβ8.Figure 14
Effects of the linear dampingμ2.Figure 15
Effects of the natural frequencyω2.Figure 16
Effects of the external excitationf2.Figure 17
Effects of the detuning parameterσ1.Figure 18
Effects of the external excitation forcef1.Figure 19
Effects of the nonlinear parameters (α1, β1, and γ1 coupling terms).Figure 20
Effects of the nonlinear parameters (α2, β2, and γ2 coupling terms).Figure 21
Effects of the nonlinear parameters (α8, β8, and γ8 coupling terms).Figure 22
Effects of the nonlinear parameters (α7, β7, and γ7 coupling terms).Figure 23
Effects of the nonlinear parameters (α3, β3, and γ3 coupling terms).Figure 24
Effects of the nonlinear parameters (α5, β5, and γ5 coupling terms).Figure 25
Effects of the linear damping coefficients.


#### 4.2.1. Response Curve of Case1

To check the accuracy of the analytical solution derived by the multiple time scale in predicting the amplitude of the first mode, we compare the amplitude of the first mode obtained from frequency response equation of Case1 with values obtained from numerical integration of (4a). Figure 6 shows a comparison of these outputs for the first mode. The effects of the detuning parameter σ1 on the steady-state amplitude of the first mode for the stability first case, where a1≠0, a2=0, and a3=0, for the parameters μ1=0.2, α7=0.01, ω1=2.3, f1=4,  f11=0.1, f12=0.2, f14=0.5, Ω1=1, Ω2=1.2, and Ω4=1.4, as shown in Figure 6.Figures7–11, show the effects of the damping coefficient μ1, the first mode natural frequency ω1, the external excitation amplitude f1, the nonlinear spring stiffness α7, and the parametric excitation f14. Figures 7 and 8 show that the steady-state amplitude a1 is inversely proportional to μ1 and ω1, and also for decreasing μ1 or ω1 the curve is bending to the left.  It is clear from Figure 9 that the steady-state amplitude a1 is increasing for increasing value of external excitation force f1, and the zone of instability is increased. Figure 10 shows that as the nonlinear spring stiffness α7 is increased; the continuous curve is moved downwards. Also, the negative and positive values of α7 produce either hard or soft spring, respectively, as the curve is either bent to the right or to the left, leading to the appearance of the jump phenomenon. The region of stability is increased for increasing value of α7. From Figure 11, we observe that for increasing value of parametric excitation amplitude f14, the steady-state amplitude of the first mode is increased, and the curve is shifted to the left.


#### 4.2.2. Response Curve of Case2

Figures12–16, show the frequency response curves for the stability of the second case, where a2≠0, a1=0, and a3=0. Figure 12 shows a comparison of these outputs for the second mode. The effects of the detuning parameter σ2 on the steady-state amplitude of the second mode for the stability second case, where a2≠0, a1=0, a3=0, for the parameters: μ2=0.2, β8=0.01, ω2=4, f2=4, f21=0.1, f22=0.2, f24=0.5, Ω1=1, Ω2=1.2, and Ω4=1.4, as shown in Figure 12. It can be seen from the figure that maximum steady-state amplitude occurs when ω2≅ω1. Figure 13 shows that as the nonlinear spring stiffness β8 is increased, the continuous curve is moved downwards. Also, the positive and negative values of β8 produce either soft or hard spring, respectively, as the curve is either bent to the left or to the right, leading to the appearance of the jump phenomenon. Figures 14, 15, and 16 show that the steady-state amplitude a2 is inversely proportional to μ2, ω2 and directly proportional to the external excitation f2. Also, for decreasing μ2, ω2 the curve is bending to the left.


#### 4.2.3. Response Curve of Case6

Figures17, 18, 19, 20, 21, 22, 23, 24, and 25 show that the frequency response curves for practical case stability, where a1≠0, a2≠0, and a3≠0. Figure 17 shows that the effects of the detuning parameter σ1 on the amplitudes of the three modes. From this figure, we observe that these modes intersect, and for positive value of the detuning parameter σ1 the amplitudes are stable. For negative value of the detuning parameter σ1 down to and including −0.5 the system becomes unstable.Figure18 shows that the steady-state amplitudes of the three modes a1, a2, and a3 are directly proportional to the external excitation force  f1. Also, form this figure we show that the stability region is decreased for increasing  f1. Figures 19–21 show that the steady-state amplitudes and stability of the three modes a1, a2, and a3 are inversely proportional to the nonlinear parameters (α1,α2,α8), (β1,β2,β8), and (γ1,γ2,γ8), respectively.Figures22–24, show the effects of the nonlinear parameters (α7,β7,γ7), (α3,β3,γ3), and (α5,β5,  γ5) on the steady-state amplitudes of the three modes. It is clear that the stability regions are increased for increasing these nonlinear parameters. For increasing value of linear viscous damping coefficients, we note that the steady-state amplitudes are increasing or decreasing for the first, second, and third modes, respectively, as shown in Figure 25. The region of stability system is increased for decreasing value of damping coefficients.


### 4.3. Comparison Study

In the previous work [11], the chaotic dynamics of a six-dimensional nonlinear system which represents the averaged equation of a composite laminated piezoelectric rectangular plate subjected to the transverse, in-plane excitations and the excitation loaded by piezoelectric layers are analyzed. The case of 1 : 2 : 4 internal resonances is considered.In our study, the nonlinear analysis and stability of a composite laminated piezoelectric rectangular thin plate under simultaneous external and parametric excitation forces are investigated. The second-order approximation is obtained to consider the influence of the cubic terms on nonlinear dynamic characteristics of the composite laminated piezoelectric rectangular plate using the multiple scale method. All possible resonance cases are extracted at this approximation order. The case of 1 : 1 : 3 internal resonance and primary resonance is considered. The stability of the system and the effects of different parameters on system behavior have been studied using phase plane and frequency response curves. The analytical results given by the method of multiple time scale are verified by comparison with results from numerical integration of the modal equations. Reliability of the obtained results is verified by comparison between the finite difference method (FDM) and Runge-Kutta method (RKM). Variation of the some parameters leads to multivalued amplitudes and hence to jump phenomena.  It is quite clear that some of the simultaneous resonance cases are undesirable in the design of such system. Such cases should be avoided as working conditions for the system.


## 4.1. FDM with Approximation O(c2)

The infinite equations (2a)–(2e) will be solved via a finite difference method (FDM). We briefly describe the procedure here. More details are available in [23, 26]. The infinite dimensional equations (2a)–(2e) can be reduced to the finite dimensional one via the finite difference method with second-order approximation O(c2). Namely, at each mesh node the following system of ordinary differential equations is obtained:
(31)L1,c(Nxx,Nxy)=I0(u¨0)i,j+J1(ϕ¨x)i,j-23h2cI3((w¨0)i+1,j-(w¨0)i-1,j),L2,c(Nxy,Nyy)=I0(v¨0)i,j+J1(ϕ¨y)i,j-23h2cI3((w¨0)i,j+1-(w¨0)i,j-1),L3,c(Q-x,Q-y,Nxx,w0,Nxy)=I0(w¨0)i,j-169h4c2I6((w¨0)i+1,j-2(w¨0)i,j+(w¨0)i+1,j+(w¨0)i,j+1-2(w¨0)i,j+(w¨0)i,j+1)+43h2c2I3((u¨0)i+1,j-(u¨0)i-1,j+(v¨0)i,j+1-(v¨0)i,j-1)+J42c((ϕ¨x)i+1,j-(ϕ¨x)i-1,j+(ϕ¨y)i,j+1-(ϕ¨y)i,j-1),L4,c(M-xx,M-xy)=J1(u¨0)i,j+k2(ϕ¨x)i,j-2J43h2c((w¨0)i+1,j-(w¨0)i-1,j),L5,c(M-xy,M-yy)=J1(v¨0)i,j+k2(ϕ¨y)i,j-2J43h2c((w¨0)i,j+1-(w¨0)i,j-1),(i=1,2,…,n),(j=1,2,…,n),
where n denotes the partition number regarding a spatial coordinate, c is the computational step regarding spatial coordinate, and L1,c(·), L2,c(·), L3,c(·), L4,c(·), and L5,c(·) are the difference operators as follows:(32)L1,c(Nxx,Nxy)=12c((Nxx)i+1,j-(Nxx)i-1,j+(Nxy)i,j+1-(Nxy)i,j-1),L2,c(Nxy,Nyy)=12c((Nxy)i+1,j-(Nxy)i-1,j+(Nyy)i,j+1-(Nyy)i,j-1),L3,c(Q-x,Q-y,Nxx,w0,Nxy)=12c((Q-x)i+1,j-(Q-x)i-1,j+(Q-y)i,j+1-(Q-y)i,j-1)+Nxx((w0)i-1,j-2(w0)i,j+(w0)i+1,jc2)+((w0)i+1,j-(w0)i-1,j2c)((Nxx)i+1,j-(Nxx)i-1,j2c)+Nxy(((w0)j+1,i+1-2(w0)j+1,i-1)((w0)j-1,i+1-2(w0)j-1,i-1)4c2)+((w0)i,j+1-(w0)i,j-12c)((Nxy)i+1,j-(Nxy)i-1,j2c)+Nxy(((w0)i+1,j+1-(w0)i+1,j-1)-((w0)i-1,j+1-(w0)i-1,j-1)8c3)+((w0)i+1,j-(w0)i-1,j2c)((Nyy)i,j+1-(Nyy)i,j-12c)+Nyy((w0)i,j-1-2(w0)i,j+(w0)i,j+1c2)+((w0)i,j+1-(w0)i,j-12c)((Nyy)i,j+1-(Nyy)i,j-12c)+43h2((Pxx)i-1,j-2(Pxx)i,j+(Pxx)i+1,jc2)+43h2(2(((Pxy)i+1,j+1-(Pxy)i-1,j+1)-((Pxy)i+1,j-1-(Pxy)i-1,j-1)8c3)+(Pyy)i,j-1-2(Pyy)i,j+(Pyy)i,j+1c2)+q,L4,c(M-xx,M-xy)=12c((M-xx)i+1,j-(M-xx)i-1,j+(M-xy)i,j+1-(M-xy)i,j-1)-Q-x,L5,c(M-xy,M-yy)=12c((M-xy)i+1,j-(M-xy)i-1,j+(M-yy)i,j+1-(M-yy)i,j-1)-Q-y.The obtained system of (31) with the supplemented boundary conditions equation and the initial conditions equation is solved by the fourth-order Runge-Kutta method. Figure 5 shows a comparison between the time responses of the system equations (4a), (4b), and (4c) using the Runge-Kutta of fourth-order method and the time response of the problems (31), using the finite difference method at the same values of the parameters shown in Figure 2.Figure 5
Time response of composite laminated rectangular plate using RKM and FDM.


## 4.2. Frequency Response Curves

When the amplitude achieves a constant nontrivial value, a steady-state vibration exists. Using the frequency response equations we can assess the influence of the damping coefficients, the nonlinear parameters, and the excitation amplitude on the steady-state amplitudes. The frequency response equations (24)–(29) are nonlinear equations in a1, a2, and a3 which are solved numerically. The numerical results are shown in Figures 6–25, and in all figures the region of stability of the nonlinear solutions is determined by applying the Routh-Hurwitz criterion. The solid lines stand for the stable solution, and the dotted lines stand for the unstable solution. From the geometry of the figures, we observe that each curve is continuous and has stable and unstable solutions.Figure 6
Comparison between analytical prediction using multiple time scale and numerical integration of the first mode.Figure 7
Effects of the linear dampingμ1.Figure 8
Effects of the natural frequencyω1.Figure 9
Effects of the external excitationf1.Figure 10
Effects of the nonlinear parameterα7.Figure 11
Effects of the parametric excitationf14.Figure 12
Comparison between analytical prediction using multiple time scale and numerical integration of the second mode.Figure 13
Effects of the nonlinear parameterβ8.Figure 14
Effects of the linear dampingμ2.Figure 15
Effects of the natural frequencyω2.Figure 16
Effects of the external excitationf2.Figure 17
Effects of the detuning parameterσ1.Figure 18
Effects of the external excitation forcef1.Figure 19
Effects of the nonlinear parameters (α1, β1, and γ1 coupling terms).Figure 20
Effects of the nonlinear parameters (α2, β2, and γ2 coupling terms).Figure 21
Effects of the nonlinear parameters (α8, β8, and γ8 coupling terms).Figure 22
Effects of the nonlinear parameters (α7, β7, and γ7 coupling terms).Figure 23
Effects of the nonlinear parameters (α3, β3, and γ3 coupling terms).Figure 24
Effects of the nonlinear parameters (α5, β5, and γ5 coupling terms).Figure 25
Effects of the linear damping coefficients.


### 4.2.1. Response Curve of Case1

To check the accuracy of the analytical solution derived by the multiple time scale in predicting the amplitude of the first mode, we compare the amplitude of the first mode obtained from frequency response equation of Case1 with values obtained from numerical integration of (4a). Figure 6 shows a comparison of these outputs for the first mode. The effects of the detuning parameter σ1 on the steady-state amplitude of the first mode for the stability first case, where a1≠0, a2=0, and a3=0, for the parameters μ1=0.2, α7=0.01, ω1=2.3, f1=4,  f11=0.1, f12=0.2, f14=0.5, Ω1=1, Ω2=1.2, and Ω4=1.4, as shown in Figure 6.Figures7–11, show the effects of the damping coefficient μ1, the first mode natural frequency ω1, the external excitation amplitude f1, the nonlinear spring stiffness α7, and the parametric excitation f14. Figures 7 and 8 show that the steady-state amplitude a1 is inversely proportional to μ1 and ω1, and also for decreasing μ1 or ω1 the curve is bending to the left.  It is clear from Figure 9 that the steady-state amplitude a1 is increasing for increasing value of external excitation force f1, and the zone of instability is increased. Figure 10 shows that as the nonlinear spring stiffness α7 is increased; the continuous curve is moved downwards. Also, the negative and positive values of α7 produce either hard or soft spring, respectively, as the curve is either bent to the right or to the left, leading to the appearance of the jump phenomenon. The region of stability is increased for increasing value of α7. From Figure 11, we observe that for increasing value of parametric excitation amplitude f14, the steady-state amplitude of the first mode is increased, and the curve is shifted to the left.


### 4.2.2. Response Curve of Case2

Figures12–16, show the frequency response curves for the stability of the second case, where a2≠0, a1=0, and a3=0. Figure 12 shows a comparison of these outputs for the second mode. The effects of the detuning parameter σ2 on the steady-state amplitude of the second mode for the stability second case, where a2≠0, a1=0, a3=0, for the parameters: μ2=0.2, β8=0.01, ω2=4, f2=4, f21=0.1, f22=0.2, f24=0.5, Ω1=1, Ω2=1.2, and Ω4=1.4, as shown in Figure 12. It can be seen from the figure that maximum steady-state amplitude occurs when ω2≅ω1. Figure 13 shows that as the nonlinear spring stiffness β8 is increased, the continuous curve is moved downwards. Also, the positive and negative values of β8 produce either soft or hard spring, respectively, as the curve is either bent to the left or to the right, leading to the appearance of the jump phenomenon. Figures 14, 15, and 16 show that the steady-state amplitude a2 is inversely proportional to μ2, ω2 and directly proportional to the external excitation f2. Also, for decreasing μ2, ω2 the curve is bending to the left.


### 4.2.3. Response Curve of Case6

Figures17, 18, 19, 20, 21, 22, 23, 24, and 25 show that the frequency response curves for practical case stability, where a1≠0, a2≠0, and a3≠0. Figure 17 shows that the effects of the detuning parameter σ1 on the amplitudes of the three modes. From this figure, we observe that these modes intersect, and for positive value of the detuning parameter σ1 the amplitudes are stable. For negative value of the detuning parameter σ1 down to and including −0.5 the system becomes unstable.Figure18 shows that the steady-state amplitudes of the three modes a1, a2, and a3 are directly proportional to the external excitation force  f1. Also, form this figure we show that the stability region is decreased for increasing  f1. Figures 19–21 show that the steady-state amplitudes and stability of the three modes a1, a2, and a3 are inversely proportional to the nonlinear parameters (α1,α2,α8), (β1,β2,β8), and (γ1,γ2,γ8), respectively.Figures22–24, show the effects of the nonlinear parameters (α7,β7,γ7), (α3,β3,γ3), and (α5,β5,  γ5) on the steady-state amplitudes of the three modes. It is clear that the stability regions are increased for increasing these nonlinear parameters. For increasing value of linear viscous damping coefficients, we note that the steady-state amplitudes are increasing or decreasing for the first, second, and third modes, respectively, as shown in Figure 25. The region of stability system is increased for decreasing value of damping coefficients.


## 4.2.1. Response Curve of Case1

To check the accuracy of the analytical solution derived by the multiple time scale in predicting the amplitude of the first mode, we compare the amplitude of the first mode obtained from frequency response equation of Case1 with values obtained from numerical integration of (4a). Figure 6 shows a comparison of these outputs for the first mode. The effects of the detuning parameter σ1 on the steady-state amplitude of the first mode for the stability first case, where a1≠0, a2=0, and a3=0, for the parameters μ1=0.2, α7=0.01, ω1=2.3, f1=4,  f11=0.1, f12=0.2, f14=0.5, Ω1=1, Ω2=1.2, and Ω4=1.4, as shown in Figure 6.Figures7–11, show the effects of the damping coefficient μ1, the first mode natural frequency ω1, the external excitation amplitude f1, the nonlinear spring stiffness α7, and the parametric excitation f14. Figures 7 and 8 show that the steady-state amplitude a1 is inversely proportional to μ1 and ω1, and also for decreasing μ1 or ω1 the curve is bending to the left.  It is clear from Figure 9 that the steady-state amplitude a1 is increasing for increasing value of external excitation force f1, and the zone of instability is increased. Figure 10 shows that as the nonlinear spring stiffness α7 is increased; the continuous curve is moved downwards. Also, the negative and positive values of α7 produce either hard or soft spring, respectively, as the curve is either bent to the right or to the left, leading to the appearance of the jump phenomenon. The region of stability is increased for increasing value of α7. From Figure 11, we observe that for increasing value of parametric excitation amplitude f14, the steady-state amplitude of the first mode is increased, and the curve is shifted to the left.


## 4.2.2. Response Curve of Case2

Figures12–16, show the frequency response curves for the stability of the second case, where a2≠0, a1=0, and a3=0. Figure 12 shows a comparison of these outputs for the second mode. The effects of the detuning parameter σ2 on the steady-state amplitude of the second mode for the stability second case, where a2≠0, a1=0, a3=0, for the parameters: μ2=0.2, β8=0.01, ω2=4, f2=4, f21=0.1, f22=0.2, f24=0.5, Ω1=1, Ω2=1.2, and Ω4=1.4, as shown in Figure 12. It can be seen from the figure that maximum steady-state amplitude occurs when ω2≅ω1. Figure 13 shows that as the nonlinear spring stiffness β8 is increased, the continuous curve is moved downwards. Also, the positive and negative values of β8 produce either soft or hard spring, respectively, as the curve is either bent to the left or to the right, leading to the appearance of the jump phenomenon. Figures 14, 15, and 16 show that the steady-state amplitude a2 is inversely proportional to μ2, ω2 and directly proportional to the external excitation f2. Also, for decreasing μ2, ω2 the curve is bending to the left.


## 4.2.3. Response Curve of Case6

Figures17, 18, 19, 20, 21, 22, 23, 24, and 25 show that the frequency response curves for practical case stability, where a1≠0, a2≠0, and a3≠0. Figure 17 shows that the effects of the detuning parameter σ1 on the amplitudes of the three modes. From this figure, we observe that these modes intersect, and for positive value of the detuning parameter σ1 the amplitudes are stable. For negative value of the detuning parameter σ1 down to and including −0.5 the system becomes unstable.Figure18 shows that the steady-state amplitudes of the three modes a1, a2, and a3 are directly proportional to the external excitation force  f1. Also, form this figure we show that the stability region is decreased for increasing  f1. Figures 19–21 show that the steady-state amplitudes and stability of the three modes a1, a2, and a3 are inversely proportional to the nonlinear parameters (α1,α2,α8), (β1,β2,β8), and (γ1,γ2,γ8), respectively.Figures22–24, show the effects of the nonlinear parameters (α7,β7,γ7), (α3,β3,γ3), and (α5,β5,  γ5) on the steady-state amplitudes of the three modes. It is clear that the stability regions are increased for increasing these nonlinear parameters. For increasing value of linear viscous damping coefficients, we note that the steady-state amplitudes are increasing or decreasing for the first, second, and third modes, respectively, as shown in Figure 25. The region of stability system is increased for decreasing value of damping coefficients.


## 4.3. Comparison Study

In the previous work [11], the chaotic dynamics of a six-dimensional nonlinear system which represents the averaged equation of a composite laminated piezoelectric rectangular plate subjected to the transverse, in-plane excitations and the excitation loaded by piezoelectric layers are analyzed. The case of 1 : 2 : 4 internal resonances is considered.In our study, the nonlinear analysis and stability of a composite laminated piezoelectric rectangular thin plate under simultaneous external and parametric excitation forces are investigated. The second-order approximation is obtained to consider the influence of the cubic terms on nonlinear dynamic characteristics of the composite laminated piezoelectric rectangular plate using the multiple scale method. All possible resonance cases are extracted at this approximation order. The case of 1 : 1 : 3 internal resonance and primary resonance is considered. The stability of the system and the effects of different parameters on system behavior have been studied using phase plane and frequency response curves. The analytical results given by the method of multiple time scale are verified by comparison with results from numerical integration of the modal equations. Reliability of the obtained results is verified by comparison between the finite difference method (FDM) and Runge-Kutta method (RKM). Variation of the some parameters leads to multivalued amplitudes and hence to jump phenomena.  It is quite clear that some of the simultaneous resonance cases are undesirable in the design of such system. Such cases should be avoided as working conditions for the system.


## 5. Conclusions

Multiple time scale perturbation method is applied to determine second-order approximate solutions for rectangular symmetric cross-ply laminated composite thin plate subjected to external and parametric excitations. Second-order approximate solutions are obtained to study the influence of the cubic terms on nonlinear dynamic characteristics of the composite laminated piezoelectric rectangular plate. All possible resonance cases are extracted at this approximation order. The study is focused on the case of 1 : 1 : 3 primary resonance and internal resonance, whereΩ3≅ω1, ω2≅ω1, and ω3≅3ω1. The analytical results given by the method of multiple time scale are verified by comparison with results from numerical integration of the modal equations. Reliability of the obtained results is verified by comparison between the finite difference method (FDM) and Runge-Kutta method (RKM). The stability of a composite laminated piezoelectric rectangular thin plate is investigated. The phase-plane method and frequency response curves are applied to study the stability of the system. From the previous study the following may be concluded.(1)
The simultaneous resonance caseΩ3≅ω1,ω2≅ω1, and ω3≅3ω1 is one of the worst case, and they should be avoided in design of such system. Of course, the excitation frequency Ω3 is out of control. But this case should be avoided through having ω2≠ω1 or ω3≠3ω1.(2)
A comparison between the solutions obtained numerically with that prediction from the multiple scales shows an excellent agreement.(3)
Reliability of the obtained results is verified by comparison between the finite difference method (FDM) and Runge-Kutta method (RKM).(4)
Variation of the parametersμ1, μ2, α7, β8, ω1, ω2, f1, and f2 leads to multivalued amplitudes and hence to jump phenomena.(5)
For the first and second modes, the steady-state amplitudesa1 and a2 are directly proportional to the excitation amplitude f1 and f2 and inversely proportional to the linear damping μ1 and μ2, respectively.(6)
The negative and positive values of nonlinear stiffnessα7, β8 produce either hard or soft spring, respectively, as the curve is either bent to the right or to the left.(7)
The region of instability increase, which is undesirable, for increasing excitation amplitudesf1, f2 and for negative values of nonlinear stiffness α7, β8.(8)
The multivalued solutions are disappeared for increasing linear damping coefficientsμ1, μ2.(9)
The region of stability increases, which is desirable for decreasing excitation amplitudef1, nonlinear parameters (α1,β1,γ1), (α2,β2,γ2), and for increasing nonlinear stiffness (α7,β7,γ7).(10)
The steady-state amplitudes of the three modesa1, a2, and a3 are directly proportional to the excitation amplitude f1 and inversely proportional to the nonlinear parameters (α1,α2,α8), (β1,β2,β8), and (γ1,γ2,γ8), respectively.For further work, we intend to extend our work to explore the modeling formulation by investigating the role static loading (in addition to the dynamic component).



---

*Source: 418374-2013-06-10.xml*