# Chaotic Dynamics of Piezoelectric Composite Laminated Beams

Authors: Liangqiang Zhou; Wen Wang
Journal: Shock and Vibration
 (2021)

---

## 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.


---

## Body

## 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 Figure1. 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:(1a)p=p0+p1cosΩ1t,(1b)q=q0cosΩ2t+q1cosΩ3t.Figure 1
Dynamic model of the piezoelectric laminated beam.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]:(2a)x˙1=−12c3x1+14a13p1−12a13G0e31−12σ1x2−964a11x2x12+x22−2716a11x2x32+x42,(2b)x˙2=−12c3x2+14a13p1−12a13G0e31+12σ1x1+964a11x1x12+x22+2716a11x1x32+x42−2a12q0,(2c)x˙3=−12c3x3−118σ2x4−8164a11x4x32+x42−316a11x4x12+x22,(2d)x˙4=−12c3x4+118σ2x3+8164a11x3x32+x42+316a11x3x12+x22−23a12q1,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 x3=Icosγandx4=Isinγ, system (2a–2d) can be rewritten as follows [26]:(3a)u˙1=u2,(3b)u˙2=c1u1−c2u2−α1u13−β1u1I2,(3c)I˙=−c3I−f12sinγ,(3d)Iγ˙=σ¯2I−α2I3−β1u12I−f12cosγ,where c1=−μ¯12+σ¯11−σ¯1,c2=2μ¯1,c3=μ¯2,α1=3a14a11/2a23,β1=2a11,andα2=6a21.Introduce the following scale transformation:(4)c2⟶εc2,c3⟶εc3,f12⟶εf12.The normal form with perturbation terms is obtained as follows:(5a)u˙1=∂H∂u2+εgu1=u2,(5b)u˙2=−∂H∂u1+εgu2=c1u1−α1u13−β1u1I2−εc2u2,(5c)I˙=∂H∂γ+εgI−εf12sinγ=−εc3I−εf12sinγ,(5d)Iγ˙=−∂H∂I+εgγ−εf12cosγ=σ¯2I−α2I3−β1u12I−εf12cosγ.The Hamilton function of system (5a–5d) has the following form:(6)Hu1,u2,I,γ=12u22−12c1u12+14α1u14+12β1u12I2−12σ¯2I2+14α2I4,and the disturbance terms are as follows:(7)gu1=0,gu2=−c2u2,gI=−c3I,gγ=0.


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

### 3.1. Dynamics of Unperturbed Systems

Whenε=0, system (5a–5d) is a decoupled nonlinear system. We analyze the first two equations:(8)u˙1=u2,u˙2=c1u1−α1u13−β1u1I2.The Hamilton function of system (8) is as follows:(9)H1u1,u2=12u22−12c1u12+14α1u14+12β1u12I2.Ifα1<0, all possible fixed points of the system are(10)q0I=0,0,q±I=B,0=±c1−β1I2α1,0.Whenc1−β1I2>0, system (8) has a unique equilibrium point u1,u2=0,0, which is a saddle point. When c1−β1I2<0, system (8) has three equilibrium points, i.e., the center q0I, saddle points q±I, and heteroclinic bifurcation will occur.Ifα1>0, when c1−β1I2<0, it is easy to see that system (8) has a unique center q0I. When c1−β1I2>0, system (8) has one saddle q0I and two centers q±I. Homoclinic bifurcation will occur. Here, we only discuss the case of homoclinic bifurcation. The discussion of heteroclinic bifurcation is similar.Define the critical curvec1=β1I2, namely,(11)I=±−μ¯12+σ¯11−σ¯1β1.Along this critical line, pitchfork bifurcation will occur. SinceI and γ represent the amplitude and phase of vibration, respectively, it is obvious that I≥0, so equation (11) can be written as follows:(12)I1=0,I2=−μ¯12+σ¯11−σ¯1β1.For allI∈I1,I2, system (8) has a pair of homoclinic orbits u±hT1,I connected to hyperbolic saddle point 0,0, namely, limT1⟶±∞u±hT1,I=q0I. Therefore, there exists a two-dimensional invariant manifold in the four-dimensional phase space as follows:(13)M=u,I,γ|u=q0I,I1≤I≤I2,0≤γ≤2π.According to the research results of [28–31], it is concluded that the four-dimensional normal invariant manifold M has three-dimensional stable manifold WsM and unstable manifold WuM. Under small disturbance, the three-dimensional stable manifold WsM and the unstable manifold WuM intersect nontransversely along the three-dimensional heteroclinic manifold Γ with the following expression:(14)Γ=u,I,γ|u=u±hT1,I,I1<I<I2,γ=∫0T1DIHu±hT1,I,Ids+γ0.Lettingε1=c1−β1I2andδ1=α1, system (8) is simplified as(15)u˙1=u2,u˙2=ε1u1−δ1u13,and the Hamilton function is(16)H1u1,u2=12u22−12ε1u12+14δ1u14.According to Hamilton principle, the energy at saddle point is equal to that on homoclinic orbit, namely,(17)H1u1,u2=H1q0I=0.From equations (15)–(17), one can obtain that(18)u2=du1dt=±εu12−12δ1u14.Separating the variables and integrating yields(19)1ε1arctanhε1ε1−δ1/2u12=∓t.From (19), one can obtain a pair of homoclinic orbits with the following expressions:(20)u1T1=±2ε1δ1sechε1T1,u2T1=∓2δ1ε1tanhε1T1sechε1T1.The undisturbed system restricted toM is(21)I˙=0,Iγ˙=σ¯2I−α2I3−β1u12I=DIHq0I,I,where I1≤I≤I2. According to Haller and Wiggins [15], if DIHq0I,I≠0, I is a constant, and it is a periodic orbit; if DIHq0I,I=0, I is a constant and system (21) is a singularity ring. For any I∈I1,I2, The value of I satisfying DIHq0I,I=0 is called the resonant value, denoted as Ir, i.e.,(22)DIHq0I,I=σ¯2I−α2I3=0.It is easy to obtain that the resonance value isIr=σ¯2/α2.The geometric structure of the stable manifoldWsM and the unstable manifold WuM of the undisturbed system (5a–5d) in the four-dimensional phase space is shown in Figure 2.Figure 2
Geometric structure diagram of manifolds (M), WsM, and WuM in four-dimensional phase space.Next, phase drift is calculated. Substituting (20) into (5d), one can obtain(23)γ˙=σ¯2−α2I2−β12ε1δ1⋅sech2ε1T1.Integrating equation (23) yields(24)γT1=wrT1−2β1ε1δ1tanhε1T1+γ0,where wr=σ¯2−α2I2. At I=Ir, wr=0. Therefore, the phase shift is expressed as follows:(25)ΔγT1=−4β1ε1δ1I=Ir=−4β1α1−μ¯12+σ¯11−σ¯1−β1Ir2.


### 3.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(26)Mε=u1ε,u2ε,I,γ|u1εI,γ=0+εu1′I,γ+Oε2,u2εI,γ=0+εu2′I,γ+Oε2,I1≤I≤I2,0≤γ≤2π.System (5a–5d) restricted to the manifold Mε is(27)I˙=−εc3I−εf12sinγ+Oε2,γ˙=σ¯2−α2I2−εf12Icosγ+Oε2.We study the dynamical behavior of manifoldMε near the resonance value I=Ir. Introduce the following scale transformation:(28)I=Ir+εh,T1⟶T1ε.Substituting the above transformation into equation (27) and expanding it into the Taylor series of ε yield(29)h˙=−c3Ir−f12sinγ−εc3h+Oε,γ˙=−2α2Irh−εα2h2+f12Ircosγ+Oε.Whenε=0, the undisturbed part of system (29) is as follows:(30)h˙=−c3Ir−f12sinγ,γ˙=−2α2Irh.It is easy to know that system (30) is a Hamilton system with the following Hamilton function:(31)H^Dh,γ=−c3Irγ+f12cosγ+a2Irh2.The equilibrium points of system (30) is(32)p0=0,γc=0,−arcsinc3Irf12,p1=0,γs=0,π+arcsinc3Irf12.It is obvious that the undisturbed system has a centerp0 and a saddle point p1. In the case of small disturbance, the saddle point p1 remains a saddle point after being disturbed, but p0 changes from the center to the focus p0ε. The phase diagram of the disturbance system (29) is shown in Figure 3.Figure 3
(a) Phase diagram of the undisturbed system (30) and (b) the disturbed system (29).
(a)(b)


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

#### 3.3.1. The Case ofc2=c3=0

Whenc2=c3=0, system (5a–5d) can be written as follows:(33)u˙1=∂H∂u2=u2,u˙2=−∂H∂u1=c1u1−α1u13−β1u1I2,I˙=∂H∂γ−εf12sinγ=−εf12sinγ,Iγ˙=−∂H∂I−εf12cosγ=σ¯2I−α2I3−β1u12I−εf12cosγ.The Hamilton functionHu1,u2,I,γ 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:(34)h˙=−f12sinγ+Oε,γ˙=−2Irα2h−εα2h2+f12Ircosγ+Oε.Whenε=0, the unperturbed Hamilton system is as follows:(35)h˙=−f12sinγ=∂H¯D∂γ,γ˙=−2Irα2h=−∂H¯D∂h,(36)H¯Dh,γ=f12cosγ+α2Irh2.According to the energy-phase method proposed by Haller and Wiggins [16–19], the n-order energy difference function is(37)ΔnH¯Dh,γ=H¯Dh,γ+nΔγ−H¯Dh,γ.According to equation (36), the n-order energy difference function is expressed as follows:(38)ΔnH¯Dh,γ=f12cosγ+nΔγ−cosγ.The transversal zero set of then-order order energy difference function is(39)Z−n=h,γ|ΔnH¯Dh,γ=0,DγΔnH¯Dh,γ≠0.According to (36), for any n, it is satisfied that(40)nΔγ≠2kπ,k∈Z.So, the transversal zeros ofΔnH¯Dh,γ in γ∈−π,π are(41)γ−1n=π−nΔγ2mod2π,γ−2n=π−π+nΔγ2mod2π.Introduce the following angle transformation:(42)γ+1n=γ−1n+nΔγmod2π,γ+2n=γ−2n+nΔγmod2π.Then, the transversal set ofΔnH¯Dh,γ is(43)Z−n=h,γ|γ∈γ−1n,γ−2n,Z+n=h,γ|γ∈γ+1n,γ+2n.According to (40), all the periodic orbits outside S0 intersect with the cross section of Z−1, and it can be obtained that the number of pulses is 1. The periodic orbits in S0 can be classified according to the number of pulses in the orbit. The internal energy sequence can be defined as follows:(44)h∞=H¯D0,0=f12,hn=minH¯D0,γ−1n,H¯D0,γ−2n,h0=H¯D0,π=−f12.Define the open set of inner orbits as follows:(45)A0=ϕ,An=h,γ∈S0|H¯Dh,γ>hn,n≥1,where hn is the energy of the inner boundary of set An. Define the pulse sequence as follows:(46)N1=1,Nk=minn∈Z+|n>Nk−1,hn<hNk−1,k≥2.It is deduced that as the closer the periodic orbit ofS0 to the center, the pulse number Nk gradually increases, but the corresponding energy sequence hn gradually decreases, i.e.,(47)AN1⊂AN2⊂⋯⊂ANk⊂⋯.Define the layer sequence inside the track as follows:(48)LNk=IntANk\ANk−1.It is easy to know thatLNk is an open set, the inner boundary is a periodic orbit in S0, and the energy of the orbit is hNk. According to Figure 4, the layer radius of LNk can be defined as follows:(49)rNk=minγc−γ−1Nk,γc−γ−2Nk.Figure 4
Layer sequence structure diagram of the undamped system.


#### 3.3.2. The Case ofc2≠,c3≠0

Whenc2≠,c3≠0, the n-order energy difference function with dissipation terms is(50)ΔnH⌢Dγ=H⌢Dh,γ+nΔγ−H⌢Dh,γ−n∫Addu1gu1u1,u2,I,γ+ddu2gu2u1,u2,I,γdu1du2−n∫∂AgIdγ,where A represents the area enclosed by homoclinic orbits and ∂A is the boundary of A, i.e.,(51)H⌢Dh,γ+nΔγ−H⌢Dh,γ=f12cosγ+nΔγ−cosγ−nc3IrΔγ.From (7), (20), and (25), the third term of (50) can be simplified as follows:(52)n∫Addu1gu1u1,u2,I,γ+ddu2gu2u1,u2,I,γdu1du2=−nc2∫Adu1du2=2nc2ε1Δγ3β1.Substituting (7) into the fourth item of (50), one can obtain(53)−n∫∂AgIdγ=nc3Ir∫∂Adγ=nc3IrΔγ.According to (51)–(53), the expression of the energy difference function can be obtained as follows:(54)ΔnH⌢Dγ=f12cosγ+nΔγ−cosγ−2nc2ε1Δγ3β1=−2f12sinγ+nΔγ2sinnΔγ2−2nc2ε1Δγ3β1.Next, define the dissipation factord=c/f12c2=c3=c. From (54), the zeros of ΔnH⌢Dγ satisfy(55)sinγ+nΔγ2=−ndε1Δγ3β1sinnΔγ/2.It is easy to know the upper bound of the dissipation factor is(56)d<dmax=3β1sinnΔγ/2nε1Δγ.Whend<1, we can get the upper bound of the maximum number of pulses as follows:(57)n<nmax=3β1dε1Δγ.Therefore, it can be concluded that the upper bound of the maximum number of pulses is inversely proportional to the dissipation factord.Define the transversal zero set of the energy difference function as follows:(58)zn=h,γ|ΔnH⌢Dγ=0,DγΔnH⌢Dγ≠0.For any integern satisfying nΔγ≠2kπ,k∈Z, the transversal zero of the n-order energy difference function in the interval −π,π is as follows:(59)γ^−1n=π−nΔγ2+αmod2π,γ^−2n=π−π+nΔγ2−αmod2π,where α=arcsin−ndε1Δγ/3β1sinnΔγ/2.Introducing the following angle transformation:(60)γ^+1n=γ^−1n+nΔγmod2π,γ^+2n=γ^−2n+nΔγmod2π,one can obtain the following transversal zeros with dissipation term ΔnH⌢Dγ:(61)Z^−n=h,γ|γ∈γ^−1n,γ^−2n,Z^+n=h,γ|γ∈γ^+1n,γ^+2n.Since each multipulse orbit has only one energy function, the energy function corresponding to the center point is denoted byh¯∞, 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. The following energy function sequence can be obtained:(62)h¯∞=H⌢D0,γc=c3IrarcsindIr+f121−dIr2,(63)h¯0=H⌢D0,γs=−c3Irπ+arcsindIr−f121−dIr2,(64)h¯n=minH⌢D0,γ^−1n,H⌢D0,γ^−2n.The set sequence of the inner orbit is denoted asA^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 γ^⊂L^N^k, the number of pulses N^γ^=N^k, and the layer radius is defined as follows, and the sketch is given in Figure 5.(65)r^Nk=minγ^c−γ^−1Nk,γ^c−γ^−2Nk.Figure 5
Layer sequence structure diagram of the system with damping.According to equation (32), the perturbation of center p0 becomes saddle focus p0ε. Next, we need to study the existence of multipulse orbits that reside at saddle focus p0ε. For the number of pulses n=3, the corresponding Silnikov type 3-pulse trajectory is shown in Figure 6.Figure 6
Silnikov type 3-pulse homoclinic orbit.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(66)p0=0,γc=0,−arcsindIr.IfΔnH⌢Dγ has a cross-sectional zero point with respect to d,Δγ,N, then(67)f12cos−arcsindIr+NΔγ−cos−arcsindIr−2Nε1cΔγ3β1=0.Equation (67) can be simplified as follows:(68)1−dIr21−cosNΔγ=dIrsinNΔγ−2Nε1Δγ3β1.Consequently, one can solve the dissipation factord value as follows:(69)d=1−cosNΔγIrsinNΔγ−2Nε1Δγ/3β12+Ir21−cosNΔγ2.When (69) is established, it only needs to satisfy(70)Δγ≠2kπN,k∈Z.Next, we need to verify whetherp0 meets the nondegenerate condition:(71)DdΔNH⌢Dp0=Ddf12cos−arcsindIr+NΔγ−cos−arcsindIr−2Nε1cΔγ3β1≠0.After simplification, the opposite condition of (71) is obtained as follows:(72)dIr2cosNΔγ−1=1−dIr2IrsinNΔγ−2Nε1Δγ3β1.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 theN-pulse orbit from the convergence point of the slow manifold finally falls into the attraction region of the convergence point. We define the point γ∗N on the interval −π,π, and the distance from the regression point γc+NΔγ is about 2kπ long. We have(73)γ∗N=γs+γc+NΔγ−γsmod2π,where(74)γc=0,−arcsindIr,γs=0,π+arcsindIr.Then, it is necessary to find the regression point closest to saddle pointγs in −π,π. If γ∗N>γs, γ∗N is redefined as γ∗N−2π. If the energy at γ∗N is less than the energy at saddle point γs, that is,(75)H⌢D0,γ∗N<H⌢D0,γs,then(76)cosγ∗N−cosγs<dIrγ∗N−γs.Therefore, 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.


## 3.1. Dynamics of Unperturbed Systems

Whenε=0, system (5a–5d) is a decoupled nonlinear system. We analyze the first two equations:(8)u˙1=u2,u˙2=c1u1−α1u13−β1u1I2.The Hamilton function of system (8) is as follows:(9)H1u1,u2=12u22−12c1u12+14α1u14+12β1u12I2.Ifα1<0, all possible fixed points of the system are(10)q0I=0,0,q±I=B,0=±c1−β1I2α1,0.Whenc1−β1I2>0, system (8) has a unique equilibrium point u1,u2=0,0, which is a saddle point. When c1−β1I2<0, system (8) has three equilibrium points, i.e., the center q0I, saddle points q±I, and heteroclinic bifurcation will occur.Ifα1>0, when c1−β1I2<0, it is easy to see that system (8) has a unique center q0I. When c1−β1I2>0, system (8) has one saddle q0I and two centers q±I. Homoclinic bifurcation will occur. Here, we only discuss the case of homoclinic bifurcation. The discussion of heteroclinic bifurcation is similar.Define the critical curvec1=β1I2, namely,(11)I=±−μ¯12+σ¯11−σ¯1β1.Along this critical line, pitchfork bifurcation will occur. SinceI and γ represent the amplitude and phase of vibration, respectively, it is obvious that I≥0, so equation (11) can be written as follows:(12)I1=0,I2=−μ¯12+σ¯11−σ¯1β1.For allI∈I1,I2, system (8) has a pair of homoclinic orbits u±hT1,I connected to hyperbolic saddle point 0,0, namely, limT1⟶±∞u±hT1,I=q0I. Therefore, there exists a two-dimensional invariant manifold in the four-dimensional phase space as follows:(13)M=u,I,γ|u=q0I,I1≤I≤I2,0≤γ≤2π.According to the research results of [28–31], it is concluded that the four-dimensional normal invariant manifold M has three-dimensional stable manifold WsM and unstable manifold WuM. Under small disturbance, the three-dimensional stable manifold WsM and the unstable manifold WuM intersect nontransversely along the three-dimensional heteroclinic manifold Γ with the following expression:(14)Γ=u,I,γ|u=u±hT1,I,I1<I<I2,γ=∫0T1DIHu±hT1,I,Ids+γ0.Lettingε1=c1−β1I2andδ1=α1, system (8) is simplified as(15)u˙1=u2,u˙2=ε1u1−δ1u13,and the Hamilton function is(16)H1u1,u2=12u22−12ε1u12+14δ1u14.According to Hamilton principle, the energy at saddle point is equal to that on homoclinic orbit, namely,(17)H1u1,u2=H1q0I=0.From equations (15)–(17), one can obtain that(18)u2=du1dt=±εu12−12δ1u14.Separating the variables and integrating yields(19)1ε1arctanhε1ε1−δ1/2u12=∓t.From (19), one can obtain a pair of homoclinic orbits with the following expressions:(20)u1T1=±2ε1δ1sechε1T1,u2T1=∓2δ1ε1tanhε1T1sechε1T1.The undisturbed system restricted toM is(21)I˙=0,Iγ˙=σ¯2I−α2I3−β1u12I=DIHq0I,I,where I1≤I≤I2. According to Haller and Wiggins [15], if DIHq0I,I≠0, I is a constant, and it is a periodic orbit; if DIHq0I,I=0, I is a constant and system (21) is a singularity ring. For any I∈I1,I2, The value of I satisfying DIHq0I,I=0 is called the resonant value, denoted as Ir, i.e.,(22)DIHq0I,I=σ¯2I−α2I3=0.It is easy to obtain that the resonance value isIr=σ¯2/α2.The geometric structure of the stable manifoldWsM and the unstable manifold WuM of the undisturbed system (5a–5d) in the four-dimensional phase space is shown in Figure 2.Figure 2
Geometric structure diagram of manifolds (M), WsM, and WuM in four-dimensional phase space.Next, phase drift is calculated. Substituting (20) into (5d), one can obtain(23)γ˙=σ¯2−α2I2−β12ε1δ1⋅sech2ε1T1.Integrating equation (23) yields(24)γT1=wrT1−2β1ε1δ1tanhε1T1+γ0,where wr=σ¯2−α2I2. At I=Ir, wr=0. Therefore, the phase shift is expressed as follows:(25)ΔγT1=−4β1ε1δ1I=Ir=−4β1α1−μ¯12+σ¯11−σ¯1−β1Ir2.


## 3.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(26)Mε=u1ε,u2ε,I,γ|u1εI,γ=0+εu1′I,γ+Oε2,u2εI,γ=0+εu2′I,γ+Oε2,I1≤I≤I2,0≤γ≤2π.System (5a–5d) restricted to the manifold Mε is(27)I˙=−εc3I−εf12sinγ+Oε2,γ˙=σ¯2−α2I2−εf12Icosγ+Oε2.We study the dynamical behavior of manifoldMε near the resonance value I=Ir. Introduce the following scale transformation:(28)I=Ir+εh,T1⟶T1ε.Substituting the above transformation into equation (27) and expanding it into the Taylor series of ε yield(29)h˙=−c3Ir−f12sinγ−εc3h+Oε,γ˙=−2α2Irh−εα2h2+f12Ircosγ+Oε.Whenε=0, the undisturbed part of system (29) is as follows:(30)h˙=−c3Ir−f12sinγ,γ˙=−2α2Irh.It is easy to know that system (30) is a Hamilton system with the following Hamilton function:(31)H^Dh,γ=−c3Irγ+f12cosγ+a2Irh2.The equilibrium points of system (30) is(32)p0=0,γc=0,−arcsinc3Irf12,p1=0,γs=0,π+arcsinc3Irf12.It is obvious that the undisturbed system has a centerp0 and a saddle point p1. In the case of small disturbance, the saddle point p1 remains a saddle point after being disturbed, but p0 changes from the center to the focus p0ε. The phase diagram of the disturbance system (29) is shown in Figure 3.Figure 3
(a) Phase diagram of the undisturbed system (30) and (b) the disturbed system (29).
(a)(b)


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

### 3.3.1. The Case ofc2=c3=0

Whenc2=c3=0, system (5a–5d) can be written as follows:(33)u˙1=∂H∂u2=u2,u˙2=−∂H∂u1=c1u1−α1u13−β1u1I2,I˙=∂H∂γ−εf12sinγ=−εf12sinγ,Iγ˙=−∂H∂I−εf12cosγ=σ¯2I−α2I3−β1u12I−εf12cosγ.The Hamilton functionHu1,u2,I,γ 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:(34)h˙=−f12sinγ+Oε,γ˙=−2Irα2h−εα2h2+f12Ircosγ+Oε.Whenε=0, the unperturbed Hamilton system is as follows:(35)h˙=−f12sinγ=∂H¯D∂γ,γ˙=−2Irα2h=−∂H¯D∂h,(36)H¯Dh,γ=f12cosγ+α2Irh2.According to the energy-phase method proposed by Haller and Wiggins [16–19], the n-order energy difference function is(37)ΔnH¯Dh,γ=H¯Dh,γ+nΔγ−H¯Dh,γ.According to equation (36), the n-order energy difference function is expressed as follows:(38)ΔnH¯Dh,γ=f12cosγ+nΔγ−cosγ.The transversal zero set of then-order order energy difference function is(39)Z−n=h,γ|ΔnH¯Dh,γ=0,DγΔnH¯Dh,γ≠0.According to (36), for any n, it is satisfied that(40)nΔγ≠2kπ,k∈Z.So, the transversal zeros ofΔnH¯Dh,γ in γ∈−π,π are(41)γ−1n=π−nΔγ2mod2π,γ−2n=π−π+nΔγ2mod2π.Introduce the following angle transformation:(42)γ+1n=γ−1n+nΔγmod2π,γ+2n=γ−2n+nΔγmod2π.Then, the transversal set ofΔnH¯Dh,γ is(43)Z−n=h,γ|γ∈γ−1n,γ−2n,Z+n=h,γ|γ∈γ+1n,γ+2n.According to (40), all the periodic orbits outside S0 intersect with the cross section of Z−1, and it can be obtained that the number of pulses is 1. The periodic orbits in S0 can be classified according to the number of pulses in the orbit. The internal energy sequence can be defined as follows:(44)h∞=H¯D0,0=f12,hn=minH¯D0,γ−1n,H¯D0,γ−2n,h0=H¯D0,π=−f12.Define the open set of inner orbits as follows:(45)A0=ϕ,An=h,γ∈S0|H¯Dh,γ>hn,n≥1,where hn is the energy of the inner boundary of set An. Define the pulse sequence as follows:(46)N1=1,Nk=minn∈Z+|n>Nk−1,hn<hNk−1,k≥2.It is deduced that as the closer the periodic orbit ofS0 to the center, the pulse number Nk gradually increases, but the corresponding energy sequence hn gradually decreases, i.e.,(47)AN1⊂AN2⊂⋯⊂ANk⊂⋯.Define the layer sequence inside the track as follows:(48)LNk=IntANk\ANk−1.It is easy to know thatLNk is an open set, the inner boundary is a periodic orbit in S0, and the energy of the orbit is hNk. According to Figure 4, the layer radius of LNk can be defined as follows:(49)rNk=minγc−γ−1Nk,γc−γ−2Nk.Figure 4
Layer sequence structure diagram of the undamped system.


### 3.3.2. The Case ofc2≠,c3≠0

Whenc2≠,c3≠0, the n-order energy difference function with dissipation terms is(50)ΔnH⌢Dγ=H⌢Dh,γ+nΔγ−H⌢Dh,γ−n∫Addu1gu1u1,u2,I,γ+ddu2gu2u1,u2,I,γdu1du2−n∫∂AgIdγ,where A represents the area enclosed by homoclinic orbits and ∂A is the boundary of A, i.e.,(51)H⌢Dh,γ+nΔγ−H⌢Dh,γ=f12cosγ+nΔγ−cosγ−nc3IrΔγ.From (7), (20), and (25), the third term of (50) can be simplified as follows:(52)n∫Addu1gu1u1,u2,I,γ+ddu2gu2u1,u2,I,γdu1du2=−nc2∫Adu1du2=2nc2ε1Δγ3β1.Substituting (7) into the fourth item of (50), one can obtain(53)−n∫∂AgIdγ=nc3Ir∫∂Adγ=nc3IrΔγ.According to (51)–(53), the expression of the energy difference function can be obtained as follows:(54)ΔnH⌢Dγ=f12cosγ+nΔγ−cosγ−2nc2ε1Δγ3β1=−2f12sinγ+nΔγ2sinnΔγ2−2nc2ε1Δγ3β1.Next, define the dissipation factord=c/f12c2=c3=c. From (54), the zeros of ΔnH⌢Dγ satisfy(55)sinγ+nΔγ2=−ndε1Δγ3β1sinnΔγ/2.It is easy to know the upper bound of the dissipation factor is(56)d<dmax=3β1sinnΔγ/2nε1Δγ.Whend<1, we can get the upper bound of the maximum number of pulses as follows:(57)n<nmax=3β1dε1Δγ.Therefore, it can be concluded that the upper bound of the maximum number of pulses is inversely proportional to the dissipation factord.Define the transversal zero set of the energy difference function as follows:(58)zn=h,γ|ΔnH⌢Dγ=0,DγΔnH⌢Dγ≠0.For any integern satisfying nΔγ≠2kπ,k∈Z, the transversal zero of the n-order energy difference function in the interval −π,π is as follows:(59)γ^−1n=π−nΔγ2+αmod2π,γ^−2n=π−π+nΔγ2−αmod2π,where α=arcsin−ndε1Δγ/3β1sinnΔγ/2.Introducing the following angle transformation:(60)γ^+1n=γ^−1n+nΔγmod2π,γ^+2n=γ^−2n+nΔγmod2π,one can obtain the following transversal zeros with dissipation term ΔnH⌢Dγ:(61)Z^−n=h,γ|γ∈γ^−1n,γ^−2n,Z^+n=h,γ|γ∈γ^+1n,γ^+2n.Since each multipulse orbit has only one energy function, the energy function corresponding to the center point is denoted byh¯∞, 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. The following energy function sequence can be obtained:(62)h¯∞=H⌢D0,γc=c3IrarcsindIr+f121−dIr2,(63)h¯0=H⌢D0,γs=−c3Irπ+arcsindIr−f121−dIr2,(64)h¯n=minH⌢D0,γ^−1n,H⌢D0,γ^−2n.The set sequence of the inner orbit is denoted asA^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 γ^⊂L^N^k, the number of pulses N^γ^=N^k, and the layer radius is defined as follows, and the sketch is given in Figure 5.(65)r^Nk=minγ^c−γ^−1Nk,γ^c−γ^−2Nk.Figure 5
Layer sequence structure diagram of the system with damping.According to equation (32), the perturbation of center p0 becomes saddle focus p0ε. Next, we need to study the existence of multipulse orbits that reside at saddle focus p0ε. For the number of pulses n=3, the corresponding Silnikov type 3-pulse trajectory is shown in Figure 6.Figure 6
Silnikov type 3-pulse homoclinic orbit.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(66)p0=0,γc=0,−arcsindIr.IfΔnH⌢Dγ has a cross-sectional zero point with respect to d,Δγ,N, then(67)f12cos−arcsindIr+NΔγ−cos−arcsindIr−2Nε1cΔγ3β1=0.Equation (67) can be simplified as follows:(68)1−dIr21−cosNΔγ=dIrsinNΔγ−2Nε1Δγ3β1.Consequently, one can solve the dissipation factord value as follows:(69)d=1−cosNΔγIrsinNΔγ−2Nε1Δγ/3β12+Ir21−cosNΔγ2.When (69) is established, it only needs to satisfy(70)Δγ≠2kπN,k∈Z.Next, we need to verify whetherp0 meets the nondegenerate condition:(71)DdΔNH⌢Dp0=Ddf12cos−arcsindIr+NΔγ−cos−arcsindIr−2Nε1cΔγ3β1≠0.After simplification, the opposite condition of (71) is obtained as follows:(72)dIr2cosNΔγ−1=1−dIr2IrsinNΔγ−2Nε1Δγ3β1.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 theN-pulse orbit from the convergence point of the slow manifold finally falls into the attraction region of the convergence point. We define the point γ∗N on the interval −π,π, and the distance from the regression point γc+NΔγ is about 2kπ long. We have(73)γ∗N=γs+γc+NΔγ−γsmod2π,where(74)γc=0,−arcsindIr,γs=0,π+arcsindIr.Then, it is necessary to find the regression point closest to saddle pointγs in −π,π. If γ∗N>γs, γ∗N is redefined as γ∗N−2π. If the energy at γ∗N is less than the energy at saddle point γs, that is,(75)H⌢D0,γ∗N<H⌢D0,γs,then(76)cosγ∗N−cosγs<dIrγ∗N−γs.Therefore, 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.


## 3.3.1. The Case ofc2=c3=0

Whenc2=c3=0, system (5a–5d) can be written as follows:(33)u˙1=∂H∂u2=u2,u˙2=−∂H∂u1=c1u1−α1u13−β1u1I2,I˙=∂H∂γ−εf12sinγ=−εf12sinγ,Iγ˙=−∂H∂I−εf12cosγ=σ¯2I−α2I3−β1u12I−εf12cosγ.The Hamilton functionHu1,u2,I,γ 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:(34)h˙=−f12sinγ+Oε,γ˙=−2Irα2h−εα2h2+f12Ircosγ+Oε.Whenε=0, the unperturbed Hamilton system is as follows:(35)h˙=−f12sinγ=∂H¯D∂γ,γ˙=−2Irα2h=−∂H¯D∂h,(36)H¯Dh,γ=f12cosγ+α2Irh2.According to the energy-phase method proposed by Haller and Wiggins [16–19], the n-order energy difference function is(37)ΔnH¯Dh,γ=H¯Dh,γ+nΔγ−H¯Dh,γ.According to equation (36), the n-order energy difference function is expressed as follows:(38)ΔnH¯Dh,γ=f12cosγ+nΔγ−cosγ.The transversal zero set of then-order order energy difference function is(39)Z−n=h,γ|ΔnH¯Dh,γ=0,DγΔnH¯Dh,γ≠0.According to (36), for any n, it is satisfied that(40)nΔγ≠2kπ,k∈Z.So, the transversal zeros ofΔnH¯Dh,γ in γ∈−π,π are(41)γ−1n=π−nΔγ2mod2π,γ−2n=π−π+nΔγ2mod2π.Introduce the following angle transformation:(42)γ+1n=γ−1n+nΔγmod2π,γ+2n=γ−2n+nΔγmod2π.Then, the transversal set ofΔnH¯Dh,γ is(43)Z−n=h,γ|γ∈γ−1n,γ−2n,Z+n=h,γ|γ∈γ+1n,γ+2n.According to (40), all the periodic orbits outside S0 intersect with the cross section of Z−1, and it can be obtained that the number of pulses is 1. The periodic orbits in S0 can be classified according to the number of pulses in the orbit. The internal energy sequence can be defined as follows:(44)h∞=H¯D0,0=f12,hn=minH¯D0,γ−1n,H¯D0,γ−2n,h0=H¯D0,π=−f12.Define the open set of inner orbits as follows:(45)A0=ϕ,An=h,γ∈S0|H¯Dh,γ>hn,n≥1,where hn is the energy of the inner boundary of set An. Define the pulse sequence as follows:(46)N1=1,Nk=minn∈Z+|n>Nk−1,hn<hNk−1,k≥2.It is deduced that as the closer the periodic orbit ofS0 to the center, the pulse number Nk gradually increases, but the corresponding energy sequence hn gradually decreases, i.e.,(47)AN1⊂AN2⊂⋯⊂ANk⊂⋯.Define the layer sequence inside the track as follows:(48)LNk=IntANk\ANk−1.It is easy to know thatLNk is an open set, the inner boundary is a periodic orbit in S0, and the energy of the orbit is hNk. According to Figure 4, the layer radius of LNk can be defined as follows:(49)rNk=minγc−γ−1Nk,γc−γ−2Nk.Figure 4
Layer sequence structure diagram of the undamped system.


## 3.3.2. The Case ofc2≠,c3≠0

Whenc2≠,c3≠0, the n-order energy difference function with dissipation terms is(50)ΔnH⌢Dγ=H⌢Dh,γ+nΔγ−H⌢Dh,γ−n∫Addu1gu1u1,u2,I,γ+ddu2gu2u1,u2,I,γdu1du2−n∫∂AgIdγ,where A represents the area enclosed by homoclinic orbits and ∂A is the boundary of A, i.e.,(51)H⌢Dh,γ+nΔγ−H⌢Dh,γ=f12cosγ+nΔγ−cosγ−nc3IrΔγ.From (7), (20), and (25), the third term of (50) can be simplified as follows:(52)n∫Addu1gu1u1,u2,I,γ+ddu2gu2u1,u2,I,γdu1du2=−nc2∫Adu1du2=2nc2ε1Δγ3β1.Substituting (7) into the fourth item of (50), one can obtain(53)−n∫∂AgIdγ=nc3Ir∫∂Adγ=nc3IrΔγ.According to (51)–(53), the expression of the energy difference function can be obtained as follows:(54)ΔnH⌢Dγ=f12cosγ+nΔγ−cosγ−2nc2ε1Δγ3β1=−2f12sinγ+nΔγ2sinnΔγ2−2nc2ε1Δγ3β1.Next, define the dissipation factord=c/f12c2=c3=c. From (54), the zeros of ΔnH⌢Dγ satisfy(55)sinγ+nΔγ2=−ndε1Δγ3β1sinnΔγ/2.It is easy to know the upper bound of the dissipation factor is(56)d<dmax=3β1sinnΔγ/2nε1Δγ.Whend<1, we can get the upper bound of the maximum number of pulses as follows:(57)n<nmax=3β1dε1Δγ.Therefore, it can be concluded that the upper bound of the maximum number of pulses is inversely proportional to the dissipation factord.Define the transversal zero set of the energy difference function as follows:(58)zn=h,γ|ΔnH⌢Dγ=0,DγΔnH⌢Dγ≠0.For any integern satisfying nΔγ≠2kπ,k∈Z, the transversal zero of the n-order energy difference function in the interval −π,π is as follows:(59)γ^−1n=π−nΔγ2+αmod2π,γ^−2n=π−π+nΔγ2−αmod2π,where α=arcsin−ndε1Δγ/3β1sinnΔγ/2.Introducing the following angle transformation:(60)γ^+1n=γ^−1n+nΔγmod2π,γ^+2n=γ^−2n+nΔγmod2π,one can obtain the following transversal zeros with dissipation term ΔnH⌢Dγ:(61)Z^−n=h,γ|γ∈γ^−1n,γ^−2n,Z^+n=h,γ|γ∈γ^+1n,γ^+2n.Since each multipulse orbit has only one energy function, the energy function corresponding to the center point is denoted byh¯∞, 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. The following energy function sequence can be obtained:(62)h¯∞=H⌢D0,γc=c3IrarcsindIr+f121−dIr2,(63)h¯0=H⌢D0,γs=−c3Irπ+arcsindIr−f121−dIr2,(64)h¯n=minH⌢D0,γ^−1n,H⌢D0,γ^−2n.The set sequence of the inner orbit is denoted asA^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 γ^⊂L^N^k, the number of pulses N^γ^=N^k, and the layer radius is defined as follows, and the sketch is given in Figure 5.(65)r^Nk=minγ^c−γ^−1Nk,γ^c−γ^−2Nk.Figure 5
Layer sequence structure diagram of the system with damping.According to equation (32), the perturbation of center p0 becomes saddle focus p0ε. Next, we need to study the existence of multipulse orbits that reside at saddle focus p0ε. For the number of pulses n=3, the corresponding Silnikov type 3-pulse trajectory is shown in Figure 6.Figure 6
Silnikov type 3-pulse homoclinic orbit.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(66)p0=0,γc=0,−arcsindIr.IfΔnH⌢Dγ has a cross-sectional zero point with respect to d,Δγ,N, then(67)f12cos−arcsindIr+NΔγ−cos−arcsindIr−2Nε1cΔγ3β1=0.Equation (67) can be simplified as follows:(68)1−dIr21−cosNΔγ=dIrsinNΔγ−2Nε1Δγ3β1.Consequently, one can solve the dissipation factord value as follows:(69)d=1−cosNΔγIrsinNΔγ−2Nε1Δγ/3β12+Ir21−cosNΔγ2.When (69) is established, it only needs to satisfy(70)Δγ≠2kπN,k∈Z.Next, we need to verify whetherp0 meets the nondegenerate condition:(71)DdΔNH⌢Dp0=Ddf12cos−arcsindIr+NΔγ−cos−arcsindIr−2Nε1cΔγ3β1≠0.After simplification, the opposite condition of (71) is obtained as follows:(72)dIr2cosNΔγ−1=1−dIr2IrsinNΔγ−2Nε1Δγ3β1.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 theN-pulse orbit from the convergence point of the slow manifold finally falls into the attraction region of the convergence point. We define the point γ∗N on the interval −π,π, and the distance from the regression point γc+NΔγ is about 2kπ long. We have(73)γ∗N=γs+γc+NΔγ−γsmod2π,where(74)γc=0,−arcsindIr,γs=0,π+arcsindIr.Then, it is necessary to find the regression point closest to saddle pointγs in −π,π. If γ∗N>γs, γ∗N is redefined as γ∗N−2π. If the energy at γ∗N is less than the energy at saddle point γs, that is,(75)H⌢D0,γ∗N<H⌢D0,γs,then(76)cosγ∗N−cosγs<dIrγ∗N−γs.Therefore, 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.


## 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 numberNk and phase difference Δγ and the relationship between layer radius rNk and phase difference Δγ in the case of no perturbation for Nk=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 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.Figure 7
Relationship between pulse sequenceNk and phase difference Δγ: (a) Nk=10; (b) Nk=20; (c) Nk=50; (d) Nk=100.
(a)(b)(c)(d)Figure 8
Relationship between layer radiusrNk and phase difference Δγ: (a) Nk=10; (b) Nk=20; (c) Nk=50; (d) Nk=100.
(a)(b)(c)(d)In the case of disturbance, select the parameter valueε1/3β1=1,Ir=0.005,f12=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.Figure 9
Relationship between pulse sequenceNk and phase difference Δγ (damped): (a) d = 0; (b) d = 10−5; (c) d = 0.01; (d) d = 0.5.
(a)(b)(c)(d)Figure 10
Relationship between layer radiusrNk and phase difference Δγ (damped): (a) d = 0; (b) d = 10−5; (c) d = 0.01; (d) d = 0.5.
(a)(b)(c)(d)According to conditions (69) and (76), with the square column geometrical shape, choose the system parameters as follows:a11=−0.028,a12=0.028,q0=10.94,σ2=2,c3=0.006,a13=0.0193,p1=163,G0=7,ande31=0.56, and the initial condition of the system x1,x2,x3,x4=1.4,−1.3,1.2,−1.6; 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.Figure 11
Phase portraits and time histories of the system. (a) Phase portrait in spacex1,x2,x3. (b) Phase portrait in space x1,x3,x4. (c) Phase portrait on the plane x1,x2. (d) Phase portrait on the plane x3,x4. (e) The waveforms of x1. (f) The waveforms of x3.
(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.



---

*Source: 5582201-2021-04-27.xml*