# Multipulse Homoclinic Orbits and Chaotic Dynamics of a Reinforced Composite Plate with Carbon Nanotubes

Authors: Fengxian An; Fangqi Chen; Xiaoxia Bian; Li Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The multipulse homoclinic orbits and chaotic dynamics of a reinforced composite plate with the carbon nanotubes (CNTs) under combined in-plane and transverse excitations are studied in the case of 1 : 1 internal resonance. The method of multiple scales is adopted to derive the averaged equations. From the averaged equations, the normal form theory is applied to reduce the equations to a simpler normal form associated with a double zero and a pair of pure imaginary eigenvalues. The energy-phase method proposed by Haller and Wiggins is utilized to examine the global bifurcations and chaotic dynamics of the CNT-reinforced composite plate. The analytical results demonstrate that the multipulse Shilnikov-type homoclinic orbits and chaotic motions exist in the system. Homoclinic trees are constructed to illustrate the repeated bifurcations of multipulse solutions. In order to verify the theoretical results, numerical simulations are given to show the multipulse Shilnikov-type chaotic motions in the CNT-reinforced composite plate. The results obtained here imply that the motion is chaotic in the sense of the Smale horseshoes for the CNT-reinforced composite plate.


---

## Body

## 1. Introduction

Carbon nanotubes (CNTs), as a new type of advanced materials, have attracted a lot of attention of researchers. This is because CNTs possess high strength and stiffness with high aspect ratio and low density. Due to these properties of CNTs, a number of studies on the nonlinear vibrations and dynamic responses of CNT-reinforced composite plates have been carried out by many researchers in recent years. By employing an equivalent continuum model, Formica et al. [1] investigated the vibration behaviors of CNT-reinforced composite plates. Zhu et al. [2] presented the bending and free vibration analysis of CNT-reinforced composite plates using the finite element method based on the first-order shear deformation plate theory. In their work, the authors showed the effects of the volume fractions of CNTs and the edge-to-thickness ratios on the bending responses, natural frequencies, and mode shapes of the plates. In addition, Wang and Shen [3] examined the nonlinear dynamic response of CNT-reinforced composite plates resting on elastic foundations in thermal environments. The motion equations were derived based on a higher-order shear deformation theory with a von Kármán-type of kinematic nonlinearity. In [4], a mixed Navier-layerwise differential quadrature method was employed by Malekzadeh and Heydarpour for the free vibration and static response analysis of functionally graded carbon nanotube- (FG-CNT-) reinforced composite laminated plates. Based on the conventional Ritz method accompanied with the Lagrangian multiplier technique, Kiani [5] analyzed the free vibration characteristics of FG-CNT-reinforced composite plates located on point supports. Rafiee et al. [6] applied Galerkin’s method to deal with the nonlinear dynamic stability of initially imperfect piezoelectric FG-CNT-reinforced composite plates under combined thermal and electrical loadings. Using the Fourier series expansion and state-space technique, Alibeigloo and Liew [7] investigated the bending behavior of FG-CNT-reinforced composite plates with simply supported edges subjected to thermomechanical loads. Subsequently, this work was extended by Alibeigloo and Emtehani [8] for various boundary conditions by employing the differential quadrature method. Recently, Zhang et al. [9] carried out the analysis of geometrically nonlinear large deformation of triangular FG-CNT-reinforced composite plates using the element-free improved moving least-squares Ritz (IMLS-Ritz) method.In addition, with the increasing applications of functionally graded materials (FGM) in modern technology, the buckling and vibration analysis of FGM structures began to attract a widespread attention. Thus, many studies on the buckling and vibration analysis of FGM shells have been published in the literature. Zhang and Li [10] discussed the dynamic buckling of FGM truncated conical shells subjected to normal impact loads. Bagherizadeh et al. [11] investigated the mechanical buckling of simply-supported FGM cylindrical shells using third-order shear deformation shell theory. They found that system parameters have great influence on the buckling characteristics of FGM shells. Buckling behavior analysis based on Reddy’s high-order shear deformation theory has also been performed by Sun et al. [12] for FGM cylindrical shells subjected to an axial compression in the thermal environment. In recent years, there have been some numerical methods for analyzing the plates and shells. The multiquadric radial basis function (MQ) method was developed and applied by Ferreira [13] to discuss the effects of system parameters on the laminated composite plates. The method of discrete singular convolution (DSC) gives a fast and accurate solution of the mathematical physics and engineering problems. Then, Civalek and his coworkers made a number of remarkable studies using the DSC method; such investigations involved the analysis of composite conical shells and panels [14–16], plates on elastic foundations [17, 18], and so on.However, in the course of our study, we found that there are only few studies on the global bifurcations and multipulse chaotic dynamics for the CNT-reinforced composite plate. In order to eliminate or suppress large nonlinear vibrations and chaotic motions of the CNT-reinforced composite plate, we should deepen and complete the theoretical analysis on the CNT-reinforced composite plate model, discuss the complex dynamic behaviors, explore the existence conditions of the multipulse Shilnikov-type orbits, and analyze the impact of parameters on the system, so as to ensure the stability and controllability of the CNT-reinforced composite plate. The present work is therefore motivated to examine the global bifurcations and multipulse chaotic dynamics of the CNT-reinforced composite plate considered by Guo and Zhang [19]. The energy-phase method developed by Haller and Wiggins [20] is applied to study the multipulse homoclinic bifurcations and chaotic dynamics for the CNT-reinforced composite plate. In recent decades, two theories that are usually employed to detect the global bifurcations and multipulse Shilnikov-type chaotic dynamics in high-dimensional systems are the energy-phase method [20] and the extended Melnikov method [21]. Also, some researchers have utilized the developed theories to engineering applications, see, for example, [22–29].Based on the research by Guo and Zhang [19], the averaged equations are obtained for the case of 1 : 1 internal resonance, principal parametric resonance, and 1/2 subharmonic resonance of the CNT-reinforced composite plate. From the averaged equations, the explicit expressions of the 3-order normal form are derived employing the normal form theory. The energy-phase method is utilized to detect the existence of the multipulse Shilnikov-type homoclinic orbits in the CNT-reinforced composite plate. Homoclinic trees are constructed to illustrate the repeated bifurcations of multipulse solutions. In addition, numerical simulations are given to verify the analytical predictions.The organization of this paper is as follows. In Section2, we derive the averaged equations and normal form for the CNT-reinforced composite plate. According to the normal form derived in Section 2, the dynamics of the unperturbed system is discussed in Section 3. In Section 4, a detailed dynamic analysis of the perturbed system is presented, and the existence of multipulse homoclinic orbits is confirmed. Numerical simulations are given in Section 5 to show that the chaotic motions can occur in the CNT-reinforced composite plate. A short conclusion is reported in Section 6.


## 2. Formulation of the Problem

### 2.1. Equations of Motion and Perturbation Analysis

This paper focuses on studies of the multipulse homoclinic orbits and chaotic dynamics of a four-edge simply-supported CNT-reinforced composite rectangular plate subjected to the in-plane and transverse excitations. The model is shown in Figure1 [19]. The edge length and width of the plate in the x and y directions are, respectively, a and b, and the thickness is h˜. A Cartesian coordinate system is located in the middle surface of the CNT-reinforced composite rectangular plate. The displacements of an arbitrary point within the plate are u, v, and w in the x, y, and z directions, respectively. The in-plane excitation along the y direction at x=a is given by P=P˜0−P˜1cosΩ2t. The transverse excitation subjecting to the CNT-reinforced composite rectangular plate is of the form F=F0−F1cosΩ1t.Figure 1
The model of a CNT-reinforced composite rectangular plate.According to the research by Guo and Zhang [19], the axial and transverse displacement fields at any point for the CNT-reinforced composite plate are given as(1a)ux,y,t=u0x,y,t+zφxx,y,t−z343h˜2φx+∂w0∂x,(1b)vx,y,t=v0x,y,t+zφyx,y,t−z343h˜2φy+∂w0∂y,(1c)wx,y,t=w0x,y,t,where u0, v0, and w0 represent the displacements at the midplane of the plate in the x, y, and z directions. φx and φy denote the rotations of the transverse normal at the midplane about the x and y axes. Employing the von Kármán-type plate theory and in terms of the displacements, the strains εii=xx,yy and the curvatures γii=xy,yz,zx in the midplane can be expressed as(2)εxx=∂u∂x+12∂w∂x2,εyy=∂v∂y+12∂w∂y2,γxy=12∂u∂x+∂v∂y+∂w∂x∂w∂y,γyz=12∂v∂z+∂w∂y,γzx=12∂u∂z+∂w∂x.The force and moment resultants associated with the strains and curvatures of the plate constitutive equations are(3a)NxNyNxy=A11A12A16A12A22A26A16A26A66εx0εy0γxy0,(3b)MxMyMxy=D11D12D16D12D22D26D16D26D66κx0κy0κxy0+F11F12F16F12F22F26F16F26F66κx2κy2κxy2,(3c)PxPyPxy=F11F12F16F12F22F26F16F26F66κx0κy0κxy0+H11H12H16H21H22H26H16H26H66κx2κy2κxy2,(3d)QyQx=A44A45A45A55γyz0γxz0,(3e)RyRx=D44D45D45D55γyz0γxz0+F44F45F54F55γyz2γxz2,where the coupling matrix Bij and Eij is equal to zero in the stiffness matrix. Then, the nonlinear governing equations of motion for the CNT-reinforced composite plate are obtained as follows [19]:(4a)Nx,x+Nxy,y=I0u¨0+I˜1−43h˜2I˜3φ¨x−43h˜2I˜3∂w¨0∂x,(4b)Ny,y+Nxy,x=I0v¨0+I˜1−43h˜2I˜3φ¨y−43h˜2I˜3∂w¨0∂y,(4c)Ny,y∂w0∂y+Ny∂2w0∂y2+Nx,x∂w0∂x+Nx∂2w0∂x2+Nxy,x∂w0∂y+Nxy,y∂w0∂x+2Nxy∂2w0∂y∂x+43h˜2Px,xx+2Pxy,xy+Py,yy−4h˜2Rx,x+Ry,y−μw˙0+Qx,x+Qy,y+F=I0w¨0+43h˜2I˜3∂u¨0∂x+∂v¨0∂y−43h˜22I˜6∂2w¨0∂x2+∂2w¨0∂y2+43h˜2I˜4−43h˜2I˜6∂φ¨x∂x+∂φ¨y∂y,(4d)Mx,x+Mxy,y−43h˜2Px,x−43h˜2Pxy,y−Qx−13h˜2Rx=I˜1−43h˜2I˜3u¨0+I˜2−83h˜2I˜4+43h˜22I˜6φ¨x−43h˜2I˜4−43h˜2I˜6∂w¨0∂x,(4e)My,y+Mxy,x−43h˜2Py,y−43h˜2Pxy,x−Qy−13h˜2Ry=I˜1−43h˜2I˜3v¨0+I˜2−83h˜2I˜4+43h˜22I˜6φ¨y−43h˜2I˜4−43h˜2I˜6∂w¨0∂y.Substituting the stress resultants of equations (1a)–(1c), (2), and (3a)–(3e) into equations (4a)–(4e), the governing equations of motion for the CNT-reinforced composite plate in terms of generalized displacements are derived as follows [19]:(5a)A11∂2u0∂x2+A66∂2u0∂y2+A26∂2v0∂y2+A16∂2v0∂x2+A12+A66∂2u0∂x∂y+2A16∂2u0∂x∂y+A11∂w0∂x∂2w0∂x2+A66∂w0∂x∂2w0∂y2+A16∂w0∂y∂2w0∂x2+A26∂w0∂y∂2w0∂y2+2A16∂w0∂x∂2w0∂x∂y+A12+A66∂w0∂y∂2w0∂x∂y=I0u¨0+I˜1−43h˜2I˜3φ¨x−43h˜2I˜3∂w¨0∂x,(5b)A16∂2u0∂x2+A26∂2u0∂y2+A22∂2v0∂y2+A66∂2v0∂x2+A12+A66∂2u0∂x∂y+2A26∂2v0∂x∂y+A16∂w0∂x∂2w0∂x2+A26∂w0∂x∂2w0∂y2+A66∂w0∂y∂2w0∂x2+A22∂w0∂y∂2w0∂y2+2A26∂w0∂y∂2w0∂x∂y+A12+A66∂w0∂x∂2w0∂x∂y=I0v¨0+I˜1−43h˜2I˜3φ¨y−43h˜2I˜3∂w¨0∂y,(5c)−c12H11∂4w0∂x4−c12H22∂4w0∂y4+A66+32A11∂2w0∂y2+32A22∂2w0∂x2+3A26∂2w0∂x∂y∂w0∂x2−2c12H12+2H66∂4w0∂x2∂y2+A66+12A11∂2w0∂x2+32A22∂2w0∂y2+3A26∂2w0∂x∂y∂w0∂y2−4c12H26∂4w0∂x∂y3−4c12H16∂4w0∂x3∂y+3A16∂w0∂x∂w0∂y∂2w0∂x2+3A26∂w0∂x∂w0∂y∂2w0∂y2+2A12+4A66∂2w0∂x∂y∂w0∂x∂w0∂y+A55−c22F55+2c2D55∂2w0∂x2+A44−c22F44+2c2D44∂2w0∂y2+2A45−c22F45+c2D45∂2w0∂x∂y+A26∂u0∂y+A12∂u0∂x+A22∂v0∂y+A26∂v0∂x∂2w0∂y2+A66∂u0∂y+A16∂u0∂x+A26∂v0∂y+A66∂v0∂x∂2w0∂x∂y+A11∂u0∂x+A16∂u0∂y+A12∂v0∂y+A16∂v0∂x∂2w0∂x2+A11+A16∂2u0∂x2+A26+A22∂2v0∂y2+A16∂2v0∂x2+A66∂2u0∂y2+A12+A66∂2v0∂x∂y+A16+A66∂2u0∂x∂y∂w0∂x+A16+A66∂2u0∂x2+A22+A26∂2v0∂y2+A66∂2v0∂x2+A26∂2u0∂y2A26∂2v0∂x∂y+A12+A66∂2u0∂x∂y∂w0∂y+c1F11−c12H11∂3φx∂x3+c1F26−c12H26∂3φx∂y3+c13F16−c1H16∂3φx∂x2∂y+c1F21+2F66−c1H21−c1H66∂3φx∂x∂y2+A45−2c1D45+c12F45∂φx∂y+A55−2c1D55+c12F55∂φx∂x+c1F22−c1H22∂3φy∂y3+c1F16−c1H16∂3φy∂x3+c1F12+2F66−c1H21−c1H66∂3φy∂x2∂y+c13F26−c1H26∂3φy∂x∂y2+A45−c2D45+c22F45∂φy∂x+A44−c2D44+c22F44∂φy∂y+FcosΩ1t−μw˙0=I0∂2w¨0∂t2−c12I˜6∂2w¨0∂x2+∂2w¨0∂y2+c1I˜3∂u¨0∂x+∂v¨0∂y,(5d)D11−2c1F11+c12H11∂2φx∂x2+D66−2c1F66+c12H66∂2φx∂y2+D16−2c1F16+c12H16∂2φy∂x2+D22−2c1F22+c12H22∂2φy∂y2+D21+D16−2c1F16−2c12H12∂2φx∂x∂y+−A55+2c2D55−c22F55φx+−A45+2c2D45−c22F45φy+c12H11−c1F11∂3w0∂x3+c12H26−c1F26∂3w0∂y3−c12H16+c1F16∂2w0∂x2∂y+4c12H66−c12H12−2c1F66−c1F12∂3w0∂x∂y2+−A55+2c2D55−c22F55∂w0∂x+−A45+2c2D45−c22F45∂w0∂y=I˜2−43h˜2I˜4u¨0−43h˜2I˜5−43h˜2I˜7∂w¨0∂x+I˜3−83h˜2I˜5+169h˜4I˜7φ¨x,(5e)D16−2c1F16+c12H16∂2φx∂x2+D26−2c1F26+c12H26∂2φx∂y2+D66−2c1F66+c12H66∂2φy∂x2+D22−2c1F22+c12H22∂2φy∂y2+D22+D26−4c1F16+2c12H12∂2φy∂x∂y+D21+D66−2c1F21−2c1F66+c12H21+c12H66∂2φx∂x∂y−A45−2c2D45+c22F45φx−A55−2c2D55+c22F55φy+c12H16−c1F16∂3w0∂x3+c12H22−c1F22∂3w0∂y3+2c12H66−c12H21−c1F66−c1F21∂3w0∂x2∂y−c12H26+c1F26∂3w0∂x∂y2−A45−c2D45+16h˜4F45∂w0∂x−A55−8h˜2D55+16h˜4F55∂w0∂y=I˜2−c1I˜4v¨0−c1I˜5−c1I˜7∂w¨0∂y+I˜3−2c1I˜5+c12I˜7φ¨y,where μ is the damping coefficient, and all the other coefficients can be found in [19]. The associated boundary conditions for the simply-supported CNT composite plate can be written as(6a)v=w=0,φy=0,Mxx=Nxy=0,x=0,a,(6b)v=w=0,φx=0,Myy=Nxy=0,y=0,b,(6c)∫0h˜Nxxx=0dz=−∫0h˜F0−F1cosΩ1tdz.Also, the following nondimensional parameters are introduced:(7)u¯0=u0a,v¯0=v0b,w¯0=w0h˜,φ¯x=φx,φ¯y=φy,x¯=xa,y¯=yb,I¯i=1Li+1ρI˜i,F¯=ab7/2Eh˜7F,P¯=b2El3P,μ¯=1ρEa2b2π2h˜4μ,T¯=1LEρT,A¯ij=ab1/2Eh˜2Aij,D¯ij=ab11/2Eh˜4Dij,F¯ij=ab1/2Eh˜4Dij,H¯ij=ab1/2Eh˜4Hij,i=1,2,3,4,5,6;j=1,2,3,4,5,6.In this paper, our research is focused on the multipulse global bifurcations and chaotic dynamics of the CNT composite plate in its first two modes. Hence, we expressw in the following form:(8)w0=w1sin3πxasinπyb+w2sinπxasin3πyb,where w1 and w2 are the amplitudes of two modes. The other variables and transverse excitation are given as(9a)u0=u1sin3πxacosπyb+u2sinπxacos3πyb,(9b)v0=v1cos3πxasinπyb+v2cosπxasin3πyb,(9c)φx=φ1cos3πxasinπyb+φ2cosπxasin3πyb,(9d)φy=φ3sin3πxacosπyb+φ4sinπxacos3πyb,(9e)F=F1sin3πxasinπyb+F2sinπxasin3πyb.Then, substituting these expressions into equations (5a)–(5e) and applying the Galerkin integration procedure, we can obtain the displacements u0, v0, φx, and φy with respect to w0. Consequently, the nonlinear ordinary differential equations of this system in terms of transverse displacements are derived. According to the results of convergence studies given by Hao et al. [30], the dimensionless equations of motion for the CNT-reinforced composite plate in the first two modes are derived as [19](10a)w¨1+μw˙1+ω12w1+a1w1cosΩ2t+a2w13+a3w23+a4w12w2+a5w1w22+a6w1w2+a7w12+a8w22=f1cosΩ1t,(10b)w¨2+μw˙2+ω22w2+b1w2cosΩ2t+b2w13+b3w23+b4w12w2+b5w1w22+b6w1w2+b7w12+b8w22=f2cosΩ1t,where all the coefficients can be found in [19]. We consider the case of 1 : 1 internal resonance, principal parametric resonance, and 1/2 subharmonic resonance of the CNT-reinforced composite plate, for which the resonant relations are listed as follows:(11)ω1=12Ω1+ε˜σ1,ω2=12Ω2+ε˜σ2,Ω1=Ω2=Ω,where σ1 and σ2 are two detuning parameters. For convenience, we let Ω=2 in the following analysis. The uniform solutions of equations (10a) and (10b) take the following form:(12)wnt,ε˜=xn0T0,T1+ε˜xn1T0,T1+⋯,n=1,2,where T0=t and T1=ε˜t. We employ the method of multiple scales to equations (10a) and (10b) and obtain the averaged equations as follows:(13a)x˙1=−12μx1−σ1−14α1x2−32α2x2x12+x22−32α3x4x32+x42+12α4x4x12−x22−α4x4x12+x22+12α5x2x32−x42−α5x2x32+x42−α4x1x2x3−α5x1x3x4,(13b)x˙2=−12μx2+σ1+14α1x1+32α2x1x12+x22+32α3x3x32+x42+12α4x3x12−x22+α4x3x12+x22+12α5x1x32−x42+α5x1x32+x42+α4x1x2x4+α5x2x3x4,(13c)x˙3=−12μx3−σ2−f0x4−32β2x2x12+x22−32β3x4x32+x42+12β4x4x12−x22−β4x4x12+x22+12β5x2x32−x42−β5x2x32+x42−β4x1x2x3−β5x1x3x4,(13d)x˙4=−12μx4+σ2+f0x3+32β2x1x12+x22+32β3x3x32+x42+12β4x3x12−x22+β4x3x12+x22+12β5x1x32−x42+β5x1x32+x42+β4x1x2x4+β5x2x3x4,where f0=1/4β1.


### 2.2. Computation of Normal Form

In order to study the multipulse chaotic dynamics of the CNT-reinforced composite plate, the first is to reduce equations (13a)–(13d) to a simpler normal form. It is clear that equations (13a)–(13d) possess a trivial zero solution x1,x2,x3,x4=0,0,0,0 at which the Jacobian matrix can be written as(14)J=−12μ−σ1+14α100σ1+14α1−12μ0000−12μ−σ2−f000σ2+f0−12μ.The characteristic polynomial of Jacobian matrix (14) is(15)λI−J=λ2+μλ+14μ2+σ12−116α12λ2+μλ+14μ2+σ22−f02.It is easy to see that whenμ=0, α1=−4σ1, and σ22−f02>0 are satisfied simultaneously, systems (13a)–(13d) have a pair of pure imaginary eigenvalues and one non-semisimple double zero:(16)λ1,2=0,λ3,4=±iσ22−f02.Settingα1=2 and σ=σ1+1/4α1 and considering μ, σ, and f0 as the perturbation parameters, then equations (13a)–(13d) without the perturbation parameters become(17a)x˙1=x2−32α2x2x12+x22−32α3x4x32+x42+12α4x4x12−x22−α4x4x12+x22+12α5x2x32−x42−α5x2x32+x42−α4x1x2x3−α5x1x3x4,(17b)x˙2=32α2x1x12+x22+32α3x3x32+x42+12α4x3x12−x22+α4x3x12+x22+12α5x1x32−x42+α5x1x32+x42+α4x1x2x4+α5x2x3x4,(17c)x˙3=−σ2x4−32β2x2x12+x22−32β3x4x32+x42+12β4x4x12−x22−β4x4x12+x22+12β5x2x32−x42−β5x2x32+x42−β4x1x2x3−β5x1x3x4,(17d)x˙4=σ2x3+32β2x1x12+x22+32β3x3x32+x42+12β4x3x12−x22+β4x3x12+x22+12β5x1x32−x42+β5x1x32+x42+β4x1x2x4+β5x2x3x4.For the trivial zero solutionx1,x2,x3,x4=0,0,0,0, the Jacobian matrix of equations (17a)–(17d) is evaluated as(18)A=01000000000−σ200σ20.Executing the Maple program designed by Zhang et al. [31] leads to the following 3-order normal form of equations (17a)–(17d):(19a)y˙1=y2,(19b)y˙2=32α2y13+α5y1y32+y42,(19c)y˙3=−σ2y4−β4y12y4−32β3y4y32+y42,(19d)y˙4=σ2y3+β4y12y3+32β3y3y32+y42.Accordingly, the normal form of equations (13a)–(13d) with the perturbation parameters is obtained as(20a)y˙1=−12μy1+1−σy2,(20b)y˙2=σy1−12μy2+32α2y13+α5y1y32+y42,(20c)y˙3=−12μy3−σ2−f0y4−β4y12y4−32β3y4y32+y42,(20d)y˙4=σ2+f0y3−12μy4+β4y12y3+32β3y3y32+y42.Letting(21)y3=Icosϕ,y4=Isinϕ,we have(22a)y˙1=−12μy1+1−σy2,(22b)y˙2=σy1−12μy2+32α2y13+α5y1I2,(22c)I˙=−12μI+f0Isin2ϕ,(22d)Iϕ˙=σ2I+β4y12I+32β3I3+f0Icos2ϕ.We then construct the following transformation to find the unfolding of equations (22a)–(22d):(23)y1y2=α5β41−σ012μ1z1z2.Substituting equation (23) into equations (22a)–(22d) and canceling the nonlinear terms including parameter σ, the unfolding of equations (22a)–(22d) can be derived as(24a)z˙1=z2,(24b)z˙2=−μ¯z1−μu2+ηz13+α5z1I2,(24c)I˙=−12μI+f0Isin2ϕ,(24d)Iϕ˙=σ2I+α5z12I+32β3I3+f0Icos2ϕ,where μ¯=1/4μ2−σ1−σ and η=3α2α5/2β4. Introduce the scale transformations(25)μ⟶εμ,f0⟶εf0.Then, equation (24a)–(24d) can be expressed as the Hamilton form with perturbations:(26a)z˙1=∂H∂z2+εgz1=z2,(26b)z˙2=−∂H∂z1+εgz2=−μ¯z1+ηz13+α5z1I2−εμz2,(26c)I˙=∂H∂ϕ+εgI+εf0Isin2ϕ=−12εμI+εf0Isin2ϕ,(26d)Iϕ˙=−∂H∂I+εgϕ+εf0Icos2ϕ=σ2I+α5z12I+32β3I3+εf0Icos2ϕ.The Hamiltonian functionH takes the form(27)H=12z22+12μ¯z12−14ηz14−12α5z12I2−12σ2I2−38β3I4,gz1=0,gz2=−μz2,gI=−12μI,gϕ=0.


## 2.1. Equations of Motion and Perturbation Analysis

This paper focuses on studies of the multipulse homoclinic orbits and chaotic dynamics of a four-edge simply-supported CNT-reinforced composite rectangular plate subjected to the in-plane and transverse excitations. The model is shown in Figure1 [19]. The edge length and width of the plate in the x and y directions are, respectively, a and b, and the thickness is h˜. A Cartesian coordinate system is located in the middle surface of the CNT-reinforced composite rectangular plate. The displacements of an arbitrary point within the plate are u, v, and w in the x, y, and z directions, respectively. The in-plane excitation along the y direction at x=a is given by P=P˜0−P˜1cosΩ2t. The transverse excitation subjecting to the CNT-reinforced composite rectangular plate is of the form F=F0−F1cosΩ1t.Figure 1
The model of a CNT-reinforced composite rectangular plate.According to the research by Guo and Zhang [19], the axial and transverse displacement fields at any point for the CNT-reinforced composite plate are given as(1a)ux,y,t=u0x,y,t+zφxx,y,t−z343h˜2φx+∂w0∂x,(1b)vx,y,t=v0x,y,t+zφyx,y,t−z343h˜2φy+∂w0∂y,(1c)wx,y,t=w0x,y,t,where u0, v0, and w0 represent the displacements at the midplane of the plate in the x, y, and z directions. φx and φy denote the rotations of the transverse normal at the midplane about the x and y axes. Employing the von Kármán-type plate theory and in terms of the displacements, the strains εii=xx,yy and the curvatures γii=xy,yz,zx in the midplane can be expressed as(2)εxx=∂u∂x+12∂w∂x2,εyy=∂v∂y+12∂w∂y2,γxy=12∂u∂x+∂v∂y+∂w∂x∂w∂y,γyz=12∂v∂z+∂w∂y,γzx=12∂u∂z+∂w∂x.The force and moment resultants associated with the strains and curvatures of the plate constitutive equations are(3a)NxNyNxy=A11A12A16A12A22A26A16A26A66εx0εy0γxy0,(3b)MxMyMxy=D11D12D16D12D22D26D16D26D66κx0κy0κxy0+F11F12F16F12F22F26F16F26F66κx2κy2κxy2,(3c)PxPyPxy=F11F12F16F12F22F26F16F26F66κx0κy0κxy0+H11H12H16H21H22H26H16H26H66κx2κy2κxy2,(3d)QyQx=A44A45A45A55γyz0γxz0,(3e)RyRx=D44D45D45D55γyz0γxz0+F44F45F54F55γyz2γxz2,where the coupling matrix Bij and Eij is equal to zero in the stiffness matrix. Then, the nonlinear governing equations of motion for the CNT-reinforced composite plate are obtained as follows [19]:(4a)Nx,x+Nxy,y=I0u¨0+I˜1−43h˜2I˜3φ¨x−43h˜2I˜3∂w¨0∂x,(4b)Ny,y+Nxy,x=I0v¨0+I˜1−43h˜2I˜3φ¨y−43h˜2I˜3∂w¨0∂y,(4c)Ny,y∂w0∂y+Ny∂2w0∂y2+Nx,x∂w0∂x+Nx∂2w0∂x2+Nxy,x∂w0∂y+Nxy,y∂w0∂x+2Nxy∂2w0∂y∂x+43h˜2Px,xx+2Pxy,xy+Py,yy−4h˜2Rx,x+Ry,y−μw˙0+Qx,x+Qy,y+F=I0w¨0+43h˜2I˜3∂u¨0∂x+∂v¨0∂y−43h˜22I˜6∂2w¨0∂x2+∂2w¨0∂y2+43h˜2I˜4−43h˜2I˜6∂φ¨x∂x+∂φ¨y∂y,(4d)Mx,x+Mxy,y−43h˜2Px,x−43h˜2Pxy,y−Qx−13h˜2Rx=I˜1−43h˜2I˜3u¨0+I˜2−83h˜2I˜4+43h˜22I˜6φ¨x−43h˜2I˜4−43h˜2I˜6∂w¨0∂x,(4e)My,y+Mxy,x−43h˜2Py,y−43h˜2Pxy,x−Qy−13h˜2Ry=I˜1−43h˜2I˜3v¨0+I˜2−83h˜2I˜4+43h˜22I˜6φ¨y−43h˜2I˜4−43h˜2I˜6∂w¨0∂y.Substituting the stress resultants of equations (1a)–(1c), (2), and (3a)–(3e) into equations (4a)–(4e), the governing equations of motion for the CNT-reinforced composite plate in terms of generalized displacements are derived as follows [19]:(5a)A11∂2u0∂x2+A66∂2u0∂y2+A26∂2v0∂y2+A16∂2v0∂x2+A12+A66∂2u0∂x∂y+2A16∂2u0∂x∂y+A11∂w0∂x∂2w0∂x2+A66∂w0∂x∂2w0∂y2+A16∂w0∂y∂2w0∂x2+A26∂w0∂y∂2w0∂y2+2A16∂w0∂x∂2w0∂x∂y+A12+A66∂w0∂y∂2w0∂x∂y=I0u¨0+I˜1−43h˜2I˜3φ¨x−43h˜2I˜3∂w¨0∂x,(5b)A16∂2u0∂x2+A26∂2u0∂y2+A22∂2v0∂y2+A66∂2v0∂x2+A12+A66∂2u0∂x∂y+2A26∂2v0∂x∂y+A16∂w0∂x∂2w0∂x2+A26∂w0∂x∂2w0∂y2+A66∂w0∂y∂2w0∂x2+A22∂w0∂y∂2w0∂y2+2A26∂w0∂y∂2w0∂x∂y+A12+A66∂w0∂x∂2w0∂x∂y=I0v¨0+I˜1−43h˜2I˜3φ¨y−43h˜2I˜3∂w¨0∂y,(5c)−c12H11∂4w0∂x4−c12H22∂4w0∂y4+A66+32A11∂2w0∂y2+32A22∂2w0∂x2+3A26∂2w0∂x∂y∂w0∂x2−2c12H12+2H66∂4w0∂x2∂y2+A66+12A11∂2w0∂x2+32A22∂2w0∂y2+3A26∂2w0∂x∂y∂w0∂y2−4c12H26∂4w0∂x∂y3−4c12H16∂4w0∂x3∂y+3A16∂w0∂x∂w0∂y∂2w0∂x2+3A26∂w0∂x∂w0∂y∂2w0∂y2+2A12+4A66∂2w0∂x∂y∂w0∂x∂w0∂y+A55−c22F55+2c2D55∂2w0∂x2+A44−c22F44+2c2D44∂2w0∂y2+2A45−c22F45+c2D45∂2w0∂x∂y+A26∂u0∂y+A12∂u0∂x+A22∂v0∂y+A26∂v0∂x∂2w0∂y2+A66∂u0∂y+A16∂u0∂x+A26∂v0∂y+A66∂v0∂x∂2w0∂x∂y+A11∂u0∂x+A16∂u0∂y+A12∂v0∂y+A16∂v0∂x∂2w0∂x2+A11+A16∂2u0∂x2+A26+A22∂2v0∂y2+A16∂2v0∂x2+A66∂2u0∂y2+A12+A66∂2v0∂x∂y+A16+A66∂2u0∂x∂y∂w0∂x+A16+A66∂2u0∂x2+A22+A26∂2v0∂y2+A66∂2v0∂x2+A26∂2u0∂y2A26∂2v0∂x∂y+A12+A66∂2u0∂x∂y∂w0∂y+c1F11−c12H11∂3φx∂x3+c1F26−c12H26∂3φx∂y3+c13F16−c1H16∂3φx∂x2∂y+c1F21+2F66−c1H21−c1H66∂3φx∂x∂y2+A45−2c1D45+c12F45∂φx∂y+A55−2c1D55+c12F55∂φx∂x+c1F22−c1H22∂3φy∂y3+c1F16−c1H16∂3φy∂x3+c1F12+2F66−c1H21−c1H66∂3φy∂x2∂y+c13F26−c1H26∂3φy∂x∂y2+A45−c2D45+c22F45∂φy∂x+A44−c2D44+c22F44∂φy∂y+FcosΩ1t−μw˙0=I0∂2w¨0∂t2−c12I˜6∂2w¨0∂x2+∂2w¨0∂y2+c1I˜3∂u¨0∂x+∂v¨0∂y,(5d)D11−2c1F11+c12H11∂2φx∂x2+D66−2c1F66+c12H66∂2φx∂y2+D16−2c1F16+c12H16∂2φy∂x2+D22−2c1F22+c12H22∂2φy∂y2+D21+D16−2c1F16−2c12H12∂2φx∂x∂y+−A55+2c2D55−c22F55φx+−A45+2c2D45−c22F45φy+c12H11−c1F11∂3w0∂x3+c12H26−c1F26∂3w0∂y3−c12H16+c1F16∂2w0∂x2∂y+4c12H66−c12H12−2c1F66−c1F12∂3w0∂x∂y2+−A55+2c2D55−c22F55∂w0∂x+−A45+2c2D45−c22F45∂w0∂y=I˜2−43h˜2I˜4u¨0−43h˜2I˜5−43h˜2I˜7∂w¨0∂x+I˜3−83h˜2I˜5+169h˜4I˜7φ¨x,(5e)D16−2c1F16+c12H16∂2φx∂x2+D26−2c1F26+c12H26∂2φx∂y2+D66−2c1F66+c12H66∂2φy∂x2+D22−2c1F22+c12H22∂2φy∂y2+D22+D26−4c1F16+2c12H12∂2φy∂x∂y+D21+D66−2c1F21−2c1F66+c12H21+c12H66∂2φx∂x∂y−A45−2c2D45+c22F45φx−A55−2c2D55+c22F55φy+c12H16−c1F16∂3w0∂x3+c12H22−c1F22∂3w0∂y3+2c12H66−c12H21−c1F66−c1F21∂3w0∂x2∂y−c12H26+c1F26∂3w0∂x∂y2−A45−c2D45+16h˜4F45∂w0∂x−A55−8h˜2D55+16h˜4F55∂w0∂y=I˜2−c1I˜4v¨0−c1I˜5−c1I˜7∂w¨0∂y+I˜3−2c1I˜5+c12I˜7φ¨y,where μ is the damping coefficient, and all the other coefficients can be found in [19]. The associated boundary conditions for the simply-supported CNT composite plate can be written as(6a)v=w=0,φy=0,Mxx=Nxy=0,x=0,a,(6b)v=w=0,φx=0,Myy=Nxy=0,y=0,b,(6c)∫0h˜Nxxx=0dz=−∫0h˜F0−F1cosΩ1tdz.Also, the following nondimensional parameters are introduced:(7)u¯0=u0a,v¯0=v0b,w¯0=w0h˜,φ¯x=φx,φ¯y=φy,x¯=xa,y¯=yb,I¯i=1Li+1ρI˜i,F¯=ab7/2Eh˜7F,P¯=b2El3P,μ¯=1ρEa2b2π2h˜4μ,T¯=1LEρT,A¯ij=ab1/2Eh˜2Aij,D¯ij=ab11/2Eh˜4Dij,F¯ij=ab1/2Eh˜4Dij,H¯ij=ab1/2Eh˜4Hij,i=1,2,3,4,5,6;j=1,2,3,4,5,6.In this paper, our research is focused on the multipulse global bifurcations and chaotic dynamics of the CNT composite plate in its first two modes. Hence, we expressw in the following form:(8)w0=w1sin3πxasinπyb+w2sinπxasin3πyb,where w1 and w2 are the amplitudes of two modes. The other variables and transverse excitation are given as(9a)u0=u1sin3πxacosπyb+u2sinπxacos3πyb,(9b)v0=v1cos3πxasinπyb+v2cosπxasin3πyb,(9c)φx=φ1cos3πxasinπyb+φ2cosπxasin3πyb,(9d)φy=φ3sin3πxacosπyb+φ4sinπxacos3πyb,(9e)F=F1sin3πxasinπyb+F2sinπxasin3πyb.Then, substituting these expressions into equations (5a)–(5e) and applying the Galerkin integration procedure, we can obtain the displacements u0, v0, φx, and φy with respect to w0. Consequently, the nonlinear ordinary differential equations of this system in terms of transverse displacements are derived. According to the results of convergence studies given by Hao et al. [30], the dimensionless equations of motion for the CNT-reinforced composite plate in the first two modes are derived as [19](10a)w¨1+μw˙1+ω12w1+a1w1cosΩ2t+a2w13+a3w23+a4w12w2+a5w1w22+a6w1w2+a7w12+a8w22=f1cosΩ1t,(10b)w¨2+μw˙2+ω22w2+b1w2cosΩ2t+b2w13+b3w23+b4w12w2+b5w1w22+b6w1w2+b7w12+b8w22=f2cosΩ1t,where all the coefficients can be found in [19]. We consider the case of 1 : 1 internal resonance, principal parametric resonance, and 1/2 subharmonic resonance of the CNT-reinforced composite plate, for which the resonant relations are listed as follows:(11)ω1=12Ω1+ε˜σ1,ω2=12Ω2+ε˜σ2,Ω1=Ω2=Ω,where σ1 and σ2 are two detuning parameters. For convenience, we let Ω=2 in the following analysis. The uniform solutions of equations (10a) and (10b) take the following form:(12)wnt,ε˜=xn0T0,T1+ε˜xn1T0,T1+⋯,n=1,2,where T0=t and T1=ε˜t. We employ the method of multiple scales to equations (10a) and (10b) and obtain the averaged equations as follows:(13a)x˙1=−12μx1−σ1−14α1x2−32α2x2x12+x22−32α3x4x32+x42+12α4x4x12−x22−α4x4x12+x22+12α5x2x32−x42−α5x2x32+x42−α4x1x2x3−α5x1x3x4,(13b)x˙2=−12μx2+σ1+14α1x1+32α2x1x12+x22+32α3x3x32+x42+12α4x3x12−x22+α4x3x12+x22+12α5x1x32−x42+α5x1x32+x42+α4x1x2x4+α5x2x3x4,(13c)x˙3=−12μx3−σ2−f0x4−32β2x2x12+x22−32β3x4x32+x42+12β4x4x12−x22−β4x4x12+x22+12β5x2x32−x42−β5x2x32+x42−β4x1x2x3−β5x1x3x4,(13d)x˙4=−12μx4+σ2+f0x3+32β2x1x12+x22+32β3x3x32+x42+12β4x3x12−x22+β4x3x12+x22+12β5x1x32−x42+β5x1x32+x42+β4x1x2x4+β5x2x3x4,where f0=1/4β1.


## 2.2. Computation of Normal Form

In order to study the multipulse chaotic dynamics of the CNT-reinforced composite plate, the first is to reduce equations (13a)–(13d) to a simpler normal form. It is clear that equations (13a)–(13d) possess a trivial zero solution x1,x2,x3,x4=0,0,0,0 at which the Jacobian matrix can be written as(14)J=−12μ−σ1+14α100σ1+14α1−12μ0000−12μ−σ2−f000σ2+f0−12μ.The characteristic polynomial of Jacobian matrix (14) is(15)λI−J=λ2+μλ+14μ2+σ12−116α12λ2+μλ+14μ2+σ22−f02.It is easy to see that whenμ=0, α1=−4σ1, and σ22−f02>0 are satisfied simultaneously, systems (13a)–(13d) have a pair of pure imaginary eigenvalues and one non-semisimple double zero:(16)λ1,2=0,λ3,4=±iσ22−f02.Settingα1=2 and σ=σ1+1/4α1 and considering μ, σ, and f0 as the perturbation parameters, then equations (13a)–(13d) without the perturbation parameters become(17a)x˙1=x2−32α2x2x12+x22−32α3x4x32+x42+12α4x4x12−x22−α4x4x12+x22+12α5x2x32−x42−α5x2x32+x42−α4x1x2x3−α5x1x3x4,(17b)x˙2=32α2x1x12+x22+32α3x3x32+x42+12α4x3x12−x22+α4x3x12+x22+12α5x1x32−x42+α5x1x32+x42+α4x1x2x4+α5x2x3x4,(17c)x˙3=−σ2x4−32β2x2x12+x22−32β3x4x32+x42+12β4x4x12−x22−β4x4x12+x22+12β5x2x32−x42−β5x2x32+x42−β4x1x2x3−β5x1x3x4,(17d)x˙4=σ2x3+32β2x1x12+x22+32β3x3x32+x42+12β4x3x12−x22+β4x3x12+x22+12β5x1x32−x42+β5x1x32+x42+β4x1x2x4+β5x2x3x4.For the trivial zero solutionx1,x2,x3,x4=0,0,0,0, the Jacobian matrix of equations (17a)–(17d) is evaluated as(18)A=01000000000−σ200σ20.Executing the Maple program designed by Zhang et al. [31] leads to the following 3-order normal form of equations (17a)–(17d):(19a)y˙1=y2,(19b)y˙2=32α2y13+α5y1y32+y42,(19c)y˙3=−σ2y4−β4y12y4−32β3y4y32+y42,(19d)y˙4=σ2y3+β4y12y3+32β3y3y32+y42.Accordingly, the normal form of equations (13a)–(13d) with the perturbation parameters is obtained as(20a)y˙1=−12μy1+1−σy2,(20b)y˙2=σy1−12μy2+32α2y13+α5y1y32+y42,(20c)y˙3=−12μy3−σ2−f0y4−β4y12y4−32β3y4y32+y42,(20d)y˙4=σ2+f0y3−12μy4+β4y12y3+32β3y3y32+y42.Letting(21)y3=Icosϕ,y4=Isinϕ,we have(22a)y˙1=−12μy1+1−σy2,(22b)y˙2=σy1−12μy2+32α2y13+α5y1I2,(22c)I˙=−12μI+f0Isin2ϕ,(22d)Iϕ˙=σ2I+β4y12I+32β3I3+f0Icos2ϕ.We then construct the following transformation to find the unfolding of equations (22a)–(22d):(23)y1y2=α5β41−σ012μ1z1z2.Substituting equation (23) into equations (22a)–(22d) and canceling the nonlinear terms including parameter σ, the unfolding of equations (22a)–(22d) can be derived as(24a)z˙1=z2,(24b)z˙2=−μ¯z1−μu2+ηz13+α5z1I2,(24c)I˙=−12μI+f0Isin2ϕ,(24d)Iϕ˙=σ2I+α5z12I+32β3I3+f0Icos2ϕ,where μ¯=1/4μ2−σ1−σ and η=3α2α5/2β4. Introduce the scale transformations(25)μ⟶εμ,f0⟶εf0.Then, equation (24a)–(24d) can be expressed as the Hamilton form with perturbations:(26a)z˙1=∂H∂z2+εgz1=z2,(26b)z˙2=−∂H∂z1+εgz2=−μ¯z1+ηz13+α5z1I2−εμz2,(26c)I˙=∂H∂ϕ+εgI+εf0Isin2ϕ=−12εμI+εf0Isin2ϕ,(26d)Iϕ˙=−∂H∂I+εgϕ+εf0Icos2ϕ=σ2I+α5z12I+32β3I3+εf0Icos2ϕ.The Hamiltonian functionH takes the form(27)H=12z22+12μ¯z12−14ηz14−12α5z12I2−12σ2I2−38β3I4,gz1=0,gz2=−μz2,gI=−12μI,gϕ=0.


## 3. Dynamics of the Unperturbed System

Settingε=0 in equations (26a)–(26d) results in the completely integrable equations defined as the unperturbed system. Hence, we will now study the nonlinear dynamics of the following unperturbed system:(28a)z˙1=z2,(28b)z˙2=−μ¯z1+ηz13+α5z1I2,(28c)I˙=0,(28d)Iϕ˙=σ2I+α5z12I+32β3I3.Consider the first two equations of systems (28a)–(28d):(29a)z˙1=z2,(29b)z˙2=−μ¯z1+ηz13+α5z1I2.Note that homoclinic bifurcations occur in systems (29a) and (29b) for η<0. Also, when μ¯−α5I2>0, we can obtain that systems (29a) and (29b) have only the trivial zero solution z1,z2=0,0 being a center. On a curve expressed by μ¯=α5I2, that is,(30)I1=1/4μ2−σ1−σα51/2,three solutions bifurcate from the trivial zero solution through a pitchfork bifurcation denoted by q0I=0,0 and q±I=B,0, respectively, where(31)B=±1/4μ2−σ1−σ−α5I2η1/2.For allI∈I1,+∞, systems (29a) and (29b) possess one hyperbolic saddle point q0I connected to itself by a pair of homoclinic orbits, z±hT1,I, that is, limT1⟶±∞z±hT1,I=q0I. Therefore, the set defined by(32)M=z,I,ϕz=q0I,I>I1,0≤ϕ≤2π,is a two-dimensional invariant manifold in the full four-dimensional phase space. Based on the analysis given by Haller and Wiggins [20], we know that the two-dimensional invariant manifold M is normally hyperbolic and contains three-dimensional stable and unstable manifolds. Let WsM and WuM, respectively, denote the stable and unstable manifolds. The existence of the homoclinic orbit in systems (29a) and (29b) to q0I=0,0 means that WsM and WuM intersect nontransversally along a three-dimensional homoclinic manifold Γ described as follows:(33)Γ=z,I,ϕz=z±hT1,I,I>I1,ϕ=1I∫0T1DIHz±hs,I,Ids+ϕ0,where DIHz,I=−∂H/∂I. Letting ξ1=−μ¯+α5I2, δ1=−η, systems (29a) and (29b) can be rewritten as(34a)z˙1=z2,(34b)z˙2=ξ1z1−δ1z13.The Hamiltonian for Hamilton systems (34a) and (34b) is(35)H¯z1,z2=12z22−12ξ1z12+14δ1z14.WhenH¯=0, there is a homoclinic loop Γ0 formed by one saddle point q0I and a pair of homoclinic orbits z±T1. Using equations (34a) and (34b) and (35), the analytical expressions for the homoclinic orbits of (34a) and (34b) are then obtained as(36a)z1T1=±2ξ1δ1sechξ1T1,(36b)z2T1=∓2ξ1δ1tanhξ1T1sechξ1T1.We now consider the dynamics of the unperturbed system of (26a) and (26d) restricted to M given by(37a)I˙=0,(37b)Iϕ˙=σ2I+α5Iq02I+32β3I3=DIHq0I,I,I>I1.If the conditionDIHq0I,I≠0 holds, I= constant is a periodic orbit, and if DIHq0I,I=0, I= constant is a circle of fixed points. A value of I∈I1,+∞ for which DIHq0I,I=0 is referred to as a resonant value I, and these fixed points are identified as resonant fixed points. A resonant value is denoted by Ir such that(38)DIHq0I,I=σ2I+32β3I3=0.Thus, the resonant valueIr is derived as(39)Ir=−2σ23β3.ForI=Ir, the phase shift Δϕ can be defined as(40)Δϕ=ϕ+∞,Ir−ϕ−∞,Ir,which will be used in subsequent studies to determine the condition under which the multipulse Shilnikov-type homoclinic orbits may exist. Substituting the first equations of (36a) and (36b) into the fourth equations of (28a)–(28d) leads to(41)ϕ˙=σ2+32β3I2+2α5ξ1δ1sech2ξ1T1.After integration, we obtain(42)ϕT1=σ2+32β3I2T1+2α5ξ1δ1tanhξ1T1+ϕ0.Accordingly, the phase shift is also obtained as(43)△ϕ=4α5ξ1δ1I=Ir=−4α5η−14μ2+σ1−σ+α5Ir2.


## 4. Dynamics of the Perturbed System

After obtaining the detailed properties of the subspacez1,z2 for systems (28a)–(28d), the next step is to examine the effects of small perturbation terms 0<ε≪1 on systems (28a)–(28d). The energy-phase method developed by Haller and Wiggins [20] is utilized to determine the existence of multipulse orbits and chaotic dynamics in the CNT-reinforced composite rectangular plate.


### 4.1. The Case of Dissipative Perturbations

We start by considering the influence of such small perturbations on manifoldM. Based on the research by Haller and Wiggins [20], we know that the manifold M is invariant for small nonzero ε. Therefore, the perturbed annulus Mε is taken to be the same as M. That is,(44)Mε=z,I,ϕz=q0I,I>I1,0≤ϕ≤2π.In order to study the dynamic behaviors of the perturbed vector field restricted toMε near the resonance I=Ir, we now introduce the following scale transformations:(45)I=Ir+εh,τ=εT1.The last two equations of (26a)–(26d) then become(46a)h′=−12μIr+f0Irsin2ϕ+ε−12μh+f0hsin2ϕ,(46b)ϕ′=3β3Irh+ε32β3h2+f0cos2ϕ,where the prime represents the differentiation with respect to τ. When ε=0, systems (46a) and (46b) become(47a)h′=−12μIr+f0Irsin2ϕ,(47b)ϕ′=3β3Irh.Systems (47a) and (47b) are a Hamiltonian system with the Hamiltonian(48)H^D=−12μIrϕ−12f0Ircos2ϕ−32β3Irh2.A simple calculation indicates that systems (47a) and (47b) have two singular points expressed by(49a)P0=0,ϕc=0,12arcsinμ2f0,(49b)Q0=0,ϕs=0,12π−12arcsinμ2f0.It is known that if6β3f0Ir2cos2ϕc<0, the singular point P0 is a center point. And if 6β3f0Ir2cos2ϕs>0, the singular point Q0 is a saddle connected to itself by a homoclinic orbit. Figure 2(a) shows the phase portrait of systems (47a) and (47b). It is observed that, for sufficiently small parameter ε, the singular point Q0 still remains a hyperbolic saddle point Qε, while for small perturbations, the singular point P0 changes to a hyperbolic sink Pε, see Figure 2(b) for the phase portrait of perturbed systems (46a) and (46b).Figure 2
Dynamics on the normally hyperbolic manifold: (a) the unperturbed case; (b) the perturbed case.
(a)
(b)On the basis of equation (48), at h=0, the basin of the attractor for ϕmin can be estimated as follows:(50)−12μIrϕmin−12f0Ircos2ϕmin=−12μIrϕs−12f0Ircos2ϕs.Substitutingϕs in (49a) and (49b) into (50) yields(51)ϕmin+f0μcos2ϕmin=12π−12arcsinμ2f0−4f02−μ22μ.The neighborhood ofIr is defined as(52)Aε=z,I,ϕz=0,I−Ir<εC,0≤ϕ≤2π,where the constant C is sufficiently large so that the unperturbed homoclinic orbit is enclosed within the annulus.


### 4.2. The Existence of Multipulse Homoclinic Orbits

Based on the results obtained by Haller and Wiggins [20], the dissipative energy-difference function is of the form(53)ΔnH^Dϕ=H^Dh,ϕ+n△ϕ−H^Dh,ϕ−∑i=1n∫−∞+∞DH,gz1T1,z2T1dT1,where(54)H^Dh,ϕ+nΔϕ−H^Dh,ϕ=−12μIrnΔϕ−12f0Ircos2ϕ+2nΔϕ−cos2ϕ,and A represents the region contained between a pair of homoclinic orbits in the z1,z2 plane, ∂Al is the boundary of A, and Δϕ denotes the phase difference given in (43). The third term in (53) can be calculated as follows:(55)∑i=1n∫−∞+∞DH,gz1T1,z2T1dT1=−2nμξ13α5Δϕ−12μIr2nΔϕ.Thus, the dissipative energy-difference function becomes(56)ΔnH^Dϕ=−12f0Ircos2ϕ+2nΔϕ−cos2ϕ+2nμξ13α5Δϕ+12nμIrIr−1Δϕ=f0Irsin2ϕ+nΔϕsinnΔϕ+2nμξ1˜3α5Δϕ,where ξ1˜=ξ1+3/4α5IrIr−1. A dissipation factor is defined as d=μ/f0. Equation (56) can then be rewritten in the following form:(57)ΔnH^Dϕ=f0Irsin2ϕ+nΔϕsinnΔϕ+2nξ1˜d3α5Δϕ.The zeros ofΔnH^Dϕ can be achieved by solving the following equation:(58)sin2ϕ+nΔϕ=−2nξ1˜dΔϕ3α5IrsinnΔϕ.The upper bound of the dissipation factord is derived from equation (58):(59)d<dmax=3α5Ir2nξ1˜sinn△ϕΔϕ.According to (59), it is found that, for small dissipative effects d<1, the following upper bound on the pulse numbers is obtained:(60)n<nmax=3α5Ir2dξ1˜Δϕ,which means that the homoclinic tree is finite even for arbitrarily small dissipation. The result in equation (60) also indicates that the upper bound nmax is inversely proportional to the dissipation factor d. We now identify zeros of the energy-difference function ΔnH^Dϕ in the interval −π/2,π/2. From (57), we find that, for any n satisfying(61)nΔϕ≠lπ,l∈Z,the zeros of ΔnH^Dϕ in −π/2,π/2 are(62a)ϕ−,1n=π2−12nΔϕ+φmodπ,(62b)ϕ−,2n=π2−12π+12nΔϕ−φmodπ,where φ=−1/2arcsin2nξ1˜dΔϕ/3α5IrsinnΔϕ. Define a set that contains all transverse zeros of ΔnH^Dϕ as(63)Z−n=h,ϕΔnH^Dϕ=0,DϕΔnH^Dϕ≠0.It is easy to know that both zeros are transverse under condition (61). We now need to introduce the two additional angles(64a)ϕ+,1n=ϕ−,1n+nΔϕmod2π,(64b)ϕ+,2n=ϕ−,2n+nΔϕmod2π,and utilize them to construct the two sets of transverse zeros of ΔnH^Dϕ as(65a)Z−n=h,ϕϕ∈ϕ−,1n,ϕ−,2n,n≥1,(65b)Z+n=h,ϕϕ∈ϕ+,1n,ϕ+,2n,n≥1.Consider a domainS0∈Aε enclosed inside the homoclinic orbit of (47a) and (47b) located in the interval ϕ∈−π/2,π/2. The periodic orbits in S0 are classified according to their pulse numbers, see Figure 3. Note that, under condition (61), all internal orbits outside S0 intersect Z−1 transversally. Thus, for any periodic or homoclinic orbit γ outside S0, the pulse number is Nγ=1. Define the energy sequence(66a)h0=H^D0,ϕs=−12μIrπ2−12arcsind2+14f0Ir4−d2,(66b)hn=minH^D0,ϕ−,1n,H^D0,ϕ−,2n,such that hn provides the energy level related to an orbit closer to the center. Define the open set of internal orbits in S0 as(67)A0=∅,An=h,ϕ∈S0H^Dh,ϕ>hn,n≥1.Figure 3
The construction of the layer sequence for (47a) and (47b).The pulse sequence is then defined as(68)N1=1,Nk=minn∈Z+n>Nk−1,hn<hNk−1,k≥2.Since the energy of the periodic orbits inS0 decreases monotonically with the orbits shrinking to the center, we necessarily have(69)AN1⊂AN2⊂⋯ANk⊂⋯.The layer sequence is defined as(70)LNk=IntANkANk−1,where Int(⋅) denotes the interior of a set. And the construction of the layer sequence is shown in Figure 3. We can easily see that all these sequences defined above are finite by (60). For any periodic orbit γ∈S0, the pulse number is Nγ=Nk. As illustrated in Figure 3, the layer radii are defined as(71)rNk=minϕc−ϕ−,1Nk,ϕc−ϕ−,2Nk.A recursive algorithm for the calculation of the distribution of pulse numbers and the layer radii is implemented on a computer. For different values of the dissipation factord, we obtain the corresponding pulse diagrams and layer radius diagrams as a function of the phase shift Δϕ for Nk<20, illustrated in Figures 4 and 5. As can be seen in Figure 4, the horizontal line segments at each level Nk identify that an infinity of Nk-pulse orbits exist for all values of the phase shift in the interval below that line. Besides, the diagrams in Figure 5 exhibit a gradual breakup of the homoclinic tree as d is increased. This observation means that the system parameters have a vital effect on the distribution of the pulse numbers and the layer radii.Figure 4
The pulse sequence as a function of the phase shiftξ1˜/α5Ir=2/15. (a) d=0. (b) d=10−5. (c) d=0.01. (d) d=0.05. (e) d=0.1. (f) d=0.4.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
The layer radius sequence as a function of the phase shiftξ1˜/α5Ir=2/15. (a) d=0. (b) d=10−5. (c) d=0.01. (d) d=0.05. (e) d=0.1. (f) d=0.4.
(a)
(b)
(c)
(d)
(e)
(f)We now detect the existence of the multipulse Shilnikov-type orbits based on Theorem 4.5 in [20]. Note that the sinks of (46a) and (46b) become the centers of (47a) and (47b). Therefore, systems (47a) and (47b) have a nondegenerate equilibrium(72)P0≡hc,ϕc=0,ϕc=0,12arcsind2.The values ofd,Δϕ,N for which the center P0 falls in the zero set Z−n given in (63) must be determined, that is,(73)−12Ircosarcsind2+2NΔϕ−cosarcsind2+2Nξ1˜dΔϕ3α5=0,which yields(74)4−d2cos2NΔϕ−1=dsin2NΔϕ+8Nξ1˜Δϕ3α5Ir.Solving (74), the dissipation factor d is obtained as(75)d=21−cos2NΔϕ1−cos2NΔϕ2+sin2NΔϕ+8Nξ1˜Δϕ/3α5Ir2.The result is only valid when(76)Δϕ≠mπN,m∈Z.It is also necessary to prove the nondegeneracy condition(77)Dd−12Ircosarcsind2+2NΔϕ−cosarcsind2+2Nξ1˜dΔϕ3α5≠0,whenever (74) and (76) hold. Carrying out the differentiation in expression (74), we find that (77) fails to be satisfied when(78)−dcos2NΔϕ−1=4−d2sin2NΔϕ+8Nξ1˜Δϕ3α5Ir.It is easy to see that equations (74) and (78) cannot hold simultaneously under condition (76). Thus, nondegeneracy condition (77) holds. We now have to guarantee that the landing point of any N-pulse orbit taking off from a slow sink lies in the domain of attraction of a sink on Aε. To verify it, we consider the point in the interval −π/2,π/2 and is kπ apart from the approximate landing point ϕc+NΔϕ. The ϕ-coordinate of this point is calculated as(79)ϕ∗N=ϕs+ϕc+NΔϕ−ϕsmodπ,in which we recall that(80a)ϕc=12arcsind2,(80b)ϕs=π2−12arcsind2.Ifϕ∗N>ϕs, we redefine ϕ∗N by subtracting π. The main purpose here is to achieve kπ translation of the landing point which is the closest to the saddle point ϕs in the interval −π/2,π/2. We require that(81)H^D0,ϕ∗N<H^D0,ϕs,or, equivalently, that(82)cos2ϕ∗N−cos2ϕs>dϕs−ϕ∗N,to make sure that the projection of landing points falls in the domain of attraction of one of the sinks.Based on the above analysis, the results in this section give rise to the following theorem.Theorem 1.
For any integerN≥1, there exists a positive number ε0N>0 and a finite union CN of codimension-one surfaces in the ξ1˜,α5,μ,f0,ε∈CN parameter space near the set satisfying 0<d<1, (75), and (82) such that, for any ξ1˜,α5,μ,f0,ε∈CN and 0<ε<ε0N, the following conclusions hold:(1)
If the integer(83)Q=INT12+NΔϕ+arcsind/2π,is even, then each of the saddle-focus-type equilibria contained in the slow manifoldAε admits two generalized Shilnikov-type homoclinic orbits. If Q is odd, then there exist two cycles of Shilnikov-type heteroclinic orbits connecting the two saddle-foci to each other. In both cases, the N-pulse orbits form pairs which are symmetric with respect to the subspace z1,z2=0,0.(2)
There exists an open set of parameters containingCN for which system (26) admits Smale horseshoes in its dynamics.


## 4.1. The Case of Dissipative Perturbations

We start by considering the influence of such small perturbations on manifoldM. Based on the research by Haller and Wiggins [20], we know that the manifold M is invariant for small nonzero ε. Therefore, the perturbed annulus Mε is taken to be the same as M. That is,(44)Mε=z,I,ϕz=q0I,I>I1,0≤ϕ≤2π.In order to study the dynamic behaviors of the perturbed vector field restricted toMε near the resonance I=Ir, we now introduce the following scale transformations:(45)I=Ir+εh,τ=εT1.The last two equations of (26a)–(26d) then become(46a)h′=−12μIr+f0Irsin2ϕ+ε−12μh+f0hsin2ϕ,(46b)ϕ′=3β3Irh+ε32β3h2+f0cos2ϕ,where the prime represents the differentiation with respect to τ. When ε=0, systems (46a) and (46b) become(47a)h′=−12μIr+f0Irsin2ϕ,(47b)ϕ′=3β3Irh.Systems (47a) and (47b) are a Hamiltonian system with the Hamiltonian(48)H^D=−12μIrϕ−12f0Ircos2ϕ−32β3Irh2.A simple calculation indicates that systems (47a) and (47b) have two singular points expressed by(49a)P0=0,ϕc=0,12arcsinμ2f0,(49b)Q0=0,ϕs=0,12π−12arcsinμ2f0.It is known that if6β3f0Ir2cos2ϕc<0, the singular point P0 is a center point. And if 6β3f0Ir2cos2ϕs>0, the singular point Q0 is a saddle connected to itself by a homoclinic orbit. Figure 2(a) shows the phase portrait of systems (47a) and (47b). It is observed that, for sufficiently small parameter ε, the singular point Q0 still remains a hyperbolic saddle point Qε, while for small perturbations, the singular point P0 changes to a hyperbolic sink Pε, see Figure 2(b) for the phase portrait of perturbed systems (46a) and (46b).Figure 2
Dynamics on the normally hyperbolic manifold: (a) the unperturbed case; (b) the perturbed case.
(a)
(b)On the basis of equation (48), at h=0, the basin of the attractor for ϕmin can be estimated as follows:(50)−12μIrϕmin−12f0Ircos2ϕmin=−12μIrϕs−12f0Ircos2ϕs.Substitutingϕs in (49a) and (49b) into (50) yields(51)ϕmin+f0μcos2ϕmin=12π−12arcsinμ2f0−4f02−μ22μ.The neighborhood ofIr is defined as(52)Aε=z,I,ϕz=0,I−Ir<εC,0≤ϕ≤2π,where the constant C is sufficiently large so that the unperturbed homoclinic orbit is enclosed within the annulus.


## 4.2. The Existence of Multipulse Homoclinic Orbits

Based on the results obtained by Haller and Wiggins [20], the dissipative energy-difference function is of the form(53)ΔnH^Dϕ=H^Dh,ϕ+n△ϕ−H^Dh,ϕ−∑i=1n∫−∞+∞DH,gz1T1,z2T1dT1,where(54)H^Dh,ϕ+nΔϕ−H^Dh,ϕ=−12μIrnΔϕ−12f0Ircos2ϕ+2nΔϕ−cos2ϕ,and A represents the region contained between a pair of homoclinic orbits in the z1,z2 plane, ∂Al is the boundary of A, and Δϕ denotes the phase difference given in (43). The third term in (53) can be calculated as follows:(55)∑i=1n∫−∞+∞DH,gz1T1,z2T1dT1=−2nμξ13α5Δϕ−12μIr2nΔϕ.Thus, the dissipative energy-difference function becomes(56)ΔnH^Dϕ=−12f0Ircos2ϕ+2nΔϕ−cos2ϕ+2nμξ13α5Δϕ+12nμIrIr−1Δϕ=f0Irsin2ϕ+nΔϕsinnΔϕ+2nμξ1˜3α5Δϕ,where ξ1˜=ξ1+3/4α5IrIr−1. A dissipation factor is defined as d=μ/f0. Equation (56) can then be rewritten in the following form:(57)ΔnH^Dϕ=f0Irsin2ϕ+nΔϕsinnΔϕ+2nξ1˜d3α5Δϕ.The zeros ofΔnH^Dϕ can be achieved by solving the following equation:(58)sin2ϕ+nΔϕ=−2nξ1˜dΔϕ3α5IrsinnΔϕ.The upper bound of the dissipation factord is derived from equation (58):(59)d<dmax=3α5Ir2nξ1˜sinn△ϕΔϕ.According to (59), it is found that, for small dissipative effects d<1, the following upper bound on the pulse numbers is obtained:(60)n<nmax=3α5Ir2dξ1˜Δϕ,which means that the homoclinic tree is finite even for arbitrarily small dissipation. The result in equation (60) also indicates that the upper bound nmax is inversely proportional to the dissipation factor d. We now identify zeros of the energy-difference function ΔnH^Dϕ in the interval −π/2,π/2. From (57), we find that, for any n satisfying(61)nΔϕ≠lπ,l∈Z,the zeros of ΔnH^Dϕ in −π/2,π/2 are(62a)ϕ−,1n=π2−12nΔϕ+φmodπ,(62b)ϕ−,2n=π2−12π+12nΔϕ−φmodπ,where φ=−1/2arcsin2nξ1˜dΔϕ/3α5IrsinnΔϕ. Define a set that contains all transverse zeros of ΔnH^Dϕ as(63)Z−n=h,ϕΔnH^Dϕ=0,DϕΔnH^Dϕ≠0.It is easy to know that both zeros are transverse under condition (61). We now need to introduce the two additional angles(64a)ϕ+,1n=ϕ−,1n+nΔϕmod2π,(64b)ϕ+,2n=ϕ−,2n+nΔϕmod2π,and utilize them to construct the two sets of transverse zeros of ΔnH^Dϕ as(65a)Z−n=h,ϕϕ∈ϕ−,1n,ϕ−,2n,n≥1,(65b)Z+n=h,ϕϕ∈ϕ+,1n,ϕ+,2n,n≥1.Consider a domainS0∈Aε enclosed inside the homoclinic orbit of (47a) and (47b) located in the interval ϕ∈−π/2,π/2. The periodic orbits in S0 are classified according to their pulse numbers, see Figure 3. Note that, under condition (61), all internal orbits outside S0 intersect Z−1 transversally. Thus, for any periodic or homoclinic orbit γ outside S0, the pulse number is Nγ=1. Define the energy sequence(66a)h0=H^D0,ϕs=−12μIrπ2−12arcsind2+14f0Ir4−d2,(66b)hn=minH^D0,ϕ−,1n,H^D0,ϕ−,2n,such that hn provides the energy level related to an orbit closer to the center. Define the open set of internal orbits in S0 as(67)A0=∅,An=h,ϕ∈S0H^Dh,ϕ>hn,n≥1.Figure 3
The construction of the layer sequence for (47a) and (47b).The pulse sequence is then defined as(68)N1=1,Nk=minn∈Z+n>Nk−1,hn<hNk−1,k≥2.Since the energy of the periodic orbits inS0 decreases monotonically with the orbits shrinking to the center, we necessarily have(69)AN1⊂AN2⊂⋯ANk⊂⋯.The layer sequence is defined as(70)LNk=IntANkANk−1,where Int(⋅) denotes the interior of a set. And the construction of the layer sequence is shown in Figure 3. We can easily see that all these sequences defined above are finite by (60). For any periodic orbit γ∈S0, the pulse number is Nγ=Nk. As illustrated in Figure 3, the layer radii are defined as(71)rNk=minϕc−ϕ−,1Nk,ϕc−ϕ−,2Nk.A recursive algorithm for the calculation of the distribution of pulse numbers and the layer radii is implemented on a computer. For different values of the dissipation factord, we obtain the corresponding pulse diagrams and layer radius diagrams as a function of the phase shift Δϕ for Nk<20, illustrated in Figures 4 and 5. As can be seen in Figure 4, the horizontal line segments at each level Nk identify that an infinity of Nk-pulse orbits exist for all values of the phase shift in the interval below that line. Besides, the diagrams in Figure 5 exhibit a gradual breakup of the homoclinic tree as d is increased. This observation means that the system parameters have a vital effect on the distribution of the pulse numbers and the layer radii.Figure 4
The pulse sequence as a function of the phase shiftξ1˜/α5Ir=2/15. (a) d=0. (b) d=10−5. (c) d=0.01. (d) d=0.05. (e) d=0.1. (f) d=0.4.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
The layer radius sequence as a function of the phase shiftξ1˜/α5Ir=2/15. (a) d=0. (b) d=10−5. (c) d=0.01. (d) d=0.05. (e) d=0.1. (f) d=0.4.
(a)
(b)
(c)
(d)
(e)
(f)We now detect the existence of the multipulse Shilnikov-type orbits based on Theorem 4.5 in [20]. Note that the sinks of (46a) and (46b) become the centers of (47a) and (47b). Therefore, systems (47a) and (47b) have a nondegenerate equilibrium(72)P0≡hc,ϕc=0,ϕc=0,12arcsind2.The values ofd,Δϕ,N for which the center P0 falls in the zero set Z−n given in (63) must be determined, that is,(73)−12Ircosarcsind2+2NΔϕ−cosarcsind2+2Nξ1˜dΔϕ3α5=0,which yields(74)4−d2cos2NΔϕ−1=dsin2NΔϕ+8Nξ1˜Δϕ3α5Ir.Solving (74), the dissipation factor d is obtained as(75)d=21−cos2NΔϕ1−cos2NΔϕ2+sin2NΔϕ+8Nξ1˜Δϕ/3α5Ir2.The result is only valid when(76)Δϕ≠mπN,m∈Z.It is also necessary to prove the nondegeneracy condition(77)Dd−12Ircosarcsind2+2NΔϕ−cosarcsind2+2Nξ1˜dΔϕ3α5≠0,whenever (74) and (76) hold. Carrying out the differentiation in expression (74), we find that (77) fails to be satisfied when(78)−dcos2NΔϕ−1=4−d2sin2NΔϕ+8Nξ1˜Δϕ3α5Ir.It is easy to see that equations (74) and (78) cannot hold simultaneously under condition (76). Thus, nondegeneracy condition (77) holds. We now have to guarantee that the landing point of any N-pulse orbit taking off from a slow sink lies in the domain of attraction of a sink on Aε. To verify it, we consider the point in the interval −π/2,π/2 and is kπ apart from the approximate landing point ϕc+NΔϕ. The ϕ-coordinate of this point is calculated as(79)ϕ∗N=ϕs+ϕc+NΔϕ−ϕsmodπ,in which we recall that(80a)ϕc=12arcsind2,(80b)ϕs=π2−12arcsind2.Ifϕ∗N>ϕs, we redefine ϕ∗N by subtracting π. The main purpose here is to achieve kπ translation of the landing point which is the closest to the saddle point ϕs in the interval −π/2,π/2. We require that(81)H^D0,ϕ∗N<H^D0,ϕs,or, equivalently, that(82)cos2ϕ∗N−cos2ϕs>dϕs−ϕ∗N,to make sure that the projection of landing points falls in the domain of attraction of one of the sinks.Based on the above analysis, the results in this section give rise to the following theorem.Theorem 1.
For any integerN≥1, there exists a positive number ε0N>0 and a finite union CN of codimension-one surfaces in the ξ1˜,α5,μ,f0,ε∈CN parameter space near the set satisfying 0<d<1, (75), and (82) such that, for any ξ1˜,α5,μ,f0,ε∈CN and 0<ε<ε0N, the following conclusions hold:(1)
If the integer(83)Q=INT12+NΔϕ+arcsind/2π,is even, then each of the saddle-focus-type equilibria contained in the slow manifoldAε admits two generalized Shilnikov-type homoclinic orbits. If Q is odd, then there exist two cycles of Shilnikov-type heteroclinic orbits connecting the two saddle-foci to each other. In both cases, the N-pulse orbits form pairs which are symmetric with respect to the subspace z1,z2=0,0.(2)
There exists an open set of parameters containingCN for which system (26) admits Smale horseshoes in its dynamics.


## 5. Numerical Simulations

In order to verify the analytical predictions, we choose systems (13a)–(13d) to perform numerical simulations. The fourth-order Runge–Kutta algorithm is employed to indicate the existence of the multipulse Shilnikov-type homoclinic orbits and chaotic motions in the CNT-reinforced composite plate. Figure 6 is obtained to show that there exist multipulse chaotic motions for the CNT-reinforced composite plate. In Figure 6, the parameters and initial conditions are selected as μ=0.05, σ1=−1.86, σ2=1.5, α1=7.2, α2=−10.8, α3=−3.9, α4=11.6, α5=15, β2=−4.78, β3=−4, β4=5, β5=−7.18, f0=12.235, x1,x2,x3,x4=0.105,−0.101,−0.013,−0.506, and the phase shift at I=Ir=−2σ2/3β3=1/2 is Δϕ=2.37. We choose the pulse number N=4; then, (75) gives d=μ/f0=0.00409. Using equations (79) and (80a) and (80b), we have ϕ∗N=3.199, ϕs=1.569, cos2ϕ∗N−cos2ϕs=1.993, and dϕs−ϕ∗N=−0.067. Obviously, condition (82) holds, which implies that the chaotic motion presented in Figure 6 is multipulse Shilnikov-type chaotic motion.Figure 6
The multipulse chaotic motion obtained based on equation (13): (a) the phase portrait on the plane x1,x2, (b) the phase portrait on the plane x3,x4, (c) the waveform on the plane t,x1, (d) the waveform on the plane t,x3, (e) the phase portraits in the three-dimensional space x1,x2,x3, and (f) the phase portraits in the three-dimensional space x2,x3,x4.
(a)
(b)
(c)
(d)
(e)
(f)


## 6. Conclusions

In this paper, the multipulse homoclinic orbits and chaotic dynamics of a CNT-reinforced composite plate under combined in-plane and transverse excitations are studied. The method of multiple scales is adopted to acquire the averaged equations in the case of 1 : 1 internal resonance, principal parametric resonance, and 1/2 subharmonic resonance. On the basis of the averaged equations obtained, the normal form theory is employed to derive the expressions of normal form associated with a double zero and a pair of pure imaginary eigenvalues. The energy-phase method proposed by Haller and Wiggins [20] is utilized to detect the presence of the multipulse homoclinic orbits and chaotic dynamics for the resonant case. As the trajectory of motion comes close to the sink point Pε, every Shilnikov-type orbit takes off again and repeats this similar motion in the full four-dimensional phase space and eventually leads to the multipulse Shilnikov-type orbits. Our analysis in Section 4 also indicates that the multipulse Shilnikov-type orbits depend on the system parameters and dissipative perturbations. It is known that the existence of multipulse Shilnikov-type orbits implies the existence of chaos in the sense of Smale horseshoes. Homoclinic trees are presented to describe the repeated bifurcations of multipulse solutions. From the diagrams, we can see a gradual breakup of the homoclinic tree with the increase of the dissipation factor. This observation denotes that the damping coefficient and system parameters affect the distribution of the pulse numbers and the layer radii. The analytical results obtained here are extensions of those appearing in [19].To confirm the theoretical results, numerical simulations are applied to examine the chaotic dynamics of the CNT-reinforced composite plate. The numerical results demonstrate that there exist multipulse homoclinic orbits and chaotic dynamics in systems (13a)–(13d) when the conditions are satisfied. As we all know, under certain conditions, the multipulse chaotic motions in the averaged equations can result in the multipulse amplitude-modulated chaotic vibrations in the original system, which indicates the existence of multipulse chaotic dynamics in the CNT-reinforced composite plate. The research in this paper provides a detailed explanation of the multipulse jumping behaviors for the CNT-reinforced composite plate under combined in-plane and transverse excitations. It is established simultaneously that the damping coefficient and system parameters have a significant influence on the nonlinear dynamics of the CNT-reinforced composite plate. Therefore, as mentioned above, the nonlinear dynamic behaviors of the CNT-reinforced composite plate can be controlled by varying the system parameters.



---

*Source: 7310187-2020-04-29.xml*