# Modeling and Chaotic Dynamics of the Laminated Composite Piezoelectric Rectangular Plate

Authors: Minghui Yao; Wei Zhang; D. M. Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper investigates the multipulse heteroclinic bifurcations and chaotic dynamics of a laminated composite piezoelectric rectangular plate by using an extended Melnikov method in the resonant case. According to 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. The method of multiple scales and Galerkin’s approach are applied to the partial differential governing equation. Then, the four-dimensional averaged equation is obtained for the case of 1 : 3 internal resonance and primary parametric resonance. The extended Melnikov method is used to study the Shilnikov type multipulse heteroclinic bifurcations and chaotic dynamics of the laminated composite piezoelectric rectangular plate. The necessary conditions of the existence for the Shilnikov type multipulse chaotic dynamics are analytically obtained. From the investigation, the geometric structure of the multipulse orbits is described in the four-dimensional phase space. Numerical simulations show that the Shilnikov type multipulse chaotic motions can occur. To sum up, both theoretical and numerical studies suggest that chaos for the Smale horseshoe sense in motion exists for the laminated composite piezoelectric rectangular plate.


---

## Body

## 1. Introduction

The need for high-speed, light-weight, and energy-saving structures in the aerospace and aviation industry has led to the composite materials instead of traditional materials. Additional requirements for multifunctionality, active vibration, shape control, vibration suppression, and acoustic control have made the development of smart and intelligent structures. A piezoelectric composite laminate is composed of piezoelectric layers which are embedded in laminated composite structures or are boned on the surface of structures. The direct and converse piezoelectric effects are used to suppress the transient vibration and to control the deformation, shape, and buckling of the structures. Such lightweight flexible structures generate large deformations, geometrical nonlinearity, and structural instability when piezoelectric composite laminates are subjected to the coupling between the mechanical and electrical loads. Therefore, it is necessary to study geometrically nonlinear effects on dynamic characteristics of structures in order to accurately design and effectively control vibrations of piezoelectric composite laminate structures. It is very important to investigate the large amplitude nonlinear vibrations of smart structures with piezoelectric materials in order to achieve and predict the desired performance of the systems.Recently, the studies on dynamics of composite structures with piezoelectric materials have made some progress. Tzou et al. [1] used spatially distributed orthogonal piezoelectric actuators to perform the distributed structural control of elastic shell. They utilized a gain factor and a spatially distributed mode actuator function to describe modal feedback functions. Purekar et al. [2] presented phased array filters with piezoelectric sensors to detect damage in isotropic plates and adopted wave propagation to describe plate dynamics. Ishihara and Noda [3] took into account the effect of transverse shear to analyze the dynamic behavior of the laminate composed of fiber-reinforced laminae and piezoelectric layers constituting a symmetric cross-ply laminate rectangular plate with simply supported edges. Oh [4] considered snap-through thermopiezo-elastic behaviors to examine the buckling bifurcation and sling-shot buckling of active piezo-laminated plates. Lee et al. [5] employed third-order shear deformation theory and nonlinear finite element to canvass deflection suppression characteristics of laminated composite shell structures with smart material laminae. Panda and Ray [6] exploited the first-order shear deformation theory and the three-dimensional finite element method to delve into the open-loop and closed-loop nonlinear dynamics of functionally graded plates with the piezoelectric fiber-reinforced composite material under the thermal environment. Dumir et al. [7] used the extended Hamilton’s principle to derive the coupled nonlinear equations of motion and the boundary conditions for buckling and vibration of symmetrically laminated hybrid angle-ply piezoelectric panels under in-plane electrothermomechanical loading. Yao and Zhang [8] employed the third-order shear deformation plate theory to explore the bifurcations and chaotic dynamics of the four-edge simply supported laminated composite piezoelectric rectangular plate in the case of the 1 : 2 internal resonances.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 Shilnikov type chaotic dynamics of high-dimensional nonlinear systems. One is Shilnikov type single-pulse chaotic dynamics and the other is Shilnikov type multipulse chaotic dynamics. Most researchers focused on Shilnikov type single-pulse chaotic dynamics of high-dimensional nonlinear systems. Much research in this field has concentrated on Shilnikov type single-pulse chaotic dynamics of thin plate structures. Feng and Sethna [9] utilized the global perturbation method to study the global bifurcations and chaotic dynamics of the thin plate under parametric excitation and obtained the conditions in which the Shilnikov type homoclinic orbits and chaos can occur. Tien et al. [10] applied the Melnikov method to investigate the global bifurcation and chaos for the Smale horseshoe sense of a two-degree-of-freedom shallow arch subjected to simple harmonic excitation for the case of 1 : 2 internal resonance. Malhotra and Sri Namachchivaya [11] employed the averaging method and Melnikov technique to canvass the local, global bifurcations and chaotic motions of a two-degree-of-freedom shallow arch subjected to simple harmonic excitation for the case of 1 : 1 internal resonance. The global bifurcations and chaotic dynamics were investigated by Zhang [12] for the simply supported rectangular thin plates subjected to the parametrical-external excitation and the parametrical excitation. Yeo and Lee [13] made use of the global perturbation technique to examine the global dynamics of an imperfect circular plate for the case of 1 : 1 internal resonance and obtained the criteria for chaotic motions of homoclinic orbits and heteroclinic orbits. Yu and Chen [14] adopted the global perturbation method to explore the global bifurcations of a simply supported rectangular metallic plate subjected to a transverse harmonic excitation for the case of 1 : 1 internal resonance.While most of studies are on the Shilnikov type single-pulse global bifurcations and chaotic dynamics of high-dimensional nonlinear systems, 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 extended Melnikov method and the other theory is the energy phase method. Much achievement is made in the former theory of high-dimensional nonlinear systems. In 1996, Kovačič and Wettergren [15] used a modified 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č [16] studied the existence of several classes of the multibump orbits homoclinic to resonance bands for completely integral Hamiltonian systems subjected to small amplitude Hamiltonian and damped perturbations. Camassa et al. [17] presented a new Melnikov method which is called the extended Melnikov method to explore the multipulse jumping of homoclinic and heteroclinic orbits in a class of perturbed Hamiltonian systems. Until recently, Zhang and Yao [18] 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.The study on the second theory of the Shilnikov type multipulse chaotic dynamics was stated by Haller and Wiggins [19]. They presented the energy phase method to investigate the existence of the multipulse jumping homoclinic and heteroclinic orbits in perturbed Hamiltonian systems. Up to now, few researchers have made use of the energy phase method to study the Shilnikov type multipulse homoclinic and heteroclinic bifurcations and chaotic dynamics of high-dimensional nonlinear systems in engineering applications. Malhotra et al. [20] used the energy-phase method to investigate multipulse homoclinic orbits and chaotic dynamics for the motion of flexible spinning discs. Yu and Chen [21] made use of the energy-phase method to examine the Shilnikov type multipulse homoclinic orbits of a harmonically excited circular plate.This paper focuses on the Shilnikov type multipulse orbits and chaotic dynamics for a simply supported laminated composite piezoelectric rectangular plate under combined parametric excitations and transverse load. Based on the von Karman type equations and Reddy’s third-order shear deformation plate theory, 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 Galerkin’s approach and the method of multiple scales to the partial differential governing equations to obtain the four-dimensional averaged equation for the case of 1 : 3 internal resonance and primary parametric resonance. From the averaged equation, the theory of normal form is used to find the explicit formulas of normal form. We study the heteroclinic bifurcations of the unperturbed system and the characteristic of the hyperbolic dynamics of the dissipative system, respectively. Finally, we employ the extended Melnikov method to analyze the Shilnikov type multipulse orbits and chaotic dynamics in the laminated composite piezoelectric plate. In this paper, the extended Melnikov function can be simplified in the resonant case and does not depend on the perturbation parameter. We have used the extended Melnikov method to investigate heteroclinic bifurcations and multipulse chaotic dynamics of the laminated composite piezoelectric plate under the case of 1 : 3 internal resonances. The analysis indicates that there exist the Shilnikov type multipulse jumping orbits in the perturbed phase space for the averaged equations. We present the geometric structure of the multipulse orbits in the four-dimensional phase space. 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. This paper demonstrates how to employ the extended Melnikov method to analyze the Shilnikov type multipulse heteroclinic bifurcations and chaotic dynamics of high-dimensional nonlinear systems in engineering applications.The laminated composite piezoelectric rectangular plates are widely applied in space stations, satellite solar panels, sensors, and actuators for the active control of structures and so on. In this paper, we have investigated the multipulse global bifurcations and chaotic dynamics of a laminated composite piezoelectric rectangular plate by using an extended Melnikov method and numerical simulations in detail. We have understood nonlinear vibration characteristics of a laminated composite piezoelectric rectangular plate. Our theoretical results can be used to solve some engineering problems. Since these smart structures are generally light weight and relatively large structural flexibility, laminated composite piezoelectric rectangular plates can induce large vibration deformation during the rapid deployment. In order to eliminate or suppress large vibration and chaotic motion, theoretical results can help optimize the design of the structural parameters of laminated composite piezoelectric rectangular plates. Therefore, the theoretical studies on the multipulse global bifurcations and chaotic dynamics of laminated composite piezoelectric rectangular plates play a very important role in applications in aerospace and mechanical engineering.


## 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 laminated coordinate axes. Some of 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 O
x
y
z is located in the middle surface of the composite laminated piezoelectric rectangular plate. Assume that (
w
,
v
,
u
) and (
w
0
,
v
0
,
u
0
) describe the displacements of an arbitrary point and a point in 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 q
0
+
q
x
cos
⁡
Ω
1
t and q
1
+
q
y
cos
⁡
Ω
2
t, respectively. Transverse excitation loaded to the composite laminated piezoelectric rectangular plate is expressed as q
=
q
3
cos
⁡
Ω
3
t. The dynamic electrical loading is represented by E
z
=
E
z
cos
⁡
Ω
4
t.Figure 1
The model of a laminated composite piezoelectric rectangular plate is given.Considering Reddy’s third-order shear deformation description of the displacement field, we have(1a)
u
(
x
,
y
,
z
,
t
)
=
u
0
(
x
,
y
,
t
)
+
z
ϕ
x
(
x
,
y
,
t
)
-
z
3
4
3
h
2
(
ϕ
x
+
∂
w
0
∂
x
)
,
(1b)
v
(
x
,
y
,
z
,
t
)
=
v
0
(
x
,
y
,
t
)
+
z
ϕ
y
(
x
,
y
,
t
)
-
z
3
4
3
h
2
(
ϕ
y
+
∂
w
0
∂
y
)
,
(1c)
w
(
x
,
y
,
z
,
t
)
=
w
0
(
x
,
y
,
t
)
,where (u
0
,
v
0
,
w
0) 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 represent 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)
ε
x
x
=
∂
u
∂
x
+
1
2
(
∂
w
∂
x
)
2
,
ε
x
z
=
1
2
(
∂
u
∂
z
+
∂
w
∂
x
)
,
ε
x
y
=
1
2
(
∂
u
∂
y
+
∂
v
∂
x
+
∂
w
∂
x
∂
w
∂
y
)
ε
y
y
=
∂
v
∂
y
+
1
2
(
∂
w
∂
y
)
2
,
ε
y
z
=
1
2
(
∂
v
∂
z
+
∂
w
∂
y
)
,
ε
z
z
=
∂
w
∂
z
.Stress constitutive relations are presented as follows:(3)
σ
i
j
=
σ
i
j
k
l
s
ε
k
l
-
e
i
j
k
E
k
,
(
i
,
j
,
k
,
l
=
x
,
y
,
z
)
,
where σ
i
j and ε
k
l denote the mechanical stresses and strains in extended vector notation, σ
i
j
k
l
srepresents the elastic stiffness tensor, E
k stands for the electric field vector, and e
i
j is the piezoelectric tensor.According to Hamilton’s principle, the nonlinear governing equations of motion in terms of generalized displacements(
u
0
,
v
0
,
w
0
,
ϕ
x
,
ϕ
y
) for the composite laminated piezoelectric rectangular plate are given in the previous studies as follows [8]:
(4a)
A
11
∂
2
u
0
∂
x
2
+
A
66
∂
2
u
0
∂
y
2
+
(
A
12
+
A
66
)
∂
2
v
0
∂
x
∂
y
+
A
11
∂
w
0
∂
x
∂
2
w
0
∂
x
2
+
A
66
∂
w
0
∂
x
∂
2
w
0
∂
y
2
+
(
A
12
+
A
66
)
∂
w
0
∂
y
∂
2
w
0
∂
x
∂
y
=
I
0
u
¨
0
+
J
1
ϕ
¨
x
-
c
1
I
3
∂
w
¨
0
∂
x
,
(4b)
A
66
∂
2
v
0
∂
x
2
+
A
22
∂
2
v
0
∂
y
2
+
(
A
21
+
A
66
)
∂
2
u
0
∂
x
∂
y
+
A
66
∂
w
0
∂
y
∂
2
w
0
∂
x
2
+
A
22
∂
w
0
∂
y
∂
2
w
0
∂
y
2
+
(
A
21
+
A
66
)
∂
w
0
∂
x
∂
2
w
0
∂
x
∂
y
=
I
0
v
¨
0
+
J
1
ϕ
¨
y
-
c
1
I
3
∂
w
¨
0
∂
y
,
(4c)
A
66
∂
w
0
∂
x
∂
2
u
0
∂
y
2
-
H
22
c
1
2
∂
4
w
0
∂
y
4
+
c
1
(
2
F
66
+
F
12
-
2
H
66
c
1
-
H
12
c
1
)
∂
3
ϕ
y
∂
y
∂
x
2
+
c
1
(
F
22
-
H
22
c
1
)
∂
3
ϕ
y
∂
y
3
-
H
11
c
1
2
∂
4
w
0
∂
x
4
+
A
11
∂
w
0
∂
x
∂
2
u
0
∂
x
2
+
(
F
44
c
2
2
-
2
D
44
c
2
+
A
44
)
∂
ϕ
y
∂
y
+
c
1
(
F
21
+
2
F
66
-
H
21
c
1
-
2
H
66
c
1
)
∂
3
ϕ
x
∂
y
2
∂
x
-
c
1
2
(
H
21
+
4
H
66
+
H
12
)
∂
4
w
0
∂
y
2
∂
x
2
+
(
A
44
-
N
y
P
cos
⁡
⁡
(
Ω
4
t
)
+
F
44
c
2
2
-
2
D
44
c
2
)
∂
2
w
0
∂
y
2
-
∂
N
y
P
∂
y
cos
⁡
⁡
(
Ω
2
t
)
∂
w
0
∂
y
+
(
A
21
+
4
A
66
+
A
12
)
×
∂
w
0
∂
x
∂
w
0
∂
y
∂
2
w
0
∂
y
∂
x
+
c
1
(
F
11
-
H
11
c
1
)
∂
3
ϕ
x
∂
x
3
+
(
A
21
+
A
66
)
∂
w
0
∂
y
∂
2
u
0
∂
y
∂
x
+
A
21
∂
u
0
∂
x
∂
2
w
0
∂
y
2
×
A
66
∂
w
0
∂
y
∂
2
v
0
∂
x
2
+
A
22
∂
w
0
∂
y
∂
2
v
0
∂
y
2
+
1
2
(
A
12
+
2
A
66
)
(
∂
w
0
∂
y
)
2
∂
2
w
0
∂
x
2
+
A
22
∂
2
w
0
∂
y
2
∂
v
0
∂
y
+
(
A
12
+
A
66
)
∂
w
0
∂
x
∂
2
v
0
∂
y
∂
x
+
1
2
(
A
21
+
2
A
66
)
∂
2
w
0
∂
y
2
×
(
∂
w
0
∂
x
)
2
+
3
2
A
11
(
∂
w
0
∂
x
)
2
∂
2
w
0
∂
x
2
+
A
11
∂
2
w
0
∂
x
2
∂
u
0
∂
x
+
A
12
∂
2
w
0
∂
x
2
∂
v
0
∂
y
+
2
A
66
∂
2
w
0
∂
y
∂
x
∂
v
0
∂
x
+
2
A
66
∂
2
w
0
∂
y
∂
x
∂
u
0
∂
y
+
3
2
A
22
(
∂
w
0
∂
y
)
2
∂
2
w
0
∂
y
2
+
(
A
55
+
q
x
cos
⁡
⁡
(
Ω
1
t
)
-
N
x
P
cos
⁡
⁡
(
Ω
4
t
)
+
F
55
c
2
2
-
2
D
55
c
2
)
∂
2
w
0
∂
x
2
-
∂
N
x
P
∂
x
cos
⁡
⁡
(
Ω
3
t
)
∂
w
∂
x
+
(
F
55
c
2
2
-
2
D
55
c
2
+
A
55
)
∂
ϕ
x
∂
x
-
q
cos
⁡
(
Ω
3
t
)
+
k
f
∂
w
0
∂
t
=
I
0
w
¨
0
-
c
1
2
I
6
(
∂
2
w
¨
0
∂
x
2
+
∂
2
w
¨
0
∂
y
2
)
+
c
1
I
3
(
∂
u
¨
0
∂
x
+
∂
v
¨
0
∂
y
)
+
c
1
J
4
(
∂
ϕ
¨
x
∂
x
+
∂
ϕ
¨
y
∂
y
)
,
(4d)
(
D
11
-
2
F
11
c
1
+
H
11
c
1
2
)
∂
2
ϕ
x
∂
x
2
+
(
D
66
-
2
F
66
c
1
+
H
66
c
1
2
)
∂
2
ϕ
x
∂
y
2
-
c
1
(
F
11
-
H
11
c
1
)
∂
3
w
0
∂
x
3
-
(
F
55
c
2
2
-
2
D
55
c
2
+
A
55
)
∂
w
0
∂
x
+
(
D
12
+
D
66
+
H
66
c
1
2
-
2
F
66
c
1
+
H
12
c
1
2
-
2
F
12
c
1
)
×
∂
2
ϕ
y
∂
y
∂
x
-
c
1
(
2
F
66
+
F
12
-
2
H
66
c
1
-
H
12
c
1
)
×
∂
3
w
0
∂
y
2
∂
x
+
(
2
D
55
c
2
-
A
55
-
F
55
c
2
2
)
ϕ
x
=
J
1
u
¨
0
+
K
2
ϕ
¨
x
-
c
1
J
4
∂
w
¨
0
∂
x
,
(4e)
(
D
66
-
2
F
66
c
1
+
H
66
c
1
2
)
∂
2
ϕ
y
∂
x
2
-
c
1
(
F
21
+
2
F
66
-
H
21
c
1
-
2
H
66
c
1
)
∂
3
w
0
∂
y
∂
x
2
+
(
H
21
c
1
2
+
D
66
+
D
21
-
2
F
21
c
1
+
H
66
c
1
2
-
2
F
66
c
1
)
×
∂
2
ϕ
x
∂
y
∂
x
+
(
H
22
c
1
2
+
D
22
-
2
F
22
c
1
)
∂
2
ϕ
y
∂
y
2
-
c
1
(
F
22
-
H
22
c
1
)
∂
3
w
0
∂
y
3
-
(
F
44
c
2
2
-
2
D
44
c
2
+
A
44
)
×
∂
w
0
∂
y
+
(
2
D
44
c
2
-
F
44
c
2
2
-
A
44
)
ϕ
y
=
J
1
v
¨
0
+
K
2
ϕ
¨
y
-
c
1
J
4
∂
w
¨
0
∂
y
.The simply supported boundary conditions of the composite laminated piezoelectric rectangular plate can be represented as follows [8, 22]:
(5a)
x
=
0
:
w
=
0
,
ϕ
y
=
0
,
N
x
y
=
0
,
M
x
x
=
0
,
(5b)
x
=
a
:
w
=
0
,
ϕ
y
=
0
,
N
x
y
=
0
,
M
x
x
=
0
,
(5c)
y
=
0
:
w
=
0
,
ϕ
x
=
0
,
N
x
y
=
0
,
M
y
y
=
0
,
(5d)
y
=
b
:
w
=
0
,
ϕ
x
=
0
,
N
x
y
=
0
,
M
y
y
=
0
,
(5e)
∫
0
h
N
x
x
|
x
=
0
d
z
=
-
∫
0
h
(
q
0
+
q
x
cos
⁡
Ω
1
t
)
d
z
,
(5f)
∫
0
h
N
y
y
|
y
=
0
d
z
=
-
∫
0
h
(
q
1
+
q
y
cos
⁡
Ω
2
t
)
d
z
.The boundary condition (5f) 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 u
0,  v
0,  w
0,  ϕ
x, and ϕ
y. It is desirable that we select an appropriate mode function to satisfy the boundary condition. Thus, we can rewrite u
0,  v
0,  w
0,  ϕ
x, and ϕ
y in the following forms:
(6a)
u
0
=
u
01
(
t
)
cos
⁡
π
x
2
a
cos
⁡
π
y
2
b
+
u
02
(
t
)
cos
⁡
3
π
x
2
a
cos
⁡
π
y
2
b
,
(6b)
v
0
=
v
1
(
t
)
cos
⁡
π
y
2
b
cos
⁡
π
x
2
a
+
v
2
(
t
)
cos
⁡
π
y
2
b
cos
⁡
3
π
x
2
a
,
(6c)
w
0
=
w
1
(
t
)
sin
π
x
a
sin
π
y
b
+
w
2
(
t
)
sin
3
π
x
a
sin
π
y
b
,
(6d)
ϕ
x
=
ϕ
1
(
t
)
cos
⁡
π
x
a
sin
π
y
b
+
ϕ
2
(
t
)
cos
⁡
3
π
x
a
sin
π
y
b
,
(6e)
ϕ
y
=
ϕ
3
(
t
)
cos
⁡
π
y
b
sin
π
x
a
+
ϕ
4
(
t
)
cos
⁡
π
y
b
sin
3
π
x
a
.By means of the Galerkin method, substituting (6a), (6b), (6c), (6d), (6e) into (4a), (4b), (4c), (4d), (4e), integrating, and neglecting all inertia terms in (4a), (4b), (4d), and (4e), we obtain the expressions of u
01,  u
02,  v
1,  v
2,  ϕ
1,  ϕ
2,  ϕ
3, and ϕ
4 via w
1 and w
2 as follows:
(7a)
u
01
=
k
1
w
1
2
+
k
2
w
2
2
+
k
3
w
1
w
2
,
(7b)
u
02
=
k
4
w
1
2
+
k
5
w
2
2
+
k
6
w
1
w
2
,
(7c)
v
1
=
k
7
w
1
2
+
k
8
w
2
2
+
k
9
w
1
w
2
,
(7d)
v
2
=
k
10
w
1
2
+
k
11
w
2
2
+
k
12
w
1
w
2
,
(7e)
ϕ
1
=
k
19
w
1
,
ϕ
2
=
k
20
w
2
,
(7f)
ϕ
3
=
k
21
w
1
,
ϕ
4
=
k
22
w
2
,where the coefficients presented in (7a), (7b), (7c), (7d), (7e), (7f) can be found in the previous studies [8].In order to obtain the dimensionless governing equations of motion, we introduce the transformations of the variables and parameters(8)
u
-
=
u
0
a
,
v
-
=
v
0
b
,
w
-
=
w
0
h
,
ϕ
-
x
=
ϕ
x
,
ϕ
-
y
=
ϕ
y
,
x
-
=
x
a
,
y
-
=
y
b
,
q
-
=
b
2
E
h
3
q
,
q
-
x
=
b
2
E
h
3
q
x
,
q
-
y
=
b
2
E
h
3
q
y
,
t
-
=
π
2
(
E
a
b
ρ
)
1
/
2
t
,
Ω
-
i
=
1
π
2
(
a
b
ρ
E
)
1
/
2
Ω
i
(
i
=
1,2
)
,
A
-
i
j
=
(
a
b
)
1
/
2
E
h
2
A
i
j
,
B
-
i
j
=
(
a
b
)
1
/
2
E
h
3
B
i
j
,
D
-
i
j
=
(
a
b
)
1
/
2
E
h
4
D
i
j
,
E
-
i
j
=
(
a
b
)
1
/
2
E
h
5
E
i
j
,
F
-
i
j
=
(
a
b
)
1
/
2
E
h
6
F
i
j
,
H
-
i
j
=
(
a
b
)
1
/
2
E
h
8
H
i
j
,
I
-
i
=
1
(
a
b
)
(
i
+
1
)
/
2
ρ
I
i
.For simplicity, we drop the overbar in the following analysis. Substituting (5a), (5b), (5c), (5d), (5e), (5f)–(8) into (4c) and applying the Galerkin procedure, we obtain the governing equations of motion of the composite laminated piezoelectric rectangular plate for the dimensionless as follows:
(9a)
w
¨
1
+
μ
1
w
˙
1
+
ω
1
2
w
1
+
(
a
2
cos
⁡
⁡
Ω
1
t
+
a
3
cos
⁡
⁡
Ω
2
t
-
a
4
cos
⁡
⁡
Ω
4
t
)
w
1
+
a
5
w
1
2
w
2
+
a
6
w
2
2
w
1
+
a
7
w
1
3
+
a
8
w
2
3
=
f
1
cos
⁡
⁡
Ω
3
t
,
(9b)
w
¨
2
+
μ
2
w
˙
2
+
ω
2
2
w
2
+
(
b
2
cos
⁡
⁡
Ω
1
t
+
b
3
cos
⁡
⁡
Ω
2
t
+
b
4
cos
⁡
⁡
Ω
4
t
)
w
2
+
b
5
w
2
2
w
1
+
b
6
w
1
2
w
2
+
b
7
w
2
3
+
b
8
w
1
3
=
f
2
cos
⁡
⁡
Ω
3
t
,where the coefficients presented in (9a), (9b) are given in the previous studies [8].The above equations include the cubic terms, in-plane excitation, transverse excitation, and piezoelectric excitation. Equation (9a), (9b) can describe the nonlinear transverse oscillations of the composite laminated piezoelectric rectangular plate. We only study the case of primary parametric resonance and 1 : 3 internal resonances. In this resonant case, there are the following resonant relations:
(10)
ω
1
2
=
ω
2
9
+
ε
σ
1
,
ω
2
2
=
ω
2
+
ε
σ
2
,
Ω
3
=
ω
,
Ω
1
=
Ω
2
=
Ω
4
=
2
ω
3
,
ω
2
≈
3
ω
1
,
where σ
1 and σ
2 are two detuning parameters.The method of multiple scales [23] is employed to (9a), (9b) to find the uniform solutions in the following form:
(11a)
w
1
(
t
,
ε
)
=
x
10
(
T
0
,
T
1
)
+
ε
x
11
(
T
0
,
T
1
)
+
⋯
,
(11b)
w
2
(
t
,
ε
)
=
x
20
(
T
0
,
T
1
)
+
ε
x
21
(
T
0
,
T
1
)
+
⋯
,where T
0
=
t
,
T
1
=
ε
t.Substituting (10) and (11a), (11b) into (9a), (9b) and balancing the coefficients of 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:
(12a)
x
˙
1
=
-
1
2
μ
1
x
1
-
1
2
σ
1
x
2
+
1
4
(
a
2
+
a
3
-
a
4
)
x
2
-
3
2
a
7
x
2
(
x
1
2
+
x
2
2
)
-
1
2
a
5
x
4
(
x
1
2
-
x
2
2
)
-
a
6
x
2
(
x
3
2
+
x
4
2
)
+
a
5
x
1
x
2
x
3
,
(12b)
x
˙
2
=
-
1
2
μ
1
x
2
+
1
2
σ
1
x
1
+
1
4
(
a
2
+
a
3
-
a
4
)
x
1
+
3
2
a
7
x
1
(
x
1
2
+
x
2
2
)
+
1
2
a
5
x
3
(
x
1
2
-
x
2
2
)
+
a
6
x
1
(
x
3
2
+
x
4
2
)
+
a
5
x
1
x
2
x
4
,
(12c)
x
˙
3
=
-
1
2
μ
2
x
3
-
1
6
σ
2
x
4
-
1
3
b
6
x
4
(
x
1
2
+
x
2
2
)
-
1
2
b
7
x
4
(
x
3
2
+
x
4
2
)
-
1
6
b
8
x
2
(
3
x
1
2
-
x
2
2
)
,
(12d)
x
˙
4
=
-
1
2
μ
2
x
4
+
1
6
σ
2
x
3
+
1
3
b
6
x
3
(
x
1
2
+
x
2
2
)
+
1
2
b
7
x
3
(
x
3
2
+
x
4
2
)
+
1
6
b
8
x
1
(
x
1
2
-
3
x
2
2
)
-
1
12
f
2
.


## 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 equation (12a), (12b), (12c), (12d) to a simpler normal form. It is found that there are Z
2
⊕
Z
2 and D
4 symmetries in the averaged equation (12a), (12b), (12c), (12d) without the parameters. Therefore, these symmetries are also held in normal form.We take into account the excitation amplitudef
2 as a perturbation parameter. Amplitude f
2 can be considered as an unfolding parameter when the Shilnikov type multipulse orbits are investigated. Obviously, when we do not consider the perturbation parameter, (12a), (12b), (12c), (12d) become(13a)
x
˙
1
=
-
1
2
μ
1
x
1
+
(
f
0
-
1
2
σ
1
)
x
2
-
3
2
a
7
x
2
(
x
1
2
+
x
2
2
)
-
1
2
a
5
x
4
(
x
1
2
-
x
2
2
)
-
a
6
x
2
(
x
3
2
+
x
4
2
)
+
a
5
x
1
x
2
x
3
,
(13b)
x
˙
2
=
-
1
2
μ
1
x
2
+
(
f
0
+
1
2
σ
1
)
x
1
+
3
2
a
7
x
1
(
x
1
2
+
x
2
2
)
+
1
2
a
5
x
3
(
x
1
2
-
x
2
2
)
+
a
6
x
1
(
x
3
2
+
x
4
2
)
+
a
5
x
1
x
2
x
4
,
(13c)
x
˙
3
=
-
1
2
μ
2
x
3
-
1
6
σ
2
x
4
-
1
3
b
6
x
4
(
x
1
2
+
x
2
2
)
-
1
2
b
7
x
4
(
x
3
2
+
x
4
2
)
-
1
6
b
8
x
2
(
3
x
1
2
-
x
2
2
)
,
(13d)
x
˙
4
=
-
1
2
μ
2
x
4
+
1
6
σ
2
x
3
+
1
3
b
6
x
3
(
x
1
2
+
x
2
2
)
+
1
2
b
7
x
3
(
x
3
2
+
x
4
2
)
+
1
6
b
8
x
1
(
x
1
2
-
3
x
2
2
)
,where f
0
=
(
1
/
4
)
(
a
2
+
a
3
-
a
4
).Executing the Maple program given by Zhang et al. [24], the nonlinear transformation used here is given as follows:
(14a)
x
1
=
y
1
-
1
4
a
7
y
1
3
+
3
a
5
(
σ
2
-
6
)
σ
2
2
y
1
2
y
3
-
a
6
y
1
y
3
2
-
a
6
y
1
y
4
2
-
3
a
5
(
σ
2
3
-
18
σ
2
2
+
216
σ
2
-
1296
)
σ
2
4
y
2
2
y
3
+
a
5
(
6
σ
2
2
-
72
σ
2
+
432
)
σ
2
3
y
1
y
2
y
4
,
(14b)
x
2
=
y
2
+
3
2
a
7
y
2
3
+
3
4
a
7
y
1
2
y
2
+
3
a
5
σ
2
y
1
2
y
4
-
3
a
5
(
σ
2
2
-
12
σ
2
+
72
)
σ
2
3
y
2
2
y
4
-
6
a
5
(
σ
2
-
6
)
σ
2
2
y
1
y
2
y
3
,
(14c)
x
3
=
y
3
-
b
8
σ
2
y
1
3
+
3
b
8
(
σ
2
2
-
12
σ
2
+
72
)
σ
2
3
y
1
y
2
2
-
1
3
b
6
y
1
y
2
y
4
,
(14d)
x
4
=
y
4
+
b
8
(
σ
2
3
-
18
σ
2
2
+
216
σ
2
-
1296
)
σ
2
4
y
2
3
-
3
b
8
(
σ
2
-
6
)
σ
2
2
y
1
2
y
2
+
1
3
b
6
y
1
y
2
y
3
.Substituting (14a), (14b), (14c), (14d) into (13a), (13b), (13c), (13d) yields a simpler 3rd-order normal form with the parameters for averaged equation (12a), (12b), (12c), (12d) as follows:
(15a)
y
˙
1
=
-
μ
-
1
y
1
+
(
1
-
σ
-
1
)
y
2
,
(15b)
y
˙
2
=
σ
-
1
y
1
-
μ
-
1
y
2
+
a
6
y
1
(
y
3
2
+
y
4
2
)
+
3
2
a
7
y
1
3
,
(15c)
y
˙
3
=
-
μ
-
2
y
3
-
σ
-
2
y
4
-
1
3
b
6
y
1
2
y
4
-
1
2
b
7
y
4
(
y
3
2
+
y
4
2
)
,
(15d)
y
˙
4
=
σ
-
2
y
3
-
μ
-
2
y
4
+
1
3
b
6
y
1
2
y
3
+
1
2
b
7
y
3
(
y
3
2
+
y
4
2
)
-
f
-
2
,where the coefficients are μ
-
1
=
(
1
/
2
)
μ
1, μ
-
2
=
(
1
/
2
)
μ
2, σ
-
2
=
(
1
/
6
)
σ
2, and f
-
2
=
(
1
/
12
)
f
2, respectively.Further, let(16)
y
3
=
I
cos
⁡
γ
,
y
4
=
I
sin
γ
.
Substituting (16) into (15a), (15b), (15c), (15d) yields(17a)
y
˙
1
=
-
μ
-
1
y
1
+
(
1
-
σ
-
1
)
y
2
,
(17b)
y
˙
2
=
σ
-
1
y
1
-
μ
-
1
y
2
+
a
6
y
1
I
2
+
3
2
a
7
y
1
3
,
(17c)
I
˙
=
-
μ
-
2
I
-
f
-
2
sin
γ
,
(17d)
I
γ
˙
=
σ
-
2
I
+
1
3
b
6
y
1
2
I
+
1
2
b
7
I
3
-
f
-
2
cos
⁡
γ
.In order to get the unfolding of (17a), (17b), (17c), (17d) a linear transformation is introduced:
(18)
[
y
1
y
2
]
=
3
|
a
6
|
|
b
6
|
[
1
-
σ
-
1
0
μ
-
1
1
]
[
u
1
u
2
]
.Substituting (18) into (17a), (17b), (17c), (17d) and canceling nonlinear terms including the parameter σ
-
1 yield the unfolding as follows:
(19a)
u
˙
1
=
u
2
,
(19b)
u
˙
2
=
-
μ
u
1
-
μ
3
u
2
+
η
1
u
1
3
+
a
6
u
1
I
2
,
(19c)
I
˙
=
-
μ
-
2
I
-
f
-
2
sin
γ
,
(19d)
I
γ
˙
=
σ
-
2
I
+
a
6
u
1
2
I
+
α
2
I
3
-
f
-
2
cos
⁡
γ
,where μ
=
μ
-
1
2
-
σ
-
1
(
1
-
σ
-
1
),  μ
3
=
2
μ
-
1,  η
1
=
9
a
6
a
7
/
2
b
6 and α
2
=
(
1
/
2
)
b
7.The scale transformations to be introduced into (19a), (19b), (19c), (19d) are
(20)
μ
-
2
⟶
ε
μ
-
2
,
μ
3
⟶
ε
μ
3
,
f
-
2
⟶
ε
f
-
2
,
η
1
⟶
η
1
,
α
2
⟶
α
2
,
a
6
⟶
a
6
.Then, normal form (19a), (19b), (19c), (19d) can be rewritten in the form with the perturbations(21a)
u
˙
1
=
∂
H
∂
u
2
+
ε
g
u
1
=
u
2
,
(21b)
u
˙
2
=
-
∂
H
∂
u
1
+
ε
g
u
2
=
-
μ
u
1
+
η
1
u
1
3
+
a
6
u
1
I
2
-
ε
μ
3
u
2
,
(21c)
I
˙
=
∂
H
∂
γ
+
ε
g
I
=
-
ε
μ
-
2
I
-
ε
f
-
2
sin
γ
,
(21d)
I
γ
˙
=
-
∂
H
∂
I
+
ε
g
γ
=
σ
-
2
I
+
α
2
I
3
+
a
6
I
u
1
2
-
ε
f
-
2
cos
⁡
γ
,where the Hamiltonian function H is of the form
(22)
H
(
u
1
,
u
2
,
I
,
γ
)
=
1
2
u
2
2
+
1
2
μ
u
1
2
-
1
4
η
1
u
1
4
-
1
2
a
6
I
2
u
1
2
-
1
2
σ
-
2
I
2
-
1
4
α
2
I
4
,
and g
u
1, g
u
2, g
I, and g
γ are the perturbation terms induced by the dissipative effects
(23)
g
u
1
=
0
,
g
u
2
=
-
μ
3
u
2
,
g
I
=
-
μ
-
2
I
-
f
-
2
sin
γ
,
g
γ
=
-
f
-
2
cos
⁡
γ
.


## 4. Heteroclinic Bifurcations of Unperturbed System

In this section, we focus on studying the nonlinear dynamics characteristic of the unperturbed system. Whenε
=
0, it can be known that system from (21a), (21b), (21c), (21d) is an uncoupled two-degree-of-freedom nonlinear system. The variable I appears in the subspace (
u
1
,
u
2
) of (21a), (21b), (21c), (21d) as a parameter since I
˙
=
0. Consider the first two decoupled equations of (21a), (21b), (21c), (21d),(24a)
u
˙
1
=
u
2
,
(24b)
u
˙
2
=
-
μ
u
1
+
η
1
u
1
3
+
a
6
u
1
I
2
.Sinceη
1
>
0, (24a), (24b) can exhibit the heteroclinic bifurcations. It is obvious from (24a), (24b) that when μ
-
a
6
I
2
<
0, the only solution to (24a), (24b) is the trivial zero solution, (
u
1
,
u
2
)
=
(
0
,
0
), which is the saddle point. On the curve defined by μ
=
a
6
I
2, that is,
(25)
I
1,2
=
±
[
μ
-
1
2
-
σ
-
1
(
1
-
σ
-
1
)
a
6
]
1
/
2
,
the trivial zero solution bifurcates into three solutions through a pitchfork bifurcation, which are given by q
0
=
(
0,0
) and q
±
(
I
)
=
(
B
,
0
), respectively, where
(26)
B
=
±
{
1
η
1
[
μ
-
1
2
-
σ
-
1
(
1
-
σ
-
1
)
-
a
6
I
2
]
}
1
/
2
.From the Jacobian matrix evaluated at the nonzero solutions, it can be found that the singular pointq
0 is the center point and the singular points q
±
(
I
) are saddle points. It is observed that I and γ actually represent the amplitude and phase of vibrations. Therefore, we assume that I
≥
0 and (25) becomes
(27)
I
1
=
0
,
I
2
=
[
μ
-
1
2
-
σ
-
1
(
1
-
σ
-
1
)
a
6
]
1
/
2
,
such that for all I
∈
[
I
1
,
I
2
], (24a), (24b) have two hyperbolic saddle points, q
±
(
I
), which are connected by a pair of heteroclinic orbits, u
±
h
(
T
1
,
I
); that is, lim
⁡
T
1
→
±
∞
⁡
u
±
h
(
T
1
,
I
)
=
q
±
(
I
). Thus, in the full four-dimensional phase space, the set defined by(28)
M
=
{
(
u
,
I
,
γ
)
∣
u
=
q
±
(
I
)
,
I
1
<
I
<
I
2
,
0
≤
γ
<
2
π
}
is a two-dimensional invariant manifold.From the results obtained by Feng et al. [9–11], it is known that 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 W
s
(
M
) and W
u
(
M
), respectively. The existence of the heteroclinic orbit of (24a), (24b) to q
±
(
I
)
=
(
B
,
0
) indicates that W
s
(
M
) and W
u
(
M
) intersect nontransversally along a three-dimensional manifold denoted by Γ, which can be written as
(29)
Γ
=
{
∫
0
T
1
D
I
H
(
u
±
h
(
T
1
,
I
)
,
I
)
d
s
+
γ
0
(
u
,
I
,
γ
)
∣
u
=
u
±
h
(
T
1
,
I
)
,
I
1
<
I
<
I
2
,
γ
=
∫
0
T
1
D
I
H
(
u
±
h
(
T
1
,
I
)
,
I
)
d
s
+
γ
0
}
.We analyze the dynamics of the unperturbed system of (21a), (21b), (21c), (21d) restricted to M. Considering the unperturbed system of (21a), (21b), (21c), (21d) restricted to M yields(30a)
I
˙
=
0
,
(30b)
I
γ
˙
=
D
I
H
(
q
±
(
I
)
,
I
)
,
I
1
<
I
<
I
2
,where
(31)
D
I
H
(
q
±
(
I
)
,
I
)
=
-
∂
H
(
q
±
(
I
)
,
I
)
∂
I
=
σ
-
2
I
+
α
2
I
3
+
a
6
I
u
1
2
.From the results obtained by Feng et al. [9–11], it is known that if D
I
H
(
q
±
(
I
)
,
I
)
≠
0, I = constant is called a periodic orbit, and if D
I
H
(
q
±
,
I
)
=
0, I = constant is known as a circle of the singular points. Any value of I
∈
[
I
1
,
I
2
] at which D
I
H
(
q
±
,
I
)
=
0 is a resonant value I and these singular points are resonant singular points. We denote a resonant value by I
r such that
(32)
D
I
H
(
q
±
,
I
)
=
σ
-
2
I
r
+
α
2
I
r
3
+
a
6
I
r
u
1
2
=
0
.Then, we obtain(33)
I
r
=
±
{
σ
-
2
η
1
+
a
6
[
μ
-
1
2
-
σ
-
1
(
1
-
σ
-
1
)
]
a
6
2
-
α
2
η
1
}
1
/
2
.The geometric structure of the stable and unstable manifolds ofM in the full four-dimensional phase space for the unperturbed system of (21a), (21b), (21c), (21d) is given in Figure 2. Since γ represents the phase of the oscillations, when I
=
I
r, the phase shift Δ
γ of oscillations is defined by
(34)
Δ
γ
=
γ
(
+
∞
,
I
r
)
-
γ
(
-
∞
,
I
r
)
.The geometric structure of manifoldsM, W
s
(
M
), and W
u
(
M
) is 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(
u
1
,
u
2
), there exists a pair of heteroclinic orbits connecting to saddle points. Therefore, the homoclinic orbit in the subspace (
I
,
γ
) is, in fact, a heteroclinic connecting in the full four-dimensional space (
u
1
,
u
2
,
I
,
γ
). The phase shift denotes the difference of the value γ when a trajectory leaves and returns to the basin of attraction of M. We will use the phase shift in 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 (24a), (24b). Let ε
1
=
μ
-
a
6
I
2 and μ
3
=
ε
2, (24a), (24b) can be rewritten as(35a)
u
˙
1
=
u
2
,
(35b)
u
˙
2
=
-
ε
1
u
1
+
η
1
u
1
3
-
ε
ε
2
u
2
.Setε
=
0; (35a), (35b) is a system with the Hamiltonian function
(36)
H
-
(
u
1
,
u
2
)
=
1
2
u
2
2
+
1
2
ε
1
u
1
2
-
1
4
η
1
u
1
4
.WhenH
-
=
0, there is a heteroclinic loop Γ
0 which consists of the two hyperbolic saddle points q
±
(
I
) and a pair of heteroclinic orbits u
±
(
T
1
). In order to calculate the phase shift and the extended Melnikov function, it is necessary to obtain the equations of a pair of heteroclinic orbits, which are given as follows:
(37a)
u
1
(
T
1
)
=
±
ε
1
η
1
tanh
⁡
(
2
ε
1
2
T
1
)
,
(37b)
u
2
(
T
1
)
=
±
ε
1
2
η
1
sech
2
(
2
ε
1
2
T
1
)
.We turn our attention to the computation of the phase shift. Substituting the first equation of (37a), (37b) into the fourth equation of the unperturbed system of (21a), (21b), (21c), (21d) and integrating yield
(38)
γ
(
T
1
)
=
ω
r
T
1
-
a
6
2
ε
1
η
1
tanh
⁡
(
2
ε
1
2
T
1
)
+
γ
0
,
where ω
r
=
σ
-
2
+
α
2
I
2
+
a
6
ε
1
/
η
1.AtI
=
I
r, there is ω
r
≡
0. Therefore, the phase shift may be expressed as
(39)
Δ
γ
=
[
-
2
a
6
2
ε
1
η
1
]
I
=
I
r
=
-
2
a
6
η
1
2
[
μ
-
1
2
-
σ
-
1
(
1
-
σ
-
1
)
-
a
6
I
r
2
]
.


## 5. Existence of Multipulse Orbits

After obtaining detailed information on the nonlinear dynamic characteristics of the subspace(
u
1
,
u
2
) for the unperturbed system from (21a), (21b), (21c), (21d), the next step is to examine the effects of small perturbation terms (
0
<
ε
≪
1
) on the unperturbed system from (21a), (21b), (21c), (21d). The extended Melnikov method developed by Kovačič et al. [15–17] is utilized to discover the existence of the multipulse orbits and chaotic dynamics of the nonlinear vibration for the laminated composite piezoelectric rectangular plate. We start by studying the influence of such small perturbations on the manifold M. The objective of the research is to identify the parameter regions where the existence of the multipulse orbits is possible in the perturbed phase space. The main aim is to verify whether these parameters satisfy the transversality condition of multipulse chaotic dynamics. It will be shown that these multipulse orbits can occur in the Hamilton system with dissipative perturbations if the parameters meet the transversality condition. The existence of such multipulse orbits provides a robust mechanism for the existence of the complicated dynamics in the perturbed system. In this section, the emphasis is put on the application aspects of the extended Melnikov method to (21a), (21b), (21c), (21d).


### 5.1. Dissipative Perturbations of the Homoclinic Loop

We analyze dynamics of the perturbed system and the influence of small perturbations onM. Based on the analysis by Kovačič et al. [15–17], we know that M along with its stable and unstable manifolds is invariant under small, sufficiently differentiable perturbations. It is noticed that q
±
(
I
) in (24a), (24b) maintains the characteristic of the hyperbolic singular point under small perturbations, in particular, M
→
M
ε. Therefore, we obtain
(40)
M
=
M
ε
=
{
(
u
,
I
,
γ
)
∣
u
=
q
±
(
I
)
,
I
1
<
I
<
I
2
,
0
≤
γ
<
2
π
}
.Considering the last two equations of (21a), (21b), (21c), (21d) yields(41a)
I
˙
=
-
μ
-
2
I
-
f
-
2
sin
⁡
γ
,
(41b)
γ
˙
=
σ
-
2
+
α
2
I
2
+
a
6
u
1
2
-
f
-
2
cos
⁡
⁡
γ
I
.It is known from the above analysis that the last two equations of (21a), (21b), (21c), (21d) are of a pair of pure imaginary eigenvalues. Therefore, the resonance can occur in (41a), (41b). Also introduce the scale transformations
(42)
μ
-
2
⟶
ε
μ
-
2
,
I
=
I
r
+
ε
h
,
f
-
2
⟶
ε
f
-
2
,
T
1
⟶
T
1
ε
.Substituting the above transformations into (41a), (41b) yields(43a)
h
˙
=
-
μ
-
2
I
r
-
f
-
2
sin
⁡
γ
-
ε
μ
-
2
h
,
(43b)
γ
˙
=
-
2
δ
η
1
I
r
h
-
ε
(
δ
η
1
h
2
+
f
-
2
I
r
cos
⁡
⁡
γ
)
,where δ
=
a
6
2
-
α
2
η
1.Whenε
=
0, (43a), (43b) become(44a)
h
˙
=
-
μ
-
2
I
r
-
f
-
2
sin
⁡
γ
,
(44b)
γ
˙
=
-
2
δ
η
1
I
r
h
.The unperturbed system from (44a), (44b) is a Hamilton system with the function
(45)
H
^
D
(
h
,
γ
)
=
-
μ
-
2
I
r
γ
+
f
-
2
cos
⁡
⁡
γ
+
δ
η
1
I
r
h
2
.The singular points of (44a), (44b) are given as
(46)
P
0
=
(
0
,
γ
c
)
=
(
0
,
-
arcsin
⁡
(
μ
-
2
I
r
f
-
2
)
)
,
Q
0
=
(
0
,
γ
s
)
=
(
0
,
π
+
arcsin
⁡
(
μ
-
2
I
r
f
-
2
)
)
.Based on the characteristic equations evaluated at the two singular pointsP
0 and Q
0, we can know the stabilities of these singular points. Therefore, it is known that the singular point P
0 is a center point. The singular point Q
0 is a saddle which is connected to itself by a homoclinic orbit. The phase portrait of system for (44a), (44b) 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 Q
0 remains a hyperbolic singular point Q
ε of the saddle stability type. For small perturbations, the singular point P
0 becomes a hyperbolic sink P
ε. The phase portrait of the perturbed system from (43a), (43b) is depicted in Figure 3(b).Using the function (45), at h
=
0, and substituting γ
s in (46) into (45), the estimate of the basin of the attractor for γ
min
⁡ is obtained as
(47)
γ
min
⁡
-
f
-
2
μ
-
2
I
r
cos
⁡
⁡
γ
min
⁡
=
π
+
arcsin
⁡
μ
-
2
I
r
f
-
2
+
f
-
2
2
-
μ
-
2
2
I
r
2
μ
-
2
I
r
.Define an annulusA
ε near I
=
I
r as
(48)
A
ε
=
{
(
u
1
,
u
2
,
I
,
γ
)
∣
u
1
=
B
,
u
2
=
0
,
|
I
-
I
r
|
<
ε
C
,
γ
∈
T
L
}
,
where C is a constant and is sufficiently large so that the unperturbed orbit is enclosed within the annulus.It is noticed that the three-dimensional stable and unstable manifolds ofA
ε, denoted as W
s
(
A
ε
) and W
u
(
A
ε
), are the subsets of the manifolds W
s
(
M
ε
) and W
u
(
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 orbits are obtained.


### 5.2. Thek-Pulse Melnikov Function

Most researchers focused on Shilnikov type single-pulse chaotic dynamics of the high-dimensional nonlinear systems from the thin plate structures in the past. There exist multipulse chaotic dynamics in the practical engineering systems. The extended Melnikov method is a kind of theory which can be used to investigate the multipulse jumping orbits in the high-dimensional nonlinear systems. Since the theory on multipulse chaotic dynamics is very esoteric and abstract, it is difficult to be extended to solve the engineering problems. Up to now, few researchers have made use of the extended Melnikov method to study the Shilnikov type multipulse homoclinic and heteroclinic bifurcations and chaotic dynamics of high-dimensional nonlinear systems in engineering applications.The extended Melnikov method was first presented by Kovačič et al. [15–17], which is an extension of the global perturbation method developed by Feng et al. [9–11]. Camassa et al. [17] gave the detailed procedure of mathematical proof on the extended Melnikov method, which unifies several disjoint perturbation theoretical methods. This method can be also utilized to detect the Shilnikov type multipulse homoclinic or heteroclinic orbits to the slow manifolds of four-dimensional, near-integrable Hamilton systems or higher-dimensional, nonlinear systems. The extended Melnikov function is different from the usual Melnikov function and describes slow dynamics of the multipulse orbits on the hyperbolic manifold.The key of the extended Melnikov method is how to calculate the extended Melnikov method. The extended Melnikov function is computed by a recursion procedure from the usual 1-pulse Melnikov function and depends on the small perturbation parameterε through a logarithmic function which calculates the asymptotic in the particularly delicate small ε limit. In this paper, the extended Melnikov function can be simplified in the resonant case and does not depend on the perturbation parameter. We have used the extended Melnikov method to investigate heteroclinic bifurcations and multipulse chaotic dynamics of the laminated composite piezoelectric rectangular plate.We use the extended Melnikov method described by Kovačič et al. [15–17] to find the Shilnikov type multipulse orbits for nonlinear vibration for the laminated composite piezoelectric rectangular plate. We search for the multipulse excursions to find the nondegenerate zeroes of the extended Melnikov function M
k
(
ε
,
I
,
γ
0
,
μ
-
2
) with the certain combination of parameters ε
,
I
,
γ
0, and μ
-
2, which we name the k-pulse Melnikov function.It is important to obtain the detailed expression of thek-pulse Melnikov function. We compute the 1-pulse Melnikov function based on the formula obtained by Kovačič et al. [15–17] at the resonant case I
=
I
r. The 1-pulse Melnikov function M
1
(
ε
,
I
r
,
γ
0
,
μ
-
2
) coincides with the standard Melnikov function M
(
I
r
,
γ
0
,
μ
-
2
). The 1-pulse Melnikov function M
(
I
r
,
γ
0
,
μ
-
2
) on both heteroclinic manifolds W
s
(
M
) and W
u
(
M
) is given as follows:
(49)
M
(
I
r
,
γ
0
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
∫
-
∞
+
∞
〈
n
(
p
h
(
t
)
)
,
g
(
p
h
(
t
)
,
μ
-
2
,
0
)
〉
d
T
1
=
∫
-
∞
+
∞
(
∂
H
∂
u
1
g
u
1
+
∂
H
∂
u
2
g
u
2
+
∂
H
∂
I
g
I
+
∂
H
∂
γ
g
γ
)
d
T
1
=
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
)
]
.Based on the results given by Kovačič et al. [15–17], it is known that the k-pulse Melnikov function M
k
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
(
k
=
1,2
,
…
) is defined as
(50)
M
k
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
=
∑
j
=
0
k
-
1
M
(
I
r
,
j
Δ
γ
(
I
r
)
+
Γ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
+
γ
0
,
μ
-
2
)
,
where
(51)
Γ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
=
Ω
(
x
-
0
(
I
r
)
,
I
r
)
λ
(
I
r
)
∑
r
=
1
j
log
⁡
⁡
|
ς
(
I
r
)
ε
M
r
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
|
,
for j
=
1
,
…
,
k
-
1 and Γ
0
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
=
0.It is noticed that the angle functionΓ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
) is the complex formula where M
k
(
ε
,
I
r
,
γ
0
,
μ
-
2
) appears as the argument of a logarithm. When resonance occurs, the periodic orbit corresponding to the value I
r degenerates into a circle of equilibria. Under this case, there exists Γ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
)  =  0,  (
j
=
0,1
,
…
,
k
-
1
). Based on the expression obtained by Kovačič et al. [15–17], the k-pulse Melnikov function can be written as follows:
(52)
M
k
(
I
r
,
γ
0
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
∑
j
=
0
k
-
1
M
(
I
r
,
γ
0
+
j
Δ
γ
(
I
r
)
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
)
]
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
-
2
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
-
2
a
6
2
ε
1
η
1
)
]
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
+
⋯
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
-
2
(
k
-
1
)
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
-
2
(
k
-
1
)
a
6
2
ε
1
η
1
)
]
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
=
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
-
2
(
k
-
1
)
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
)
]
-
2
2
k
μ
3
3
η
1
ε
1
3
/
2
-
2
2
k
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
.If we setΔ
γ
=
-
2
a
6
(
2
ε
1
/
η
1
) and γ
k
-
1
=
γ
0
+
(
k
-
1
)
(
Δ
γ
/
2
), (52) can be rewritten as follows:
(53)
M
k
(
I
r
,
γ
0
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
M
k
(
I
r
,
γ
k
-
1
-
(
k
-
1
)
Δ
γ
2
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
k
-
1
+
1
2
k
Δ
γ
)
-
cos
⁡
⁡
(
γ
k
-
1
-
1
2
k
Δ
γ
)
]
+
k
μ
3
ε
1
3
a
6
Δ
γ
+
μ
-
2
I
r
2
(
k
Δ
γ
)
=
2
f
-
2
I
r
sin
⁡
γ
k
-
1
sin
⁡
(
1
2
k
Δ
γ
)
+
μ
3
ε
1
3
a
6
(
k
Δ
γ
)
+
μ
-
2
I
r
2
(
k
Δ
γ
)
.Based on Proposition 3.1 given by Kovačič et al. [15–17], the nonfolding condition is always satisfied in the resonant case. We obtain the following two conditions:
(54)
|
(
1
/
2
)
k
Δ
γ
sin
⁡
(
(
1
/
2
)
k
Δ
γ
)
(
μ
3
ε
1
+
3
a
6
μ
-
2
I
r
2
)
3
a
6
f
-
2
I
r
|
<
1
,
a
a
a
a
a
a
a
1
2
k
Δ
γ
≠
n
π
,
n
=
0
,
±
1
,
±
2
,
…
.The main aim of the following analysis focuses on identifying simple zeroes of thek-pulse Melnikov function. Define a set that contains all such simple zeroes to be
(55)
Z
-
n
=
{
(
I
r
,
γ
k
-
1
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
∣
M
k
=
0
,
D
γ
0
M
k
≠
0
}
.
The k-pulse Melnikov function has two simple zeroes in the interval γ
k
-
1
∈
[
0
,
π
]
(56)
γ
-
k
-
1
,
1
=
-
arcsin
⁡
(
1
/
2
)
k
Δ
γ
sin
⁡
(
(
1
/
2
)
k
Δ
γ
)
(
μ
3
ε
1
+
3
a
6
μ
-
2
I
r
2
)
(
3
a
6
f
-
2
I
r
)
,
γ
-
k
-
1
,
2
=
π
-
γ
-
k
-
1
,
1
.


### 5.3. Geometric Structure of Multipulse Orbits

Based on the aforementioned analysis, we obtain the following conclusions. When the parameters ofk, μ
3, ε
1, μ
-
2, a
6, and f
-
2 satisfy condition (54), the k-pulse Melnikov function (53) has simple zeroes at γ
k
-
1
=
γ
-
k
-
1
,
1 and γ
k
-
1
=
γ
-
k
-
1
,
2
=
π
-
γ
-
k
-
1
,
1. For i
=
1 or i
=
2, when the j-pulse Melnikov function M
j
(
I
r
,
γ
-
0
,
i
,
μ
-
2
,
η
1
,
a
6
,
ε
1
) has no simple zeroes, the stable and unstable manifolds W
s
(
M
ε
) and W
u
(
M
ε
) intersect transversely along a symmetric pair of the two-dimensional, k-pulse surfaces ∑
±
,
ε
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
i
). This signifies that the presence of the Shilnikov type  n-pulse orbits leads to chaotic dynamics in the sense of the Smale horseshoes for the nonlinear motion for the laminated composite piezoelectric rectangular plate. In the phase space of the unperturbed system from (21a), (21b), (21c), (21d), this symmetric pair of the two-dimensional, k-pulse surfaces breaks down smoothly onto a pair of limiting k-pulse surfaces, ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
i
), parametrized by (37a), (37b), and (38) with I
=
I
r,  γ
0
=
γ
-
k
-
1
,
i
-
(
k
-
1
)
(
Δ
γ
/
2
)
+
j
Δ
γ, and an arbitrary h. The sign in (49) is determined by the sign of the corresponding j-pulse Melnikov function M
j
(
I
r
,
γ
-
0
,
i
,
μ
-
2
,
η
1
,
a
6
,
ε
1
).From the discussion given by Kovačič et al. [15–17], it is easily found that for γ
-
0
,
i
=
γ
-
k
-
1
,
i
-
(
k
-
1
)
(
Δ
γ
/
2
)
+
j
Δ
γ (i
=
1 or i
=
2), the values of the j-pulse Melnikov functions M
j
(
I
r
,
γ
-
0
,
i
,
μ
-
2
,
η
1
,
a
6
,
ε
1
) are not zero for all j
=
1
,
…
,
k
-
1, and all j have the same sign. It is known that this sign is negative for γ
-
0
,
1 and positive for γ
-
0
,
2. Therefore, the k-pulse heteroclinic surfaces ∑
±
,
ε
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) and ∑
±
,
ε
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
) indeed exist, and the limiting k-pulse surfaces ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) and ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
) also exist when ε
=
0. Since the regions enclosed by the stable and unstable heteroclinic manifolds W
s
(
M
) and W
u
(
M
) are both convex, and the normal vector
(57)
n
=
(
(
-
μ
u
1
+
η
1
u
1
3
+
a
6
I
2
u
1
)
,
-
u
2
,
0,0
)
is known to point out of these manifolds, it demonstrates that the orbits forming each of the surfaces ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) are parametrized by (37a), (37b), and (38) with the alternating signs, and the orbits forming each of the surfaces ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
) are parametrized by (37a), (37b), and (38) with the same signs.For the parameterμ
-
2
=
μ
-, there exist N
-
1 orbit segments O
i
(
μ
-
)
(
i
=
2
,
…
,
N
) on the annulus M, where the end points of the segments O
i
(
μ
-
) are d
i
(
μ
-
) and c
i
(
μ
-
), respectively. The trajectories of (44a), (44b) on the segments O
i
(
μ
-
) travel from the end points d
i
(
μ
-
) to c
i
(
μ
-
) in forward time. Therefore, the end points d
i
(
μ
-
) and c
i
(
μ
-
) are, respectively, referred to as the departure and landing points of the heteroclinic jumping Γ
i. In addition, the line γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
-
Δ
γ
-
(
I
r
) transversely intersects the segments O
i
(
μ
-
) at the end point c
i
(
μ
-
) for i
=
2
,
…
,
N, while the line γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
+
Δ
γ
+
(
I
r
) transversely intersects the segments O
i
+
1
(
μ
-
) at the end point d
i
+
1
(
μ
-
) when i
=
1
,
…
,
N
-
1. For all i
=
2
,
…
,
N
-
1, the difference in the coordinates h of two end points c
i
(
μ
-
) and d
i
+
1
(
μ
-
) is zero, namely,
(58)
h
(
c
i
(
μ
-
)
)
-
h
(
d
i
+
1
(
μ
-
)
)
=
0
.
For each i
=
2
,
…
,
N
-
1, one of the heteroclinic orbits represented by Γ
i and contained in the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
i
) at the value μ
=
μ
-, connects two intersection points c
i
(
μ
-
) and d
i
+
1
(
μ
-
). Therefore, a heteroclinic orbit Γ
1 on the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
1
) connects the certain point c
1
(
μ
-
) on the annulus M to the end point d
2
(
μ
-
) on the segment O
2
(
μ
-
). It is also known that a heteroclinic orbit Γ
N on the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
N
) connects the end point c
N
(
μ
-
) on the segments O
N
(
μ
-
) to the certain point d
N
+
1
(
μ
-
) on the annulus M. According to the study of Kovačič et al. [15–17], there exists an n-bump singular transition orbit or a modified N-bump singular transition orbit. The 3-bump jumping orbit depicted in Figure 5 consists of the heteroclinic orbits Γ
i
(
i
=
1
,
2
,
3
) on the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
i
)
(
i
=
1
,
2
,
3
) at the parameter μ
=
μ
- and the orbit segments O
1
(
μ
-
) and O
2
(
μ
-
) of (44a), (44b). It is known from the above analysis that the orbit segments O
i
(
μ
-
)
(
i
=
2
,
…
,
N
) intersect transversely with the lines γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
+
Δ
γ
+
(
I
r
) and γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
-
Δ
γ
-
(
I
r
).Figure 5
The 3-bump orbit with the single-pulse is depicted.The 2-bump singular surface shown in Figure6 is composed of two single-pulse singular intersection surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) and ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
). This surface connects the singular points of (44a), (44b) that lie on the line γ
=
γ
-
0
,
1
-
Δ
γ
- to those of (44a), (44b) that lie on the line γ
=
γ
-
0
,
1
-
Δ
γ
+ on the annulus M.Figure 6
The 2-pulse singular surfaces∑
(
γ
-
0
,
1
,
γ
-
0
,
2
) are depicted.We obtain a countable infinity of the singular heteroclinic jumping orbits as follows. Each orbit starts along one branch of the manifoldW
(
Q
0
) of the saddle Q
0 on the annulus M. Then, the singular heteroclinic jumping orbit departs from the annulus M, goes along one of the singular k-pulse orbits Γ
k, and lands back at a point on the separatrix that connects the saddle Q
0 to itself on the annulus M. After traveling along the separatrix for a while, the singular heteroclinic jumping orbit takes off again along the singular l-pulse orbit Γ
l and continues such process. Eventually, the singular heteroclinic jumping orbit lands back on the separatrix.Therefore, it is concluded that the multipulse orbits of (21a), (21b), (21c), (21d) consist of several portions of the slow time scale on the hyperbolic manifold M
ε and many fast time scale heteroclinic pulses leaving from the manifold M
ε, and these multipulse heteroclinic orbits form a consecutive and recurrence process.


## 5.1. Dissipative Perturbations of the Homoclinic Loop

We analyze dynamics of the perturbed system and the influence of small perturbations onM. Based on the analysis by Kovačič et al. [15–17], we know that M along with its stable and unstable manifolds is invariant under small, sufficiently differentiable perturbations. It is noticed that q
±
(
I
) in (24a), (24b) maintains the characteristic of the hyperbolic singular point under small perturbations, in particular, M
→
M
ε. Therefore, we obtain
(40)
M
=
M
ε
=
{
(
u
,
I
,
γ
)
∣
u
=
q
±
(
I
)
,
I
1
<
I
<
I
2
,
0
≤
γ
<
2
π
}
.Considering the last two equations of (21a), (21b), (21c), (21d) yields(41a)
I
˙
=
-
μ
-
2
I
-
f
-
2
sin
⁡
γ
,
(41b)
γ
˙
=
σ
-
2
+
α
2
I
2
+
a
6
u
1
2
-
f
-
2
cos
⁡
⁡
γ
I
.It is known from the above analysis that the last two equations of (21a), (21b), (21c), (21d) are of a pair of pure imaginary eigenvalues. Therefore, the resonance can occur in (41a), (41b). Also introduce the scale transformations
(42)
μ
-
2
⟶
ε
μ
-
2
,
I
=
I
r
+
ε
h
,
f
-
2
⟶
ε
f
-
2
,
T
1
⟶
T
1
ε
.Substituting the above transformations into (41a), (41b) yields(43a)
h
˙
=
-
μ
-
2
I
r
-
f
-
2
sin
⁡
γ
-
ε
μ
-
2
h
,
(43b)
γ
˙
=
-
2
δ
η
1
I
r
h
-
ε
(
δ
η
1
h
2
+
f
-
2
I
r
cos
⁡
⁡
γ
)
,where δ
=
a
6
2
-
α
2
η
1.Whenε
=
0, (43a), (43b) become(44a)
h
˙
=
-
μ
-
2
I
r
-
f
-
2
sin
⁡
γ
,
(44b)
γ
˙
=
-
2
δ
η
1
I
r
h
.The unperturbed system from (44a), (44b) is a Hamilton system with the function
(45)
H
^
D
(
h
,
γ
)
=
-
μ
-
2
I
r
γ
+
f
-
2
cos
⁡
⁡
γ
+
δ
η
1
I
r
h
2
.The singular points of (44a), (44b) are given as
(46)
P
0
=
(
0
,
γ
c
)
=
(
0
,
-
arcsin
⁡
(
μ
-
2
I
r
f
-
2
)
)
,
Q
0
=
(
0
,
γ
s
)
=
(
0
,
π
+
arcsin
⁡
(
μ
-
2
I
r
f
-
2
)
)
.Based on the characteristic equations evaluated at the two singular pointsP
0 and Q
0, we can know the stabilities of these singular points. Therefore, it is known that the singular point P
0 is a center point. The singular point Q
0 is a saddle which is connected to itself by a homoclinic orbit. The phase portrait of system for (44a), (44b) 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 Q
0 remains a hyperbolic singular point Q
ε of the saddle stability type. For small perturbations, the singular point P
0 becomes a hyperbolic sink P
ε. The phase portrait of the perturbed system from (43a), (43b) is depicted in Figure 3(b).Using the function (45), at h
=
0, and substituting γ
s in (46) into (45), the estimate of the basin of the attractor for γ
min
⁡ is obtained as
(47)
γ
min
⁡
-
f
-
2
μ
-
2
I
r
cos
⁡
⁡
γ
min
⁡
=
π
+
arcsin
⁡
μ
-
2
I
r
f
-
2
+
f
-
2
2
-
μ
-
2
2
I
r
2
μ
-
2
I
r
.Define an annulusA
ε near I
=
I
r as
(48)
A
ε
=
{
(
u
1
,
u
2
,
I
,
γ
)
∣
u
1
=
B
,
u
2
=
0
,
|
I
-
I
r
|
<
ε
C
,
γ
∈
T
L
}
,
where C is a constant and is sufficiently large so that the unperturbed orbit is enclosed within the annulus.It is noticed that the three-dimensional stable and unstable manifolds ofA
ε, denoted as W
s
(
A
ε
) and W
u
(
A
ε
), are the subsets of the manifolds W
s
(
M
ε
) and W
u
(
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 orbits are obtained.


## 5.2. Thek-Pulse Melnikov Function

Most researchers focused on Shilnikov type single-pulse chaotic dynamics of the high-dimensional nonlinear systems from the thin plate structures in the past. There exist multipulse chaotic dynamics in the practical engineering systems. The extended Melnikov method is a kind of theory which can be used to investigate the multipulse jumping orbits in the high-dimensional nonlinear systems. Since the theory on multipulse chaotic dynamics is very esoteric and abstract, it is difficult to be extended to solve the engineering problems. Up to now, few researchers have made use of the extended Melnikov method to study the Shilnikov type multipulse homoclinic and heteroclinic bifurcations and chaotic dynamics of high-dimensional nonlinear systems in engineering applications.The extended Melnikov method was first presented by Kovačič et al. [15–17], which is an extension of the global perturbation method developed by Feng et al. [9–11]. Camassa et al. [17] gave the detailed procedure of mathematical proof on the extended Melnikov method, which unifies several disjoint perturbation theoretical methods. This method can be also utilized to detect the Shilnikov type multipulse homoclinic or heteroclinic orbits to the slow manifolds of four-dimensional, near-integrable Hamilton systems or higher-dimensional, nonlinear systems. The extended Melnikov function is different from the usual Melnikov function and describes slow dynamics of the multipulse orbits on the hyperbolic manifold.The key of the extended Melnikov method is how to calculate the extended Melnikov method. The extended Melnikov function is computed by a recursion procedure from the usual 1-pulse Melnikov function and depends on the small perturbation parameterε through a logarithmic function which calculates the asymptotic in the particularly delicate small ε limit. In this paper, the extended Melnikov function can be simplified in the resonant case and does not depend on the perturbation parameter. We have used the extended Melnikov method to investigate heteroclinic bifurcations and multipulse chaotic dynamics of the laminated composite piezoelectric rectangular plate.We use the extended Melnikov method described by Kovačič et al. [15–17] to find the Shilnikov type multipulse orbits for nonlinear vibration for the laminated composite piezoelectric rectangular plate. We search for the multipulse excursions to find the nondegenerate zeroes of the extended Melnikov function M
k
(
ε
,
I
,
γ
0
,
μ
-
2
) with the certain combination of parameters ε
,
I
,
γ
0, and μ
-
2, which we name the k-pulse Melnikov function.It is important to obtain the detailed expression of thek-pulse Melnikov function. We compute the 1-pulse Melnikov function based on the formula obtained by Kovačič et al. [15–17] at the resonant case I
=
I
r. The 1-pulse Melnikov function M
1
(
ε
,
I
r
,
γ
0
,
μ
-
2
) coincides with the standard Melnikov function M
(
I
r
,
γ
0
,
μ
-
2
). The 1-pulse Melnikov function M
(
I
r
,
γ
0
,
μ
-
2
) on both heteroclinic manifolds W
s
(
M
) and W
u
(
M
) is given as follows:
(49)
M
(
I
r
,
γ
0
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
∫
-
∞
+
∞
〈
n
(
p
h
(
t
)
)
,
g
(
p
h
(
t
)
,
μ
-
2
,
0
)
〉
d
T
1
=
∫
-
∞
+
∞
(
∂
H
∂
u
1
g
u
1
+
∂
H
∂
u
2
g
u
2
+
∂
H
∂
I
g
I
+
∂
H
∂
γ
g
γ
)
d
T
1
=
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
)
]
.Based on the results given by Kovačič et al. [15–17], it is known that the k-pulse Melnikov function M
k
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
(
k
=
1,2
,
…
) is defined as
(50)
M
k
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
=
∑
j
=
0
k
-
1
M
(
I
r
,
j
Δ
γ
(
I
r
)
+
Γ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
+
γ
0
,
μ
-
2
)
,
where
(51)
Γ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
=
Ω
(
x
-
0
(
I
r
)
,
I
r
)
λ
(
I
r
)
∑
r
=
1
j
log
⁡
⁡
|
ς
(
I
r
)
ε
M
r
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
|
,
for j
=
1
,
…
,
k
-
1 and Γ
0
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
=
0.It is noticed that the angle functionΓ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
) is the complex formula where M
k
(
ε
,
I
r
,
γ
0
,
μ
-
2
) appears as the argument of a logarithm. When resonance occurs, the periodic orbit corresponding to the value I
r degenerates into a circle of equilibria. Under this case, there exists Γ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
)  =  0,  (
j
=
0,1
,
…
,
k
-
1
). Based on the expression obtained by Kovačič et al. [15–17], the k-pulse Melnikov function can be written as follows:
(52)
M
k
(
I
r
,
γ
0
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
∑
j
=
0
k
-
1
M
(
I
r
,
γ
0
+
j
Δ
γ
(
I
r
)
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
)
]
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
-
2
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
-
2
a
6
2
ε
1
η
1
)
]
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
+
⋯
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
-
2
(
k
-
1
)
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
-
2
(
k
-
1
)
a
6
2
ε
1
η
1
)
]
-
2
2
μ
3
3
η
1
ε
1
3
/
2
-
2
2
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
=
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
0
-
a
6
2
ε
1
η
1
-
2
(
k
-
1
)
a
6
2
ε
1
η
1
)
-
cos
⁡
⁡
(
γ
0
+
a
6
2
ε
1
η
1
)
]
-
2
2
k
μ
3
3
η
1
ε
1
3
/
2
-
2
2
k
a
6
μ
-
2
I
r
2
ε
1
1
/
2
η
1
.If we setΔ
γ
=
-
2
a
6
(
2
ε
1
/
η
1
) and γ
k
-
1
=
γ
0
+
(
k
-
1
)
(
Δ
γ
/
2
), (52) can be rewritten as follows:
(53)
M
k
(
I
r
,
γ
0
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
M
k
(
I
r
,
γ
k
-
1
-
(
k
-
1
)
Δ
γ
2
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
=
-
f
-
2
I
r
[
cos
⁡
⁡
(
γ
k
-
1
+
1
2
k
Δ
γ
)
-
cos
⁡
⁡
(
γ
k
-
1
-
1
2
k
Δ
γ
)
]
+
k
μ
3
ε
1
3
a
6
Δ
γ
+
μ
-
2
I
r
2
(
k
Δ
γ
)
=
2
f
-
2
I
r
sin
⁡
γ
k
-
1
sin
⁡
(
1
2
k
Δ
γ
)
+
μ
3
ε
1
3
a
6
(
k
Δ
γ
)
+
μ
-
2
I
r
2
(
k
Δ
γ
)
.Based on Proposition 3.1 given by Kovačič et al. [15–17], the nonfolding condition is always satisfied in the resonant case. We obtain the following two conditions:
(54)
|
(
1
/
2
)
k
Δ
γ
sin
⁡
(
(
1
/
2
)
k
Δ
γ
)
(
μ
3
ε
1
+
3
a
6
μ
-
2
I
r
2
)
3
a
6
f
-
2
I
r
|
<
1
,
a
a
a
a
a
a
a
1
2
k
Δ
γ
≠
n
π
,
n
=
0
,
±
1
,
±
2
,
…
.The main aim of the following analysis focuses on identifying simple zeroes of thek-pulse Melnikov function. Define a set that contains all such simple zeroes to be
(55)
Z
-
n
=
{
(
I
r
,
γ
k
-
1
,
μ
-
2
,
η
1
,
a
6
,
ε
1
)
∣
M
k
=
0
,
D
γ
0
M
k
≠
0
}
.
The k-pulse Melnikov function has two simple zeroes in the interval γ
k
-
1
∈
[
0
,
π
]
(56)
γ
-
k
-
1
,
1
=
-
arcsin
⁡
(
1
/
2
)
k
Δ
γ
sin
⁡
(
(
1
/
2
)
k
Δ
γ
)
(
μ
3
ε
1
+
3
a
6
μ
-
2
I
r
2
)
(
3
a
6
f
-
2
I
r
)
,
γ
-
k
-
1
,
2
=
π
-
γ
-
k
-
1
,
1
.


## 5.3. Geometric Structure of Multipulse Orbits

Based on the aforementioned analysis, we obtain the following conclusions. When the parameters ofk, μ
3, ε
1, μ
-
2, a
6, and f
-
2 satisfy condition (54), the k-pulse Melnikov function (53) has simple zeroes at γ
k
-
1
=
γ
-
k
-
1
,
1 and γ
k
-
1
=
γ
-
k
-
1
,
2
=
π
-
γ
-
k
-
1
,
1. For i
=
1 or i
=
2, when the j-pulse Melnikov function M
j
(
I
r
,
γ
-
0
,
i
,
μ
-
2
,
η
1
,
a
6
,
ε
1
) has no simple zeroes, the stable and unstable manifolds W
s
(
M
ε
) and W
u
(
M
ε
) intersect transversely along a symmetric pair of the two-dimensional, k-pulse surfaces ∑
±
,
ε
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
i
). This signifies that the presence of the Shilnikov type  n-pulse orbits leads to chaotic dynamics in the sense of the Smale horseshoes for the nonlinear motion for the laminated composite piezoelectric rectangular plate. In the phase space of the unperturbed system from (21a), (21b), (21c), (21d), this symmetric pair of the two-dimensional, k-pulse surfaces breaks down smoothly onto a pair of limiting k-pulse surfaces, ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
i
), parametrized by (37a), (37b), and (38) with I
=
I
r,  γ
0
=
γ
-
k
-
1
,
i
-
(
k
-
1
)
(
Δ
γ
/
2
)
+
j
Δ
γ, and an arbitrary h. The sign in (49) is determined by the sign of the corresponding j-pulse Melnikov function M
j
(
I
r
,
γ
-
0
,
i
,
μ
-
2
,
η
1
,
a
6
,
ε
1
).From the discussion given by Kovačič et al. [15–17], it is easily found that for γ
-
0
,
i
=
γ
-
k
-
1
,
i
-
(
k
-
1
)
(
Δ
γ
/
2
)
+
j
Δ
γ (i
=
1 or i
=
2), the values of the j-pulse Melnikov functions M
j
(
I
r
,
γ
-
0
,
i
,
μ
-
2
,
η
1
,
a
6
,
ε
1
) are not zero for all j
=
1
,
…
,
k
-
1, and all j have the same sign. It is known that this sign is negative for γ
-
0
,
1 and positive for γ
-
0
,
2. Therefore, the k-pulse heteroclinic surfaces ∑
±
,
ε
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) and ∑
±
,
ε
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
) indeed exist, and the limiting k-pulse surfaces ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) and ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
) also exist when ε
=
0. Since the regions enclosed by the stable and unstable heteroclinic manifolds W
s
(
M
) and W
u
(
M
) are both convex, and the normal vector
(57)
n
=
(
(
-
μ
u
1
+
η
1
u
1
3
+
a
6
I
2
u
1
)
,
-
u
2
,
0,0
)
is known to point out of these manifolds, it demonstrates that the orbits forming each of the surfaces ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) are parametrized by (37a), (37b), and (38) with the alternating signs, and the orbits forming each of the surfaces ∑
±
,
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
) are parametrized by (37a), (37b), and (38) with the same signs.For the parameterμ
-
2
=
μ
-, there exist N
-
1 orbit segments O
i
(
μ
-
)
(
i
=
2
,
…
,
N
) on the annulus M, where the end points of the segments O
i
(
μ
-
) are d
i
(
μ
-
) and c
i
(
μ
-
), respectively. The trajectories of (44a), (44b) on the segments O
i
(
μ
-
) travel from the end points d
i
(
μ
-
) to c
i
(
μ
-
) in forward time. Therefore, the end points d
i
(
μ
-
) and c
i
(
μ
-
) are, respectively, referred to as the departure and landing points of the heteroclinic jumping Γ
i. In addition, the line γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
-
Δ
γ
-
(
I
r
) transversely intersects the segments O
i
(
μ
-
) at the end point c
i
(
μ
-
) for i
=
2
,
…
,
N, while the line γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
+
Δ
γ
+
(
I
r
) transversely intersects the segments O
i
+
1
(
μ
-
) at the end point d
i
+
1
(
μ
-
) when i
=
1
,
…
,
N
-
1. For all i
=
2
,
…
,
N
-
1, the difference in the coordinates h of two end points c
i
(
μ
-
) and d
i
+
1
(
μ
-
) is zero, namely,
(58)
h
(
c
i
(
μ
-
)
)
-
h
(
d
i
+
1
(
μ
-
)
)
=
0
.
For each i
=
2
,
…
,
N
-
1, one of the heteroclinic orbits represented by Γ
i and contained in the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
i
) at the value μ
=
μ
-, connects two intersection points c
i
(
μ
-
) and d
i
+
1
(
μ
-
). Therefore, a heteroclinic orbit Γ
1 on the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
1
) connects the certain point c
1
(
μ
-
) on the annulus M to the end point d
2
(
μ
-
) on the segment O
2
(
μ
-
). It is also known that a heteroclinic orbit Γ
N on the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
N
) connects the end point c
N
(
μ
-
) on the segments O
N
(
μ
-
) to the certain point d
N
+
1
(
μ
-
) on the annulus M. According to the study of Kovačič et al. [15–17], there exists an n-bump singular transition orbit or a modified N-bump singular transition orbit. The 3-bump jumping orbit depicted in Figure 5 consists of the heteroclinic orbits Γ
i
(
i
=
1
,
2
,
3
) on the limiting surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
0
,
i
)
(
i
=
1
,
2
,
3
) at the parameter μ
=
μ
- and the orbit segments O
1
(
μ
-
) and O
2
(
μ
-
) of (44a), (44b). It is known from the above analysis that the orbit segments O
i
(
μ
-
)
(
i
=
2
,
…
,
N
) intersect transversely with the lines γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
+
Δ
γ
+
(
I
r
) and γ
=
γ
-
0
,
i
(
I
r
,
μ
-
)
-
Δ
γ
-
(
I
r
).Figure 5
The 3-bump orbit with the single-pulse is depicted.The 2-bump singular surface shown in Figure6 is composed of two single-pulse singular intersection surfaces ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
1
) and ∑
0
μ
-
2
,
η
1
,
a
6
,
ε
1
(
γ
-
k
-
1
,
2
). This surface connects the singular points of (44a), (44b) that lie on the line γ
=
γ
-
0
,
1
-
Δ
γ
- to those of (44a), (44b) that lie on the line γ
=
γ
-
0
,
1
-
Δ
γ
+ on the annulus M.Figure 6
The 2-pulse singular surfaces∑
(
γ
-
0
,
1
,
γ
-
0
,
2
) are depicted.We obtain a countable infinity of the singular heteroclinic jumping orbits as follows. Each orbit starts along one branch of the manifoldW
(
Q
0
) of the saddle Q
0 on the annulus M. Then, the singular heteroclinic jumping orbit departs from the annulus M, goes along one of the singular k-pulse orbits Γ
k, and lands back at a point on the separatrix that connects the saddle Q
0 to itself on the annulus M. After traveling along the separatrix for a while, the singular heteroclinic jumping orbit takes off again along the singular l-pulse orbit Γ
l and continues such process. Eventually, the singular heteroclinic jumping orbit lands back on the separatrix.Therefore, it is concluded that the multipulse orbits of (21a), (21b), (21c), (21d) consist of several portions of the slow time scale on the hyperbolic manifold M
ε and many fast time scale heteroclinic pulses leaving from the manifold M
ε, and these multipulse heteroclinic orbits form a consecutive and recurrence process.


## 6. 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 (12a), (12b), (12c), (12d) appear when η
1
>
0. Therefore, the above theoretical analysis is focused on the situation which there exist heteroclinic bifurcations in (12a), (12b), (12c), (12d). The parameter η
1 is the combination of the parameters a
6, a
7, and b
6, where η
1
=
(
9
a
6
a
7
)
/
2
b
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 a
6, a
7, and b
6 are chosen to satisfy η
1
>
0.We choose the averaged equation (12a), (12b), (12c), (12d) to conduct numerical simulations. A numerical approach through the computer software Matlabis 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
,
μ
2 and transverse excitation f
2 play an important role in the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate. In addition, the parameters a
2 and a
3 are related to the in-plane excitation in the x-direction and the in-plane excitation in the y-direction, respectively. The parameter a
4 is the piezoelectric excitation which reflects the characteristics of the piezoelectric material. Hence, the parameters μ
1,  a
2,  a
4, and f
2 are selected as the controlling parameters to discover the law for complicated nonlinear dynamics of the laminated composite piezoelectric rectangular plate.We begin to draw bifurcation diagrams of the parametersf
2, μ
1, a
2, and a
4. Bifurcation diagrams describe the vibration law of the modal displacements x
1 and x
3, respectively, when the parameters f
2,  μ
1,  a
2, and a
4 change in a certain region. We draw bifurcation diagrams according to the rules of the Runge-Kutta algorithm and the Poincaré map theory. For the periodic motions, Poincaré map is of several separate points. For a chaotic motion, the Poincaré map consists of a number of points on the limited Poincaré section. Therefore, it can be observed that chaotic motion and periodic motion of nonlinear system appear from bifurcation diagrams. The chaotic and periodic responses can be identified by several conventional criteria such as phase portraits and Poincaré map. Based on the response law of bifurcation diagrams, phase portraits and Poincaré map are utilized to further verify the existence of the chaotic motions and the periodic motions. In order to compare the influence of these parameters f
2,  μ
1,  a
2, and a
4 on nonlinear vibration in the laminated composite piezoelectric rectangular plate, we choose the same initial conditions to carry out numerical simulation.Figure7 illustrates the bifurcation diagram of the laminated composite piezoelectric rectangular plate when the excitation f
2 varies in the interval f
2
=
2
~
200. Other parameters and initial conditions are chosen as σ
1
=
1.83,  σ
2
=
1.97,  μ
1
=
0.2,  μ
2
=
0.2
,
a
2
=
23.0,  a
3
=
12.0,  a
4
=
13.0,  a
5
=
-
1.01,  a
6
=
-
2.03,  a
7
=
-
2.05,  b
6
=
4.07,  b
7
=
-
3.08,  b
8
=
1.09,  x
10
=
-
0.01,  x
20
=
-
0.05,  x
30
=
-
0.01,  x
40
=
-
0.01. Figures 7(a) and 7(b) represent the bifurcation diagram on the plane (
x
1
,
f
2
) and (
x
3
,
f
2
), respectively. It is observed from Figure 7 that the excitation f
2 is an important parameter that influences on the nonlinear dynamic responses of the laminated composite piezoelectric rectangular plate. Figure 7 shows that the chaotic motion of the laminated composite piezoelectric rectangular plate appears first, followed by a periodic motion of that. With the increase of excitation f
2, Figure 7 presents the following law: chaotic motion→multi-period motion.The bifurcation diagram is obtained for the excitationf
2
=
2~200, and initial conditions x
10
=
-
0.01,  x
20
=
-
0.05,  x
30
=
-
0.01,  x
40
=
-
0.01; (a) the bifurcation diagram on the plane (
x
1
,
f
2
); (b) the bifurcation diagram on the plane (
x
3
,
f
2
).
(a)
(b)We study the impact of the damping parameter on the nonlinear dynamic responses of the laminated composite piezoelectric rectangular plate. Figure8 is the bifurcation diagram of the laminated composite piezoelectric rectangular plate with the damping coefficient μ
1. The figure demonstrates that system is beginning to enter into the region of the chaotic motion then appears the periodic motion window and finally comes into the region of the chaotic motion again as the damping coefficient μ
1 varies in the interval μ
1
=
0.01
~
0.7. Other parameters and initial conditions are the same as those in Figure 7 when excitation is chosen as f
2
=
82.7. Figures 8(a) and 8(b) describe the nonlinear motion of the laminated composite piezoelectric rectangular plate on the planes  (
x
1
,
μ
1
) and (
x
3
,
μ
1
), respectively, as well as the impact of the damping coefficient μ
1 on the system.The bifurcation diagram is obtained for the damping coefficientμ
1
=
0.01~0.7, the excitation f
2
=
82.7, and initial conditions x
10
=
-
0.01,  x
20
=
-
0.05,  x
30
=
-
0.01,  x
40
=
-
0.01; (a) the bifurcation diagram on the plane (
x
1
,
μ
1
); (b) the bifurcation diagram on the plane (
x
3
,
μ
1
).
(a)
(b)Figure9 portraysthe bifurcation diagram for the laminated composite piezoelectric rectangular plate when the in-plane excitation a
2 in the x-direction varies in the interval a
2
=
2~65. Other parameters and initial conditions remain the same as those in Figure 8 when the damping coefficient μ
1 is selected as μ
1
=
0.2. Figures 9(a) and 9(b) display the bifurcation diagram on the plane (
x
1
,
a
2
) and (
x
3
,
a
2
), respectively. Figure 9 presents that the beginning movement of the system is the periodic motion; then the system appears the chaotic motion. With the increase of the excitation a
2, Figure 9 shows the following evolution law: periodic motion→chaotic motion.The bifurcation diagram is obtained for the in-plane excitationa
2
=
2~65, the damping coefficients μ
1
=
0.2 and μ
2
=
0.2, the excitation f
2
=
82.7, and initial conditions x
10
=
-
0.01,  x
20
=
-
0.05,  x
30
=
-
0.01,  x
40
=
-
0.01; (a) the bifurcation diagram on the plane (
x
1
,
a
2
); (b) the bifurcation diagram on the plane (
x
3
,
a
2
).
(a)
(b)Figure10 indicates the bifurcation diagram for the laminated composite piezoelectric rectangular plate when the piezoelectric excitation a
4 varies from a
4
=
2 to a
4
=
120. Other parameters and initial conditions remain the same as those in Figure 9 when the in-plane excitation a
2 is chosen as a
2
=
23. Figures 10(a) and 10(b) demonstrate the bifurcation diagram on the planes  (
x
1
,
a
4
) and (
x
3
,
a
4
), respectively. It is observed from Figure 10 that the piezoelectric excitation a
4 has a significant influence on the complicated nonlinear dynamic behaviors of the laminated composite piezoelectric rectangular plate. As the piezoelectric excitation a
4 increases, Figure 10 reveals the following law: chaotic motion→multiperiod motion→one-period motion→multiperiod motion→chaotic motion.The bifurcation diagram is obtained for the piezoelectric excitationa
4
=
2~120, the excitation f
2
=
82.7, the damping coefficients μ
1
=
0.2 and μ
2
=
0.2, the in-plane excitations a
2
=
23 and a
3
=
12.0, and initial conditions x
10
=
-
0.01,  x
20
=
-
0.05,  x
30
=
-
0.01,  x
40
=
-
0.01; (a) the bifurcation diagram on the plane (
x
1
,
a
4
); (b) the bifurcation diagram on the plane (
x
3
,
a
4
).
(a)
(b)Based on the above bifurcation diagram, the excitationsf
2, a
2, a
4, and the damping coefficient μ
1 are selected as specific values in order to find the multipulse chaotic motions of the laminated composite piezoelectric rectangular plate. Figure 11 indicates existence of the multipulse chaotic motion of the laminated composite piezoelectric rectangular plate when the excitation f
2 is 82.7. In this case, the chosen parameters and initial conditions are the same as those in Figure 7. Figures 11(a) and 11(b) are the three-dimensional phase portrait in the space (
x
1
,
x
2
,
x
3
) and the Poincaré map on the plane (
x
1
,
x
2
), respectively. Figure 11 shows that the excitation f
2 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 whenf
2
=
82.7 and μ
1
=
0.2; (a) the phase portrait in the three-dimensional space (
x
1
,
x
2
,
x
3
); (b) Poincaré map on the plane (
x
1
,
x
2
).
(a)
(b)Besides the excitationsf
2, a
2, a
4 and the damping coefficient μ
1, 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
=
14.37,  σ
2
=
11.42,  μ
1
=
0.2,  μ
2
=
0.2,  a
2
=
30.0,  a
3
=
75.0,  a
4
=
45.0,  a
5
=
-
11.66,  a
6
=
12.27,  f
2
=
122.7,  a
7
=
-
2.68,  b
6
=
-
2.2,  b
7
=
-
9.69,  b
8
=
-
22.32,  x
10
=
-
1.08,  x
20
=
0.5,  x
30
=
-
0.01,  x
40
=
9.16. Comparing with Figures 11 and 12, it is found that there are differences in the phase portrait and the Poincaré map, respectively. From the three-dimensional phase portrait in Figure 12, we can see that there exists obvious multipulse jumping phenomenon. The three-dimensional phase portrait is composed of the four regions. The different regions are connected by the multipulse orbits.The multipulse chaotic motion is obtained whenσ
1
=
14.37,  σ
2
=
11.42,  a
2
=
30.0,  a
3
=
75.0,  a
4
=
45.0,  a
5
=
-
11.66,  a
6
=
12.27,  f
2
=
122.7,  a
7
=
-
2.68,  b
6
=
-
2.2,  b
7
=
-
9.69,  b
8
=
-
22.32,  x
10
=
-
1.08,  x
20
=
0.5,  x
30
=
-
0.01,  x
40
=
9.16; (a) the phase portrait in the three-dimensional space (
x
1
,
x
2
,
x
3
); (b) Poincaré map on the plane (
x
1
,
x
2
).
(a)
(b)In the following numerical simulations, several different sets of parameters and initial conditions 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 f
2
=
92.38. Some parameters and initial conditions are chosen as σ
1
=
3.61,  σ
2
=
3.13,  a
5
=
-
15.01,  b
8
=
4.09,  x
10
=
-
0.01,  x
20
=
-
0.09,  x
30
=
-
0.05,  x
40
=
-
0.05. In this case, other parameters are the same as those in Figure 7. 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. From the three-dimensional phase portrait in Figure 13, it is found that multipulse jumping phenomenon is more prominent.The multipulse chaotic motion is obtained whenσ
1
=
3.61,  σ
2
=
3.13,  a
2
=
23.0,  a
3
=
12.0,  a
4
=
13.0,  a
5
=
-
15.01,  a
6
=
-
2.03,  f
2
=
92.38,  a
7
=
-
2.05,  b
6
=
4.07,  b
7
=
-
3.08,  b
8
=
4.09,  x
10
=
-
0.01,  x
20
=
-
0.09,  x
30
=
-
0.05,  x
40
=
-
0.05; (a) the phase portrait in the three-dimensional space (
x
1
,
x
2
,
x
3
); (b) Poincaré map on the plane (
x
1
,
x
2
).
(a)
(b)


## 7. 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 extended Melnikov method for the case where the averaged equations have one nonsemisimple double zero and a pair of pure imaginary eigenvalues. The extended Melnikov method can be applied to study the Shilnikov type multipulse heteroclinic bifurcations and chaotic dynamics of high-dimensional nonlinear systems in engineering applications. Analysis of the multipulse heteroclinic orbits in the laminated composite piezoelectric rectangular plate demonstrates that such an analysis is a typical singular perturbation problem in which there are two different time scales. Dynamics on the hyperbolic manifoldM
ε are of the slow time scale and the multipulse jumping orbits taking off from this manifold are of the fast time scale. It is shown that the transfer of energy between the two different modes occurs through the multipulse jumping orbits. The studies have led to the following conclusions.(
1
) There exist the Shilnikov type multipulse chaotic motions in nonlinear vibration of the laminated composite piezoelectric rectangular plate. The geometric interpretation of the k-pulse Melnikov function is a signed distance measured along the normal to a heteroclinic manifold, which gives the more delicate local estimates near the hyperbolic manifold. In the resonant case, the k-pulse extended Melnikov function M
k
(
ε
,
I
,
γ
0
,
μ
-
2
) does not depend on the small perturbation parameter 0
<
ε
≪
1, and the nonfolding condition is automatically satisfied, resulting in the angle function Γ
j
(
ε
,
I
r
,
γ
0
,
μ
-
2
)
(
j
=
0
,
1
,
…
,
k
-
1
) being zero. Therefore, the computing procedure of the extended Melnikov function can be simplified.(
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 excitation f
2, the in-plane excitation a
2, the piezoelectric excitation a
4, and the damping coefficient μ
1 are chosen as several different kinds of control parameters. Based on the bifurcation diagrams, the nonlinear complicated dynamic behavior of the laminated composite piezoelectric rectangular plate is controlled by varying the excitations f
2,  a
2,  a
4 and the damping coefficient μ
1, respectively. Therefore, the excitations f
2,  a
2,  a
4 and the damping coefficient μ
1 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 numerical simulations that the shapes of the chaotic motions are completely different. From the three-dimensional phase portraits in Figures 12 and 13, it is found that there exist obvious multipulse jumping phenomena. Therefore, parameters and initial conditions impact 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 summary, both theoretical and numerical studies suggest that chaos for the Smale horseshoe sense in nonlinear motion of the simply supported laminated composite piezoelectric rectangular plate exists.



---

*Source: 345072-2014-03-02.xml*