# Multipulse Heteroclinic Orbits and Chaotic Dynamics of the Laminated Composite Piezoelectric Rectangular Plate

Authors: Minghui Yao; Wei Zhang
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

This paper investigates the multipulse global bifurcations and chaotic dynamics for the nonlinear oscillations of the laminated composite piezoelectric rectangular plate by using an energy phase method in the resonant case. Using the von Karman type equations, Reddy’s third-order shear deformation plate theory, and Hamilton’s principle, the equations of motion are derived for the laminated composite piezoelectric rectangular plate with combined parametric excitations and transverse excitation. Applying the method of multiple scales and Galerkin’s approach to the partial differential governing equation, the four-dimensional averaged equation is obtained for the case of 1 : 2 internal resonance and primary parametric resonance. The energy phase method is used for the first time to investigate the Shilnikov type multipulse heteroclinic bifurcations and chaotic dynamics of the laminated composite piezoelectric rectangular plate. The paper demonstrates how to employ the energy phase method to analyze the Shilnikov type multipulse heteroclinic bifurcations and chaotic dynamics of high-dimensional nonlinear systems in engineering applications. Numerical simulations show that for the nonlinear oscillations of the laminated composite piezoelectric rectangular plate, the Shilnikov type multipulse chaotic motions can occur. Overall, both theoretical and numerical studies suggest that chaos for the Smale horseshoe sense in motion exists.


---

## Body

## 1. Introduction

A piezoelectric material subjected to the mechanical force produces an electrical charge which is called the direct piezoelectric effect. Conversely, the material under an electrical field can generate mechanical stress which is called converse piezoelectric effect. In the development of intelligent structures, piezoelectric materials are widely designed as sensors and actuators for the active control of structures. Piezoelectric materials are usually clung to structural laminates used as devices to control deformation, shape, and vibration. There has been the model proposed for analysis of laminated composite plates containing active and passive piezoelectric layers. In aerospace applications, these smart structures are generally lightweight and have relatively large structural flexibility, and their flexibility can induce large deformation. Due to small material damping or the lack of environmental damping in space, fast motion or high-speed operation of laminated composite piezoelectric plates often generates the large amplitude vibration or geometrical nonlinearity. The complex nonlinear motions are strongly incited under the certain ranges of exciting frequency. Therefore, it is important to investigate large deformation and geometrically nonlinear effects of laminated composite piezoelectric plates in order to accurately design and effectively control structures.Most of the studies on piezoelectric structures are based on linear piezoelectricity and linear elasticity theories, while research on the nonlinear dynamics of the piezoelectric structures is less. Lee et al. [1] incorporated the piezoelectric effect into the classical laminate plate theory and derived the distributed sensors and actuators which are capable of sensing and controlling the modal vibration of the cantilever plate. Tzou and Bao [2] formulated a new theory on thick anisotropic composite piezoelectric shell transducer laminates. This theory is applicable to a variety of composite piezoelectric structures. Reddy and Mitchell [3] developed geometrically nonlinear theories of laminated composite plates with piezoelectric laminae and used Hamilton's principle to derive equations of motion and boundary conditions of these theories. Hagood IV and McFarland [4] employed the Rayleigh-Ritz mode energy method to study the distributed piezoceramics and the traveling wave dynamics of the stator. Sadri et al. [5] presented the theoretical vibration model of the piezoelectric plate and used the Rayleigh-Ritz method to obtain the governing nonlinear equations of the piezoelectric plate. Vel and Batra [6] utilized Eshelby-Stroh formalism to analyze the quasistatic deformations of linear piezoelectric laminated plates. Yu and Hodges [7] took advantage of the variational-asymptotic method to construct a Reissner-Mindlin model for the laminated composite piezoelectric plates subjected to mechanical, thermal, and electric loads. Huang and Shen [8] exploited the higher order shear deformation plate theory and von Karman-type equation to investigate the nonlinear vibration and dynamic response of the laminated piezoelectric plates subjected to mechanical, electrical, and thermal loads. Della and Shu [9] gave a review on the various analytical models and numerical analyses for free vibration of delaminated composites including composite piezoelectric laminates under axial loads. Zhang et al. [10] made use of Reddy’s third-order shear deformation plate theory to analyze the bifurcations and chaotic dynamics of the four-edge simply supported, laminated composite piezoelectric rectangular plate in the case of 1 : 2 internal resonances. Sarangi and Ray [11] analyzed geometrically nonlinear transient vibrations of doubly curved laminated composite piezoelectric shells.The global bifurcations and chaotic dynamics of high-dimensional nonlinear systems have been at the forefront of nonlinear dynamics for the past two decades. There are two ways of solutions on the Shilnikov type chaotic dynamics of high-dimensional nonlinear systems. One is Shilnikov type single-pulse chaotic dynamics and the other is the Shilnikov type multipulse chaotic dynamics. Most researchers focused on the Shilnikov type single-pulse chaotic dynamics of high-dimensional nonlinear systems. Wiggins [12] used the Melnikov method to investigate the global bifurcations and chaotic dynamics of three types high-dimensional perturbed Hamiltonian systems. With the aid of new global perturbation technique, Kovačič [13, 14] investigated the existence of the orbits homoclinic to resonance bands for the Hamiltonian systems and dissipative systems. Yagasaki [15] used the extended subharmonic Melnikov method and the modified homoclinic Melnikov method to examine periodic orbits and homoclinic motions in the coupled oscillators. Feng and Liew [16] canvassed the existence of the Shilnikov type single-pulse homoclinic orbits in the averaged equation which represents the modal interactions between two modes and influence of the fast mode on the slow mode. The global bifurcations and chaotic dynamics were investigated by Zhang et al. [17, 18] for the simply supported rectangular thin plates subjected to the parametrical-external excitation and the parametrical excitation. Yeo and Lee [19] employed the global perturbation technique to probe into the global dynamics of an imperfect circular plate with one-to-one internal resonance and obtained the criteria for chaotic motions of homoclinic orbits. Vakakis [20] adopted subharmonic and homoclinic Melnikov theory to study the strong nonlinear dynamic in a lattice and a lightweight of a nonlinear oscillator with nonlinearity stiffness.Most research is on the Shilnikov type single-pulse global bifurcations and chaotic dynamics of high-dimensional nonlinear systems, but there are researchers investigating the Shilnikov type multipulse homoclinic and heteroclinic bifurcations and chaotic dynamics. So far, there are two theories of the Shilnikov type multipulse chaotic dynamics. One is the energy phase method and the other theory is the extended Melnikov method. Much achievement is made in the former theory of high-dimensional nonlinear systems. Haller and Wiggins [21] first established a simple energy-phase criterion which combined geometric singular perturbation theory, higher-dimensional Melnikov method, and transversality theory. Haller [22] derived a normal form for weak-strong resonance junctions in n-degree-of-freedom, nearly integrable Hamiltonian systems. Haller [23] proposed a new unified theory of orbits homoclinic to resonance bands in a class of near-integrable dissipative systems. Subsequently, Haller and Wiggins [24, 25] further developed the energy phase method to examine the existence of the multipulse homoclinic orbits in the damped-forced nonlinear Schrodinger equation and perturbed the Hamiltonian systems. Haller and Wiggins [26] proved the existence of the Shilnikov type multipulse homoclinic orbits to invariant 3 spheres and utilized the energy-phase method to investigate chaotic dynamics near resonant equilibria in three-degree-of-freedom Hamiltonian systems. Haller [27] derived a universal homoclinic tree describing the bifurcations of the multipulse homoclinic orbits near the intersection of a weaker and a stronger resonance in n-degree-of-freedom, nearly integrable Hamiltonian systems. Haller [28, 29] developed the energy-phase method to detect the existence of the Shilnikov-type multipulse orbits that repeatedly leave and come back to an invariant manifold with two different time scales in the perturbed nonlinear Schrodinger equation. Haller et al. [30] verified the existence of the Shilnikov-type multipulse homoclinic orbits to a spatially independent invariant torus in two coupled nonlinear Schrodinger equations with the damping and the quasiperiodical force. In book [31] published by Haller in 1999, they summarized the energy-phase method and presented a detailed procedure of the application to several problems in mechanics, which include the Shilnikov type multipulse homoclinic and heteroclinic orbits and chaotic dynamics.Up to present, few researchers have made use of the energy phase method to study the Shilnikov type multipulse homoclinic and heteroclinic orbits and chaotic dynamics of high-dimensional nonlinear systems in engineering applications. Malhotra et al. [32] used the energy-phase method to investigate multipulse homoclinic orbits and chaotic dynamics for the motion of flexible spinning discs. McDonald and Namachchivaya [33] applied the Melnikov method, the Shilnikov method, and the energy-phase method to detect the presence of chaotic dynamics of parametrically excited pipes conveying fluid near a zero-to-one resonance. Yao and Zhang [34, 35] utilized the energy-phase method to analyze the Shilnikov type multipulse heteroclinic or homoclinic orbits and chaotic dynamics in a parametrically and externally excited rectangular thin plate and a laminated composite piezoelectric rectangular plate. Yu and Chen [36, 37] made use of the energy-phase method to examine the Shilnikov type multipulse homoclinic orbits of a nonlinear cyclic system and a harmonically excited circular plate. Zhang et al. [38] extended the energy-phase method to six-dimensional system from four-dimensional system. They adopted the energy-phase method for six-dimensional system to delve into multipulse chaotic dynamics of a composite laminated piezoelectric rectangular plate subjected to the transverse, in-plane excitations and the piezoelectric excitation.The study on the second theory of the Shilnikov type multipulse chaotic dynamics was stated by Kovačič and Wettergren [39]. They presented the extended Melnikov method to investigate the existence of the multipulse jumping of homoclinic orbits and chaotic dynamics in resonantly forced coupled pendula. Furthermore, Kaper and Kovačič [40] studied the existence of several classes of the multibump orbits homoclinic to resonance bands for completely integrable Hamiltonian systems subjected to small Hamiltonian amplitude and damped perturbations. Camassa et al. [41] applied the extended Melnikov method to analyze the multipulse jumping of homoclinic and heteroclinic orbits in a class of perturbed Hamiltonian systems. Until recently, Zhang and Yao [42, 43] introduced the extended Melnikov method to the engineering field. They came up with a simplification of the extended Melnikov method in the resonant case and utilized it to analyze the Shilnikov type multipulse homoclinic bifurcations and chaotic dynamics for the nonlinear nonplanar oscillations of the cantilever beam subjected to a harmonic axial excitation and two transverse excitations at the free end. Yao et al. [44] made use of the extended Melnikov method and numerical method to investigate multipulse chaotic dynamics in nonplanar motion of parametrically excited viscoelastic moving belt.Despite extensive research of nonlinear dynamics in the laminated composite piezoelectric plate, multipulse chaotic dynamics has been studied rarely. Few researchers utilized the energy-phase method to investigate the Shilnikov type multipulse chaotic dynamics of high-dimensional nonlinear systems in engineering applications in the past several years. We have previously studied homoclinic bifurcations and multipulse chaotic dynamics of the laminated composite piezoelectric plate under the case of 1 : 3 internal resonances by applying the energy-phase method. In this paper, we have used the energy-phase method to investigate heteroclinic bifurcations and multipulse chaotic dynamics of the laminated composite piezoelectric plate under the case of 1 : 2 internal resonances.It is the purpose of this paper to fill the research gap by investigating the Shilnikov type multipulse chaotic dynamics in the complex motion of the laminated composite piezoelectric plate. Based on the von Karman type equations and the Reddy’s third-order shear deformation plate theory, the Hamilton’s principle is employed to obtain the governing nonlinear equations of the laminated composite piezoelectric rectangular plate with combined parametric excitation and transverse load. We apply the method of multiple scales and Galerkin’s approach to the partial differential governing equations to obtain the four-dimensional averaged equation for the case of 1 : 2 internal resonances and primary parametric resonance. From the averaged equation, the theory of normal form is used to find the explicit formulas of normal form. Finally, we employ the energy-phase method to analyze the Shilnikov type multipulse orbits and chaotic dynamics in the laminated composite piezoelectric plate. The analysis indicates that there exist the Shilnikov type multipulse jumping orbits in the perturbed phase space for the averaged equations. The results from numerical simulation also show that the chaotic motion can occur in the motion of the laminated composite piezoelectric plate, which verifies the analytical prediction. The Shilnikov type multipulse orbits are discovered from the results of numerical simulation. In summary, both theoretical and numerical studies demonstrate that chaos for the Smale horseshoe sense in the motion exists.


## 2. Equations of Motion and Perturbation Analysis

We consider a four-edge simply supported laminated composite piezoelectric rectangular plate, where the length, the width, and the thickness are denoted bya, b, and h, respectively. The laminated composite piezoelectric rectangular plate is subjected to in-plane excitation, transverse excitation, and piezoelectric excitation, as shown in Figure 1. We consider the laminated composite piezoelectric rectangular plate as regular symmetric cross-ply laminates with n layers with respect to principal material coordinates alternatively oriented at 0∘ and 90∘ to the laminate coordinate axes. Some of the layers are made of the piezoelectric materials as actuators, and the other layers are made of fiber-reinforced composite materials. It is assumed that different layers of the symmetric cross-ply composite laminated piezoelectric rectangular plate are perfectly clung to each other, and piezoelectric actuator layers are embedded in the plate. The fiber direction of odd-numbered layers is the x-direction of the laminate. The fiber direction of even-numbered layers is the y-direction of the laminate. Simply supported plate with immovable edges satisfies the symmetry requirement that eliminates the coupling between bending and extension. However, the displacement of x is free to move at the edge of y=0, and the displacement of y is free to move at the edge of x=0. Therefore, the membrane stress is smaller and there exists the coupling between bending and extension. A Cartesian coordinate system Oxyz is located on the middle surface of the composite laminated piezoelectric rectangular plate. Assume that (w,v,u) and (w0,v0,u0) describe the displacements of an arbitrary point and a point on the middle surface of the composite laminated piezoelectric rectangular plate in the x, y, and z directions, respectively. It is also assumed that in-plane excitations of the composite laminated piezoelectric rectangular plate are loaded along the y-direction at x=0 and the x-direction at y=0 with the form of q0+qxcosΩ1t and q1+qycosΩ2t, respectively. Transverse excitation loaded to the composite laminated piezoelectric rectangular plate is expressed as q=q3cosΩ3t. The dynamic electrical loading is represented by Ez=EzcosΩ4t.Figure 1
The model of a laminated composite piezoelectric rectangular plate is given.In this paper, Reddy’s third-order shear deformation description of the displacement field is adopted:(1a)u(x,y,z,t)=u0(x,y,t)+zϕx(x,y,t)-z343h2(ϕx+∂w0∂x),(1b)v(x,y,z,t)=v0(x,y,t)+zϕy(x,y,t)-z343h2(ϕy+∂w0∂y),(1c)w(x,y,z,t)=w0(x,y,t),where (u0,v0,w0) are the deflection of a point on the middle surface, (u,v,w) are the displacement components along the (x,y,z) coordinate directions, and ϕx and ϕy denote the rotation components of normal to the middle surface about the y and x axes, respectively.The nonlinear strain-displacement relations are assumed to have the following form:(2)εxx=∂u∂x+12(∂w∂x)2,εxz=12(∂u∂z+∂w∂x),εxy=12(∂u∂y+∂v∂x+∂w∂x∂w∂y)εyy=∂v∂y+12(∂w∂y)2,εyz=12(∂v∂z+∂w∂y),εzz=∂w∂z.Stress constitutive relations are presented as follows:(3)σij=σijklsεkl-eijkEk,(i,j,k,l=x,y,z),
where σij and εkl denote the mechanical stresses and strains in extended vector notation, σijklsrepresents the elastic stiffness tensor, Ek stands for the electric field vector, and eij is the piezoelectric tensor.According to Hamilton’s principle, the nonlinear governing equations of motion for the composite laminated piezoelectric rectangular plate are given in the previou studies as follows [10, 35]:(4a)∂Nxx∂x+∂Nxy∂y=I0u¨0+J1ϕ¨x-c1I3∂w¨0∂x,(4b)∂Nxy∂x+∂Nyy∂y=I0v¨0+J1ϕ¨y-c1I3∂w¨0∂y,(4c)∂Q-x∂x+∂Q-y∂y+∂∂x(Nxx∂w0∂x+Nxy∂w0∂y)+∂∂y(Nxy∂w0∂x+Nyy∂w0∂y)+c1(∂2Pxx∂x2+2∂2Pxy∂x∂y+∂2Pyy∂y2)-γw˙0+q=I0w¨0-c12I6(∂2w¨0∂x2+∂2w¨0∂y2)+c1[I3(∂u¨0∂x+∂v¨0∂y)+I4(∂ϕ¨x∂x+∂ϕ¨y∂y)],(4d)∂M-xx∂x+∂M-xy∂y-Q-x=J1u¨0+k2ϕ¨x-c1J4∂w¨0∂x,(4e)∂M-xy∂x+∂M-yy∂y-Q-y=J1v¨0+k2ϕ¨y-c1J4∂w¨0∂y,where the dot represents the partial differentiation with respect to time t, the comma denotes the partial differentiation with respect to a specified coordinate, γ is the damping coefficient, and all kinds of inertias in (4a), (4b), (4c), (4d), and (4e) are calculated by
(5)Ii=∑k=1N∫zkzk+1ρkzidz,Ji=Ii-c1Ii+2,K2=I2-2c1I4.In addition, the stress resultants are represented as follows:(6a)Nxx=A11[∂u0∂x+12(∂w0∂x)2]+A12[∂v0∂y+12(∂w0∂y)2]-Nxxp,(6b)Nyy=A21[∂u0∂x+12(∂w0∂x)2]+A22[∂v0∂y+12(∂w0∂y)2]-Nyyp,(6c)Nxy=A66(∂v0∂x+∂w0∂x∂w0∂y+∂u0∂y)-Nxyp,(6d)Mxx=(D11-c1F11)∂ϕx∂x+(D12-c1F12)∂ϕy∂y-c1F11∂2w0∂x2-c1F12∂2w0∂y2-Mxxp,(6e)Myy=(D21-c1F21)∂ϕx∂x+(D22-c1F22)∂ϕy∂y-c1F21∂2w0∂x2-c1F22∂2w0∂y2-Myyp,(6f)Mxy=(D66-c1F66)(∂ϕx∂x+∂ϕy∂y)-2c1F66∂2w0∂x∂y-Mxyp,(6g)Pxx=(F11-c1H11)∂ϕx∂x+(F12-c1H12)∂ϕy∂y-c1H11∂2w0∂x2-c1H12∂2w0∂y2,(6h)Pyy=(F21-c1vH21)(∂ϕx∂x+∂ϕy∂y)-c1H21∂2w0∂x2-c1H22∂2w0∂y2,(6i)Pxy=(P66-c1H66)(∂ϕx∂x+∂ϕy∂y)-2c1H66∂2w0∂x∂y,(6j)Q-x=(A44-c2D44)ϕx+(A44-c2D44)∂w0∂y,(6k)Q-y=(A55-c2D55)ϕx+(A55-c2D55)∂w0∂x,(6l)NxP=∑k=1N∫zkzk+1Q-11ke31kEzdz,(6m)NyP=∑k=1N∫zkzk+1Q-22ke32kEzdz,where NiP=NiPcos(Ω4t) (i=x,y) represent the piezoelectric stress resultants, e31 and e32 are piezoelectric constants, Ez denotes electric field, and Aij, Bij, Dij, Eij, Fij, and Hij are, respectively, the stiffness elements of the laminated composite piezoelectric rectangular plate, which are denoted as(7a)(Aij,Bij,Dij,Eij,Fij,Hij)=∑k=1N∫zkzk+1Q-ijk(1,z,z2,z3,z4,z6)dz,(i,j=1,2,6),(7b)(Aij,Dij,Fij)=∑k=1N∫zkzk+1Q-ijk(1,z2,z4)dz,(i,j=4,5).Substituting (6a), (6b), (6c), (6d), (6e), (6f), (6g), (6h), (6i), (6j), (6k), (6l), and (6m) into (4a), (4b), (4c), (4d), and (4e) yields governing equations of motion in terms of generalized displacements (u0,v0,w0,ϕx,ϕy) as follows:(8a)A11∂2u0∂x2+A66∂2u0∂y2+(A12+A66)∂2v0∂x∂y+A11∂w0∂x∂2w0∂x2+A66∂w0∂x∂2w0∂y2+(A12+A66)∂w0∂y∂2w0∂x∂y=I0u¨0+J1ϕ¨x-c1I3∂w¨0∂x,(8b)A66∂2v0∂x2+A22∂2v0∂y2+(A21+A66)∂2u0∂x∂y+A66∂w0∂y∂2w0∂x2+A22∂w0∂y∂2w0∂y2+(A21+A66)∂w0∂x∂2w0∂x∂y=I0v¨0+J1ϕ¨y-c1I3∂w¨0∂y,(8c)A66∂w0∂x∂2u0∂y2-H22c12∂4w0∂y4+c1(2F66+F12-2H66c1-H12c1)∂3ϕy∂y∂x2+c1(F22-H22c1)∂3ϕy∂y3-H11c12∂4w0∂x4+A11∂w0∂x∂2u0∂x2+(F44c22-2D44c2+A44)∂ϕy∂y+c1(F21+2F66-H21c1-2H66c1)∂3ϕx∂y2∂x-c12(H21+4H66+H12)∂4w0∂y2∂x2+(A44-NyPcos(Ω4t)+F44c22-2D44c2)∂2w0∂y2-∂NyP∂ycos(Ω2t)∂w0∂y+(A21+4A66+A12)∂w0∂x∂w0∂y∂2w0∂y∂x+c1(F11-H11c1)∂3ϕx∂x3+(A21+A66)∂w0∂y∂2u0∂y∂x+A21∂u0∂x∂2w0∂y2A66∂w0∂y∂2v0∂x2+A22∂w0∂y∂2v0∂y2+12(A12+2A66)(∂w0∂y)2∂2w0∂x2+A22∂2w0∂y2∂v0∂y+(A12+A66)∂w0∂x∂2v0∂y∂x+12(A21+2A66)∂2w0∂y2(∂w0∂x)2+32A11(∂w0∂x)2∂2w0∂x2+A11∂2w0∂x2∂u0∂x+A12∂2w0∂x2∂v0∂y+2A66∂2w0∂y∂x∂v0∂x+2A66∂2w0∂y∂x∂u0∂y+32A22(∂w0∂y)2∂2w0∂y2+(A55+qxcos(Ω1t)-NxPcos(Ω4t)+F55c22-2D55c2)∂2w0∂x2-∂NxP∂xcos(Ω3t)∂w∂x+(F55c22-2D55c2+A55)∂ϕx∂x-qcos(Ω3t)+kf∂w0∂t=I0w¨0-c12I6(∂2w¨0∂x2+∂2w¨0∂y2)+c1I3(∂u¨0∂x+∂v¨0∂y)+c1J4(∂ϕ¨x∂x+∂ϕ¨y∂y),(8d)(D11-2F11c1+H11c12)∂2ϕx∂x2+(D66-2F66c1+H66c12)∂2ϕx∂y2-c1(F11-H11c1)∂3w0∂x3-(F55c22-2D55c2+A55)∂w0∂x+(D12+D66+H66c12-2F66c1+H12c12-2F12c1)∂2ϕy∂y∂x-c1(2F66+F12-2H66c1-H12c1)∂3w0∂y2∂x+(2D55c2-A55-F55c22)ϕx=J1u¨0+K2ϕ¨x-c1J4∂w¨0∂x,(8e)(D66-2F66c1+H66c12)∂2ϕy∂x2-c1(F21+2F66-H21c1-2H66c1)∂3w0∂y∂x2+(H21c12+D66+D21-2F21c1+H66c12-2F66c1)∂2ϕx∂y∂x+(H22c12+D22-2F22c1)∂2ϕy∂y2-c1(F22-H22c1)∂3w0∂y3-(F44c22-2D44c2+A44)∂w0∂y+(2D44c2-F44c22-A44)ϕy=J1v¨0+K2ϕ¨y-c1J4∂w¨0∂y.The simply supported boundary conditions of the composite laminated piezoelectric rectangular plate can be represented as follows [10, 35, 45]: (9a)x=0:w=0,ϕy=0,Nxy=0,Mxx=0,(9b)x=a:w=0,ϕy=0,Nxy=0,Mxx=0,(9c)y=0:w=0,ϕx=0,Nxy=0,Myy=0,(9d)y=b:w=0,ϕx=0,Nxy=0,Myy=0,(9e)∫0hNxx|x=0dz=-∫0h(q0+qxcosΩ1t)dz,(9f)∫0hNyy|y=0dz=-∫0h(q1+qycosΩ2t)dz.The boundary condition (9f) includes the influence of the in-plane load. We consider complicated nonlinear dynamics of the composite laminated piezoelectric rectangular plate in the first two modes of u0, v0, w0, ϕx, and ϕy. It is desirable that we select an appropriate mode function to satisfy the boundary condition. Thus, we can write u0, v0, w0, ϕx, and ϕy in the following forms:(10a)u0=u1(t)cosπx2acosπy2b+u2(t)cos3πx2acosπy2b,(10b)v0=v1(t)cosπy2bcosπx2a+v2(t)cosπy2bcos3πx2a,(10c)w0=w1(t)sinπxasinπyb+w2(t)sin3πxasinπyb,(10d)ϕx=ϕ1(t)cosπxasinπyb+ϕ2(t)cos3πxasinπyb,(10e)ϕy=ϕ3(t)cosπybsinπxa+ϕ4(t)cosπybsin3πxa.Based on the study given by Nosir and Reddy and Bhimaraddi [46, 47], we neglect the in-plane and rotary inertia terms in (8a), (8b), (8d), and (8e) since their influences are small compared to the transverse inertia term. By means of the Galerkin method, substituting (10a), (10b), (10c), (10d), and (10e) into (8a), (8b), (8c), (8d), and (8e), integrating and neglecting all inertia terms in (8a), (8b), (8d), and (8e), we obtain the expressions of u1, u2, v1, v2, ϕ1, ϕ2, ϕ3, and ϕ4 in terms of w1 and w2 as follows:(11a)u1=k1w12+k2w22+k3w1w2,(11b)u2=k4w12+k5w22+k6w1w2,(11c)v1=k7w12+k8w22+k9w1w2,(11d)v2=k10w12+k11w22+k12w1w2,(11e)ϕ1=k19w1,ϕ2=k20w2,(11f)ϕ3=k21w1,ϕ4=k22w2,where the coefficients presented in (11a), (11b), (11c), (11d), (11e), and (11f) can be found in Appendix A.In order to obtain the dimensionless governing equations of motion, we introduce the transformations of the variables and parameters:(12)u-=u0a,v-=v0b,w-=w0h,ϕ-x=ϕx,ϕ-y=ϕy,x-=xa,y-=yb,q-=b2Eh3q,q-x=b2Eh3qx,q-y=b2Eh3qy,t-=π2(Eabρ)1/2t,Ω-i=1π2(abρE)1/2Ωi(i=1,2),A-ij=(ab)1/2Eh2Aij,B-ij=(ab)1/2Eh3Bij,D-ij=(ab)1/2Eh4Dij,E-ij=(ab)1/2Eh5Eij,F-ij=(ab)1/2Eh6Fij,H-ij=(ab)1/2Eh8Hij,I-i=1(ab)(i+1)/2ρIi.Based on the practical work condition of laminated composite rectangular plates, and theoretical and experimental studies obtained by Reddy and Mitchell [3, 45], it is known that the nonlinear transverse vibration of laminated composite rectangular plates occupies the main aspect of the dynamical characteristics. The transverse vibration is far greater than the other additional vibrations in the u0 and v0 directions, respectively. Therefore, the equations for u0 and v0 can be neglected. We mainly consider the nonlinear transverse vibration of laminated composite rectangular plates.For simplicity, we drop the overbar in the following analysis. Substituting ((9a), (9b), (9c), (9d), (9e), and (9f))–(12) into (8c) and applying the Galerkin procedure, we obtain the governing equations of motion of the composite laminated piezoelectric rectangular plate for the dimensionless as follows:(13a)w¨1+μ1w˙1+ω12w1+(a2cosΩ1t+a3cosΩ2t-a4cosΩ4t)w1+a5w12w2+a6w22w1+a7w13+a8w23=f1cosΩ3t,(13b)w¨2+μ2w˙2+ω22w2+(b2cosΩ1t+b3cosΩ2t+b4cosΩ4t)w2+b5w22w1+b6w12w2+b7w23+b8w13=f2cosΩ3t,where the coefficients presented in (13a) and (13b) are given in Appendix B.The above equations include the cubic terms, in-plane excitation, transverse excitation, and piezoelectric excitation. Equations (13a) and (13b) can describe the nonlinear transverse oscillations of the composite laminated piezoelectric rectangular plate. We only study the case of primary parametric resonance and 1 : 2 internal resonances. In this resonant case, there are the following resonant relations:
(14)ω12=ω24+εσ1,ω22=ω2+εσ2,Ω3=ω,Ω1=Ω2=Ω4=ω,ω2≈2ω1,
where σ1 and σ2 are two detuning parameters.The method of multiple scales [48] is employed in (13a) and (13b) to find the uniform solutions in the following form:(15a)w1(t,ε)=x0(T0,T1)+εx1(T0,T1)+⋯,(15b)w2(t,ε)=y0(T0,T1)+εy1(T0,T1)+···,where T0=t, T1=εt.Substituting (14), (15a), and (15b) into (13a) and (13b) and balancing the coefficients of the corresponding powers of ε on the left-hand and right-hand sides of equations, the four-dimensional averaged equations in the Cartesian form are obtained as follows:(16a)x˙1=-12μ1x1-12σ1x2+14(a2+a3-a4)x2-a6x2(x32+x42)-32a7x2(x12+x22),(16b)x˙2=-12μ1x2+12σ1x1+14(a2+a3-a4)x1+a6x1(x32+x42)+32a7x1(x12+x22),(16c)x˙3=-12μ2x3-14σ2x4-12b6x4(x12+x22)-34b7x4(x32+x42),(16d)x˙4=-12μ2x4+14σ2x3+12b6x3(x12+x22)+34b7x3(x32+x42)-18f2.


## 3. Computation of Normal Form

In order to assist the analysis of the Shilnikov type multipulse orbits and chaotic dynamics of the laminated composite piezoelectric rectangular plate, it is necessary to reduce the averaged equations (16a), (16b), (16c), and (16d) to a simpler normal form. It is found that there are Z2⊕Z2 and D4 symmetries in the averaged equations (16a), (16b), (16c), and (16d) without the parameters. Therefore, these symmetries are also held in normal form.We take into account the excitation amplitudef2 as a perturbation parameter. Amplitude f2 can be considered as an unfolding parameter when the Shilnikov type multipulse orbits are investigated. Obviously, when we do not consider the perturbation parameter, (16a), (16b), (16c), and (16d) become(17a)x˙1=-12μ1x1+(f0-12σ1)x2-a6x2(x32+x42)-32a7x2(x12+x22),(17b)x˙2=-12μ1x2+(f0+12σ1)x1+a6x1(x32+x42)+32a7x1(x12+x22),(17c)x˙3=-12μ2x3-14σ2x4-12b6x4(x12+x22)-34b7x4(x32+x42),(17d)x˙4=-12μ2x4+14σ2x3+12b6x3(x12+x22)+34b7x3(x32+x42),where f0=(1/4)(a2+a3-a4).It is obviously known that (17a), (17b), (17c), and (17d) have a trivial zero solution (x1,x2,x3,x4)=(0,0,0,0) at which the Jacobian matrix can be written as
(18)J=DxX=[-12μ1f0-12σ100f0+12σ1-12μ10000-12μ2-14σ20014σ2-12μ2].The characteristic equation corresponding to the trivial zero solution is(19)(λ2+μ1λ+14μ12+14σ12-f02)×(λ2+μ2λ+14μ22+116σ22)=0.Let(20)Δ1=14μ12+14σ12-f02,Δ2=14μ22+116σ22.
When μ1=μ2=0, Δ1=0, and Δ2=(1/16)σ22>0 are satisfied simultaneously, (17a), (17b), (17c), and (17d) have one nonsemisimple double zero and a pair of pure imaginary eigenvalues:
(21)λ1,2=0,λ3,4=±iω-,
where ω-2=(1/16)σ22.Consideringσ-1, μ1, μ2, and f2 as the perturbation parameters, letting σ1=2(σ-1-f0), and setting f0=1/2, then, averaged equations (17a), (17b), (17c), and (17d) without the perturbation parameters become the following forms:(22a)x˙1=x2-a6x2(x32+x42)-32a7x2(x12+x22),(22b)x˙2=a6x1(x32+x42)+32a7x1(x12+x22),(22c)x˙3=-14σ2x4-12b6x4(x12+x22)-34b7x4(x32+x42),(22d)x˙4=14σ2x3+12b6x3(x12+x22)+34b7x3(x32+x42).According to (22a), (22b), (22c), and (22d), we have
(23)A=[01000000000-14σ20014σ20].Executing the Maple program given by Zhang et al. [49], the 3-order normal form of (22a), (22b), (22c), and (22d) is obtained as(24a)y˙1=y2,(24b)y˙2=a6y1(y32+y42)+32a7y13,(24c)y˙3=-14σ2y4-12b6y12y4-34b7y4(y32+y42),(24d)y˙4=14σ2y3+12b6y12y3+34b7y3(y32+y42).The nonlinear transformation used here is given as follows:(25a)x1=y1-14a7y13,(25b)x2=y2+32a7y23+34a7y12y2+a6y2y32+a6y2y42,(25c)x3=y3-12b6y1y2y4,(25d)x4=y4+12b6y1y2y3.The above nonlinear transformation is computed through the Maple program given by Zhang et al. [49], and completely agrees with those presented by Yu et al. [50]. Therefore, a simpler 3-order normal form with the parameters for the averaged equations (16a), (16b), (16c), and (16d) is obtained as follows:(26a)y˙1=-μ-1y1+(1-σ-1)y2,(26b)y˙2=σ-1y1-μ-1y2+a6y1(y32+y42)+32a7y13,(26c)y˙3=-μ-2y3-σ-2y4-12b6y12y4-34b7y4(y32+y42),(26d)y˙4=σ-2y3-μ-2y4+12b6y12y3+34b7y3(y32+y42)-f-2,where the coefficients are μ-1=(1/2)μ1, μ-2=(1/2)μ2, σ-2=(1/4)σ2, and f-2=(1/8)f2, respectively.Further, let(27)y3=Icosγ,y4=Isinγ.
Substituting (27) into (26a), (26b), (26c), and (26d) yields(28a)y˙1=-μ-1y1+(1-σ-1)y2,(28b)y˙2=σ-1y1-μ-1y2+a6y1I2+32a7y13,(28c)I˙=-μ-2I-f-2sinγ,(28d)Iγ˙=σ-2I+12b6y12I+34b7I3-f-2cosγ.In order to get the unfolding of (28a), (28b), (28c), and (28d), a linear transformation is introduced:
(29)[y1y2]=2|a6||b6|[1-σ-10μ-11][u1u2].Substituting (29) into (28a), (28b), (28c), and (28d), and canceling nonlinear terms including the parameter σ-1 yield the unfolding as follows:(30a)u˙1=u2,(30b)u˙2=-μu1-μ3u2+η1u13+a6u1I2,(30c)I˙=-μ-2I-f-2sinγ,(30d)Iγ˙=σ-2I+a6u12I+α2I3-f-2cosγ,where μ=μ-12-σ-1(1-σ-1), μ3=2μ-1, η1=3a6a7/b6, and α2=(3/4)b7.The scale transformations to be introduced into (30a), (30b), (30c), and (30d) are
(31)μ-2⟶εμ-2,μ3⟶εμ3,f-2⟶εf-2,η1⟶η1,α2⟶α2,a6⟶a6.
Then, the normal forms (30a), (30b), (30c), and (30d) can be rewritten in the form with the perturbations as follows:(32a)u˙1=∂H∂u2+εgu1=u2,(32b)u˙2=-∂H∂u1+εgu2=-μu1+η1u13+a6u1I2-εμ3u2,(32c)I˙=∂H∂γ+εgI=-εμ-2I-εf-2sinγ,(32d)Iγ˙=-∂H∂I+εgγ=σ-2I+α2I3+a6Iu12-εf-2cosγ,where the Hamiltonian function H is of the following form:
(33)H(u1,u2,I,γ)=12u22+12μu12-14η1u14-12a6I2u12-12σ-2I2-14α2I4,
and gu1, gu2, gI, and gγ are the linear perturbation terms induced by the dissipative effects:
(34)gu1=0,gu2=-μ3u2,gI=-μ-2I,gγ=0.


## 4. Heteroclinic Bifurcations

In this section, we focus on studying the nonlinear dynamics characteristic of the unperturbed system. Whenε=0, it can be seen that the system from (32a), (32b), (32c), and (32d) is an uncoupled two-degree-of-freedom nonlinear system. The variable I in the subspace (u1,u2) of (32a), (32b), (32c), and (32d) becomes a parameter since I˙=0. Consider the first two decoupled equations of (32a), (32b), (32c), and (32d),(35a)u˙1=u2,(35b)u˙2=-μu1+η1u13+a6I2u1.Sinceη1>0, (35a) and (35b) can exhibit the heteroclinic bifurcations. It is obvious from (35a) and (35b) that when μ-a6I2<0, the only solution to (35a) and (35b) is the trivial zero solution, (u1,u2)=(0,0), which is the saddle point. On the curve defined by μ=a6I2, that is,
(36)μ-12=σ-1(1-σ-1)+a6I2,
or
(37)I1,2=±[μ-12-σ-1(1-σ-1)a6]1/2,
the trivial zero solution bifurcates into three solutions through a pitchfork bifurcation, which are given by q0=(0,0) and q±(I)=(B,0), respectively, where
(38)B=±{1η1[μ-12-σ-1(1-σ-1)-a6I2]}1/2.From the Jacobian matrix evaluated at the nonzero solutions, it can be found that the singular pointsq±(I) are the saddle points. It is observed that I and γ actually represent the amplitude and the phase of vibrations, respectively. Therefore, we assume that I≥0 and (37) become
(39)I1=0,I2=[μ-12-σ-1(1-σ-1)a6]1/2,
such that for all I∈[I1,I2], (35a) and (35b) have two hyperbolic saddle points, q±(I), which are connected by a pair of heteroclinic orbits, u±h(T1,I), that is, limT1→±∞u±h(T1,I)=q±(I). Thus, in the full four-dimensional phase space, the set defined by
(40)M={(u,I,γ)∣u=q±(I),I1<I<I2,0≤γ<2π}
is a two-dimensional invariant manifold. From the results obtained by Kovačič [13, 14], it is known that the two-dimensional invariant manifold M is normally hyperbolic. The two-dimensional normally hyperbolic invariant manifold M has the three-dimensional stable and unstable manifolds represented as Ws(M) and Wu(M), respectively. The existence of the heteroclinic orbit of (35a) and (35b) to q±(I)=(B,0) indicates that Ws(M) and Wu(M) intersect nontransversally along a three-dimensional heteroclinic manifold denoted by Γ, which can be written as
(41)Γ={(u,I,γ)∣u=u±h(T1,I),I1<I<I2,∫0T1DIH(u±h(T1,I),I)ds+γ0γ=∫0T1DIH(u±h(T1,I),I)ds+γ0}.We analyze the dynamics of the unperturbed system of (32a), (32b), (32c), and (32d) restricted to M. Considering the unperturbed system of (32a), (32b), (32c), and (32d) restricted to M yields(42a)I˙=0,(42b)Iγ˙=DIH(q±(I),I),I1<I<I2,where
(43)DIH(q±(I),I)=-∂H(q±(I),I)∂I=σ-2I+α2I3+a6Iq±2(I).From the results obtained by Kovačič [13, 14], it is known that if DIH(q±(I),I)≠0; then I= constant is called a periodic orbit, and if DIH(q±(I),I)=0; then I= constant is known as a circle of the singular points. Any value of I∈[I1,I2] at which DIH(q±(I),I)=0 is a resonant value I and these singular points are resonant singular points. We denoted a resonant value by Ir such that
(44)DIH(q±(I),I)=σ-2Ir+α2Ir3+a6η1[μ-12-σ-1(1-σ-1)-a6Ir2]Ir=0.
Then, we obtain
(45)Ir=±{σ-2η1+a6[μ-12-σ-1(1-σ-1)]a62-α2η1}1/2.Figure2 shows the geometry structure of the stable and unstable manifolds of M in the full four-dimensional phase space for the unperturbed system of (32a), (32b), (32c), and (32d). Since γ represents the phase of oscillations, when I=Ir, the phase shift Δγ of  oscillations is defined by
(46)Δγ=γ(+∞,Ir)-γ(-∞,Ir).The geometric structure of manifoldsM; Ws(M) and Wu(M) are given in the full four-dimensional phase space.
(a)
(b)The physical interpretation of the phase shift is the phase difference between the two end points of the orbit. In the subspace(u1,u2), there exists a pair of the heteroclinic orbits connecting to the two saddles. Therefore, the homoclinic orbit in subspace (I,γ) represents, in fact, a heteroclinic connecting in the full four-dimensional space (u1,u2,I,γ). The phase shift denotes the difference of the value γ as a trajectory leaving and returning to the basin of the attraction of M. We will use the phase shift in the subsequent analysis to obtain the condition for the existence of the Shilnikov type multipulse orbit. The phase shift will be calculated later in the heteroclinic orbit analysis.We consider the heteroclinic orbits of (35a) and (35b). Let ε1=μ-a6I2 and let μ3=ε2; then (35a) and (35b) can be rewritten as(47a)u˙1=u2,(47b)u˙2=-ε1u1+η1u13-εε2u2.Set ε=0; then (47a) and (47b) are a system with the Hamiltonian function:
(48)H-(u1,u2)=12u22+12ε1u12-14η1u14.WhenH-=ε12/(4η1), there exists a heteroclinic loop Γ0 which consists of the two hyperbolic saddles q± and a pair of heteroclinic orbits u±(T1). In order to calculate the phase shift and the energy difference function, we obtain the equations of a pair of heteroclinic orbits given by(49a)u1(T1)=±ε1η1tanh(2ε12T1),(49b)u2(T1)=±ε12η1sech2(2ε12T1).We turn our attention to the computation of the phase shift. Substituting the first equation of (49a) and (49b) into the fourth equation of the unperturbed system of (32a), (32b), (32c), and (32d) yields
(50)γ˙=σ-2+α2I2+a6ε1η1tanh2(2ε12T1).
Integrating (50) yields
(51)γ(T1)=ωrT1-a62ε1η1tanh(2ε12T1)+γ0,
where ωr=σ-2+α2I2+(ε1a6/η1).AtI=Ir, there is ωr≡0. Therefore, the phase shift may be expressed as
(52)Δγ=[-2a62ε1η1]I=Ir=-2a6η12[μ-12-σ-1(1-σ-1)-a6Ir2].


## 5. Dissipative Perturbations of Homoclinic Loop

In this section, the effects of small perturbation terms on the unperturbed system are analyzed in detail. We now analyze dynamics of the perturbed system and the effect of small perturbations onM. Based on the analysis given by Kovačič [13, 14], it is known that M along with its stable and unstable manifolds is invariant under small, sufficiently differentiable perturbations. It is noticed that q±(I) in (35a) and (35b) maintain the characteristic of the hyperbolic singular point under small perturbations, in particular, M→Mε. Therefore, we obtain
(53)M=Mε={(u,I,γ)∣u=q±(I),I1<I<I2,0≤γ<2π}.Considering the last two equations of (32a), (32b), (32c), and (32d) yields(54a)I˙=-μ-2I-f-2sinγ,(54b)γ˙=σ-2+α2I2+a6u12-f-2cosγI.It is known from the above analysis that the last two equations of (32a), (32b), (32c), and (32d) are of a pair of pure imaginary eigenvalues. Therefore, the resonance can occur in (54a) and (54b). Also introduce the scale transformations
(55)μ-2⟶εμ-2,I=Ir+εh,f-2⟶εf-2,T1⟶T1ε.Substituting the above transformations into (54a) and (54b) yields(56a)h˙=-μ-2Ir-f-2sinγ-εμ-2h,(56b)γ˙=-2δη1Irh-ε(f-2Ircosγ+δη1h2),where δ=a62-α2η1.Whenε=0, (56a) and (56b) become(57a)h˙=-μ-2Ir-f-2sinγ,(57b)γ˙=-2δη1Irh.The unperturbed system from (57a) and (57b) is a Hamilton system with the following function:
(58)H^D(h,γ)=-μ-2Irγ+f-2cosγ+δη1Irh2.
The singular points of (57a) and (57b) are represented as
(59)P0=(0,γc)=(0,-arcsinμ-2Irf-2),Q0=(0,γs)=(0,π+arcsinμ-2Irf-2).Based on the characteristic equations evaluated at the two singular pointsP0 and Q0, we can know the stabilities of these singular points. The Jacobian matrix of (57a) and (57b) is
(60)J=[0-f-2cosγ-2δη1Ir0].
The characteristic equation corresponding to the singular point P0 is obtained as
(61)λ2-2δη1Irf-2cosγc=0.When the condition(2δ/η1)f-2Ircosγc<0 is satisfied, (57a) and (57b) have a pair of pure imaginary eigenvalues. Therefore, it is known that the singular point P0 is a center point.The characteristic equation corresponding to the singular pointQ0 is represented by
(62)λ2-2δη1f-2Ircosγs=0.
When the condition (2δ/η1)f-2Ircosγs>0 is satisfied, (57a) and (57b) have two real, unequal, and opposite sign eigenvalues. Therefore, the singular point Q0 is a saddle which is connected to itself by a homoclinic orbit. The phase portrait of the system from (57a) and (57b) is shown in Figure 3(a).Dynamics on the normally hyperbolic manifold is described; (a) the unperturbed case; (b) the perturbed case.
(a)
(b)It is found that for the sufficiently small parameterε, the singular point Q0 remains a hyperbolic singular point Qε of the saddle stability type. It is known that the Jacobian matrix of the linearization of (56a) and (56b) is of the following form:
(63)Jpε=[-εμ-2-f-2cosγc-2δη1Ir-2δη1hεεf-2Irsinγc].
or
(64)Jpε=[-εμ-2-f-2cosγc-2δη1Ir-εμ-2].Based on (64), we find that the leading order term of the trace in the linearization of (56a) and (56b) is less than zero inside the homoclinic loop. Therefore, for small perturbations, the singular point P0 becomes a hyperbolic sink Pε. The phase portrait of the perturbed system from (56a) and (56b) is also depicted in Figure 3(b).Using the function (58), at h=0, the estimate of the basin of attractor for γmin is obtained as
(65)-μ-2Irγmin+f-2cosγmin=-μ-2Irγs+f-2cosγs.
Substituting γs in (59) into (65) yields
(66)γmin-f-2μ-2Ircosγmin=π+arcsinμ-2Irf-2+f-22-μ-22Ir2μ-2Ir.Define an annulusAε near I=Ir as
(67)Aε={(u1,u2,I,γ)∣u1=B,u2=0,|I-Ir|<εC,γ∈Tl},
where C is a constant and is sufficiently large so that the unperturbed homoclinic orbit is enclosed within the annulus.It is noticed that the three-dimensional stable and unstable manifolds ofAε, denoted as Ws(Aε) and Wu(Aε), are subsets of Ws(Mε) and Wu(Mε), respectively. We will indicate that for the perturbed system, the saddle focus Pε on Aε has the multipulse orbits which come out of the annulus Aε and can return to the annulus in the full four-dimensional space. These orbits, which are asymptotic to some invariant manifolds in the slow manifold Mε, leave and enter a small neighborhood of Mε multiple times, and, finally, return and approach an invariant set in Mε asymptotically, as shown in Figure 4.  In Figure 4, this is an example of the three-pulse jumping orbit which depicts the formation mechanism of the multipulse orbits.Figure 4
The Shilnikov type three-pulse homoclinic orbits are obtained.


## 6. Energy Difference Function

In this section, the energy-phase method developed by Haller and Wiggins [21–31] is utilized to discover the existence of the Shilnikov type multipulse orbits and chaotic dynamics of nonlinear vibration for the laminated composite piezoelectric rectangular plate with combined parametric excitation and transverse load. The energy-phase method, which was first presented by Haller and Wiggins [21], is a new global perturbation method, which is different from the global perturbation technique developed by Kovačič and Wiggins [12–14]. Haller and Wiggins [21–31] gave the unperturbed system expressed by a mixture of action-angle variables, which can describe homoclinic and heteroclinic behavior. The key feature of their analysis is that the unperturbed system has a resonance with the action-angle variables. The unperturbed integrable system contains a normally hyperbolic invariant two-dimensional manifold, which has three-dimensional stable and unstable manifolds coinciding along a branch. Action-angle variables illustrate the dynamics on the normally hyperbolic two-dimensional manifold. The situation of the resonance leads to a circle of fixed points on the two-dimensional manifold. Two different points on the circle of fixed points are connected by a heteroclinic orbit. These heteroclinic orbits consist of part of a foliation of the three-dimensional stable and unstable manifolds on the normally hyperbolic two-dimensional manifold. Obviously, the perturbation can dramatically alter the dynamics near the circle of fixed points. Thus, the dynamics restricted to the normally hyperbolic two-dimensional manifold near the resonance becomes hyperbolic fixed points with periodic orbits surrounding homoclinic or heteroclinic trajectories. Each homoclinic or heteroclinic trajectory connects with the hyperbolic points. Haller and Wiggins [21–31] gave conditions for the existence of homoclinic or heteroclinic orbits created near the resonance in the full four-dimensional phase space.The energy-phase method has three steps in the process in this paper. The first step involves analyzing the perturbed dynamics restricted to the normally hyperbolic invariant two-dimensional manifold near the resonance. The used approach is based on nonlinear oscillation theory. The second step requires manifesting the existence of heteroclinic orbits to the normally hyperbolic invariant two-dimensional manifold. A higher-dimensional Melnikov theory is mainly used in the analysis process of this step. The final step is the most difficult and represents the most innovative part of the energy-phase method. Based on the Haller and Wiggins study in [21–31], we prove the existence of multipulse heteroclinic orbits to specific orbits on the two-dimensional manifold in this step. The analysis is complicated by perturbations of nontransversal intersections of manifolds. According to the investigation given by Haller and Wiggins [21–31], we use the geometric singular perturbation theory of Fenichel to analyze the existence of  heteroclinic orbits in the full four-dimensional phase space. Based on foliations of stable and unstable manifolds, Fenichel’s theory combines with energy-type arguments which are suited for the Hamiltonian systems. The analysis shows that the existence of heteroclinic orbits depends only on an energy-phase criterion which is obtained from a reduced, one-degree-of-freedom Hamiltonian system. Therefore, the energy-phase method is combined with geometric singular perturbation theory, higher-dimensional Melnikov method, and transversality theory.The energy-phase method can be utilized to detect the Shilnikov type multipulse heteroclinic orbits to slow manifolds of near-integrable four-dimensional or higher-dimensional nonlinear systems. The key to energy-phase method is the calculation of the energy difference function, which is actually the difference of the Hamiltonian function. Based on the higher-dimensional Melnikov theory and the transversality theory, Haller and Wiggins [21] derived the expression of the energy difference function. In order to illustrate the existence of the multipulse orbits, it is important to obtain the expression of the energy difference function. Using the general expression derived by Haller and Wiggins [21–31] for the class of systems, the energy difference function for the dissipative case is given as follows:
(68)ΔnH^D(γ)=H^D(h,γ+nΔγ)-H^D(h,γ)-n∫A[ddu1gu1(u1,u2,Ir,γ)+ddu2gu2(u1,u2,Ir,γ)]du1du2-n∫∂AlgIdγ,
where
(69)H^D(h,γ+nΔγ)-H^D(h,γ)=-nμ-2IrΔγ+f-2[cos(γ+nΔγ)-cosγ],
where A denotes the area enclosed by a pair of heteroclinic orbits in plane (u1,u2), ∂Al is the boundary of the area A, and Δγ is the phase difference.The coordinate transformation introduced by Haller and Wiggins [21–31] changes the topology structure of heteroclinic orbits when they used the energy-phase method to calculate the energy difference function. In the paper, we did not apply this transformation to calculate the energy difference function in order to maintain the consistency of the topology structure. We directly calculate the energy differential function on the area enclosed by a pair of heteroclinic orbits. These studies are innovations and contributions of this paper.The energy difference function gives the leading order energy difference between the taking off and landing points of ann-pulse orbit. This function contains the phase-type information from the unperturbed part and the energy-type information from the perturbed part of systems (32a), (32b), (32c), (32d).Computing (68) leads to the following expression for the dissipative energy difference function:
(70)ΔnH^D(γ)=f-2[cos(γ+nΔγ)-cosγ]-2nμ3ε1Δγ3a6.We define a dissipative factord=μ3/f-2 such that d gives the relative measure of the dissipative effect with respect to the excitation amplitude. Equation (70) can be rewritten as
(71)ΔnH^D(γ)=f-2[-2sin(γ+nΔγ2)sinnΔγ2-2ndε1Δγ3a6].The zeroes ofΔnH^D(γ) are obtained by solving the following equation:
(72)sin(γ+nΔγ2)=-ndε1Δγ3a6sin(nΔγ/2).Based on (72), the upper bound on the value of the dissipative factor is obtained as follows:
(73)|d|<dmax=3a6nε1|sin(nΔγ/2)Δγ|.
For a given value of the dissipative factor d, the multipulse orbits for all values of the pulse number n are not possible. In the presence of small dissipative effects d<1, we have an upper bound on the maximum number of pulses:
(74)n<nmax=|3a6ε1dΔγ|.
It is obviously known from (74) that the upper bound nmax is inversely proportional to the value of the dissipative factor d.


## 7. Zeroes of the Energy Difference Function

The main aim of the following research focuses on identifying the transverse zeroes of the energy difference function. Define a set that contains all such transverse zeroes:(75)Z^-n={(h,γ)∣ΔnH^D(γ)=0,DγΔnH^D(γ)≠0}.The transverse zeroes of the dissipative energy difference functionΔnH^D(γ) are given by the following equation:
(76)γ+nΔγ2=2mπ+(-1)mα,
where m∈Ζ and
(77)α=-arcsinndε1Δγ3a6sin(nΔγ/2).For anyn satisfying Δγ≠4lπ(l=0,1,2,…), there are two transverse zeroes of the dissipative energy difference function in the interval γ∈[-π/2,3π/2], that is:(78a)γ0,1n=3π2-(nΔγ2+α)mod2π,(78b)γ0,2n=3π2-(π+nΔγ2-α)mod2π.In order to classify these internal orbits, we can use the fact that each orbit has a unique energy level associated with it. The energy level associated with the homoclinic connection ish-0. The energy level of the center type singular point is represented by h-∞. Any periodic orbit enclosed inside the homoclinic connections has an energy level h-n. The energy sequence at the different energy levels is defined as(79a)h-0=H^D(0,γs)=12f-2[-kdIr(π+arcsin12kdIr)-4-k2d2Ir2],(79b)h-n=min[H^D(0,γ0,1n),H^D(0,γ0,2n)],(79c)h-∞=H^D(0,γc)=12f-2[kdIrarcsin12kdIr+4-k2d2Ir2],where μ2/μ1=k, μ-2/f-2=(1/2)kd.From (79a), (79b), and (79c), it is found that there exists the relation h-0<h-n<h-∞  in the energy sequence, which means that the energy increases monotonically as the orbits move towards the center type singular point.Due to the existence of the dissipative perturbations, the center type singular point becomes a hyperbolic sink or a saddle focus. It is known that the existence of the multipulse orbits, which are homoclinic to internal periodic orbits in the slow manifoldMε, implies the existence of the chaos for the Smale horseshoe sense in the simply supported laminated composite piezoelectric rectangular plate [24, 25].


## 8. Existence of Multipulse Shilnikov Orbits

In this section, the main interest lies in finding the existence of the multipulse Shilnikov orbits. It is known that the existence of multipulse orbits, which are homoclinic to internal periodic orbits in the slow manifoldMε, was studied by Haller and Wiggins [21–31]. There are no such internal periodic orbits that lie on the slow manifold in the case of the dissipative perturbations. Due to the existence of the dissipative perturbations, the center type singular point becomes a hyperbolic sink or a saddle focus. In the following analysis, we look for multipulse orbits which are homoclinic to the saddle focus, that is, the Shilnikov type multipulse orbits are negatively and positively asymptotic to the saddle focus itself in full four-dimensional phase space when the dissipative perturbations exist. It is indicated that the presence of N-pulse Shilnikov orbits leads to chaotic dynamics in the sense of the Smale horseshoes. In order to determine the existence of such orbits, we apply the results given by Haller and Wiggins [21–31] to the following research.First, we require the existence of nondegenerate equilibrium points forH^D. In the (h,γ) phase space, based on (59) and μ2/μ1=k, it is known that this singular point is given by
(80)Pc=(hc,γc)=(0,-arcsinkdIr2).Next, we need to compute the zeroes of the dissipative energy difference function at the saddle center(B,0,0,γc) in full four-dimensional phase space (u1,u2,h,γ). In order to obtain the zeroes of ΔnH^D(γc), we solve the following equation:
(81)ΔnH^D(γc)=-2f-2sin(γ+nΔγ2)sinnΔγ2-2nμ3ε1Δγ3a6=0.
Equation (81) leads to
(82)1-k2d2Ir24[1-cos(nΔγ)]=12d[kIrsin(nΔγ)-4nε1Δγ3a6].Based on (82), we obtain the following dissipative parameter:
(83)d=μ3f-2=2[1-cos(nΔγ)]k2Ir2[1-cos(nΔγ)]2+[kIrsin(nΔγ)-4nε1Δγ/3a6]2.
The results obtained above are only valid when the nonzero damping exists, that is,
(84)Δγ≠2mπn,m∈Ζ.In order to indicate the existence of the multipulse Shilnikov type orbits in (32a), (32b), (32c), and (32d), it is also necessary to satisfy the following two nondegeneracy conditions:
(85)DdΔnH^D(γc)≠0,DγcΔnH^D(γc)≠0,
whenever (83) and (84) hold.The first condition in (85) can be rewritten as
(86)Dd[cos(-arcsinkdIr2+nΔγ)-cos(-arcsinkdIr2)-2ndε1Δγ3a6]≠0.
Carrying out the differentiation with respect to d in (86), it is seen that the above condition is violated only if
(87)1-k2d2Ir24[kIrsin(nΔγ)-4nε1Δγ3a6]=k2dIr22[cos(nΔγ)-1].It is obviously found that (82) and (87) cannot be simultaneously satisfied if (84) is also satisfied. Therefore, the second non-degeneracy condition in (85) is satisfied by the existence of transversal zeroes of ΔnH^D. The first condition in (85) is satisfied by the inequality (73).Finally, it is necessary to verify the following condition to ensure that the landing point of anyn-pulse orbit taking off from a slow sink lies in the domain of attraction of one of the sinks. We define a point which is in the interval [0,4π], and is 2kπ apart from the approximate landing point γc+nΔγ. This point is denoted by
(88)γ*n=γs+[γc+nΔγ-γs]mod2π,
where
(89)γs=π+arcsinkdIr2,γc=-arcsinkdIr2.Ifγ*n>γs, then we redefine γ*n by subtracting 2π from it. The aim here is to find 2kπ-translation of the landing point which is closest to the saddle point γs in the region [-π/2,3π/2]. Due to the symmetry of the phase portrait γ→γ-2π, if the energy of point γ*n is more than the energy of the saddle point, that is,
(90)H^D(0,γ*n)>H^D(0,γs),
then, for ε>0 small enough, the landing point γc+nΔγ will be in the domain of attraction of one of the sinks. Calculating this energy condition yields
(91)cosγ*n-cosγs>kdIr2(γ*n-γs).The shape of such an orbit forn=3 is shown in Figure 4. Based on the aforementioned analysis, we summarize the main results in the following theorem which is similar to the theorem obtained by Haller [31].Theorem 1.
For any integern≥1, a positive number ε0(n)>0 and a finite union Cn of codimension-one surfaces in the (ε1,a6,μ3,f-2,ε) parameter space near the set satisfying 0<d<1, (83) and (91) exist, such that for any (ε1,a6,μ3,f-2,ε)∈Cn and 0<ε<ε0(n), the following conclusions hold.
(1) If the integer
(92)R=int[12+nΔγ-2arcsin(kdIr/2)2π]
is even, then each of the saddle-focus type equilibrium points contained in the slow manifold Mε admits two generalized Shilnikov orbits. If the integer R is odd, then there exist two heteroclinic cycles of generalized Shilnikov orbits connecting the saddle foci to each other. In both cases, the n-pulse orbits form pairs that are symmetric with respect to the subspace (u1,u2)=(B,0).
(2) There exists an open set of parameter values containing Cn for which the systems (32a), (32b), (32c), and (32d) admit the Smale horseshoes in its dynamics.
We give an explanation of statement(1) in Theorem 1.
The main goal is to prove that utilizing an integerR determines which region the landing point will fall in. The regions one considers are bounded by the saddle points in the unperturbed system. If one considers the domain (0,4π), then the chosen two regions are(93a)D1=(π+arcsinkdIr2,3π+arcsinkdIr2),(93b)D2=(3π+arcsinkdIr2,5π+arcsinkdIr2).
The landing pointγc+nΔγ will fall in one of these two regions or a 2π-translation of one of these regions. If the landing point falls in D1, and conditions (83) and (91) are also satisfied, then each of the saddle foci will have a homoclinic Shilnikov orbit. If, however, the landing point falls in D2, and conditions (83) and (91) are also satisfied, then there will exist heteroclinic cycles of Shilnikov orbits connecting the two saddle foci to each other. In each case, one will have pairs of orbits that are symmetric with respect to the planes (u1,u2)=(B,0), due to the existence of the symmetry (u1,u2)→(-u1,-u2).
If the landing pointγc+nΔγ is in D1 or a 2π-translation, one has
(94)γc+nΔγ∈(γs+2mπ,γs+2π+2mπ).
Equation (94) leads to
(95)γc+nΔγ-γs∈(2mπ,2mπ+2π).
Hence, the integer
(96)R=int[12+nΔγ-2arcsin(kdIr/2)2π]
is even. Reversely, if the landing point falls in D2 or a 2π-translation, one can similarly obtain
(97)γc+nΔγ∈(γs+2mπ+2π,γs+4π+2mπ).
From (97), one has
(98)γc+nΔγ-γs∈(2(m+1)π,2(m+2)π).
Therefore, it is known that the integer
(99)R=int[12+nΔγ-2arcsin(kdIr/2)2π]
is odd.The proof of the second part of Theorem1 follows from Theorem 2.8.2 presented by Haller [31].


## 9. Numerical Results of Chaotic Motions

Based on the above qualitative analysis for the multipulse orbits and chaotic dynamics of the laminated composite piezoelectric rectangular plate, the conditions of the chaotic motion under the sense of the Smale horses are obtained. The heteroclinic bifurcations of (16a), (16b), (16c), and (16d) appear when η1>0. Therefore, the above theoretical analysis is focused on the situation where there exist heteroclinic bifurcations in (16a), (16b), (16c), and (16d). The parameter η1 is the combination of the parameters a6, a7, and b6, where η1=(3a6a7)/b6. In this section, we have only performed numerical simulations of the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate under heteroclinic bifurcations in order to further verify the theoretical analysis. Consequently, the parameters a6, a7,  and b6 are chosen to satisfy η1>0.We choose the averaged equations (16a), (16b), (16c), and (16d) to conduct numerical simulations. A numerical approach through the computer software MATLAB is utilized to explore the existence of the Shilnikov type multipulse chaotic motions in the laminated composite piezoelectric rectangular plate. Based on the above qualitative analysis, it is found that the damping coefficients μ1and μ2 and the transverse excitation f2 play an important role in the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate. In addition, the parameters a2 and a3 are related to the in-plane excitation in the x-direction and the in-plane excitation in the y-direction, respectively. The parameter a4 is the piezoelectric excitation which reflects the characteristics of the piezoelectric material. Hence, the parameters μ1, μ2, a2, a3, a4, and f2 are selected as the controlling parameters to discover the law for complicated nonlinear dynamics of the laminated composite piezoelectric rectangular plate.Figure5 illustrates the bifurcation diagram of the laminated composite piezoelectric rectangular plate when the excitation f2 varies in the interval f2= 2~130. Other parameters and initial conditions are chosen as σ1=3.31, σ2=3.23, μ1=0.1, μ2=0.1, a2=11.59, a3=21.5, a4=13.7, a6=-14.37, a7=10.39, b6=-5.68, b7=-14.94, x10=-0.01, x20=-0.09, x30=-0.05, and x40=-0.6. Figures 5(a) and 5(b) represent the bifurcation diagram on the plane (x1,f2) and (x3,f2), respectively. It is observed from Figure 5 that the excitation f2 is an important parameter that influences the nonlinear dynamic responses of the laminated composite piezoelectric rectangular plate. Figure 5 shows that the chaotic motion of the laminated composite piezoelectric rectangular plate appears when the excitation f2 varies in the interval f2= 20~78, followed by a periodic motion of that. When excitation f2 continues to increase, Figure 5 presents that there is the periodic motion of the laminated composite piezoelectric rectangular plate.The bifurcation diagram is obtained for the excitationsf2= 2~130, x10=-0.01, x20=-0.09, x30=-0.05, x40=-0.6, σ1=3.31, σ2=3.23, μ1=0.1, μ2=0.1, a2=11.59, a3=21.5, a4=13.7, a6=-14.37, a7=10.39, b6=-5.68, and b7=-14.94; (a) the bifurcation diagram on the plane (x1,f2); (b) the bifurcation diagram on the plane (x3,f2).
(a)
(b)We, respectively, study the impact of the two damping parameters since the averaged equations (16a), (16b), (16c), and (16d) have two different damping coefficients. Figure 6 is the bifurcation diagram of the laminated composite piezoelectric rectangular plate with the damping coefficient μ1. Figure 6 demonstrates that the system is beginning to enter into the region of the chaotic motion and subsequently comes into the region of the periodic motion, as the damping coefficient μ1 varies in the interval μ1= 0.06~0.6. Other parameters and initial conditions are the same as those in Figure 5 when the excitation is chosen as f2=55.1. Figures 6(a) and 6(b) describe the nonlinear motion of the laminated composite piezoelectric rectangular plate on the planes (x1,μ1) and (x3,μ1), respectively, as well as the impact of the damping coefficient μ1 on the system.The bifurcation diagram is obtained for the damping coefficientsμ1= 0.06~0.6 and f2=55.1; (a) the bifurcation diagram on the plane (x1,μ1); (b) the bifurcation diagram on the plane (x3,μ1).
(a)
(b)Figure7 represents the bifurcation diagram on the damping coefficient μ2. Figure 7 shows that the laminated composite piezoelectric rectangular plate has been in a state of chaotic motion as the damping coefficient μ2 changes in the interval μ2= 0.06~0.6. Other parameters and initial conditions are the same as those in Figure 6 when the damping coefficient μ1 is only changed to μ1=0.1. Figures 7(a) and 7(b) depict the nonlinear motion of the laminated composite piezoelectric rectangular plate on the planes (x1,μ2) and (x3,μ2), respectively. In contrast to Figures 6 and 7, it can be observed that there is a great difference between them. The above numerical analysis proves that the damping coefficients μ1 and μ2 make a great impact on the emergence of the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate.The bifurcation diagram is obtained for the damping coefficientsμ2= 0.06~0.6, f2=55.1, and μ1=0.1; (a) the bifurcation diagram on the plane (x1,μ2); (b) the bifurcation diagram on the plane (x3,μ2).
(a)
(b)Figure8 portrays the bifurcation diagram for the laminated composite piezoelectric rectangular plate when the in-plane excitation along the x-direction a2 varies in the interval a2= 2~70. Other parameters and initial conditions remain the same as those in Figure 7 when the damping coefficient μ2 is selected as μ2=0.1. Figures 8(a) and 8(b) display the bifurcation diagram on the planes (x1,a2) and (x3,a2), respectively. Figure 8 presents that periodic motion is just beginning to emerge in the system, followed by a chaotic motion.The bifurcation diagram is obtained for the in-plane excitationsa2= 2~70, f2=55.1, μ1=0.1, and μ2=0.1; (a) the bifurcation diagram on the plane (x1,a2); (b) the bifurcation diagram on the plane (x3,a2).
(a)
(b)Figure9 denotes the bifurcation diagram for the laminated composite piezoelectric rectangular plate when the in-plane excitation along the y-direction a3 varies in the interval a3= 2~70. Other parameters and initial conditions are the same as those in Figure 8 when the in-plane excitation a2 is chosen as a2=11.59. Figures 9(a) and 9(b) denote the bifurcation diagram on the planes (x1,a3) and (x3,a3), respectively. Comparing with Figures 8 and 9, it is found that the evolution law of the nonlinear motion is very similar to each other. This shows that in-plane excitations a2 and a3 have the same effect on the existence of the multipulse chaotic motions in the laminated composite piezoelectric rectangular plate under this set of parameters.The bifurcation diagram is obtained for the in-plane excitationsa3= 2~70, f2=55.1, μ1=0.1, μ2=0.1, and a2=11.59; (a) the bifurcation diagram on the plane (x1,a3); (b) the bifurcation diagram on the plane (x3,a3).
(a)
(b)Figure10 indicates the bifurcation diagram for the laminated composite piezoelectric rectangular plate when the piezoelectric excitation a4 varies from a4=2 to a4=120. Other parameters and initial conditions remain the same as those in Figure 9 when the in-plane excitation a3 is chosen as a3=21.5. Figures 10(a) and 10(b) demonstrate the bifurcation diagram on the planes (x1,a4) and (x3,a4), respectively. It is observed from Figure 10 that the piezoelectric excitation a4 has a significant influence on the complicated nonlinear dynamic behaviors of the laminated composite piezoelectric rectangular plate. As the piezoelectric excitation a4 increases, Figure 10 reveals the following law: chaotic motion → periodic motion → chaotic motion.The bifurcation diagram is obtained for the piezoelectric excitationsa4= 2~120, f2=55.1, μ1=0.1, μ2=0.1, a2=11.59, and a3=21.5; (a) the bifurcation diagram on the plane (x1,a4); (b) the bifurcation diagram on the plane (x3,a4).
(a)
(b)Based on the above bifurcation diagram, the excitationsf2, a2, a3, and a4 and the damping coefficients μ1 and μ2 are selected as specific values in order to find the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate. Figure 11 indicates the existence of the multipulse chaotic motion of the laminated composite piezoelectric rectangular plate when the excitation f2 is f2=55.1. In this case, the chosen parameters and initial conditions are the same as those in Figure 5. Figures 11(a) and 11(b) represent the phase portraits on the planes (x1,x2) and (x3,x4), respectively. Figures 11(c) and 11(d) give the waveforms on the planes (t,x1) and (t,x3), respectively. Figures 11(e) and 11(f) are the three-dimensional phase portrait in the space (x1,x2,x3) and the Poincare map on the plane (x1,x2), respectively. Figure 11 shows that the excitation f2 has a noticeable effect on the existence of the multipulse chaotic motions on the laminated composite piezoelectric rectangular plate.The multipulse chaotic motion is obtained whenf2=55.1 and μ1=0.1; (a) the phase portrait on the  plane (x1,x2); (b) the phase portrait on the plane (x3,x4); (c) the waveform on the plane (t,x1); (d) the waveform on the plane (t,x3); (e) the phase portraits in the three-dimensional space (x1,x2,x3); (f) Poincare map on the plane (x1,x2).
(a)
(b)
(c)
(d)
(e)
(f)Besides the excitationsf2, a2, a3, and a4 and the damping coefficients μ1 and μ2, the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate also depend on other parameters. Figure 12 is obtained when the parameters and initial conditions are chosen as σ1=7.31, σ2=12.23, μ1=0.3, μ2=0.2, a2=20.9, a3=13.6, a4=17.3, a6=-7.37, f2=35.1, a7=13.39, b6=-3.68, b7=-11.16, x10=-3.1385, x20=-4.45, x30=-0.59, and x40=0.36. Comparing with Figures 11 and 12, it is found that there are differences in the phase portraits, the waveforms, and the Poincare map, respectively. From the three-dimensional phase portrait in Figure 12, we can see that there exists obvious multi- pulse jumping phenomenon.The multipulse chaotic motion is obtained whenσ1=7.31, σ2=12.23, μ1=0.3, μ2=0.2, a2=20.9, a3=13.6, a4=17.3, a6=-7.37, a7=13.39, f2=35.1, b6=-3.68, b7=-11.16, x10=-3.1385, x20=-4.45, x30=-0.59, and x40=0.36; (a) the phase portrait on the plane (x1,x2); (b) the phase portrait on the plane (x3,x4); (c) the waveform on the plane (t,x1); (d) the waveform on the  plane (t,x3); (e) the phase portraits in the three-dimensional space (x1,x2,x3); (f) Poincare map on the plane (x1,x2).
(a)
(b)
(c)
(d)
(e)
(f)In the following numerical simulations, different parameters are given in order to investigate the different shapes of the multipulse chaotic motion. Figure13 demonstrates the multipulse chaotic response in the laminated composite piezoelectric rectangular plate for a6=-17.37 and b6=-13.68. Other parameters and initial conditions are the same as those in Figure 12.  From Figure 13, we can see that there is another shape for the multipulse chaotic motion. It is found that the shapes of these two phenomena depicted in Figures 12 and 13 are completely different.The multipulse chaotic motion is obtained whena6=-17.37 and b6=-13.68; (a) the phase portrait on the plane (x1,x2); (b) the phase portrait on the plane (x3,x4); (c) the waveform on the plane (t,x1); (d) the waveform on the plane (t,x3); (e) the phase portraits in the three-dimensional space (x1,x2,x3); (f) Poincare map on the plane (x1,x2).
(a)
(b)
(c)
(d)
(e)
(f)


## 10. Conclusions

In this paper, the nonlinear vibrations of the laminated composite piezoelectric rectangular plate are studied by applying the theories of the global bifurcations and chaotic dynamics for high-dimensional nonlinear systems. The multipulse heteroclinic orbits and chaotic dynamics are investigated using the energy-phase method for the case of primary parametric resonance and 1 : 2 internal resonances. The energy-phase method can be applied to discover the existence of the Shilnikov type multipulse heteroclinic bifurcations and chaotic dynamics of high-dimensional nonlinear systems in engineering applications. The analysis of the multipulse heteroclinic orbits in the laminated composite piezoelectric rectangular plate demonstrates that the unperturbed integrable system contains a normally hyperbolic invariant two-dimensional manifold near the resonance. The resonance brings about a circle of fixed points on the two-dimensional manifold. Two different points on the circle of fixed points are connected by a heteroclinic orbit. The dissipative perturbation can severely change the dynamics near the circle of fixed points. The analysis shows that the existence of heteroclinic orbits depends only on an energy-phase criterion. The studies have led to the following conclusions.(1)
From the aforementioned analytical studies, it is found that the multipulse Shilnikov type orbits depend on the dissipative perturbations and excitations. As the trajectory of motion approaches to the sink pointPε, every Shilnikov type orbit takes off again and repeats the similar motion in full four-dimensional phase space, which conduces to the Shilnikov type multipulse orbits. In the Shilnikov type multipulse orbits, the final pulse orbit is similar to the one of the single-pulse. It can be conjectured that the transfer of the energy between the different two modes occurs through the Shilnikov type multipulse orbits.(2)
In order to verify the theoretical predictions, numerical simulation is used to examine the bifurcations and chaotic motions of the laminated composite piezoelectric rectangular plate. Several types of the bifurcation diagrams are obtained when the transverse excitationf2, the in-plane excitations a2 and a3, the piezoelectric excitation a4, and the damping coefficients μ1 and μ2 are chosen as several different kinds of control parameters. Based on the bifurcation diagrams, the nonlinear dynamic behavior of the laminated composite piezoelectric rectangular plate is controlled from the chaotic motion to the periodic motion by varying the excitations f2, a2, a3, and a4 and the damping coefficients μ1 and μ2, respectively. Therefore, the excitations f2, a2, a3, a4, and the damping coefficients μ1 and μ2 have important influence on the nonlinear dynamics responses of the laminated composite piezoelectric rectangular plate.(3)
There exist different shapes of the chaotic motions in the nonlinear oscillations of the laminated composite piezoelectric rectangular plate under different excitations, parameters, and initial conditions. It is found from the numerical simulations that the shapes of the chaotic motions are completely different. Therefore, the parameters and the initial conditions have impact on the shapes of the multipulse chaotic motions.(4)
There exist multipulse chaotic motions in the averaged equations. It is well known that the multipulse chaotic motions in the averaged equations can lead to the multipulse amplitude modulated chaotic vibrations in the original system under certain conditions. Therefore, the multipulse amplitude modulated chaotic motions occur in the laminated composite piezoelectric rectangular plate.In a word, we draw a conclusion through theoretical and numerical investigations that the chaos for the Smale horseshoe sense in nonlinear motion of the simply supported laminated composite piezoelectric rectangular plate exists.



---

*Source: 958219-2013-08-01.xml*