Chaotic Dynamics of Piezoelectric Composite Laminated Beams

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. +e average equation of the piezoelectric composite laminated beam is obtained by the normal form theory.+e 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. +e effect of the dissipation factor on pulse sequence and layer radius is studied in detail. +e chaotic motion of the system is verified by numerical simulations.


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. e 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]. erefore, 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. e influence of the transverse, inplane, 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 twodimensional 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 highdimensional 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. e energy-phase method [15][16][17][18][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. e existence of multipulse homoclinic orbits for undisturbed and dissipative cases is analyzed. e effect of the dissipation factor on pulse sequence and layer radius is studied in detail. Numerical simulations verify the analytical results.

Dynamic Model and Its Dynamic Analysis
Consider the piezoelectric composite material laminated beam shown in Figure 1. e 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 x 3 � I cos c and x 4 � I sin c, system (2a-2d) can be rewritten as follows [26]: q = q 0 cosΩ 2 t + q 1 cosΩ 3 t p = p 0 + p 1 cosΩ 1 t
Introduce the following scale transformation: e normal form with perturbation terms is obtained as follows: e Hamilton function of system (5a-5d) has the following form: and the disturbance terms are as follows:
If α 1 > 0, when c 1 − β 1 I 2 < 0, it is easy to see that system (8) has a unique center q 0 (I). When c 1 − β 1 I 2 > 0, system (8) has one saddle q 0 (I) and two centers q ± (I). Homoclinic bifurcation will occur. Here, we only discuss the case of homoclinic bifurcation. e discussion of heteroclinic bifurcation is similar.
Define the critical curve c 1 � β 1 I 2 , namely, Along this critical line, pitchfork bifurcation will occur. Since I and c represent the amplitude and phase of vibration, respectively, it is obvious that I ≥ 0, so equation (11) can be written as follows: For all I ∈ [I 1 , I 2 ], system (8) has a pair of homoclinic orbits u h ± (T 1 , I) connected to hyperbolic saddle point (0, 0), namely, lim T 1 ⟶ ± ∞ u h ± (T 1 , I) � q 0 (I). erefore, there exists a two-dimensional invariant manifold in the fourdimensional phase space as follows: According to the research results of [28][29][30][31], it is concluded that the four-dimensional normal invariant manifold M has three-dimensional stable manifold W s (M) and unstable manifold W u (M). Under small disturbance, the threedimensional stable manifold W s (M) and the unstable manifold W u (M) intersect nontransversely along the threedimensional heteroclinic manifold Γ with the following expression: Shock and Vibration 3 Letting and 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: where I 1 ≤ I ≤ I 2 . According to Haller and Wiggins [15], if D I H(q 0 (I), I) ≠ 0, I is a constant, and it is a periodic orbit; if D I H(q 0 (I), I) � 0, I is a constant and system (21) is a singularity ring. For any I ∈ [I 1 , I 2 ], e value of I satisfying D I H(q 0 (I), I) � 0 is called the resonant value, denoted as I r , i.e., It is easy to obtain that the resonance value is I r � ���� � σ 2 /α 2 . e geometric structure of the stable manifold W s (M) and the unstable manifold W u (M) of the undisturbed system (5a-5d) in the four-dimensional phase space is shown in Figure 2.

Dynamic Analysis of the Perturbed System.
For a sufficiently small perturbation ε > 0, the manifold M is perturbed as a locally invariant normal hyperbolic manifold M ε , where System (5a-5d) restricted to the manifold M ε is We study the dynamical behavior of manifold M ε near the resonance value I � I r . Introduce the following scale transformation: Substituting the above transformation into equation (27) and expanding it into the Taylor series of When ε � 0, 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: e equilibrium points of system (30) is It is obvious that the undisturbed system has a center p 0 and a saddle point p 1 . In the case of small disturbance, the saddle point p 1 remains a saddle point after being disturbed, but p 0 changes from the center to the focus p 0ε . e phase diagram of the disturbance system (29) is shown in Figure 3.

Energy Difference Function and Existence of Multipulse
Homoclinic Orbits When c 2 � c 3 � 0, system (5a-5d) can be written as follows: e Hamilton function H(u 1 , u 2 , I, c) 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 M ε as follows: When ε � 0, the unperturbed Hamilton system is as follows: According to the energy-phase method proposed by Haller and Wiggins [16][17][18][19], the n-order energy difference function is According to equation (36), the n-order energy difference function is expressed as follows: e 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 Introduce the following angle transformation: c n +1 � c n −1 + nΔc mod 2π, c n +2 � c n −2 + nΔc mod 2π. (42)

en, the transversal set of Δ n H D (h, c) is
According to (40), all the periodic orbits outside S 0 intersect with the cross section of Z 1 − , and it can be obtained that the number of pulses is 1. e periodic orbits in S 0 can be classified according to the number of pulses in the orbit. e internal energy sequence can be defined as follows: Shock and Vibration Define the open set of inner orbits as follows: where h n is the energy of the inner boundary of set A n . Define the pulse sequence as follows: It is deduced that as the closer the periodic orbit of S 0 to the center, the pulse number N k gradually increases, but the corresponding energy sequence h n gradually decreases, i.e., Define the layer sequence inside the track as follows: It is easy to know that L N k is an open set, the inner boundary is a periodic orbit in S 0 , and the energy of the orbit is h N k . According to Figure 4, the layer radius of L N k can be defined as follows:

3.3.2.
e Case of c 2 ≠ , c 3 ≠ 0. When c 2 ≠ , c 3 ≠ 0, the norder energy difference function with dissipation terms is where A represents the area enclosed by homoclinic orbits and zA is the boundary of A, i.e.,  (30) and (b) the disturbed system (29). Figure 4: Layer sequence structure diagram of the undamped system.

Shock and Vibration
From (7), (20), and (25), the third term of (50) can be simplified as follows: (52) Substituting (7) into the fourth item of (50), one can obtain −n zA g I dc � nc 3 I r zA dc � nc 3 I r Δc. (53) According to (51)-(53), the expression of the energy difference function can be obtained as follows: It is easy to know the upper bound of the dissipation factor is When d < 1, we can get the upper bound of the maximum number of pulses as follows: erefore, it can be concluded that the upper bound of the maximum number of pulses is inversely proportional to the dissipation factor d.
Introducing the following angle transformation: one can obtain the following transversal zeros with dissipation term Δ n H ⌢ D (c): Since each multipulse orbit has only one energy function, the energy function corresponding to the center point is denoted by h ∞ , the energy function corresponding to the homoclinic track is denoted by h 0 , and the energy function corresponding to the connecting periodic track is denoted by h n . e following energy function sequence can be obtained: e set sequence of the inner orbit is denoted as A n , the pulse sequence is denoted as N k , and the layer radius sequence is represented by L N k , and their definitions are the same as those of the undamped case. For any inner orbit c ⊂ L N k , the number of pulses N(c) � N k , and the layer radius is defined as follows, and the sketch is given in Figure 5.
According to equation (32), the perturbation of center p 0 becomes saddle focus p 0ε . Next, we need to study the existence of multipulse orbits that reside at saddle focus p 0ε . For the number of pulses n � 3, 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). e equilibrium point of system (30) is If Δ n H ⌢ D (c) has a cross-sectional zero point with respect to (d, Δc, N), then f 12 cos −arcsin dI r + NΔc − cos −arcsin dI r − 2Nε 1 cΔc Equation (67) can be simplified as follows: Shock and Vibration 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 p 0 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 c N * on the interval [−π, π], and the distance from the regression point c c + NΔc is about 2kπ long. We have where c c � 0, −arcsin dI r , en, it is necessary to find the regression point closest to saddle point c s in [−π, π]. If c N * > c s , c N * is redefined as c N * − 2π. If the energy at c N * is less than the energy at saddle point c s , that is, erefore, for a sufficiently small disturbance ε > 0, 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.

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. e relationship between pulse number N k and phase difference Δc and the relationship between layer radius r N k and phase difference Δc in the case of no perturbation for N k � 10, 30, 50, 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 Npulse orbit for the phase difference Δc 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.
In the case of disturbance, select the parameter value ε 1 /3β 1 � 1, I r � 0.005, f 12 � 0.01, 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.

Conclusions
In this paper, the multipulse chaotic motion of piezoelectric composite laminated beams excited by transverse and longitudinal excitation is studied. e 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. e 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. rough 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. erefore, we should choose appropriate parameters to avoid chaotic motion as much as possible.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this work.