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

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 themethod 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 typemultipulse 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.


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 highdimensional 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 -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-offreedom 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 -degree-of-freedom, nearly integrable Hamiltonian systems.Haller [28,29] developed the energyphase 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 energyphase 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, inplane 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 Figure 1: The model of a laminated composite piezoelectric rectangular plate is given.
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.

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 by a, 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  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 -direction of the laminate.The fiber direction of even-numbered layers is the -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  is free to move at the edge of  = 0, and the displacement of  is free to move at the edge of  = 0.
Therefore, the membrane stress is smaller and there exists the coupling between bending and extension.A Cartesian coordinate system  is located on the middle surface of the composite laminated piezoelectric rectangular plate.Assume that (, V, ) and ( 0 , V 0 ,  0 ) describe the displacements of an arbitrary point and a point on the middle surface of the composite laminated piezoelectric rectangular plate in the , , and  directions, respectively.It is also assumed that in-plane excitations of the composite laminated piezoelectric rectangular plate are loaded along the -direction at  = 0 and the -direction at  = 0 with the form of  0 +   cos Ω 1  and  1 +  cos Ω 2 , respectively.Transverse excitation loaded to the composite laminated piezoelectric rectangular plate is expressed as  =  3 cos Ω 3 .The dynamic electrical loading is represented by   =   cos Ω 4 .
In this paper, Reddy's third-order shear deformation description of the displacement field is adopted: V (, , , ) = V 0 (, , ) +   (, , ) (, , , ) =  0 (, , ) , where ( 0 , V 0 ,  0 ) are the deflection of a point on the middle surface, (, V, ) are the displacement components along the (, , ) coordinate directions, and   and   denote the rotation components of normal to the middle surface about the  and  axes, respectively.The nonlinear strain-displacement relations are assumed to have the following form: ( Stress constitutive relations are presented as follows: where   and   denote the mechanical stresses and strains in extended vector notation,    represents the elastic stiffness tensor,   stands for the electric field vector, and   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]: where the dot represents the partial differentiation with respect to time , 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 In addition, the stress resultants are represented as follows:  4 ) , (,  = 4, 5) .
In order to obtain the dimensionless governing equations of motion, we introduce the transformations of the variables and parameters: 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  0 and V 0 directions, respectively.Therefore, the equations for  0 and V 0 can be neglected.We mainly consider the nonlinear transverse vibration of laminated composite rectangular plates.
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: 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: where  0 = ,  1 = .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 fourdimensional averaged equations in the Cartesian form are obtained as follows: (16d)

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  in the subspace ( 1 ,  2 ) of (32a), (32b), (32c), and (32d) becomes a parameter since İ = 0. Consider the first two decoupled equations of (32a), (32b), (32c), and (32d), Since  1 > 0, (35a) and (35b) can exhibit the heteroclinic bifurcations.It is obvious from (35a) and (35b) that when  − a 6 I 2 < 0, the only solution to (35a) and (35b) is the trivial zero solution, ( 1 ,  2 ) = (0, 0), which is the saddle point.On the curve defined by  =  6  2 , that is, or the trivial zero solution bifurcates into three solutions through a pitchfork bifurcation, which are given by  0 = (0, 0) and  ± () = (, 0), respectively, where From the Jacobian matrix evaluated at the nonzero solutions, it can be found that the singular points  ± () are the saddle points.It is observed that  and  actually represent the amplitude and the phase of vibrations, respectively.Therefore, we assume that  ≥ 0 and (37) become such that for all  ∈ [ 1 ,  2 ], (35a) and (35b) have two hyperbolic saddle points,  ± (), which are connected by a pair of heteroclinic orbits,  ℎ ± ( 1 , ), that is, lim  1 → ±∞  ℎ ± ( 1 , ) =  ± ().Thus, in the full four-dimensional phase space, the set defined by is a two-dimensional invariant manifold.From the results obtained by Kovačič [13,14], it is known that the twodimensional invariant manifold  is normally hyperbolic.The two-dimensional normally hyperbolic invariant manifold  has the three-dimensional stable and unstable manifolds represented as   () and   (), respectively.The existence of the heteroclinic orbit of (35a) and (35b) to  ± () = (, 0) indicates that   () and   () intersect nontransversally along a three-dimensional heteroclinic manifold denoted by Γ, which can be written as We analyze the dynamics of the unperturbed system of (32a), (32b), (32c), and (32d) restricted to .Considering the unperturbed system of (32a), (32b), (32c), and (32d) restricted to  yields İ = 0, (42a) where From the results obtained by Kovačič [13,14], it is known that if   ( ± (), ) ̸ = 0; then  = constant is called a periodic orbit, and if   ( ± (), ) = 0; then  = constant is known as a circle of the singular points.Any value of  ∈ [ 1 ,  2 ] at which   ( ± (), ) = 0 is a resonant value  and these singular points are resonant singular points.We denoted a resonant value by   such that Then, we obtain Figure 2 shows the geometry structure of the stable and unstable manifolds of  in the full four-dimensional phase space for the unperturbed system of (32a), (32b), (32c), and (32d).Since  represents the phase of oscillations, when  =   , the phase shift Δ of oscillations is defined by The physical interpretation of the phase shift is the phase difference between the two end points of the orbit.In the subspace ( 1 ,  2 ), there exists a pair of the heteroclinic orbits connecting to the two saddles.Therefore, the homoclinic orbit in subspace (, ) represents, in fact, a heteroclinic connecting in the full four-dimensional space ( 1 ,  2 , , ).The phase shift denotes the difference of the value  as a trajectory leaving and returning to the basin of the attraction of .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 =  −  6  2 and let  3 =  2 ; then (35a) and (35b) can be rewritten as Set  = 0; then (47a) and (47b) are a system with the Hamiltonian function: When  =  2 1 /(4 1 ), there exists a heteroclinic loop Γ 0 which consists of the two hyperbolic saddles  ± and a pair of heteroclinic orbits  ± ( 1 ).In order to calculate the phase shift and the energy difference function, we obtain the equations of a pair of heteroclinic orbits given by 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 Integrating (50) yields where At  =   , there is   ≡ 0. Therefore, the phase shift may be expressed as

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 on .Based on the analysis given by Kovačič [13,14], it is known that  along with its stable and unstable manifolds is invariant under small, sufficiently differentiable perturbations.It is noticed that  ± () in (35a) and (35b) maintain the characteristic of the hyperbolic singular point under small perturbations, in particular,  →   .Therefore, we obtain Considering the last two equations of (32a), (32b), (32c), and (32d) yields 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 Substituting the above transformations into (54a) and (54b) yields where  =  2 6 −  2  1 .When  = 0, (56a) and (56b) become The unperturbed system from (57a) and (57b) is a Hamilton system with the following function: The singular points of (57a) and (57b) are represented as Based on the characteristic equations evaluated at the two singular points  0 and  0 , we can know the stabilities of these singular points.The Jacobian matrix of (57a) and (57b) is The characteristic equation corresponding to the singular point  0 is obtained as When the condition (2/ 1 ) 2   cos   < 0 is satisfied, (57a) and (57b) have a pair of pure imaginary eigenvalues.Therefore, it is known that the singular point  0 is a center point.
The characteristic equation corresponding to the singular point  0 is represented by When the condition (2/ 1 ) 2   cos   > 0 is satisfied, (57a) and (57b) have two real, unequal, and opposite sign eigenvalues.Therefore, the singular point  0 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).It is found that for the sufficiently small parameter , the singular point  0 remains a hyperbolic singular point   of the saddle stability type.It is known that the Jacobian matrix of the linearization of (56a) and (56b) is of the following form: or 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  0 becomes a hyperbolic sink   .The phase portrait of the perturbed system from (56a) and (56b) is also depicted in Figure 3(b).
Using the function (58), at ℎ = 0, the estimate of the basin of attractor for  min is obtained as Substituting   in (59) into (65) yields Define an annulus   near  =   as where  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 of   , denoted as   (  ) and   (  ), are subsets of   (  ) and   (  ), respectively.We will indicate that for the perturbed system, the saddle focus   on   has the multipulse orbits which come out of the annulus   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   , leave and enter a small neighborhood of   multiple times, and, finally, return and approach an invariant set in   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.
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][22][23][24][25][26][27][28][29][30][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][22][23][24][25][26][27][28][29][30][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][22][23][24][25][26][27][28][29][30][31], we use the geometric singular perturbation theory of Fenichel to analyze the existence of heteroclinic orbits in the full fourdimensional phase space.Based on foliations of stable and unstable manifolds, Fenichel's theory combines with energytype 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 higherdimensional 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][22][23][24][25][26][27][28][29][30][31] for the class of systems, the energy difference function for the dissipative case is given as follows: where where  denotes the area enclosed by a pair of heteroclinic orbits in plane ( 1 ,  2 ),   is the boundary of the area , and Δ is the phase difference.The coordinate transformation introduced by Haller and Wiggins [21][22][23][24][25][26][27][28][29][30][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 an -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: We define a dissipative factor  =  3 / 2 such that  gives the relative measure of the dissipative effect with respect to the excitation amplitude.Equation (70) can be rewritten as The zeroes of Δ  Ĥ () are obtained by solving the following equation: Based on (72), the upper bound on the value of the dissipative factor is obtained as follows: For a given value of the dissipative factor , the multipulse orbits for all values of the pulse number  are not possible.
In the presence of small dissipative effects  < 1, we have an upper bound on the maximum number of pulses: It is obviously known from (74) that the upper bound  max is inversely proportional to the value of the dissipative factor .

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: The transverse zeroes of the dissipative energy difference function Δ  Ĥ () are given by the following equation: where  ∈ Ζ and For any  satisfying Δ ̸ = 4 ( = 0, 1, 2, . ..), there are two transverse zeroes of the dissipative energy difference function in the interval  ∈ [−/2, 3/2], that is: 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 is ℎ 0 .The energy level of the center type singular point is represented by ℎ ∞ .Any periodic orbit enclosed inside the homoclinic connections has an energy level ℎ  .The energy sequence at the different energy levels is defined as where  2 / 1 = ,  2 / 2 = (1/2).
From (79a), (79b), and (79c), it is found that there exists the relation ℎ 0 < ℎ  < ℎ ∞ 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 manifold   , implies the existence of the chaos for the Smale horseshoe sense in the simply supported laminated composite piezoelectric rectangular plate [24,25].

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 manifold   , was studied by Haller and Wiggins [21][22][23][24][25][26][27][28][29][30][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 -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][22][23][24][25][26][27][28][29][30][31] to the following research.
First, we require the existence of nondegenerate equilibrium points for Ĥ .In the (ℎ, ) phase space, based on (59) and  2 / 1 = , it is known that this singular point is given by Next, we need to compute the zeroes of the dissipative energy difference function at the saddle center (, 0, 0,   ) in full four-dimensional phase space ( 1 ,  2 , ℎ, ).In order to obtain the zeroes of Δ  Ĥ (  ), we solve the following equation: Equation (81) leads to Based on (82), we obtain the following dissipative parameter: The results obtained above are only valid when the nonzero damping exists, that is, 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: whenever ( 83) and ( 84) hold.The first condition in (85) can be rewritten as Carrying out the differentiation with respect to  in (86), it is seen that the above condition is violated only if 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 Δ  Ĥ .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 any -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 2 apart from the approximate landing point   + Δ.This point is denoted by where If   * >   , then we redefine   * by subtracting 2 from it.The aim here is to find 2-translation of the landing point which is closest to the saddle point   in the region [−/2, 3/2].Due to the symmetry of the phase portrait  →  − 2, if the energy of point   * is more than the energy of the saddle point, that is, then, for  > 0 small enough, the landing point   + Δ will be in the domain of attraction of one of the sinks.Calculating this energy condition yields The shape of such an orbit for  = 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].
(1) If the integer is even, then each of the saddle-focus type equilibrium points contained in the slow manifold   admits two generalized Shilnikov orbits.If the integer  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 ( 1 ,  2 ) = (, 0).
We give an explanation of statement ( 1) in Theorem 1.
The main goal is to prove that utilizing an integer  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 The landing point   + Δ will fall in one of these two regions or a 2-translation of one of these regions.If the landing point falls in  1 , 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  2 , 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 ( 1 ,  2 ) = (, 0), due to the existence of the symmetry ( 1 ,  2 ) → (− 1 , − 2 ).
If the landing point   + Δ is in  1 or a 2-translation, one has Equation (94) leads to Hence, the integer is even.Reversely, if the landing point falls in  2 or a 2translation, one can similarly obtain From (97), one has Therefore, it is known that the integer is odd.
The proof of the second part of Theorem 1 follows from Theorem 2.8.2 presented by Haller [31].

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  6 ,  7 , and  6 , where  1 = (3 6  7 )/ 6 .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  6 ,  7 , and  6 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  1 and  2 and the transverse excitation  2 play an important role in the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate.In addition, the parameters  2 and  3 are related to the in-plane excitation in the -direction and the in-plane excitation in the -direction, respectively.The parameter  4 is the piezoelectric excitation which reflects the characteristics of the piezoelectric material.Hence, the parameters  1 ,  2 ,  2 ,  3 ,  4 , and  2 are selected as the controlling parameters to discover the law for complicated nonlinear dynamics of the laminated composite piezoelectric rectangular plate.5 that the excitation  2 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  2 varies in the interval  2 = 20∼78, followed by a periodic motion of that.When excitation  2 continues to increase, Figure 5 presents that there is the periodic motion of the laminated composite piezoelectric rectangular plate.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  2 = 55.1.Figures 6(a) and 6(b) describe the nonlinear motion of the laminated composite piezoelectric rectangular plate on the planes ( 1 ,  1 ) and ( 3 ,  1 ), respectively, as well as the impact of the damping coefficient  1 on the system.
Figure 7 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 ( 1 ,  2 ) and ( 3 ,  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.Figure 8 portrays the bifurcation diagram for the laminated composite piezoelectric rectangular plate when the inplane excitation along the -direction  2 varies in the interval  2 = 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 ( 1 ,  2 ) and ( 3 ,  2 ), respectively.Figure 8 presents that periodic motion is just beginning to emerge in the system, followed by a chaotic motion.
Figure 9 denotes the bifurcation diagram for the laminated composite piezoelectric rectangular plate when the in-plane excitation along the -direction  3 varies in the interval  3 = 2∼70.Other parameters and initial conditions are the same as those in Figure 8 when the in-plane excitation  2 is chosen as  2 = 11.59.Figures 9(a) and 9(b) denote the bifurcation diagram on the planes ( 1 ,  3 ) and ( 3 ,  3 ), 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  2 and  3 have the same effect on the existence of the multipulse chaotic motions in the laminated composite piezoelectric rectangular plate under this set of parameters.
Figure 10 indicates the bifurcation diagram for the laminated composite piezoelectric rectangular plate when the piezoelectric excitation  4 varies from  4 = 2 to  4 = 120.Other parameters and initial conditions remain the same as those in Figure 9 when the in-plane excitation  3 is chosen as 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  2 is  2 = 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 ( 1 ,  2 ) and ( 3 ,  4 ), respectively.Figures 11(c) and 11(d) give the waveforms on the planes (,  1 ) and (,  3 ), respectively.Figures 11(e) and 11(f) are the three-dimensional phase portrait in the space ( 1 ,  2 ,  3 ) and the Poincare map on the plane ( 1 ,  2 ), respectively.Figure 11 shows that the excitation  2 has a noticeable effect on the existence of the multipulse chaotic motions on the laminated composite piezoelectric rectangular plate.
Besides the excitations  2 ,  2 ,  3 , and  4 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   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.
In the following numerical simulations, different parameters are given in order to investigate the different shapes of the multipulse chaotic motion.Figure 13 demonstrates the multipulse chaotic response in the laminated composite piezoelectric rectangular plate for  6 = −17.37 and  6 = −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.

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 energyphase 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 twodimensional 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 point   , every Shilnikov type orbit takes off again and repeats the similar motion in full fourdimensional 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.
( (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.

Figure 2 :
Figure 2: The geometric structure of manifolds ;   () and   () are given in the full four-dimensional phase space.

Figure 3 :
Figure 3: Dynamics on the normally hyperbolic manifold is described; (a) the unperturbed case; (b) the perturbed case.