# Theory of Nonlinear Vibrations for Antisymmetric Cross-Ply Bistable Laminated Shells

Authors: B. Y. Zhang; W. Zhang
Journal: Shock and Vibration
 (2021)

---

## Abstract

In this paper, on the basis of taking the von Karman nonlinear factors into consideration, the constitutive equation of the antisymmetric cross-ply laminated composite was used to calculate the internal force and internal moment of the bistable structure, and the dynamic equilibrium equation and the compatible equation were constructed, respectively. The two equations were combined to establish a nonlinear dynamic model for the antisymmetric cross-ply laminated glass fiber resin bistable shell. Then, the finite element numerical simulation software ABAQUS was adopted to perform simulation modeling and numerical analysis on a series of bistable specimens, so as to study the impact of different geometric parameters on the frequencies, mode shapes, and other vibration characteristics of the antisymmetric laminated fiber resin bistable shell. Galerkin discretization was conducted on the vibration partial differential equation. Since there are only even-order partial differential terms of deflectionw with respect to x and y in the vibration partial differential equation at this time, the form of series obtained by each term is the same, which simplifies the discretization of the dynamic equilibrium equation and the compatible equation. Finally, the two equations after discretization were merged to obtain the three-degree-of-freedom nonlinear ordinary differential equation of the antisymmetric cross-ply laminated glass fiber resin bistable shell. The system averaged equation was acquired by perturbation analysis through a multiscale method, and the periodic solution of the antisymmetric laminated bistable system was studied. Moreover, the system’s nonlinear dynamic behavior characteristics such as bifurcation and chaos were explored when the main resonance Ω is close to ω1 and ω2, respectively, and the internal resonance is 1 : 2 : 3.


---

## Body

## 1. Introduction

A bistable structure is formed by laying multiple layers of fiber-reinforced materials, which is a composite laminated structure with two different stable states. Generally, under the driving of external loads such as mechanical force, smart materials, and temperature field, it can be changed from one stable state to the other, and its stable configuration can be maintained without the need for sustained energy input. Hyer [1] first discovered that asymmetric laminated cylindrical shells have bistable characteristics. In 1996, Daton-Lovett [2] found that antisymmetric laminated composite cylindrical shells also present the two stable states of a cylinder, which attracted considerable attention from many scholars, and relevant theoretical, numerical, and experimental studies are successively conducted. Innovatively, this paper found the law of deformation cloud diagrams of the antisymmetric laminated shallow shell and cylindrical shell made of T300 and studied the influence of geometric parameters (e.g., length, initial cross-section radius, initial central angle, number plies, and ply angle) on the modal frequency and mode shape of the antisymmetric laminated cylindrical shell.The theoretical study of antisymmetric bistable laminated structures is mainly focused on statics and dynamics. In statics, Iqbal [3] proposed a simple linear elastic bistable model based on the classical laminate theory and predicted the bistable behavior of antisymmetric laminated shells according to the principle of minimum potential energy. Galletly and Guest [4,5] successively proposed the beam model and the shell model regarding the theoretical model developed by Iqbal. The theoretical results are in good agreement with the finite element simulation results, but greatly different from the experimental results. Pellegrino and Guest [6] proposed the double-parameter model by referring to previous studies. This model can be used to calculate and analyze the bistable characteristics of antisymmetric and symmetric laminated isotropic cylindrical shells as well as those with prestress. The College of Engineering, University of Cambridge successively studied the bistable characteristics of multistable composite shells with surface folds [7–10] and isotropic laminated shells with prestress [11]. Concerning the studies on the nonlinear dynamics of bistable structures, Zhang and Zheng [12] adopted the extended high-dimensional Melnikov method and numerical method and found that the multipulse double-parameter chaotic motion of asymmetric bistable laminated square plates will occur under a certain coupling effect of external excitation and parametric excitation. Zhang and Liu [13] studied the dynamic jump phenomenon of asymmetric bistable laminated square plates and its nonlinear vibration under base excitation from both theoretical and experimental aspects. Fei [14] analyzed the jump behavior of a two-layer orthogonal bistable laminated plate, while considering factors such as thickness, temperature, and external excitation of the laminated plate.For the numerical analysis, the stable state transition of bistable laminated structures can be caught by the commercial finite element software, so the calculation result is more intuitive. Also, numerical analysis can be performed to perfect the theoretical model and predict the bistable characteristics of the structure, thus providing references for the experimental study. Iqbal and Pellegrinot [15] conducted numerical simulation on the bistable behavior of antisymmetric laminated shells using ABAQUS software and analyzed the stable state transition of antisymmetric laminated cylindrical shells and the stress distribution of the second stable state structure. Huang [16] introduced the bistable differential system into each cell of the grid and studied the effect of topology on the stable-state trend of the network.In terms of the experimental study, Dano and Hyer [17] measured the jump loads of the bistable laminated plate by adjusting the water load, and the experimental measurements are in good agreement with the theoretical calculations. Potter et al. [18] adopted the loading method similar to “three-point bending” to induce the jump phenomenon of the bistable laminated plate and measured the load-displacement curve of the center of the form of the laminated plate. Daynes and Tawfik [19, 20] adopted similar experimental loading schemes and obtained the jump loads of orthogonal laminated plates with different geometric dimensions and the relations of curvature variation between two stable states. Etches [21] considered the effect of humidity on the mechanical properties of a bistable laminate plate and found that the jump load of the asymmetric orthogonal bistable laminated plate changes significantly at the initial stage of the experiment and then gradually levels off. Dano and Hyer [22] realized the intelligent control of asymmetric bistable laminated plates by controlling the temperature of the alloy wires to generate different driving forces. Schultz et al. [23] utilized a bistable laminated plate driven by piezoelectric materials to realize the stable state transition, but there is a big error between experimental and theoretical results. Tawfik et al. [24] realized the mutual transition between the two stable states by applying voltage, and the error between experimental and finite element simulation results is less than 10%. Giddings et al. [25] used piezoelectric patches to control the deformation of bistable laminated plates and studied the relationship between voltage and deformation, which played a guiding role in the design and control of the intelligent deformable structure. Kim et al. [26] realized the mutual conversion of bistable cantilever beams in two stable states by using MFC and SMA.In the research of practical engineering application, the University of Cambridge [15] studied the bistable thin-walled tube similar to a tape measure and realized the stable state transition by controlling the energy accumulated in the folding process, which enables the bistable tap structure to be used for automatic hinge locks, detectors, expandable structures, and other aerospace devices. Schlecht [27] applied the bistable torsional assembly structure to the wing configuration of unmanned aerial vehicles, which in turn, can improve the sensitivity of aircraft during rotation because it requires only a small amount of energy to drive the wing to deformation. Diaconu and Weaver [28] used a bistable structure to design a wing model that can meet the actual conditions and analyzed the application prospect of such wing as a variant structure. Mattioni [29] made variable wing structures using bistable structures and used intelligent components to control their deformation to change the wing structure. The bistable structure can serve as both a deformable structure and as an energy harvesting structure. By utilizing the nonlinear vibration in the stable-state jump process, the bistable structure can be shaped into a broadband energy harvester [30–32]. Three bistable structures have been applied to the piezoelectric energy capture, which is formed by using the principle of homopolar repulsion of magnetic force [33–38], buckling of the beam, and thermal stress [39, 40], respectively, and the conversion of mechanical energy to electrical energy is achieved during vibration. Yao [41–45] separated the end of the upper piezoelectric layer from the base layer in the traditional piezoelectric cantilever beam and conducted an experimental study on broadband energy collection and dynamic response of the L-shaped piezoelectric cantilever beam. Chen [46, 47] studied a micro/nanoscale bistable plate electrical-thermal-mechanical coupling system for energy harvesting and analyzed the voltage-modal frequency response of a nonlinear system.


## 2. Dynamic Modeling

The antisymmetric cross-ply laminated structure is formed by alternately laying the orthotropic single-layer structure at an angle of0∘ and 90∘ between the main direction of the material and the coordinate axis. The stiffness coefficient of the antisymmetric cross-ply laminated structure has the following relationship: Q110°=Q2290°, Q220°=Q1190°, and Q16=Q16=0. Thus, A11=A22, D11=D22, and A16=A26=D16=D26=B16=B26=0, and it can be proved that B12=B66 and B22=−B11. In this study, four layers 0°/90°/0°/90° antisymmetric cross-ply laminated bistable thin shell was selected as the model, as shown in Figure 1, and the specific laying mode is shown in Figure 2. Glass fiber and epoxy resin substrate were used as the materials for the shell. The thickness of a single layer is 0.185 mm. Other material parameters are displayed in Table 1.Figure 1
Model of the bistable structure under the action of exciter. (a) First natural equilibrium position. (b) Second natural equilibrium position.
(a)(b)Figure 2
Schematic diagram of antisymmetrically cross-ply by0°/90°/0°/90°.Table 1
S1002 Single-layer material parameters.
E11GPaE22GPaE12GPaE12GPaE12GPav12ΡFiber (%)398.44.24.24.20.261.645Based on the classical nonlinear shallow shell theory and the von Karman large deformation theory, the nonlinear response of the antisymmetric cross-ply laminated bistable shallow shell was studied. The nonlinear geometric equation of the middle plane of the Donnell cylindrical shell was applied:(1a)εx0=∂u0∂x+12∂w∂x2+wRx,(1b)εy0=∂v0∂y+12∂w∂x2+wRy,(1c)γxy0=∂u0∂y+∂v0∂x+∂w∂x∂w∂y,where u0, v0, and w represent the displacement of the point on the middle plane of the cylindrical shell in the directions of x, y, and z and Rx and Ry represent the initial radius of curvature of the cross-section in the x and y directions of the cylindrical shell.Since the main direction of vibration of the shallow shell is transverse deflection, the inertia in the plane is ignored. Because of the coupling between the in-plane displacement and the transverse displacement in the geometric equation, it is necessary to introduce a compatible equation to form a complete set of equations. The compatible equation of the double-curved shallow shell under plane stress is expressed as(2)∂2εx0∂y2+∂2εy0∂x2−∂2γxy0∂x∂y=1Ry∂2w∂x2+∂2w∂x∂y2−∂2w∂x2∂2w∂y2.There are six independent components of internal force in the theory of shallow shell vibration:Nx, Ny, Nxy, Mx, My, and Mxy. After the shear force Qx and Qy was eliminated, three dynamic equilibrium equations can be obtained:(3a)Nx′+Nxy∗=ρhu¨,(3b)Ny∗+Nxy′=ρhv¨,(3c)Mx″+My∗∗+2M′xy∗−NxRx+NyRy+Q−ρhw¨−γw˙=0,where Q represents the external excitation received by the shell and γ represents the damping coefficient.The three components of internal forceNx, Ny, and Nxy can be connected by the Airy stress function according to equation (3) in the following way:(4a)Nx=∂2Φ∂y2,(4b)Ny=∂2Φ∂x2,(4c)Nxy=∂2Φ∂x∂y.By substituting equation (1) into equation (2), the first two dynamic equilibrium equations (3a) and (3b) can be automatically satisfied. The expressions of the internal force N and internal moment M of the antisymmetric cross-ply laminated composite material can be written as(5)NMNxNyNxyMxMyMxy=A11A120B1100A12A2200−B11000A66000B1100D11D1200−B110D12D22000000D66εx0εy0γxy0κxκyκxy.It can be seen from the above equation that the antisymmetric cross-ply laminated bistable structure has the coupling between stretching and bending. Combining equation (4a4b4c) with equation (5), the expression of the middle plane strain expressed by Pij can be obtained:(6a)εx0=P12B11∂2w∂y2+P11B11∂2w∂x2−P12∂2Φ∂y2+P11∂2Φ∂y2,(6b)εy0=−P22B11∂2w∂y2−P12B11∂2w∂x2+P22∂2Φ∂x2−P12∂2Φ∂y2a,(6c)γxy0=−P33∂2Φ∂x∂y,where P11=A22/A11A22−A122, P12=A12/A11A22−A122, P22=A11/A11A22−A122, and P33=1/A66.By substituting simplified equation (5) of the internal force and internal moment on the cross section of the antisymmetric cross-ply thin shell and midsurface strain expression (6a)–(6c) into dynamic equilibrium equation (3c), we have(7)ρhw¨+Cw˙+P22∂4Φ∂x4−B112P11−D11∂4w∂x4−B112P22−D22∂4w∂y4+B11P22−P11∂4Φ∂x2∂y2+1Rx∂2Φ∂y2+1Ry∂2Φ∂x2+P12B11∂4Φ∂x4−P12B11∂4Φ∂y4−2B112P22−D12−2D66∂4w∂x2∂y2=Q.Then, by substituting simplified equation (5) of internal force and internal moment and the expression of midsurface strain (6a6b6c) into compatible equation (2), we have(8)P22∂4Φ∂x4+P11∂4Φ∂y4+P33−2P12∂4Φ∂x2∂y2=1Rx∂2w∂y2+1Ry∂2w∂x2−∂2w∂x2∂2w∂y2+∂2w∂x∂y2+P12B11∂4w∂x4−P12B11∂4w∂y4+B11P22−P11∂4w∂x2∂y2.In the discussion of the mode superposition method, the solution of a continuous system can be written as a linear combination of all mode functions. Since the boundary condition of the model is that the four sides are free, the center is fixed, the selected shape function does not need to satisfy any geometric boundary conditions, and the displacement shape function in the direction ofz is given by(9)w0x,y,t=∑i=0M∑j=0NWijtwijx,y,=∑i=0I∑j=0jWijtcosπxiLxcosπxjLy+∑i=I+12I+1∑j=J+12JWijtcosπxi−I+1Lxsinπyj+JLy+∑i=2I+13I+1∑j=2J+23JWijtsinπxi−3I+1Lxsinπyj−3J+1Ly+∑i=3I+14I+1∑j=2J+13J+1Wijtsinπxi−2I+1Lxcosπyj−2J+1Ly.where Wijt represents generalized coordinates and wijx.y is the mode function of the system, where the subscript i=0,1…,M and j=0,1…,N, and all the shape function terms are the combined operations of nonzero trigonometric functions, where M=4I+1, and M×N is the number of degrees of freedom selected in the approximation. Since the model built is a shallow shell structure and the area defined by the transverse displacement wx,y,t and stress function Φx,y,t is the same, it is assumed that they can be expanded with the same shape function wijx,y. Take the mode shape function of the transverse displacement w and stress function Φ, as shown below:(10a)wx,y,t=∑i=0N∑j=0Nwijx,yWijt,(10b)Φx,y,t=∑i=0N∑n=0Nwmnx,yFmnt,where Wij and Fmn are, respectively, the generalized coordinates of transverse displacement w and stress functions Φ. First, substituting equation (10a10b) into dynamic equilibrium equation (7) leads to(11)ρhw¨ijW¨ijt+Cw˙ijW˙ijt+1Rxwmn∗∗+1Rywmn″+P12B11wmn″″−P12B11wmn∗∗∗∗+B11P22−P11w″mn∗∗Fijt+D11−B112P11ij″″wB112P22−D22wij∗∗∗∗−2B112P22−D12−2D66″ij∗∗Wijt=Q.Then, equation (10a10b) was substituted into compatible equation (8), and we have(12)∑m=0N∑n=0NP22wmn″″+P11wmn∗∗∗∗+P33−2P12w″mn∗∗Fijt=∑i=0N∑j=0N∑p=0N∑q=0Nw′ij∗w′pq∗−wij″wpq∗∗2WijtWpqt+∑m=0N∑n=0N1Rxwij∗∗+1Rywij″+P12B11wij″″P12B11wij∗∗∗∗+B11P22−P11w″ij∗∗Wijt.


## 3. Galerkin Discretization

Then, the two equations were discretized by Galerkin so that the infinite degree of the freedom system can be discretized into a finite degree of the freedom system. Takeαi=πi/L, βj=πj/L. Both sides of equilibrium equation (11) were multiplied by wabx,y, and they were integrated over the area to obtain(13)L24ρhW¨ab+L24CW˙ab+L24D11−B112P11λi4−B112P22−D22γj4−2B112P22−D12−2D66λi2γj2Wab+L24B11P22−P11λi2γj2+P12B11λi4−P12B11γj4+1Ryλi2+1Rxγj2Fmn=Q.Both sides were divided byρhL2/4, and the equilibrium equation obtained is as follows:(14)W¨abt+2ζab,plateωab,plateW˙abt+ωab,plate2Wabt+ΓabFmnt=Qab,where ζab,plate=Cab/2ρhωab,plate, ωab,plate2=1/ρhD11−B112P11λi4−B112P22−D22γj4−2B112P22−D12−2D66λi2γj2, and Γab=1/ρhB11P22−P11λi2γj2+P12B11λi4−P12B11γj4+1/Ryλi2+1/Rxγj2. Multiply both sides of compatible equation (12) by wmnx,y and integrate over the area, and we have(15)L24P22λi4+P11γj4+P33−2P12λi2γj2Fmnt=1Ryλi2+1Rxγj2+P12B11λi4−P12B11γj4+B11P22−P11λi2γj2Wmnt,+λiγjλpγq∫∫wijwpqwmnds−λi2γj2∫∫wijwpqwmndsWijtWpqt.Let the coefficient ahead ofFmn1t be 1, and the compatible equation obtained is(16)Fmnt=Gmn−1HmnWmn+Gmn−1Tmnijpq−∏mnijpqWijtWpqt,where Tmnijpq=λiγjλpγqΘmnijpq, Πmnijpq=λi2γq2Θmnijpq, Gmn=L2/4P22λi4+P11γj4+P33−2P12λi2γj2, Θmnijpq=∫∫wijwpqwmndS, and Hmn=L2/4P12B11λi4−P12B11γj4+B11P22−P11λi2γj2+1/Ryλi2+1/Rxγj2Due to the decoupling between the generalized coordinates of the stress function and the generalized coordinates of the transverse displacement, the equation set can be decoupled and the governing motion equation of the transverse displacement can be given. By substituting equation (16) into equation (14), the final vibration differential equation can be obtained:(17)W¨abt+2ζab,plateωab,plateW˙abt+ωab2Wabt+ΓabGab−1Tabijpq−ΠabijpqWijtWpqt=Qab,where ωab2=ωab,plate2+ΓabGmn−1Hmn.Since low-order mode shapes play a major role in vibration, the higher the order, the smaller the role. To facilitate writing, the displacementwn in the deflection direction of the plate was rewritten as xn, and the load excitation Q was replaced by F. Therefore, the first three modes were selected to obtain the three-degree-of-freedom ordinary differential governing equations:(18a)x¨1+2ζ1ω1x˙1+ω12x1+M¯1x12+M¯2x22+M¯3x32+M¯4x1x2+M¯5x1x3+M¯x2x3=F1cosΩt,(18b)x¨2+2ζ2ω2x˙2+ω22x2+7x12+M¯8x22+M¯9x32+M¯10x1x2+M¯11x1x3+M¯12x2x3=F2cosΩt,(18c)x¨3+2ζ3ω3x˙3+ω32x3+M¯13x12+M¯14x22+M¯15x32+M¯16x1x2+M¯17x1x3+M¯18x2x3=F3cosΩt,where ζn is the damping parameter of the system, ωn is the one of the natural frequencies of the system, and Mn is the coefficient before each quadratic nonlinear term.


## 4. Linear System

### 4.1. Numerical Analysis of Vibration Characteristics

After obtaining the three-degree-of-freedom ordinary differential equation of the antisymmetric cross-ply bistable thin shell, Figure3 shows the mode shapes of the first 9 antisymmetric cross-ply bistable shallow shell models obtained by the finite element numerical simulation software ABAQUS. It can be seen from the figure that each is symmetric or antisymmetric with respect to the x axis and the y axis. The first mode shape is the biaxially antisymmetrical torsion, the second mode shape is the biaxially symmetrical bending vibration, and the third mode shape is the torsional vibration symmetric with respect to x-axis and antisymmetric with respect to y-axis. Starting from the fourth order, the mode shapes are all bending-torsional coupling vibrations. The fourth mode shape is antisymmetric with respect to x-axis and symmetric with respect to y-axis, the fifth mode shape is symmetric with respect to the two axes, and the sixth mode shape is antisymmetric with respect to the x axis and symmetric with respect to y axis. The seventh mode shape is antisymmetric with respect to the two axes. The eighth mode shape is symmetric with respect to the x axis and antisymmetric with respect to the y axis, and the ninth mode shape is antisymmetric with respect to the two axes. Higher-order mode shapes also satisfy this law.Figure 3
Diagram of mode shapes at the first nine orders of the antisymmetric bistable laminated shallow shell. (a)ω1=14.2 Hz. (b) ω2=30.85 Hz. (c) ω3=42.2 Hz. (d) ω4=80.4 Hz. (e) ω5=85.5 Hz. (f) ω6=94.9 Hz. (g) ω7=95.4 Hz. (h) ω8=158.3 Hz. (i) ω9=163.1 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 4.2. The Impact of Geometric Parameters on the Frequency and Mode Shapes of the Bistable Cylindrical Shell

Since the laying angle of antisymmetrical cross-ply was selected as the alternately laying mode at0∘ and 90∘, simulation modeling and analysis were performed on the first 9 groups of the test specimen in Table 2. By changing the shell length L, the initial crosssection radius R1, the initial cross-section center angle β, and the number of layers n, the impact of different geometric parameters on the frequency and mode shapes of the antisymmetric laminated bistable shell was studied.Table 2
Model parameters of the antisymmetric bistable test specimen with different shapes.
Test specimenL (mm)R1 (mm)β°nα°1300951804452350340043007558569512071508180698104301160Figure4 shows the relationship between different shell length L and the frequency value ω of the first 14 antisymmetric cross-ply laminated cylindrical shell. The black, red, and blue lines represent the three situations where L=300mm, L=350mm, and L=400mm, respectively. It can be seen from the figure that the frequency value increases with the increase of the order of mode shape and, at the same time, increase with the decrease of the shell length. The frequency values before the 7th mode shape are close, and the values of the 10th and 11th mode shapes are close to each other, while the values of other mode shapes are obviously different. The frequency values of the 8th and 9th mode shapes of each model are slightly different, and the curves in the figure tend to be flat. The frequency values of the 12th and 13th mode shapes also have the same law. The small figure gives the details of the frequency values of the first 5 models with three shell lengths.Figure 4
Influence of lengthL on modal frequency ω.Figures5 to 7 present the stress cloud of the first 9 mode shapes of the models with three shell lengths obtained by simulation. It can be seen that the main mode shape at each order of the models with three shell lengths has changed. In the figure, the distribution law of the stress on the shell in each mode shape can be observed. The magnitude of the stress increases with the changes of blue, green, yellow, and red. When there are two straight opposite sides of the mode shape, the maximum stress is distributed on the two straight sides of the shell, and the distribution of minimum values on the shell surface is close to a regular rectangle; when the four sides of the mode shape are all curved sides, the maximum values of stress are distributed on the four corner points, and the distribution of the minimum values is irregular. Only the second, fourth, and fifth mode shapes remain unchanged, while other mode shapes are different. In the first 6 mode shapes, bending vibration and torsional vibration appear alternately. Starting from the 7th order, all mode shapes are bending-torsion coupling vibrations.Figure 5
Diagram of the mode shape at the first eight orders and modal frequency (L = 300 mm). (a) ω1=14.8 Hz. (b) ω2=25.5 Hz. (c) ω3=34.6 Hz. (d) ω4=74 Hz. (e) ω5=82.7 Hz. (f) ω6=156 Hz. (g) ω7=164.2 Hz. (h) ω8=229.7 Hz. (i) ω9=230.1 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 6
Diagram of the mode shape at the first eight orders and modal frequency (L = 350 mm). (a) ω1=13 Hz. (b) ω2=25.5 Hz. (c) ω3=32.5 Hz. (d) ω4=74 Hz. (e) ω5=80.5 Hz. (f) ω6=156 Hz. (g) ω7=162 Hz. (h) ω8=194.4 Hz. (i) ω9=195.5 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 7
Diagram of the mode shape at the first eight orders and modal frequency (L = 400 mm). (a) ω1=11.6 Hz. (b) ω2=25.5 Hz. (c) ω3=31 Hz. (d) ω4=74 Hz. (e) ω5=79.1 Hz. (f) ω6=156.1 Hz. (g) ω7=160.6 Hz. (h) ω8=168.2 Hz. (i) ω9=170.9 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure8 illustrates the impact of different initial cross section radius R1 on the frequency ω of the antisymmetric cross-ply laminated cylindrical shell. The black, red, and blue lines, respectively, represent the curves of the frequency values at the first 14 orders in the three cases where R1=75mm, R1=85mm, and R1=95mm. It can be seen from the figure that the frequency value increases with the increase of the order of mode shape and, at the same time, increases with the decrease of the initial cross-section radius R1. The overall trend in the three cases is similar to that in Figure 3, except that the frequency value increases slowly at the 8th, 9th, 12th, and 13th mode shapes, and the frequency values at other orders show a significant increase. The small figure displays the detailed diagram of the frequency values of the first 5 mode shapes. Figures 9 and 10 present the stress cloud of the first 9 mode shapes when R1=75mm and R1=85mm obtained by simulation, and they were compared with the situation where R1=95mm in Figure 4. In the figure, the stress distribution on the shell in each mode shape can be found, and the stress distribution law is the same as that shown in Figures 5–7. Only the 4th, 5th, and 8th mode shapes remain unchanged, while other mode shapes are different. In the first 6 mode shapes, bending vibration and torsional vibration appear alternately. Starting from the 7th mode shape, all modes are bending-torsion coupling vibrations.Figure 8
Influence of initial radius of curvatureR1 on modal frequency ω.Figure 9
Diagram of the mode shape at the first eight orders and modal frequency (R1= 75 mm). (a) ω1=19.6 Hz. (b) ω2=41.4 Hz. (c) ω3=51.1 Hz. (d) ω4=120.4 Hz. (e) ω5=129.3 Hz. (f) ω6=254 Hz. (g) ω7=254.8 Hz. (h) ω8=257.8 Hz. (i) ω9=262.8 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 10
Diagram of the mode shape at the first eight orders and modal frequency (R1= 85 mm). (a) ω1=17 Hz. (b) ω2=32 Hz. (c) ω3=41.6 Hz. (d) ω4=93.6 Hz. (e) ω5=102.4 Hz. (f) ω6=197.6 Hz. (g) ω7=205.6 Hz. (h) ω8=240.9 Hz. (i) ω9=243.3 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure11 exhibits the relationship between different initial cross-section central angle β and the frequency value ω. The black, red, and blue lines represent the curves of frequency at the first 14 orders when the central angle β=120∘, β=150∘, and β=180∘, respectively. The small figure is the detailed diagram of the frequency values of the 6th to 9th mode shapes. It can be seen that the frequency values of the 6th to 9th mode shapes when β=150∘ are relatively close, and the curves tend to be flat. In other cases, the frequency value increases with the decrease of the central angle. Figures 12 and 13 show the stress cloud of the first 9 mode shapes under the initial cross-section center angle β=120∘ and β=150∘ obtained by simulation, respectively. The distribution rules of the maximum and minimum stress on the shell surface are the same as those described above. Only the first, second, fourth, and eighth mode shapes remain unchanged, while other mode shapes are quite different. In the first 6 mode shapes, bending vibration and torsional vibration appear alternately. Starting from the 7th mode shape, all mode shapes are bending and torsion coupling vibrations.Figure 11
Influence of initial central angleβ on modal frequency ω.Figure 12
Diagram of the mode shape at the first eight orders and modal frequency (α= 120°). (a) ω1=22.5 Hz. (b) ω2=62.7 Hz. (c) ω3=73.4 Hz. (d) ω4=112 Hz. (e) ω5=121 Hz. (f) ω6=228 Hz. (g) ω7=235.7 Hz. (h) ω8=376.5 Hz. (i) ω9=377.3 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 13
Diagram of the mode shape at the first eight orders and modal frequency (α = 150°). (a) ω1=18 Hz. (b) ω2=38.3 Hz. (c) ω3=48.2 Hz. (d) ω4=112 Hz. (e) ω5=121 Hz. (f) ω6=228 Hz. (g) ω7=231.3 Hz. (h) ω8=232.7 Hz. (i) ω9=240.9 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure14 describes the impact of different layer thicknesses n on the frequency of antisymmetric cross-ply laminated cylindrical shell. The black, red, and blue lines represent the curves of frequency at the first 14 orders when the layer thickness n=4, n=6, and n=8. It can be found from the figure that the frequency value increases with the increase of the order of mode and, at the same time, increases with the increase of the layer thickness. The overall trend in these three cases is similar to that in the previous cases, except that the frequency values of the 3rd, 5th, and 9th mode shapes are significantly increased, and the values at other orders all increase slowly. The small figure gives the detailed diagram of the frequency values of first 5 mode shapes. Figures 15 and 16 are the stress cloud of the first 9 mode shapes when n=6 and n=8 obtained by the simulation and were compared with the situation in Figure 5 where n=4. In the figure, the stress distribution on the shell in each mode shape can be seen, and the stress distribution law is the same as that shown in Figures 5–7. Only the first, second, fourth, and eighth mode shapes remain unchanged, while other mode shapes are different. Bending vibration and torsional vibration appear alternately in the first 6 mode shapes. Starting from the 7th mode shape, the mode shapes are all bending-torsion coupling vibrations.Figure 14
Influence of ply thicknessn on modal frequency ω.Figure 15
Diagram of the mode shape at the first eight orders and modal frequency (n = 6). (a) ω1=22.2 Hz. (b) ω2=39.2 Hz. (c) ω3=52.5 Hz. (d) ω4=113.6 Hz. (e) ω5=126.3 Hz. (f) ω6=239.5 Hz. (g) ω7=251.2 Hz. (h) ω8=289.5 Hz. (i) ω9=292.9 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 16
Diagram of the mode shape at the first eight orders and modal frequency (n = 8). (a) ω1=29.6 Hz. (b) ω2=52.6 Hz. (c) ω3=70.3 Hz. (d) ω4=152.6 Hz. (e) ω5=169.3 Hz. (f) ω6=321.7 Hz. (g) ω7=337 Hz. (h) ω8=342.1 Hz. (i) ω9=345.1 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 4.1. Numerical Analysis of Vibration Characteristics

After obtaining the three-degree-of-freedom ordinary differential equation of the antisymmetric cross-ply bistable thin shell, Figure3 shows the mode shapes of the first 9 antisymmetric cross-ply bistable shallow shell models obtained by the finite element numerical simulation software ABAQUS. It can be seen from the figure that each is symmetric or antisymmetric with respect to the x axis and the y axis. The first mode shape is the biaxially antisymmetrical torsion, the second mode shape is the biaxially symmetrical bending vibration, and the third mode shape is the torsional vibration symmetric with respect to x-axis and antisymmetric with respect to y-axis. Starting from the fourth order, the mode shapes are all bending-torsional coupling vibrations. The fourth mode shape is antisymmetric with respect to x-axis and symmetric with respect to y-axis, the fifth mode shape is symmetric with respect to the two axes, and the sixth mode shape is antisymmetric with respect to the x axis and symmetric with respect to y axis. The seventh mode shape is antisymmetric with respect to the two axes. The eighth mode shape is symmetric with respect to the x axis and antisymmetric with respect to the y axis, and the ninth mode shape is antisymmetric with respect to the two axes. Higher-order mode shapes also satisfy this law.Figure 3
Diagram of mode shapes at the first nine orders of the antisymmetric bistable laminated shallow shell. (a)ω1=14.2 Hz. (b) ω2=30.85 Hz. (c) ω3=42.2 Hz. (d) ω4=80.4 Hz. (e) ω5=85.5 Hz. (f) ω6=94.9 Hz. (g) ω7=95.4 Hz. (h) ω8=158.3 Hz. (i) ω9=163.1 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 4.2. The Impact of Geometric Parameters on the Frequency and Mode Shapes of the Bistable Cylindrical Shell

Since the laying angle of antisymmetrical cross-ply was selected as the alternately laying mode at0∘ and 90∘, simulation modeling and analysis were performed on the first 9 groups of the test specimen in Table 2. By changing the shell length L, the initial crosssection radius R1, the initial cross-section center angle β, and the number of layers n, the impact of different geometric parameters on the frequency and mode shapes of the antisymmetric laminated bistable shell was studied.Table 2
Model parameters of the antisymmetric bistable test specimen with different shapes.
Test specimenL (mm)R1 (mm)β°nα°1300951804452350340043007558569512071508180698104301160Figure4 shows the relationship between different shell length L and the frequency value ω of the first 14 antisymmetric cross-ply laminated cylindrical shell. The black, red, and blue lines represent the three situations where L=300mm, L=350mm, and L=400mm, respectively. It can be seen from the figure that the frequency value increases with the increase of the order of mode shape and, at the same time, increase with the decrease of the shell length. The frequency values before the 7th mode shape are close, and the values of the 10th and 11th mode shapes are close to each other, while the values of other mode shapes are obviously different. The frequency values of the 8th and 9th mode shapes of each model are slightly different, and the curves in the figure tend to be flat. The frequency values of the 12th and 13th mode shapes also have the same law. The small figure gives the details of the frequency values of the first 5 models with three shell lengths.Figure 4
Influence of lengthL on modal frequency ω.Figures5 to 7 present the stress cloud of the first 9 mode shapes of the models with three shell lengths obtained by simulation. It can be seen that the main mode shape at each order of the models with three shell lengths has changed. In the figure, the distribution law of the stress on the shell in each mode shape can be observed. The magnitude of the stress increases with the changes of blue, green, yellow, and red. When there are two straight opposite sides of the mode shape, the maximum stress is distributed on the two straight sides of the shell, and the distribution of minimum values on the shell surface is close to a regular rectangle; when the four sides of the mode shape are all curved sides, the maximum values of stress are distributed on the four corner points, and the distribution of the minimum values is irregular. Only the second, fourth, and fifth mode shapes remain unchanged, while other mode shapes are different. In the first 6 mode shapes, bending vibration and torsional vibration appear alternately. Starting from the 7th order, all mode shapes are bending-torsion coupling vibrations.Figure 5
Diagram of the mode shape at the first eight orders and modal frequency (L = 300 mm). (a) ω1=14.8 Hz. (b) ω2=25.5 Hz. (c) ω3=34.6 Hz. (d) ω4=74 Hz. (e) ω5=82.7 Hz. (f) ω6=156 Hz. (g) ω7=164.2 Hz. (h) ω8=229.7 Hz. (i) ω9=230.1 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 6
Diagram of the mode shape at the first eight orders and modal frequency (L = 350 mm). (a) ω1=13 Hz. (b) ω2=25.5 Hz. (c) ω3=32.5 Hz. (d) ω4=74 Hz. (e) ω5=80.5 Hz. (f) ω6=156 Hz. (g) ω7=162 Hz. (h) ω8=194.4 Hz. (i) ω9=195.5 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 7
Diagram of the mode shape at the first eight orders and modal frequency (L = 400 mm). (a) ω1=11.6 Hz. (b) ω2=25.5 Hz. (c) ω3=31 Hz. (d) ω4=74 Hz. (e) ω5=79.1 Hz. (f) ω6=156.1 Hz. (g) ω7=160.6 Hz. (h) ω8=168.2 Hz. (i) ω9=170.9 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure8 illustrates the impact of different initial cross section radius R1 on the frequency ω of the antisymmetric cross-ply laminated cylindrical shell. The black, red, and blue lines, respectively, represent the curves of the frequency values at the first 14 orders in the three cases where R1=75mm, R1=85mm, and R1=95mm. It can be seen from the figure that the frequency value increases with the increase of the order of mode shape and, at the same time, increases with the decrease of the initial cross-section radius R1. The overall trend in the three cases is similar to that in Figure 3, except that the frequency value increases slowly at the 8th, 9th, 12th, and 13th mode shapes, and the frequency values at other orders show a significant increase. The small figure displays the detailed diagram of the frequency values of the first 5 mode shapes. Figures 9 and 10 present the stress cloud of the first 9 mode shapes when R1=75mm and R1=85mm obtained by simulation, and they were compared with the situation where R1=95mm in Figure 4. In the figure, the stress distribution on the shell in each mode shape can be found, and the stress distribution law is the same as that shown in Figures 5–7. Only the 4th, 5th, and 8th mode shapes remain unchanged, while other mode shapes are different. In the first 6 mode shapes, bending vibration and torsional vibration appear alternately. Starting from the 7th mode shape, all modes are bending-torsion coupling vibrations.Figure 8
Influence of initial radius of curvatureR1 on modal frequency ω.Figure 9
Diagram of the mode shape at the first eight orders and modal frequency (R1= 75 mm). (a) ω1=19.6 Hz. (b) ω2=41.4 Hz. (c) ω3=51.1 Hz. (d) ω4=120.4 Hz. (e) ω5=129.3 Hz. (f) ω6=254 Hz. (g) ω7=254.8 Hz. (h) ω8=257.8 Hz. (i) ω9=262.8 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 10
Diagram of the mode shape at the first eight orders and modal frequency (R1= 85 mm). (a) ω1=17 Hz. (b) ω2=32 Hz. (c) ω3=41.6 Hz. (d) ω4=93.6 Hz. (e) ω5=102.4 Hz. (f) ω6=197.6 Hz. (g) ω7=205.6 Hz. (h) ω8=240.9 Hz. (i) ω9=243.3 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure11 exhibits the relationship between different initial cross-section central angle β and the frequency value ω. The black, red, and blue lines represent the curves of frequency at the first 14 orders when the central angle β=120∘, β=150∘, and β=180∘, respectively. The small figure is the detailed diagram of the frequency values of the 6th to 9th mode shapes. It can be seen that the frequency values of the 6th to 9th mode shapes when β=150∘ are relatively close, and the curves tend to be flat. In other cases, the frequency value increases with the decrease of the central angle. Figures 12 and 13 show the stress cloud of the first 9 mode shapes under the initial cross-section center angle β=120∘ and β=150∘ obtained by simulation, respectively. The distribution rules of the maximum and minimum stress on the shell surface are the same as those described above. Only the first, second, fourth, and eighth mode shapes remain unchanged, while other mode shapes are quite different. In the first 6 mode shapes, bending vibration and torsional vibration appear alternately. Starting from the 7th mode shape, all mode shapes are bending and torsion coupling vibrations.Figure 11
Influence of initial central angleβ on modal frequency ω.Figure 12
Diagram of the mode shape at the first eight orders and modal frequency (α= 120°). (a) ω1=22.5 Hz. (b) ω2=62.7 Hz. (c) ω3=73.4 Hz. (d) ω4=112 Hz. (e) ω5=121 Hz. (f) ω6=228 Hz. (g) ω7=235.7 Hz. (h) ω8=376.5 Hz. (i) ω9=377.3 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 13
Diagram of the mode shape at the first eight orders and modal frequency (α = 150°). (a) ω1=18 Hz. (b) ω2=38.3 Hz. (c) ω3=48.2 Hz. (d) ω4=112 Hz. (e) ω5=121 Hz. (f) ω6=228 Hz. (g) ω7=231.3 Hz. (h) ω8=232.7 Hz. (i) ω9=240.9 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure14 describes the impact of different layer thicknesses n on the frequency of antisymmetric cross-ply laminated cylindrical shell. The black, red, and blue lines represent the curves of frequency at the first 14 orders when the layer thickness n=4, n=6, and n=8. It can be found from the figure that the frequency value increases with the increase of the order of mode and, at the same time, increases with the increase of the layer thickness. The overall trend in these three cases is similar to that in the previous cases, except that the frequency values of the 3rd, 5th, and 9th mode shapes are significantly increased, and the values at other orders all increase slowly. The small figure gives the detailed diagram of the frequency values of first 5 mode shapes. Figures 15 and 16 are the stress cloud of the first 9 mode shapes when n=6 and n=8 obtained by the simulation and were compared with the situation in Figure 5 where n=4. In the figure, the stress distribution on the shell in each mode shape can be seen, and the stress distribution law is the same as that shown in Figures 5–7. Only the first, second, fourth, and eighth mode shapes remain unchanged, while other mode shapes are different. Bending vibration and torsional vibration appear alternately in the first 6 mode shapes. Starting from the 7th mode shape, the mode shapes are all bending-torsion coupling vibrations.Figure 14
Influence of ply thicknessn on modal frequency ω.Figure 15
Diagram of the mode shape at the first eight orders and modal frequency (n = 6). (a) ω1=22.2 Hz. (b) ω2=39.2 Hz. (c) ω3=52.5 Hz. (d) ω4=113.6 Hz. (e) ω5=126.3 Hz. (f) ω6=239.5 Hz. (g) ω7=251.2 Hz. (h) ω8=289.5 Hz. (i) ω9=292.9 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 16
Diagram of the mode shape at the first eight orders and modal frequency (n = 8). (a) ω1=29.6 Hz. (b) ω2=52.6 Hz. (c) ω3=70.3 Hz. (d) ω4=152.6 Hz. (e) ω5=169.3 Hz. (f) ω6=321.7 Hz. (g) ω7=337 Hz. (h) ω8=342.1 Hz. (i) ω9=345.1 Hz.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 5. Nonlinear System

### 5.1. Perturbation Analysis When the Main ResonanceΩ≈ω1

#### 5.1.1. Noninternal Resonance Situation

According to the frequency of the linear system, there is a 1 : 2 : 3 internal resonance relationship between the natural frequencies of the first three orders in the bistable system. On the premise that the main resonanceΩ is close to ω1, internal resonance situation and noninternal resonance situation were discussed separately.First, regarding noninternal resonance situation, the solvability condition (the condition for eliminating the long-term term) is(19a)2iω1Af′+μ1A1=12f1eiσ1T1,(19b)A2′+μ2A2=0,(19c)A3′+μ3A3=0,where A1′ is the derivative of A1 about T1 and σ1 is the tuning parameter.The steady-state response that can be obtained is(20a)x1=F1sinΩt−γ02εω1μ12+σ12+Οε,(20b)x2=0,(20c)x2=0,where γ0=arctanσ1/μ1. Therefore, in the case where there is no internal resonance, any nonlinear term will not produce a long-term term, and its first approximation is not affected by the nonlinear term; thus, it is actually the solution to the corresponding linear problem.


#### 5.1.2. Internal Resonance Situation

When there is a 1 : 2 : 3 internal resonance, letΩ=ω1+εσ1, ω2=2ω1+εσ2, and ω3=3ω1+εσ3, where σ1, σ2, and σ3 are the tuning parameters, and the solvability condition is(21a)−2iω1A1′+μ1A1−m4A¯1A2eiσ2T1−m6A¯2A3eiσ3−σ2T1+f12eiσ1T1=0,(21b)−2iω2A2′+μ2A2−m7A12e−iσ2T1−m11A¯2A3eiσ3−σ2T1=0,(21c)−2iω3A3′+μ3A3−m16A1A2eiσ2−σ3T1=0.At this time, the symbol was introduced, letAn=aneiθn/2, and the result was divided into a real part and an imaginary part, and we have(22a)a1′=−μ1a1−m4a1a24ω1sinγ2−m6a2a34ω1sinγ2+f12ω1sinγ3,(22b)a1θ1′=m4a1a24ω1cosγ2+m6a2a34ω1cosγ2−f12ω1cosγ3,(22c)a2′=−μ2a2+m7a124ω2sinγ2−m11a1a34ω2sinγ2,(22d)a2θ2′=m7a124ω2cosγ1+m11a1a34ω2cosγ2,(22e)a3′=−μ3a3+m16a1a24ω3sinγ1,(22f)a3θ3′=m16a1a24ω3cosγ2,where(23a)γ1=θ2−2θ1+σ2T1,(23b)γ2=θ3−θ2−θ1+σ3−σ2T1,(23c)γ3=σ1T1−θ1.Equations (22a)–(22f) shows that the phase modulation between response and excitation begins to change, which is not available in linear systems. The aforementioned phase modulation change is important since it limits the amplitude of the response to a finite value. θ1 and θ2 were eliminated to obtain the average equation under polar coordinates:(24a)D1a1=−μ1a1−m4a1a24ω1sinγ2−m6a2a34ω1sinγ2+f12ω1sinγ3,(24b)a1D1γ1=12m7a122ω2−m4a22ω1cosγ1+a32m11a12ω2−m6a22ω1a1cosγ2+f1a2ω1a1cosγ3,(24c)D1a2=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2,(24d)a2D1γ2=σ3−σ2a2−14m7a12ω2+m4a22ω1cosγ1+f1a22ω1a1cosγ3−a34m11a1ω2+m6a22ω1a1−m16a1a3ω3cosγ2,(24e)D1a3=−μ3a3+m16a1a24ω3sinγ1,(24f)a3D1γ3=σ1a1−m4a1a24ω1cosγ1−m6a2a34ω1cosγ2+f12ω1cosγ3.LetA1=x1+ix2, A2=x3+ix4, and A4=x5+ix6. The average equation under rectangular coordinates is as follows:(25a)x˙1=−μ1x1+m42ω1x2x3−x1x4+m62ω1x4x5−x3x6,(25b)x˙2=−μ1x2+m42ω1x1x3+x2x4+m62ω1x3x5+x4x6−f14ω1,(25c)x˙3=−μ2x3−m7ω2x1x2+m112ω2x1x6−x2x5,(25d)x˙4=−μ2x4−m72ω2x12−x22+m112ω2x1x5+x2x6,(25e)x˙5=−μ3x5−m162ω3x2x3+x1x4,(25f)x˙6=−μ3x6+m162ω3x1x3−x2x4.


### 5.2. Dynamic Analysis of Bifurcation and Chaos When the Main ResonanceΩ≈ω1

Based on rectangular coordinate average equations (25a)–(25f) of the bistable system under 1 : 2 : 3 internal resonance and main resonance, the nonlinear dynamic behavior characteristics of the bistable system with the change of external excitation were analyzed by numerical methods. First, the initial conditions were defined as x1=0.11, x2=0.39, x3=0.21, x4=0.19, x5=0.11, and x6=0.22, and other parameters were set to be μ1=0.1, μ2=0.12, μ3=0.5, m4=16.4, m6=14.5, m7=10, m11=9.1, and m16=7. In rectangular coordinate average equation (25a), f1 is the parameter related to the external excitation. Here, f1 was selected as the control parameter to study the nonlinear dynamic characteristics of the system affected by the external excitation.Figure17 is the bifurcation diagram of the second and third mode shapes of the antisymmetric cross-ply bistable system. The abscissa axis represents the changes of f1 related to the external excitation, and the ordinate axis represents the changes of x3 and x5 related to the second-order and third-order amplitudes of the system. As the external excitation expands, the bifurcation diagram of the system shows obvious up-and-down vibrations, indicating that the nonlinear dynamic behavior of the antisymmetric cross-ply laminated bistable cylindrical shell has the basic characteristics of the bistable state. In the interval f1∈0,2, the bistable structure is in a period-doubling state. In the interval f1∈2,80, the bistable structure is in chaotic motion, with the amplitude constantly expanding.Figure 17
Bifurcation diagram of the antisymmetric cross-ply bistable system. (a) First-order modal. (b) Second-order modal.
(a)(b)Figure18 is the maximum Lyapunov exponent diagram of the antisymmetric cross-ply laminated bistable system. The abscissa axis represents the changes of f1 related to external excitation, and the ordinate axis represents the maximum Lyapunov exponent. It can be seen from the figure that, in the interval f1∈0,2, the maximum Lyapunov exponent is less than zero. In the interval f1∈2,80, the maximum Lyapunov exponent is always greater than zero, indicating that the bistable structure has been in a state of chaos. Based on the selected value of f1, the waveform diagram, phase diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the system under different external excitation values are displayed.Figure 18
The maximum Lyapunov exponent diagram of the antisymmetric cross-ply bistable system.Figures19 to 22 are the phase diagram, time history diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the corresponding system motion when the parameter values of the external excitation f1 are different. In each figure, (a) and (b), respectively, show the three-dimensional phase diagram in the space x4,x5,x6 and the time history diagram in t,x5. (c) and (d), respectively, show the three-dimensional phase diagram in space x1,x2,x3 and the second mode shape phase diagram on the two-dimensional plane x3,x4. (e) and (f), respectively, show the Poincaré cross-section and frequency spectrum on the two-dimensional plane x1,x2.Figure 19
The periodic motion whenf1=1. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 20
The quasi-periodic motion whenf1=20. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 21
The chaotic motion whenf1=40. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum f1=40.
(a)(b)(c)(d)(e)(f)Figure 22
The chaotic motion whenf1=60. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Whenf1=1, it can be seen from Figure 19 that a periodic phenomenon occurs in the system. It can be found from the figure that a single closed loop appears in the two-dimensional and three-dimensional phase diagrams of the system, the waveform diagram shows obvious periodic characteristics, a mapping point appears in the Poincaré mapping, and there is a main frequency component in the spectrogram. Moreover, the largest Lyapunov exponent is close to zero. These characteristics indicate that the system has periodic movement at this time.Whenf1=20, it can be seen from Figure 20 that an almost periodic phenomenon occurs in the system. It can be found from the figure that the two-dimensional and three-dimensional phase diagrams of the system show dense circular orbits, dense mapping points appear in the Poincaré map, the spectrogram is a discrete spectrum, and the maximum Lyapunov exponent is close to zero. These characteristics reveal that the system has almost periodic motion at this time.Figure21 shows the motion characteristics of the system when f1=40. The phase diagram of the system is messy, dense, and clearly bounded. There are continuous spectra in the spectrogram, countless mapping points appear in the Poincaré map, and the maximum Lyapunov exponent is greater than zero. These features demonstrate that chaotic movement occurs in the system at this time.Figure22 shows the motion characteristics of the system when f1=60. The phase diagram trajectory of the system is messy, dense, and clearly bounded. The waveform diagram shows that the system motion has random characteristics. Continuous spectra appear in the spectrogram, and there are countless mapping points in the Poincaré map. Meanwhile, the maximum Lyapunov exponent is greater than zero. These features indicate that chaotic motion occurs in the system at this time.


### 5.3. Perturbation Analysis When the Main ResonanceΩ≈ω2

#### 5.3.1. Noninternal Resonance Situation

First, regarding the situation of noninternal resonance, the solvability condition (the condition for eliminating the long-term term) is(26a)2iω1A1′+μ1A1=0,(26b)2iω2A2′+μ2A2=f22eiσ1T1,(26c)2iω3A3′+μ3A3=0.The resulting solution expression is(27a)A1=12a1e−μ1T1+iθ1,(27b)A2=12a2e−μ2T1+iθ2−if2eiσ1T14ω2μ2+iσ1,(27c)A3=12a3e−μ3T1+iθ3.Whent⟶∞ and T1⟶∞, we have(28a)A1⟶0,(28b)A2⟶if2eiσ1T14ω2μ2+iσ1,(28c)A3⟶0.Therefore, in the case where there is no internal resonance, any nonlinear term will not produce a long-term term, and its first-order approximation is not affected by the nonlinear term; thus, it is actually the solution of the corresponding linear problem.


#### 5.3.2. Internal Resonance Situation

When there is a 1 : 2 : 3 internal resonance, letΩ=ω2+εσ1, ω2=2ω1+εσ2, and ω3=2ω1+εσ3, where σ1, σ2, and σ3 are the tuning parameters, and the solvability condition is(29a)−2iω1A1′+μ1A1−m4A¯1A2eiσ2T1−m6A¯2A3eiσ3−σ2T1=0,(29b)−2iω2A2′+μ2A2−m7A12e−iσ2T1−m11A¯1A3eiσ3−σ2T1+f22eiσ1T1=0,(29c)−2iω3A3′+μ3A3−m16A1A2eiσ2−σ3T1=0.At this time, the symbol was introduced, letAn=aneiθn/2, and we have(30a)An′=12an′eiθn+i2anθn′eiθn,(30b)A¯n=12ane−iθn.The result was divided into a real part and an imaginary part:(31a)a1′=−μ1a1−m4a1a24ω1sinγ1−m6a2a34ω1sinγ2,(31b)a1θ1′=m4a1a24ω1cosγ1+m6a2a34ω1cosγ2,(31c)a2′=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,(31d)a2θ2′=m7a124ω2cosγ1+m11a1a34ω2cosγ2−f22ω2cosγ3,(31e)a3′=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,=−μ3a3+m16a1a24ω3sinγ2,(31f)a3θ3′=m16a1a24ω3cosγ2,where(32a)γ1=θ2−2θ1+σ2T1,(32b)γ2=θ3−θ2−θ1+σ3−σ2T1,(32c)γ3=σ2T1−θ2.The phase modulation between response and excitation begins to change, which is not available in linear systems. The aforementioned phase modulation change is important since it limits the amplitude of the response to a finite value.θ1 and θ2 were eliminated to obtain the average equation under polar coordinates:(33a)D1a1=−μ1a1−m4a1a24ω1sinγ1−m6a2a34ω1sinγ2,(33b)D1γ1=σ2+m7a124ω2cosγ1+m11a1a34ω2cosγ2−f22ω2a2cosγ3,(33c)D1a2=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,(33d)D1γ2=σ3−σ2a2−14m4a22ω1+m7a12ω2cosγ1−14m6a22a3ω1a1+m11a1a3ω2+m16a1a2ω3cosγ2+f22ω2γ3,(33e)D1a3=−μ3a3+m16a1a24ω3sinγ2,(33f)D1γ3=σ1−m7a124ω2a2cosγ1−m11a1a34ω2a2cosγ2+f22ω2a2cosγ3.LetA1=x1+ix2, A2=x3+ix4, and A3=x5+ix6. The average equation under rectangular coordinates was obtained, as follows:(34a)x˙1=−μ1x1−m42ω1x1x4−x2x3−m62ω1x3x6−x4x5,(34b)x˙2=−μ1x2+m42ω1x1x3+x2x4+m62ω1x3x5+x4x6,(34c)x˙3=−μ2x3−m7ω2x1x2+m112ω2x1x6−x2x5,(34d)x˙4=−μ2x4+m72ω2x12−x22+m112ω2x1x5+x2x6−f24ω2,(34e)x˙5=−μ3x5−m162ω3x2x3+x1x4,(34f)x˙6=−μ3x6+m162ω3x1x3−x2x4.


### 5.4. Dynamic Analysis of Bifurcation and Chaos When the Main ResonanceΩ≈ω2

Next, based on rectangular coordinate average equation (34a34b34c34d34e34f) of the bistable system under 1 : 2 : 3 internal resonance and main resonance, numerical methods were used to analyze the nonlinear dynamic behavior characteristics of the bistable system with the change of external excitation. First, the initial conditions were defined as x1=0.11, x2=0.39, x3=0.21, x4=0.19, x5=0.11, and x6=0.22, and the values of other parameters are μ1=0.1, μ2=0.012, μ3=0.05, m4=16.4, m6=14.5, m7=10, m11=9.1, and m16=7. In rectangular coordinate average equation (34a34b34c34d34e34f), f2 is the parameter related to the external excitation. Here, f2 was selected as the control parameter to study the nonlinear dynamic characteristics of the system affected by the external excitation.Figure23(a) is the bifurcation diagram of the first mode shape of the antisymmetric cross-ply bistable system. The abscissa axis represents the changes of f2 related to external excitation, and the ordinate axis represents the changes of x1 related to the first-order amplitude of the system. As the external excitation increases, the bifurcation diagram of the system shows obvious upper and lower areas, which indicates that the nonlinear dynamic behavior of the antisymmetric laminated bistable cylindrical shell has the basic characteristics of the bistable state. In the interval f2∈0,1.2, the bistable structure has been in a monostable state and transforms between two steady states. In the interval f2∈1.2,4, the bistable structure continues to transform between two steady states, with the amplitude continuously expanding. Figure 23(b) is a bifurcation diagram of the second mode shape of the antisymmetric cross-ply laminated bistable system. The abscissa axis represents the changes of f2 related to external excitation, and the ordinate axis represents the changes of x3 related to the second-order amplitude of the system. With the increase of the external excitation, the bifurcation diagram of this system does not have obvious upper and lower areas. In the interval f2∈0,4, the bistable structure is always in the chaotic state of the two steady states up and down near the equilibrium point.Figure 23
Bifurcation diagram of the antisymmetric cross-ply bistable system. (a) First-order modal. (b) Second-order modal.
(a)(b)Figure24 is the maximum Lyapunov exponent diagram of the antisymmetric cross-ply laminated bistable system. The abscissa axis represents the changes of f2 related to external excitation, and the ordinate axis represents the maximum Lyapunov exponent. It can be seen from the figure that, in the interval f2∈0,4, the maximum Lyapunov exponent is always greater than zero, indicating that the bistable structure has been in a state of chaos. Based on the selected value of f2, the waveform diagram, phase diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the system are given below.Figure 24
The maximum Lyapunov exponent diagram of the cross-ply bistable system.Figures25 to 28 are the phase diagram, time history diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the corresponding system motion when the parameter values of external excitation f2 are different. In each figure, (a) and (b), respectively, show the three-dimensional phase diagram in the space x1,x2,x3 and the time history diagram in t,x3. (c) and (d), respectively, show the three-dimensional phase diagram in space x4,x5,x6 and the second mode shape phase diagram on the two-dimensional plane x3,x4. (e) and (f), respectively, show the Poincaré cross-section and spectrogram on a two-dimensional plane x1,x2.Figure 25
The periodic motion whenf2=0.4. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 26
The quasi-periodic motion whenf2=1.2. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 27
The chaotic motion whenf2=2. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 28
The chaotic motion whenf2=3. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Whenf2=0.4, it can be seen from Figure 25 that a periodic phenomenon occurs in the system. There are multiple closed loops in the two-dimensional and three-dimensional phase diagrams of the system, the waveform diagram has obvious periodic characteristics, multiple mapping points appear in the Poincaré map, and there are multiple main frequency components in the spectrogram. In addition, the maximum Lyapunov exponent is close to zero. These characteristics reveal that the system has periodic movement at this time.Whenf2=1.2, it can be seen from Figure 26 that an almost periodic phenomenon occurs in the system. The two-dimensional and three-dimensional phase diagrams of the system show dense circular orbits, dense mapping points appear in the Poincaré map, the spectrogram is a discrete spectrum, and the maximum Lyapunov exponent is close to zero. Based on the characteristics shown in the figure, it can be confirmed that the system has almost periodic motion at this time.Figure27 shows the motion characteristics of the system when f2=2. The phase diagram of the system is messy, dense, and clearly bounded. There are continuous spectra in the spectrogram, countless mapping points appear in the Poincaré map, and the maximum Lyapunov exponent is greater than zero. These features indicate that chaotic movement occurs in the system.Figure28 shows the motion characteristics of the system when f2=3. The phase diagram of the system is messy, dense, and clearly bounded. The waveform diagram shows that the system motion has random characteristics. Continuous spectra appear in the spectrogram, and there are countless mapping points in the Poincaré map. The maximum Lyapunov exponent exceeds zero. These characteristics indicate that chaotic motion occurs in the system at this time.


## 5.1. Perturbation Analysis When the Main ResonanceΩ≈ω1

### 5.1.1. Noninternal Resonance Situation

According to the frequency of the linear system, there is a 1 : 2 : 3 internal resonance relationship between the natural frequencies of the first three orders in the bistable system. On the premise that the main resonanceΩ is close to ω1, internal resonance situation and noninternal resonance situation were discussed separately.First, regarding noninternal resonance situation, the solvability condition (the condition for eliminating the long-term term) is(19a)2iω1Af′+μ1A1=12f1eiσ1T1,(19b)A2′+μ2A2=0,(19c)A3′+μ3A3=0,where A1′ is the derivative of A1 about T1 and σ1 is the tuning parameter.The steady-state response that can be obtained is(20a)x1=F1sinΩt−γ02εω1μ12+σ12+Οε,(20b)x2=0,(20c)x2=0,where γ0=arctanσ1/μ1. Therefore, in the case where there is no internal resonance, any nonlinear term will not produce a long-term term, and its first approximation is not affected by the nonlinear term; thus, it is actually the solution to the corresponding linear problem.


### 5.1.2. Internal Resonance Situation

When there is a 1 : 2 : 3 internal resonance, letΩ=ω1+εσ1, ω2=2ω1+εσ2, and ω3=3ω1+εσ3, where σ1, σ2, and σ3 are the tuning parameters, and the solvability condition is(21a)−2iω1A1′+μ1A1−m4A¯1A2eiσ2T1−m6A¯2A3eiσ3−σ2T1+f12eiσ1T1=0,(21b)−2iω2A2′+μ2A2−m7A12e−iσ2T1−m11A¯2A3eiσ3−σ2T1=0,(21c)−2iω3A3′+μ3A3−m16A1A2eiσ2−σ3T1=0.At this time, the symbol was introduced, letAn=aneiθn/2, and the result was divided into a real part and an imaginary part, and we have(22a)a1′=−μ1a1−m4a1a24ω1sinγ2−m6a2a34ω1sinγ2+f12ω1sinγ3,(22b)a1θ1′=m4a1a24ω1cosγ2+m6a2a34ω1cosγ2−f12ω1cosγ3,(22c)a2′=−μ2a2+m7a124ω2sinγ2−m11a1a34ω2sinγ2,(22d)a2θ2′=m7a124ω2cosγ1+m11a1a34ω2cosγ2,(22e)a3′=−μ3a3+m16a1a24ω3sinγ1,(22f)a3θ3′=m16a1a24ω3cosγ2,where(23a)γ1=θ2−2θ1+σ2T1,(23b)γ2=θ3−θ2−θ1+σ3−σ2T1,(23c)γ3=σ1T1−θ1.Equations (22a)–(22f) shows that the phase modulation between response and excitation begins to change, which is not available in linear systems. The aforementioned phase modulation change is important since it limits the amplitude of the response to a finite value. θ1 and θ2 were eliminated to obtain the average equation under polar coordinates:(24a)D1a1=−μ1a1−m4a1a24ω1sinγ2−m6a2a34ω1sinγ2+f12ω1sinγ3,(24b)a1D1γ1=12m7a122ω2−m4a22ω1cosγ1+a32m11a12ω2−m6a22ω1a1cosγ2+f1a2ω1a1cosγ3,(24c)D1a2=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2,(24d)a2D1γ2=σ3−σ2a2−14m7a12ω2+m4a22ω1cosγ1+f1a22ω1a1cosγ3−a34m11a1ω2+m6a22ω1a1−m16a1a3ω3cosγ2,(24e)D1a3=−μ3a3+m16a1a24ω3sinγ1,(24f)a3D1γ3=σ1a1−m4a1a24ω1cosγ1−m6a2a34ω1cosγ2+f12ω1cosγ3.LetA1=x1+ix2, A2=x3+ix4, and A4=x5+ix6. The average equation under rectangular coordinates is as follows:(25a)x˙1=−μ1x1+m42ω1x2x3−x1x4+m62ω1x4x5−x3x6,(25b)x˙2=−μ1x2+m42ω1x1x3+x2x4+m62ω1x3x5+x4x6−f14ω1,(25c)x˙3=−μ2x3−m7ω2x1x2+m112ω2x1x6−x2x5,(25d)x˙4=−μ2x4−m72ω2x12−x22+m112ω2x1x5+x2x6,(25e)x˙5=−μ3x5−m162ω3x2x3+x1x4,(25f)x˙6=−μ3x6+m162ω3x1x3−x2x4.


## 5.1.1. Noninternal Resonance Situation

According to the frequency of the linear system, there is a 1 : 2 : 3 internal resonance relationship between the natural frequencies of the first three orders in the bistable system. On the premise that the main resonanceΩ is close to ω1, internal resonance situation and noninternal resonance situation were discussed separately.First, regarding noninternal resonance situation, the solvability condition (the condition for eliminating the long-term term) is(19a)2iω1Af′+μ1A1=12f1eiσ1T1,(19b)A2′+μ2A2=0,(19c)A3′+μ3A3=0,where A1′ is the derivative of A1 about T1 and σ1 is the tuning parameter.The steady-state response that can be obtained is(20a)x1=F1sinΩt−γ02εω1μ12+σ12+Οε,(20b)x2=0,(20c)x2=0,where γ0=arctanσ1/μ1. Therefore, in the case where there is no internal resonance, any nonlinear term will not produce a long-term term, and its first approximation is not affected by the nonlinear term; thus, it is actually the solution to the corresponding linear problem.


## 5.1.2. Internal Resonance Situation

When there is a 1 : 2 : 3 internal resonance, letΩ=ω1+εσ1, ω2=2ω1+εσ2, and ω3=3ω1+εσ3, where σ1, σ2, and σ3 are the tuning parameters, and the solvability condition is(21a)−2iω1A1′+μ1A1−m4A¯1A2eiσ2T1−m6A¯2A3eiσ3−σ2T1+f12eiσ1T1=0,(21b)−2iω2A2′+μ2A2−m7A12e−iσ2T1−m11A¯2A3eiσ3−σ2T1=0,(21c)−2iω3A3′+μ3A3−m16A1A2eiσ2−σ3T1=0.At this time, the symbol was introduced, letAn=aneiθn/2, and the result was divided into a real part and an imaginary part, and we have(22a)a1′=−μ1a1−m4a1a24ω1sinγ2−m6a2a34ω1sinγ2+f12ω1sinγ3,(22b)a1θ1′=m4a1a24ω1cosγ2+m6a2a34ω1cosγ2−f12ω1cosγ3,(22c)a2′=−μ2a2+m7a124ω2sinγ2−m11a1a34ω2sinγ2,(22d)a2θ2′=m7a124ω2cosγ1+m11a1a34ω2cosγ2,(22e)a3′=−μ3a3+m16a1a24ω3sinγ1,(22f)a3θ3′=m16a1a24ω3cosγ2,where(23a)γ1=θ2−2θ1+σ2T1,(23b)γ2=θ3−θ2−θ1+σ3−σ2T1,(23c)γ3=σ1T1−θ1.Equations (22a)–(22f) shows that the phase modulation between response and excitation begins to change, which is not available in linear systems. The aforementioned phase modulation change is important since it limits the amplitude of the response to a finite value. θ1 and θ2 were eliminated to obtain the average equation under polar coordinates:(24a)D1a1=−μ1a1−m4a1a24ω1sinγ2−m6a2a34ω1sinγ2+f12ω1sinγ3,(24b)a1D1γ1=12m7a122ω2−m4a22ω1cosγ1+a32m11a12ω2−m6a22ω1a1cosγ2+f1a2ω1a1cosγ3,(24c)D1a2=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2,(24d)a2D1γ2=σ3−σ2a2−14m7a12ω2+m4a22ω1cosγ1+f1a22ω1a1cosγ3−a34m11a1ω2+m6a22ω1a1−m16a1a3ω3cosγ2,(24e)D1a3=−μ3a3+m16a1a24ω3sinγ1,(24f)a3D1γ3=σ1a1−m4a1a24ω1cosγ1−m6a2a34ω1cosγ2+f12ω1cosγ3.LetA1=x1+ix2, A2=x3+ix4, and A4=x5+ix6. The average equation under rectangular coordinates is as follows:(25a)x˙1=−μ1x1+m42ω1x2x3−x1x4+m62ω1x4x5−x3x6,(25b)x˙2=−μ1x2+m42ω1x1x3+x2x4+m62ω1x3x5+x4x6−f14ω1,(25c)x˙3=−μ2x3−m7ω2x1x2+m112ω2x1x6−x2x5,(25d)x˙4=−μ2x4−m72ω2x12−x22+m112ω2x1x5+x2x6,(25e)x˙5=−μ3x5−m162ω3x2x3+x1x4,(25f)x˙6=−μ3x6+m162ω3x1x3−x2x4.


## 5.2. Dynamic Analysis of Bifurcation and Chaos When the Main ResonanceΩ≈ω1

Based on rectangular coordinate average equations (25a)–(25f) of the bistable system under 1 : 2 : 3 internal resonance and main resonance, the nonlinear dynamic behavior characteristics of the bistable system with the change of external excitation were analyzed by numerical methods. First, the initial conditions were defined as x1=0.11, x2=0.39, x3=0.21, x4=0.19, x5=0.11, and x6=0.22, and other parameters were set to be μ1=0.1, μ2=0.12, μ3=0.5, m4=16.4, m6=14.5, m7=10, m11=9.1, and m16=7. In rectangular coordinate average equation (25a), f1 is the parameter related to the external excitation. Here, f1 was selected as the control parameter to study the nonlinear dynamic characteristics of the system affected by the external excitation.Figure17 is the bifurcation diagram of the second and third mode shapes of the antisymmetric cross-ply bistable system. The abscissa axis represents the changes of f1 related to the external excitation, and the ordinate axis represents the changes of x3 and x5 related to the second-order and third-order amplitudes of the system. As the external excitation expands, the bifurcation diagram of the system shows obvious up-and-down vibrations, indicating that the nonlinear dynamic behavior of the antisymmetric cross-ply laminated bistable cylindrical shell has the basic characteristics of the bistable state. In the interval f1∈0,2, the bistable structure is in a period-doubling state. In the interval f1∈2,80, the bistable structure is in chaotic motion, with the amplitude constantly expanding.Figure 17
Bifurcation diagram of the antisymmetric cross-ply bistable system. (a) First-order modal. (b) Second-order modal.
(a)(b)Figure18 is the maximum Lyapunov exponent diagram of the antisymmetric cross-ply laminated bistable system. The abscissa axis represents the changes of f1 related to external excitation, and the ordinate axis represents the maximum Lyapunov exponent. It can be seen from the figure that, in the interval f1∈0,2, the maximum Lyapunov exponent is less than zero. In the interval f1∈2,80, the maximum Lyapunov exponent is always greater than zero, indicating that the bistable structure has been in a state of chaos. Based on the selected value of f1, the waveform diagram, phase diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the system under different external excitation values are displayed.Figure 18
The maximum Lyapunov exponent diagram of the antisymmetric cross-ply bistable system.Figures19 to 22 are the phase diagram, time history diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the corresponding system motion when the parameter values of the external excitation f1 are different. In each figure, (a) and (b), respectively, show the three-dimensional phase diagram in the space x4,x5,x6 and the time history diagram in t,x5. (c) and (d), respectively, show the three-dimensional phase diagram in space x1,x2,x3 and the second mode shape phase diagram on the two-dimensional plane x3,x4. (e) and (f), respectively, show the Poincaré cross-section and frequency spectrum on the two-dimensional plane x1,x2.Figure 19
The periodic motion whenf1=1. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 20
The quasi-periodic motion whenf1=20. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 21
The chaotic motion whenf1=40. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum f1=40.
(a)(b)(c)(d)(e)(f)Figure 22
The chaotic motion whenf1=60. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Whenf1=1, it can be seen from Figure 19 that a periodic phenomenon occurs in the system. It can be found from the figure that a single closed loop appears in the two-dimensional and three-dimensional phase diagrams of the system, the waveform diagram shows obvious periodic characteristics, a mapping point appears in the Poincaré mapping, and there is a main frequency component in the spectrogram. Moreover, the largest Lyapunov exponent is close to zero. These characteristics indicate that the system has periodic movement at this time.Whenf1=20, it can be seen from Figure 20 that an almost periodic phenomenon occurs in the system. It can be found from the figure that the two-dimensional and three-dimensional phase diagrams of the system show dense circular orbits, dense mapping points appear in the Poincaré map, the spectrogram is a discrete spectrum, and the maximum Lyapunov exponent is close to zero. These characteristics reveal that the system has almost periodic motion at this time.Figure21 shows the motion characteristics of the system when f1=40. The phase diagram of the system is messy, dense, and clearly bounded. There are continuous spectra in the spectrogram, countless mapping points appear in the Poincaré map, and the maximum Lyapunov exponent is greater than zero. These features demonstrate that chaotic movement occurs in the system at this time.Figure22 shows the motion characteristics of the system when f1=60. The phase diagram trajectory of the system is messy, dense, and clearly bounded. The waveform diagram shows that the system motion has random characteristics. Continuous spectra appear in the spectrogram, and there are countless mapping points in the Poincaré map. Meanwhile, the maximum Lyapunov exponent is greater than zero. These features indicate that chaotic motion occurs in the system at this time.


## 5.3. Perturbation Analysis When the Main ResonanceΩ≈ω2

### 5.3.1. Noninternal Resonance Situation

First, regarding the situation of noninternal resonance, the solvability condition (the condition for eliminating the long-term term) is(26a)2iω1A1′+μ1A1=0,(26b)2iω2A2′+μ2A2=f22eiσ1T1,(26c)2iω3A3′+μ3A3=0.The resulting solution expression is(27a)A1=12a1e−μ1T1+iθ1,(27b)A2=12a2e−μ2T1+iθ2−if2eiσ1T14ω2μ2+iσ1,(27c)A3=12a3e−μ3T1+iθ3.Whent⟶∞ and T1⟶∞, we have(28a)A1⟶0,(28b)A2⟶if2eiσ1T14ω2μ2+iσ1,(28c)A3⟶0.Therefore, in the case where there is no internal resonance, any nonlinear term will not produce a long-term term, and its first-order approximation is not affected by the nonlinear term; thus, it is actually the solution of the corresponding linear problem.


### 5.3.2. Internal Resonance Situation

When there is a 1 : 2 : 3 internal resonance, letΩ=ω2+εσ1, ω2=2ω1+εσ2, and ω3=2ω1+εσ3, where σ1, σ2, and σ3 are the tuning parameters, and the solvability condition is(29a)−2iω1A1′+μ1A1−m4A¯1A2eiσ2T1−m6A¯2A3eiσ3−σ2T1=0,(29b)−2iω2A2′+μ2A2−m7A12e−iσ2T1−m11A¯1A3eiσ3−σ2T1+f22eiσ1T1=0,(29c)−2iω3A3′+μ3A3−m16A1A2eiσ2−σ3T1=0.At this time, the symbol was introduced, letAn=aneiθn/2, and we have(30a)An′=12an′eiθn+i2anθn′eiθn,(30b)A¯n=12ane−iθn.The result was divided into a real part and an imaginary part:(31a)a1′=−μ1a1−m4a1a24ω1sinγ1−m6a2a34ω1sinγ2,(31b)a1θ1′=m4a1a24ω1cosγ1+m6a2a34ω1cosγ2,(31c)a2′=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,(31d)a2θ2′=m7a124ω2cosγ1+m11a1a34ω2cosγ2−f22ω2cosγ3,(31e)a3′=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,=−μ3a3+m16a1a24ω3sinγ2,(31f)a3θ3′=m16a1a24ω3cosγ2,where(32a)γ1=θ2−2θ1+σ2T1,(32b)γ2=θ3−θ2−θ1+σ3−σ2T1,(32c)γ3=σ2T1−θ2.The phase modulation between response and excitation begins to change, which is not available in linear systems. The aforementioned phase modulation change is important since it limits the amplitude of the response to a finite value.θ1 and θ2 were eliminated to obtain the average equation under polar coordinates:(33a)D1a1=−μ1a1−m4a1a24ω1sinγ1−m6a2a34ω1sinγ2,(33b)D1γ1=σ2+m7a124ω2cosγ1+m11a1a34ω2cosγ2−f22ω2a2cosγ3,(33c)D1a2=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,(33d)D1γ2=σ3−σ2a2−14m4a22ω1+m7a12ω2cosγ1−14m6a22a3ω1a1+m11a1a3ω2+m16a1a2ω3cosγ2+f22ω2γ3,(33e)D1a3=−μ3a3+m16a1a24ω3sinγ2,(33f)D1γ3=σ1−m7a124ω2a2cosγ1−m11a1a34ω2a2cosγ2+f22ω2a2cosγ3.LetA1=x1+ix2, A2=x3+ix4, and A3=x5+ix6. The average equation under rectangular coordinates was obtained, as follows:(34a)x˙1=−μ1x1−m42ω1x1x4−x2x3−m62ω1x3x6−x4x5,(34b)x˙2=−μ1x2+m42ω1x1x3+x2x4+m62ω1x3x5+x4x6,(34c)x˙3=−μ2x3−m7ω2x1x2+m112ω2x1x6−x2x5,(34d)x˙4=−μ2x4+m72ω2x12−x22+m112ω2x1x5+x2x6−f24ω2,(34e)x˙5=−μ3x5−m162ω3x2x3+x1x4,(34f)x˙6=−μ3x6+m162ω3x1x3−x2x4.


## 5.3.1. Noninternal Resonance Situation

First, regarding the situation of noninternal resonance, the solvability condition (the condition for eliminating the long-term term) is(26a)2iω1A1′+μ1A1=0,(26b)2iω2A2′+μ2A2=f22eiσ1T1,(26c)2iω3A3′+μ3A3=0.The resulting solution expression is(27a)A1=12a1e−μ1T1+iθ1,(27b)A2=12a2e−μ2T1+iθ2−if2eiσ1T14ω2μ2+iσ1,(27c)A3=12a3e−μ3T1+iθ3.Whent⟶∞ and T1⟶∞, we have(28a)A1⟶0,(28b)A2⟶if2eiσ1T14ω2μ2+iσ1,(28c)A3⟶0.Therefore, in the case where there is no internal resonance, any nonlinear term will not produce a long-term term, and its first-order approximation is not affected by the nonlinear term; thus, it is actually the solution of the corresponding linear problem.


## 5.3.2. Internal Resonance Situation

When there is a 1 : 2 : 3 internal resonance, letΩ=ω2+εσ1, ω2=2ω1+εσ2, and ω3=2ω1+εσ3, where σ1, σ2, and σ3 are the tuning parameters, and the solvability condition is(29a)−2iω1A1′+μ1A1−m4A¯1A2eiσ2T1−m6A¯2A3eiσ3−σ2T1=0,(29b)−2iω2A2′+μ2A2−m7A12e−iσ2T1−m11A¯1A3eiσ3−σ2T1+f22eiσ1T1=0,(29c)−2iω3A3′+μ3A3−m16A1A2eiσ2−σ3T1=0.At this time, the symbol was introduced, letAn=aneiθn/2, and we have(30a)An′=12an′eiθn+i2anθn′eiθn,(30b)A¯n=12ane−iθn.The result was divided into a real part and an imaginary part:(31a)a1′=−μ1a1−m4a1a24ω1sinγ1−m6a2a34ω1sinγ2,(31b)a1θ1′=m4a1a24ω1cosγ1+m6a2a34ω1cosγ2,(31c)a2′=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,(31d)a2θ2′=m7a124ω2cosγ1+m11a1a34ω2cosγ2−f22ω2cosγ3,(31e)a3′=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,=−μ3a3+m16a1a24ω3sinγ2,(31f)a3θ3′=m16a1a24ω3cosγ2,where(32a)γ1=θ2−2θ1+σ2T1,(32b)γ2=θ3−θ2−θ1+σ3−σ2T1,(32c)γ3=σ2T1−θ2.The phase modulation between response and excitation begins to change, which is not available in linear systems. The aforementioned phase modulation change is important since it limits the amplitude of the response to a finite value.θ1 and θ2 were eliminated to obtain the average equation under polar coordinates:(33a)D1a1=−μ1a1−m4a1a24ω1sinγ1−m6a2a34ω1sinγ2,(33b)D1γ1=σ2+m7a124ω2cosγ1+m11a1a34ω2cosγ2−f22ω2a2cosγ3,(33c)D1a2=−μ2a2+m7a124ω2sinγ1−m11a1a34ω2sinγ2+f22ω2sinγ3,(33d)D1γ2=σ3−σ2a2−14m4a22ω1+m7a12ω2cosγ1−14m6a22a3ω1a1+m11a1a3ω2+m16a1a2ω3cosγ2+f22ω2γ3,(33e)D1a3=−μ3a3+m16a1a24ω3sinγ2,(33f)D1γ3=σ1−m7a124ω2a2cosγ1−m11a1a34ω2a2cosγ2+f22ω2a2cosγ3.LetA1=x1+ix2, A2=x3+ix4, and A3=x5+ix6. The average equation under rectangular coordinates was obtained, as follows:(34a)x˙1=−μ1x1−m42ω1x1x4−x2x3−m62ω1x3x6−x4x5,(34b)x˙2=−μ1x2+m42ω1x1x3+x2x4+m62ω1x3x5+x4x6,(34c)x˙3=−μ2x3−m7ω2x1x2+m112ω2x1x6−x2x5,(34d)x˙4=−μ2x4+m72ω2x12−x22+m112ω2x1x5+x2x6−f24ω2,(34e)x˙5=−μ3x5−m162ω3x2x3+x1x4,(34f)x˙6=−μ3x6+m162ω3x1x3−x2x4.


## 5.4. Dynamic Analysis of Bifurcation and Chaos When the Main ResonanceΩ≈ω2

Next, based on rectangular coordinate average equation (34a34b34c34d34e34f) of the bistable system under 1 : 2 : 3 internal resonance and main resonance, numerical methods were used to analyze the nonlinear dynamic behavior characteristics of the bistable system with the change of external excitation. First, the initial conditions were defined as x1=0.11, x2=0.39, x3=0.21, x4=0.19, x5=0.11, and x6=0.22, and the values of other parameters are μ1=0.1, μ2=0.012, μ3=0.05, m4=16.4, m6=14.5, m7=10, m11=9.1, and m16=7. In rectangular coordinate average equation (34a34b34c34d34e34f), f2 is the parameter related to the external excitation. Here, f2 was selected as the control parameter to study the nonlinear dynamic characteristics of the system affected by the external excitation.Figure23(a) is the bifurcation diagram of the first mode shape of the antisymmetric cross-ply bistable system. The abscissa axis represents the changes of f2 related to external excitation, and the ordinate axis represents the changes of x1 related to the first-order amplitude of the system. As the external excitation increases, the bifurcation diagram of the system shows obvious upper and lower areas, which indicates that the nonlinear dynamic behavior of the antisymmetric laminated bistable cylindrical shell has the basic characteristics of the bistable state. In the interval f2∈0,1.2, the bistable structure has been in a monostable state and transforms between two steady states. In the interval f2∈1.2,4, the bistable structure continues to transform between two steady states, with the amplitude continuously expanding. Figure 23(b) is a bifurcation diagram of the second mode shape of the antisymmetric cross-ply laminated bistable system. The abscissa axis represents the changes of f2 related to external excitation, and the ordinate axis represents the changes of x3 related to the second-order amplitude of the system. With the increase of the external excitation, the bifurcation diagram of this system does not have obvious upper and lower areas. In the interval f2∈0,4, the bistable structure is always in the chaotic state of the two steady states up and down near the equilibrium point.Figure 23
Bifurcation diagram of the antisymmetric cross-ply bistable system. (a) First-order modal. (b) Second-order modal.
(a)(b)Figure24 is the maximum Lyapunov exponent diagram of the antisymmetric cross-ply laminated bistable system. The abscissa axis represents the changes of f2 related to external excitation, and the ordinate axis represents the maximum Lyapunov exponent. It can be seen from the figure that, in the interval f2∈0,4, the maximum Lyapunov exponent is always greater than zero, indicating that the bistable structure has been in a state of chaos. Based on the selected value of f2, the waveform diagram, phase diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the system are given below.Figure 24
The maximum Lyapunov exponent diagram of the cross-ply bistable system.Figures25 to 28 are the phase diagram, time history diagram, Poincaré cross-section diagram, and frequency spectrum diagram of the corresponding system motion when the parameter values of external excitation f2 are different. In each figure, (a) and (b), respectively, show the three-dimensional phase diagram in the space x1,x2,x3 and the time history diagram in t,x3. (c) and (d), respectively, show the three-dimensional phase diagram in space x4,x5,x6 and the second mode shape phase diagram on the two-dimensional plane x3,x4. (e) and (f), respectively, show the Poincaré cross-section and spectrogram on a two-dimensional plane x1,x2.Figure 25
The periodic motion whenf2=0.4. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 26
The quasi-periodic motion whenf2=1.2. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 27
The chaotic motion whenf2=2. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Figure 28
The chaotic motion whenf2=3. (a) The 3D phase portrait. (b) The time history. (c) The 3D phase portrait. (d) The phase portrait on the plane. (e) The Poincaré map. (f) The frequency spectrum.
(a)(b)(c)(d)(e)(f)Whenf2=0.4, it can be seen from Figure 25 that a periodic phenomenon occurs in the system. There are multiple closed loops in the two-dimensional and three-dimensional phase diagrams of the system, the waveform diagram has obvious periodic characteristics, multiple mapping points appear in the Poincaré map, and there are multiple main frequency components in the spectrogram. In addition, the maximum Lyapunov exponent is close to zero. These characteristics reveal that the system has periodic movement at this time.Whenf2=1.2, it can be seen from Figure 26 that an almost periodic phenomenon occurs in the system. The two-dimensional and three-dimensional phase diagrams of the system show dense circular orbits, dense mapping points appear in the Poincaré map, the spectrogram is a discrete spectrum, and the maximum Lyapunov exponent is close to zero. Based on the characteristics shown in the figure, it can be confirmed that the system has almost periodic motion at this time.Figure27 shows the motion characteristics of the system when f2=2. The phase diagram of the system is messy, dense, and clearly bounded. There are continuous spectra in the spectrogram, countless mapping points appear in the Poincaré map, and the maximum Lyapunov exponent is greater than zero. These features indicate that chaotic movement occurs in the system.Figure28 shows the motion characteristics of the system when f2=3. The phase diagram of the system is messy, dense, and clearly bounded. The waveform diagram shows that the system motion has random characteristics. Continuous spectra appear in the spectrogram, and there are countless mapping points in the Poincaré map. The maximum Lyapunov exponent exceeds zero. These characteristics indicate that chaotic motion occurs in the system at this time.


## 6. Conclusion

First, the relevant vibration characteristics of the antisymmetric cross-ply laminated bistable shallow shell were studied. The finite element software ABAQUS was used to perform numerical simulation on the linear model to explore the influence of different geometric parameters on the vibration characteristics such as frequency and main mode shapes. It can be found that the frequency value increases with the decrease of the model size and increases with the increase of the structure thickness and the laying angle, exhibiting stepped growth overall. The main mode shape of models in each group is similar overall, with little change. The distribution of stress on the shell structure is regular. When there are two straight opposite sides of the mode shape, the maximum values of stress are distributed on the two straight sides of the shell, and the distribution of the minimum values on the shell surface is close to a regular rectangle; when all the four sides of the mode shape are curved, the maximum stress is distributed at the four corner points, and the distribution of the minimum values is irregular. Besides, bending vibration and torsional vibration appear alternately in the first 6 mode shapes, and from the 7th order, the mode shapes are all coupled bending and torsion vibrations.Secondly, the dynamic equilibrium equation and the compatible equation were simultaneously established to build a low-order nonlinear dynamic model for the antisymmetric cross-ply laminated bistable shallow shell. Perturbation analysis was conducted on the three-degree-of-freedom nonlinear ordinary differential equation of the antisymmetric laminated bistable cylindrical shallow shell under steady-state excitation through the multiscale method, and the average equations in the six-dimensional polar coordinate form and the rectangular coordinate form under 1 : 2 : 3 internal resonance and the combination of the two different main resonances were obtained. The impact of excitation amplitude on the bifurcation behavior of the system was studied, and the following conclusions were drawn: first, from the bifurcation diagram of the system, it can be seen that there are typical nonlinear phenomena such as period doubling and chaos in the bistable system. Based on the fixed parameter value, two-dimensional and three-dimensional phase diagrams, time history diagram, frequency spectrum, Poincaré mapping, and maximum Lyapunov exponent were used to analyze the specific motion state of the system and to explain that the external excitation has a significant impact on the dynamic characteristics of the bistable model. In the case where two main resonances exist and have similar parameters, there are huge differences in the variation law of the nonlinear dynamic behavior characteristics of the antisymmetric cross-ply bistable system. With the increase of damping coefficient, tuning parameters, and excitation amplitude, the system presents complex nonlinear dynamic behavior, and several typical phenomena were found, such as the alternate occurrence of period-doubling bifurcation and chaos.



---

*Source: 3303512-2021-09-27.xml*