#### Abstract

Chaos in piezoelectric composite laminated beams has significant implications in the design of this model. Some results for this model have been obtained numerically. With the energy-phase method and numerical simulations, global dynamics of piezoelectric composite laminated beams is investigated in this paper. The average equation of the piezoelectric composite laminated beam is obtained by the normal form theory. The existence of multipulse homoclinic orbits for undisturbed and dissipative cases is analyzed by the energy-phase method, and the mechanism of chaotic motion of the system is given. The effect of the dissipation factor on pulse sequence and layer radius is studied in detail. The chaotic motion of the system is verified by numerical simulations.

#### 1. Introduction

Since piezoelectric composite materials have good thermal stability and high specific strength and specific stiffness, they can be used to manufacture aircraft wings and satellite antenna shells of large launch vehicles. The composite material has the special vibration damping characteristic, so it can reduce the vibration and noise. Moreover, due to the amazing thermal, electrical, chemical, and mechanical properties of nanowires and nanobeams, they can be utilized in micro- and nano-electro-mechanical systems [1], microsensors, microactuators, and atomic force microscopes [2], and microbeams [3]. Therefore, it is of great significance to study nonlinear dynamics of piezoelectric composite laminated beams. Some achievements in this field have been obtained recently. Zhang et al. [4] studied bifurcations and chaotic dynamics of a simply supported symmetric cross-ply composite laminated piezoelectric rectangular plate which are simultaneously forced by the transverse and in-plane excitations and the excitation loaded by piezoelectric layers. The influence of the transverse, in-plane, and piezoelectric excitations on the bifurcations and chaotic behaviors of the composite laminated piezoelectric rectangular plate was investigated numerically there. Using the asymptotic perturbation method and numerical simulations, Yao et al. [5] studied nonlinear oscillations, bifurcations, and chaos of functionally graded materials plate. Abe et al. [6] studied the nonlinear dynamics characteristics of rectangular planar laminated shells with 1 : 1 internal resonance. Using the Galerkin method and multiscale method, Ma et al. [7] investigated nonlinear subharmonic resonance of an orthotropic rectangular laminated composite plate subjected to the transverse harmonic excitation. Based on the general von Karman-type equations and the Reddy third-order shear deformation plate theory, the dynamic model of simply supported laminated piezoelectric composite beam with axial and transverse loading was established by Yao et al. [8]. Bifurcations and chaotic responses were studied by the method of multiple scales and numerical simulations there. Applying an improved Fourier series method, Shi et al. [9] studied free and forced vibration characteristics of the moderately thick laminated composite rectangular plate on the elastic Winkler and Pasternak foundations. With the extended Melnikov method, Zhang and Hao [10] investigated global bifurcations and chaotic dynamics of a four-edge simply supported composite laminated piezoelectric rectangular plate under combined in-plane, transverse, and dynamic electrical excitations. Employing a novel double integral multivariable Fourier transformation method combined with discretised higher order partial differential unit step function equations, Gohari et al. proposed a new explicit solution for obtaining static deformation and optimal shape control of smart laminated cantilever piezo composite hybrid plates and beams under thermo-electro-mechanical loads using piezoelectric actuators [11], and a novel explicit analytical solution is proposed for obtaining twisting deformation and optimal shape control of smart laminated cantilever composite plates/beams using inclined piezoelectric actuators [12], respectively. By using first-order shear deformation theory, Gohari and Sharifi [13] developed an efficient two-dimensional quadratic multilayer shell element to predict the linear strain-displacement static deformation of laminated composite plates induced by macro fiber composite actuators.

Multipulse chaotic motions are very common in high-dimensional systems. With the extended Melnikov method, Zhang et al. [14] investigated multipulse jumping chaotic vibrations of the bistable asymmetric laminated composite square panel under foundation force. The energy-phase method [15–19], proposed by Haller, is a global perturbation method for analyzing multipulse chaotic motions of high-dimensional nonlinear systems. It has been used to study multiple pulse chaos in many nonlinear dynamic models recently. With this method, Zhang et al. [20] studied chaotic dynamics of piezoelectric composite laminated rectangular plates. Guo et al. [21] studied the multipulse chaotic dynamic behavior of a simply supported symmetric cross-ply composite laminated rectangular thin plate with the parametric and forcing excitations. Yao et al. [22] investigated the multipulse orbits and chaotic dynamics of a simply supported laminated composite piezoelectric rectangular plate under combined parametric excitation and transverse excitation. Li et al. [23] studied global bifurcations and Shilnikov type multipulse chaotic dynamics for a nonlinear vibration absorber. Sun et al. [24] investigated the multipulse homoclinic orbits and chaotic dynamics of an equivalent circular cylindrical shell for the circular mesh antenna in the case of 1 : 2 internal resonance.

In this paper, global dynamics of a piezoelectric composite laminated beam is investigated with the energy-phase method and numerical simulations. The existence of multipulse homoclinic orbits for undisturbed and dissipative cases is analyzed. The effect of the dissipation factor on pulse sequence and layer radius is studied in detail. Numerical simulations verify the analytical results.

#### 2. Dynamic Model and Its Dynamic Analysis

Consider the piezoelectric composite material laminated beam shown in Figure 1. The upper and lower surfaces of the beam are covered with piezoelectric composite material symmetrical about the midplane, which are subjected to the combined action of lateral and longitudinal excitation. Assume that horizontal and vertical excitations are as follows:

According to Reddy’s high-order shear deformation beam theory, von Karman theory, and Hamilton’s principle, considering the primary resonance and 1 : 9 internal resonance, the average equation can be obtained as follows with the multiscale method [25]:where the parameters in equation (2a–2d) can be seen in [26].

Using the normal form theory, combined with the Maple [27] program, and introducing the transformation , system (2a–2d) can be rewritten as follows [26]:where .

Introduce the following scale transformation:

The normal form with perturbation terms is obtained as follows:

The Hamilton function of system (5a–5d) has the following form:and the disturbance terms are as follows:

#### 3. The Existence of Multipulse Homoclinic Orbits and Chaotic Dynamics

##### 3.1. Dynamics of Unperturbed Systems

When , system (5a–5d) is a decoupled nonlinear system. We analyze the first two equations:

The Hamilton function of system (8) is as follows:

If , all possible fixed points of the system are

When , system (8) has a unique equilibrium point , which is a saddle point. When , system (8) has three equilibrium points, i.e., the center , saddle points , and heteroclinic bifurcation will occur.

If , when , it is easy to see that system (8) has a unique center . When , system (8) has one saddle and two centers . Homoclinic bifurcation will occur. Here, we only discuss the case of homoclinic bifurcation. The discussion of heteroclinic bifurcation is similar.

Define the critical curve , namely,

Along this critical line, pitchfork bifurcation will occur. Since and represent the amplitude and phase of vibration, respectively, it is obvious that , so equation (11) can be written as follows:

For all system (8) has a pair of homoclinic orbits connected to hyperbolic saddle point , namely, . Therefore, there exists a two-dimensional invariant manifold in the four-dimensional phase space as follows:

According to the research results of [28–31], it is concluded that the four-dimensional normal invariant manifold has three-dimensional stable manifold and unstable manifold . Under small disturbance, the three-dimensional stable manifold and the unstable manifold intersect nontransversely along the three-dimensional heteroclinic manifold with the following expression:

Letting , system (8) is simplified asand the Hamilton function is

According to Hamilton principle, the energy at saddle point is equal to that on homoclinic orbit, namely,

From equations (15)–(17), one can obtain that

Separating the variables and integrating yields

From (19), one can obtain a pair of homoclinic orbits with the following expressions:

The undisturbed system restricted to iswhere . According to Haller and Wiggins [15], if , is a constant, and it is a periodic orbit; if , is a constant and system (21) is a singularity ring. For any , The value of satisfying is called the resonant value, denoted as , i.e.,

It is easy to obtain that the resonance value is .

The geometric structure of the stable manifold and the unstable manifold of the undisturbed system (5a–5d) in the four-dimensional phase space is shown in Figure 2.

Next, phase drift is calculated. Substituting (20) into (5d), one can obtain

Integrating equation (23) yieldswhere . At , . Therefore, the phase shift is expressed as follows:

##### 3.2. Dynamic Analysis of the Perturbed System

For a sufficiently small perturbation , the manifold is perturbed as a locally invariant normal hyperbolic manifold , where

System (5a–5d) restricted to the manifold is

We study the dynamical behavior of manifold near the resonance value . Introduce the following scale transformation:

Substituting the above transformation into equation (27) and expanding it into the Taylor series of yield

When , the undisturbed part of system (29) is as follows:

It is easy to know that system (30) is a Hamilton system with the following Hamilton function:

The equilibrium points of system (30) is

It is obvious that the undisturbed system has a center and a saddle point . In the case of small disturbance, the saddle point remains a saddle point after being disturbed, but changes from the center to the focus . The phase diagram of the disturbance system (29) is shown in Figure 3.

**(a)**

**(b)**

##### 3.3. Energy Difference Function and Existence of Multipulse Homoclinic Orbits

###### 3.3.1. The Case of

When , system (5a–5d) can be written as follows:

The Hamilton function is shown in equation (6). Substituting the last two equations of equation (33) into the scale transformation of equation (29) and expanding them into series of , one can get the equations on as follows:

When , the unperturbed Hamilton system is as follows:

According to the energy-phase method proposed by Haller and Wiggins [16–19], the *n*-order energy difference function is

According to equation (36), the *n*-order energy difference function is expressed as follows:

The transversal zero set of the *n*-order order energy difference function is

According to (36), for any *n*, it is satisfied that

So, the transversal zeros of in are

Introduce the following angle transformation:

Then, the transversal set of is

According to (40), all the periodic orbits outside intersect with the cross section of , and it can be obtained that the number of pulses is 1. The periodic orbits in can be classified according to the number of pulses in the orbit. The internal energy sequence can be defined as follows:

Define the open set of inner orbits as follows:where is the energy of the inner boundary of set . Define the pulse sequence as follows:

It is deduced that as the closer the periodic orbit of to the center, the pulse number gradually increases, but the corresponding energy sequence gradually decreases, i.e.,

Define the layer sequence inside the track as follows:

It is easy to know that is an open set, the inner boundary is a periodic orbit in , and the energy of the orbit is . According to Figure 4, the layer radius of can be defined as follows:

###### 3.3.2. The Case of

When , the *n*-order energy difference function with dissipation terms iswhere *A* represents the area enclosed by homoclinic orbits and is the boundary of *A*, i.e.,

From (7), (20), and (25), the third term of (50) can be simplified as follows:

Substituting (7) into the fourth item of (50), one can obtain

According to (51)–(53), the expression of the energy difference function can be obtained as follows:

Next, define the dissipation factor From (54), the zeros of satisfy

It is easy to know the upper bound of the dissipation factor is

When , we can get the upper bound of the maximum number of pulses as follows:

Therefore, it can be concluded that the upper bound of the maximum number of pulses is inversely proportional to the dissipation factor *d*.

Define the transversal zero set of the energy difference function as follows:

For any integer *n* satisfying , the transversal zero of the *n*-order energy difference function in the interval is as follows:where .

Introducing the following angle transformation:one can obtain the following transversal zeros with dissipation term :

Since each multipulse orbit has only one energy function, the energy function corresponding to the center point is denoted by , the energy function corresponding to the homoclinic track is denoted by , and the energy function corresponding to the connecting periodic track is denoted by . The following energy function sequence can be obtained:

The set sequence of the inner orbit is denoted as the pulse sequence is denoted as , and the layer radius sequence is represented by , and their definitions are the same as those of the undamped case. For any inner orbit , the number of pulses , and the layer radius is defined as follows, and the sketch is given in Figure 5.

According to equation (32), the perturbation of center becomes saddle focus . Next, we need to study the existence of multipulse orbits that reside at saddle focus . For the number of pulses , the corresponding Silnikov type 3-pulse trajectory is shown in Figure 6.

Next, it needs to verify that the meeting point of equation (29) is the center of equation (30). The equilibrium point of system (30) is

If has a cross-sectional zero point with respect to , then

Equation (67) can be simplified as follows:

Consequently, one can solve the dissipation factor *d* value as follows:

When (69) is established, it only needs to satisfy

Next, we need to verify whether meets the nondegenerate condition:

After simplification, the opposite condition of (71) is obtained as follows:

When (70) is satisfied, formulas (68) and (72) cannot be established at the same time, so it is easy to know that the nondegradation condition is satisfied.

Finally, we need to verify that the regression point of the *N*-pulse orbit from the convergence point of the slow manifold finally falls into the attraction region of the convergence point. We define the point on the interval , and the distance from the regression point is about long. We havewhere

Then, it is necessary to find the regression point closest to saddle point in . If , is redefined as . If the energy at is less than the energy at saddle point , that is,then

Therefore, for a sufficiently small disturbance , the regression point of the pulse will fall in the attractive region of saddle focus, so the system has Silnikov type multipulse chaotic motion imagination.

#### 4. Numerical Simulations

Using the Mathematica software, we draw the relationship between pulse sequence and phase difference and the relationship between layer radius and phase difference. The relationship between pulse number and phase difference and the relationship between layer radius and phase difference in the case of no perturbation for , and 100 are shown in Figures 7 and 8, respectively. It can be seen that each horizontal line segment corresponds to the corresponding *N* value, and there is an *N*-pulse orbit for the phase difference in this interval. In addition, with the increase in the number of pulses, the distribution of phase difference is more dispersed and the distribution of pulse sequence is more dense. From the image of layer radius and phase difference, it can be concluded that at the bifurcation point, each layer can be bifurcated into two new layers with different pulsation numbers.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

In the case of disturbance, select the parameter value , and the relationship between pulse sequence and phase difference and the relationship between layer radius and phase difference under different dissipation values as shown in Figures 9 and 10 are plotted by Mathematica software. We can draw a conclusion from the graph: with the increase in dissipation factor, the pulse sequence distribution becomes more sparse, and when the dissipation factor reaches a certain value, the homoclinic tree begins to break.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

According to conditions (69) and (76), with the square column geometrical shape, choose the system parameters as follows:

and the initial condition of the system ; using the Runge–Kutta method, the phase portraits and time histories are shown as in Figure 11. From Figure 11, it can be seen that the system undergoes chaotic motions in this case. Since the parameters satisfy the conditions of chaos, numerical simulations agree with the analytical results.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

#### 5. Conclusions

In this paper, the multipulse chaotic motion of piezoelectric composite laminated beams excited by transverse and longitudinal excitation is studied. The energy-phase method is used to explore the existence of multipulse orbits of piezoelectric composite laminated beams residing in slow manifold under Hamilton resonance, and the existence of the homoclinic orbit will lead to chaos. Numerical simulations prove that chaos will occur in the system, and the relationship between pulse sequence, layer radius, and phase difference is obtained under different dissipation values. It is presented that there may exist multipulse orbits which are homoclinic to fixed points on the slow manifold in the resonant case for this system. The effects of the phase shift and dissipation factor on pulse sequence and layer radius are also analyzed. It is concluded that the number of pulses in multipulse orbit decreases gradually, and the homoclinic tree breaks down with the increase in dissipation factor. Through theoretical research and numerical simulation, one can obtain that the selection of system parameters and initial conditions is the main reason for the occurrence of multipulse chaos for nonlinear dynamic systems. In real life, composite laminated beams may appear chaos, which has some hidden dangers to their security and stability. Therefore, we should choose appropriate parameters to avoid chaotic motion as much as possible.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this work.

#### Acknowledgments

The research was supported by the National Natural Science Foundation of China (11772148 and 11872201).