# Nonlinear Vibrations of Laminated Cross-Ply Composite Cantilever Plate in Subsonic Air Flow

Authors: Gen Liu; Wei Zhang; An Xi
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The nonlinear vibrations and responses of a laminated composite cantilever plate under the subsonic air flow are investigated in this paper. The subsonic air flow around the three-dimensional cantilever rectangle laminated composite plate is considered to be decreasing from the wing root to the wing tip. According to the ideal incompressible fluid flow condition and the Kutta–Joukowski lift theorem, the subsonic aerodynamic lift on the three-dimensional finite length flat wing is calculated by using the Vortex Lattice (VL) method. The finite length flat wing is modeled as a laminated composite cantilever plate based on Reddy’s third-order shear deformation plate theory and the von Karman geometry nonlinearity is introduced. The nonlinear partial differential governing equations of motion for the laminated composite cantilever plate subjected to the subsonic aerodynamic force are established via Hamilton’s principle. The Galerkin method is used to separate the partial differential equations into two nonlinear ordinary differential equations, and the four-dimensional nonlinear averaged equations are obtained by the multiple scale method. Through comparing the natural frequencies of the linear system with different material and geometric parameters, the relationship of 1 : 2 internal resonance is considered. Corresponding to several selected parameters, the frequency-response curves are obtained. The hardening-spring-type behaviors and jump phenomena are exhibited. The influence of the force excitation on the bifurcations and chaotic behaviors of the laminated composite cantilever plate is investigated numerically. It is found that the system is sensitive to the exciting force according to the complicate nonlinear behaviors exhibited in this paper.


---

## Body

## 1. Introduction

The vibration of the plate and shell structures owing to air flow is a matter of interest because of its significance in design of launch vehicles and aircrafts [1]. Laminated composite structures are widely used in aerospace field due to its high strength-to-weight ratio, light weight, and long fatigue life. The dynamic behavior of the laminated composite plates in air flow has gained significant focus of attention for many researchers. However, there are few research works dealing with the complex nonlinear dynamics of the structure which is simplified as a laminated composite cantilever plate subjected to subsonic air flow. Therefore, the nonlinear dynamics of laminated composite plates in subsonic flow will be worth analyzing in this work.In the 1990s, the research on the vibration of plates has been studied comprehensively. Some literature reviews on nonlinear vibrations of plates were given by Chia [2, 3] and Mehar and Panda [4]. The nonlinear vibrations of laminated composite spherical shell panels were also entirely investigated by Mahapatra et al. [5–9]. The vibration, bending, and buckling behaviors about the functionally graded sandwich structure have been investigated by Mehar et al. [10–13] and Kar and Panda [14–16]. A third-order theory which accounted for a cubic variation of the in-plane displacements through the plate thickness was derived by Librescu and Reddy [17]. Great progress has also been achieved in the study of the nonlinear thermoelastic frequency analysis [18–22]. Zhang [23] studied the global bifurcations and chaotic dynamics of simply supported rectangular thin plates under parametrical excitation by introducing von Karman’s geometric nonlinearity plate theory.Nonlinear analysis considering the geometrical and material nonlinearity by FEM was conducted by Panda and Singh [24–26]. Onkar and Yadav [27] studied the nonlinear random vibrations of a simply supported cross-ply laminated composite plate using analytical methods, and the accuracy of response evaluation was improved in this work. Park et al. [28] investigated the nonlinear forced vibrations of skew sandwich plates subjected to multiple of dynamic loads, the influence of the skew angles, boundary conditions, and loads on the nonlinear dynamic behaviors of composite plates which were discussed. Chien and Chen [29] studied the effects of initial stresses and different parameters on the nonlinear vibrations of laminated composite plates on elastic foundations. Singh et al. [30] used the higher-order shear deformation plate theory to study the dynamic responses of a geometrically nonlinear laminated composite plate lying on different elastic foundations. Zhang et al. [31, 32] investigated nonlinear vibration, bifurcation, and chaotic dynamics of laminated composite plates. Alijani et al. [33, 34] studied nonlinear vibrations of FGM plates and shallow shells, and the bifurcations and nonlinear dynamic behaviors were analyzed in their papers. Amabili [35] used the higher-order shear deformation theory to study the nonlinear vibrations of angle-ply laminated circular cylindrical shells. Most recently, Zhang et al. [36–39] studied the nonlinear vibrations of the laminated circular cylindrical shells based on the deployed ring truss for antenna structures. The complex nonlinear dynamics of the circular cylindrical shell were studied and the bifurcations and chaotic behaviors were investigated by using the analytical methods and numerical simulations.The research on nonlinear vibrations of plates and shells excited by the aerodynamic force were also widely investigated by researchers. In the early years, Dowell [40, 41] studied the nonlinear oscillations of the fluttering plate systematically. Literature reviews on the nonlinear vibrations of circular cylindrical shells and panels with and without fluid-structure interaction were given by Amabili and Paı¨doussis [42]. The studies on the nonlinear vibration of laminated plates in air flow were mostly about the plates in the supersoninc flow. Singha and Ganapathi [43] investigated the effect of the system parameters on supersonic panel flutter behaviors of laminated composite plates. Haddadpour et al. [44] investigated the nonlinear aeroelastic behaviors of FGM. Singha and Mandal [45] used a 16-node isoparametric degenerated shell element to study the supersonic panel flutter behaviors of laminated composite plates and cylindrical panels. Kuo [46] investigated the effect of variable fiber spacing on the supersonic flutter of rectangular composite plates. Zhao and Zhang [47] presented the analysis of the nonlinear dynamics for a laminated composite cantilever rectangular plate subjected to the supersonic gas flows and the in-plane excitations. An analysis on the nonlinear dynamics of an FGM plate in hypersonic flow subjected to an external excitation and uniform temperature change was presented by Hao et al. [48]. The nonlinear dynamic behavior of an axially extendable cantilever laminated composite plate using piezoelectric materials under the combined action of aerodynamic load and piezoelectric excitation was studied by Lu et al. [49]. Yao and Li [50] investigated the nonlinear vibration of a two-dimensional laminated composite plate in subsonic air flow with simply supported boundary conditions. A simple subsonic aerodynamics model was introduced in this work, which was used to analyze two-dimensional infinite length plate.In this paper, the nonlinear dynamics of the laminated composite cantilever plate under subsonic air flow were investigated. According to the flow condition of ideal incompressible fluid and the Kutta–Joukowski lift theorem, the subsonic aerodynamic lift on the three-dimensional finite length flat wing was calculated using the Vortex Lattice (VL) method. The nonlinear partial differential governing equations of motion for the laminated composite cantilever plate subjected to the subsonic aerodynamic force were established via Hamilton’s principle. The Galerkin method was used to separate the partial differential equation into two nonlinear ordinary differential equations. The numerical method was utilized to investigate the bifurcations and periodic and chaotic motions of the composite laminated rectangular plate. The numerical results illustrate that there existed the periodic and chaotic motions of the composite laminated cantilever plate. The hardening-spring characteristics of the composite laminated cantilever plate were demonstrated by the frequency-response curves of the system.


## 2. Derivation of the Subsonic Aerodynamic Force on the Plate

### 2.1. Vortex Lattice (VL) Method

In this section, the subsonic air flow on the plate will be calculated by using the VL method. The VL method is a numerical implementation on the general 3D lifting surface problem. This method discretizes the vortex-sheet strength distributing on each lifting surface by lumping it into a collection of horseshoe vortices. Figure1 shows that a 3D lifting surface is discretized by a vortex lattice of horseshoe vortices, in which all contribute to the velocity V at any field point r.Figure 1
A 3D lifting surface discretized by a vortex lattice of horseshoe vortices.For the finite wingspan, the high-pressure air below the wing will turn over the low pressure air at the tip of the wing, which will cause the pressure on the upper surface of the wing tip to be equal to that of the lower surface. Unlike the two-dimensional flow around the airfoil, the main characteristic of the three-dimensional flow around the wing is the variation of the lift along the wingspan. In order to calculate the lift on the wing surface by using the vortex lattice method, the Biot–Savart law is used to calculate the induced velocity on the control point. The vortex strength of the vortex system is obtained, and the pressure difference on the upper and lower surface of the wing surface is deduced at the same time.The velocity induced by a vortex line with a strength ofΓn and a length of dl is calculated by the Biot–Savart law as(1)dV=Γndl×r4πr3.As shown in Figure2, the induced velocity is(2)dV=Γnsinθdl4πr2.Figure 2
The velocity induced by finite length vortex segments.LetAB represent a vortex segment in which the vortex vector points from A to B. C is a space point, and its normal distance from the AB line is rp.We can make an integral from the pointA to the point B to find out the size of the induced velocity:(3)V=Γn4πrp∫θ1θ2sinθdθ=Γn4πrpcosθ1−cosθ2.Asr0, r1,and r2, respectively, represent the AB, AC, and BC in Figure 2, the following relationships are established:(4)rp=r1×r2r0,cosθ1=r0⋅r1r0r1,cosθ2=r0⋅r2r0r2.Then, the formula calculating the induced velocity generated by the horseshoe vortex in the vortex lattice method is derived:(5)V=Γn4πr1×r2r1×r2r0⋅r1r1−r2r2.From equation (5), the velocity induced by a vortex segment AB at any point in the coordinate x,y,z can be obtained:(6a)VAB=Γn4πFac1AB⋅Fac2AB,(6b)VA∞=Γn4πz−z1nj+y−y1nkz−z1n2+y−y1n2×1.0+x−x1nx−x1n2+y−y1n2+z−z1n2,(6c)VB∞=−Γn4πz−z2nj+y1n−ykz−z1n2+y2n−y2×1.0+x−x2nx−x2n2+y−y2n2+z−z2n2,where VA∞ is the velocity induced by vortex lines from the point A to ∞ along the x-axis and VB∞ is the velocity induced by vortex lines from the point B to ∞ along the x-axis.The total velocity induced by a horseshoe vortex at a pointx,y,z representing a surface element (i.e., the nth panel element) is the sum of the various components calculated by equations (6a)–(6c). Therefore, the velocity induced by 2N vortexes is obtained to get the total induced velocity on the first m control points. It is expressed as(7)Vm=∑n=12NCm,nΓn.The contribution of all vortices to the downwash of the control point for themth panel element:(8)wm=∑n=12Nwm,n.Then, the lift on the rectangular cantilever plate with a sweep rectangular wing surface is computed according to the VL method. The vortex system is arranged on the rectangular wing surface in Figure3.Figure 3
Layout of the vortex lattice on the flat wing surface.


### 2.2. Application of VL Method on the Plate

In Figure3, the variable b is twice the span of the wing, c is the taper ratio of the wing surface, and c is equal to 1. The sweep angle is 0∘. Five panel elements are specified on the wing surface, and the vortex line is arranged at the leading edge 1/4 chord, and the control point is located at 3/4 chord. The wing surface is divided into 5 panel elements and each panel extends from the leading edge to the trailing edge.The coordinates of the5×1 vortices on the wing surface are listed in Table 1.Table 1
The coordinates of the attachment vortex and the control point in the wing surface.
Panel element number
x
m
y
m
x
1n
y
1n
x
2n
y
2n
1
0.25b
0.05b
0.833b
0b
0.833b
0.1b
2
0.25b
0.15b
0.833b
0.1b
0.833b
0.2b
3
0.25b
0.25b
0.833b
0.2b
0.833b
0.15b
4
0.25b
0.35b
0.833b
0.15b
0.833b
0.4b
5
0.25b
0.45b
0.833b
0.4b
0.833b
0.5bThe downwash velocity of each surface element induced at the control point is superimposed through the nonpenetrating condition:(9)wm=−V∞α.The solution of vortex strength is obtained by using 5 algebraic equations with unknown vortex strength:(10a)Γ1=0.100824πbV∞α,(10b)Γ2=0.131804πbV∞α,(10c)Γ3=0.139014πbV∞α,(10d)Γ4=0.124964πbV∞α,(10e)Γ5=0.098874πbV∞α.According to the boundary conditions, the air flow is tangent to the object surface at each control point to determine the strength of each vortex. The lift of the wing can be calculated by satisfying the boundary conditions. For any wing that does not have an upper counter angle, the lift is produced by the free flow that crosses the vortex line as there is no side-washing speed or postwash speed.According to a finite number of elements, the lift on the plate is expressed as(11)L0=ρ∞V∞∑n=15ΓnΔyn.As the geometric relationship of each facet isΔyn=0.1b, the lift on the flat plate can be rewritten as(12)L0=ρ∞V∞2πb2α0.238184,where ρ∞ is the flow density, V∞ is the flow velocity, b is half length of the wingspan, and α is the attack angle.The attack angle is considered to be affected by a periodic disturbanceα=α0+α1cosΩ2t. Taking the periodic perturbation into the aerodynamic force and according to equation (12), the expression of the aerodynamic force containing the perturbation term is obtained:(13)L=L0+L1cosΩ2t,where L0=ρ∞V∞2πb2α00.238184 and L0=ρ∞V∞2πb2α10.238184.


## 2.1. Vortex Lattice (VL) Method

In this section, the subsonic air flow on the plate will be calculated by using the VL method. The VL method is a numerical implementation on the general 3D lifting surface problem. This method discretizes the vortex-sheet strength distributing on each lifting surface by lumping it into a collection of horseshoe vortices. Figure1 shows that a 3D lifting surface is discretized by a vortex lattice of horseshoe vortices, in which all contribute to the velocity V at any field point r.Figure 1
A 3D lifting surface discretized by a vortex lattice of horseshoe vortices.For the finite wingspan, the high-pressure air below the wing will turn over the low pressure air at the tip of the wing, which will cause the pressure on the upper surface of the wing tip to be equal to that of the lower surface. Unlike the two-dimensional flow around the airfoil, the main characteristic of the three-dimensional flow around the wing is the variation of the lift along the wingspan. In order to calculate the lift on the wing surface by using the vortex lattice method, the Biot–Savart law is used to calculate the induced velocity on the control point. The vortex strength of the vortex system is obtained, and the pressure difference on the upper and lower surface of the wing surface is deduced at the same time.The velocity induced by a vortex line with a strength ofΓn and a length of dl is calculated by the Biot–Savart law as(1)dV=Γndl×r4πr3.As shown in Figure2, the induced velocity is(2)dV=Γnsinθdl4πr2.Figure 2
The velocity induced by finite length vortex segments.LetAB represent a vortex segment in which the vortex vector points from A to B. C is a space point, and its normal distance from the AB line is rp.We can make an integral from the pointA to the point B to find out the size of the induced velocity:(3)V=Γn4πrp∫θ1θ2sinθdθ=Γn4πrpcosθ1−cosθ2.Asr0, r1,and r2, respectively, represent the AB, AC, and BC in Figure 2, the following relationships are established:(4)rp=r1×r2r0,cosθ1=r0⋅r1r0r1,cosθ2=r0⋅r2r0r2.Then, the formula calculating the induced velocity generated by the horseshoe vortex in the vortex lattice method is derived:(5)V=Γn4πr1×r2r1×r2r0⋅r1r1−r2r2.From equation (5), the velocity induced by a vortex segment AB at any point in the coordinate x,y,z can be obtained:(6a)VAB=Γn4πFac1AB⋅Fac2AB,(6b)VA∞=Γn4πz−z1nj+y−y1nkz−z1n2+y−y1n2×1.0+x−x1nx−x1n2+y−y1n2+z−z1n2,(6c)VB∞=−Γn4πz−z2nj+y1n−ykz−z1n2+y2n−y2×1.0+x−x2nx−x2n2+y−y2n2+z−z2n2,where VA∞ is the velocity induced by vortex lines from the point A to ∞ along the x-axis and VB∞ is the velocity induced by vortex lines from the point B to ∞ along the x-axis.The total velocity induced by a horseshoe vortex at a pointx,y,z representing a surface element (i.e., the nth panel element) is the sum of the various components calculated by equations (6a)–(6c). Therefore, the velocity induced by 2N vortexes is obtained to get the total induced velocity on the first m control points. It is expressed as(7)Vm=∑n=12NCm,nΓn.The contribution of all vortices to the downwash of the control point for themth panel element:(8)wm=∑n=12Nwm,n.Then, the lift on the rectangular cantilever plate with a sweep rectangular wing surface is computed according to the VL method. The vortex system is arranged on the rectangular wing surface in Figure3.Figure 3
Layout of the vortex lattice on the flat wing surface.


## 2.2. Application of VL Method on the Plate

In Figure3, the variable b is twice the span of the wing, c is the taper ratio of the wing surface, and c is equal to 1. The sweep angle is 0∘. Five panel elements are specified on the wing surface, and the vortex line is arranged at the leading edge 1/4 chord, and the control point is located at 3/4 chord. The wing surface is divided into 5 panel elements and each panel extends from the leading edge to the trailing edge.The coordinates of the5×1 vortices on the wing surface are listed in Table 1.Table 1
The coordinates of the attachment vortex and the control point in the wing surface.
Panel element number
x
m
y
m
x
1n
y
1n
x
2n
y
2n
1
0.25b
0.05b
0.833b
0b
0.833b
0.1b
2
0.25b
0.15b
0.833b
0.1b
0.833b
0.2b
3
0.25b
0.25b
0.833b
0.2b
0.833b
0.15b
4
0.25b
0.35b
0.833b
0.15b
0.833b
0.4b
5
0.25b
0.45b
0.833b
0.4b
0.833b
0.5bThe downwash velocity of each surface element induced at the control point is superimposed through the nonpenetrating condition:(9)wm=−V∞α.The solution of vortex strength is obtained by using 5 algebraic equations with unknown vortex strength:(10a)Γ1=0.100824πbV∞α,(10b)Γ2=0.131804πbV∞α,(10c)Γ3=0.139014πbV∞α,(10d)Γ4=0.124964πbV∞α,(10e)Γ5=0.098874πbV∞α.According to the boundary conditions, the air flow is tangent to the object surface at each control point to determine the strength of each vortex. The lift of the wing can be calculated by satisfying the boundary conditions. For any wing that does not have an upper counter angle, the lift is produced by the free flow that crosses the vortex line as there is no side-washing speed or postwash speed.According to a finite number of elements, the lift on the plate is expressed as(11)L0=ρ∞V∞∑n=15ΓnΔyn.As the geometric relationship of each facet isΔyn=0.1b, the lift on the flat plate can be rewritten as(12)L0=ρ∞V∞2πb2α0.238184,where ρ∞ is the flow density, V∞ is the flow velocity, b is half length of the wingspan, and α is the attack angle.The attack angle is considered to be affected by a periodic disturbanceα=α0+α1cosΩ2t. Taking the periodic perturbation into the aerodynamic force and according to equation (12), the expression of the aerodynamic force containing the perturbation term is obtained:(13)L=L0+L1cosΩ2t,where L0=ρ∞V∞2πb2α00.238184 and L0=ρ∞V∞2πb2α10.238184.


## 3. Formulation

In this section, the dynamic equations of the laminated composite cantilever plate are derived. As shown in Figure4, the parameters x, y, and z are, respectively, the spanwise direction, the direction of the chord, and the vertical direction of the plate. The plate is clamped at the position x=0. The ply stacking sequence is 0∘/90∘S and the layer number is N. The in-plane excitation is F=F0+F1cosΩ1t, which is distributed along the chord direction of the plate. The vertical direction of the plate along the spanwise direction is subject to subsonic aerodynamic force L.Figure 4
Mechanical model of cantilever laminated composite plate.The nonlinear governing equations are established in the Cartesian coordinate system. Reddy’s third-order shear deformation plate theory is used:(14a)ux,y,t=u0xt,y,t+zϕxxt,y,t−z343h2ϕx+∂w0∂x,(14b)vx,y,t=v0xt,y,t+zϕyxt,y,t−z343h2ϕy+∂w0∂y,(14c)wx,y,t=w0xt,y,t.The von Karman nonlinear strain-displacement relation is introduced. The displacements and strain-displacement relation are given as follows:(15a)εxx=∂u0∂x+12∂w0∂x2,εyy=∂v0∂y+12∂w0∂y2,(15b)γxy=12∂u0∂y+∂v0∂x+∂w0∂x∂w0∂y,γyz=12∂v0∂z+∂w0∂y,γzx=12∂u0∂z+∂w0∂x.Substituting the expressions of the strain into the displacements, we can obtain(16a)εxxεyyγxy=εx0εy0γxy0+zεx1εy1γxy1+z3εx3εy3γxy3,(16b)γyzγzx=γyz0γzx0+z2γyz2γzx2,where(17)εx0εy0γxy0=∂u0∂x+12∂w0∂x2∂v0∂y+12∂w0∂y2∂u0∂y+∂v0∂x+∂w0∂x∂w0∂yεxx1εyy1γxy1=∂ϕx∂x∂ϕy∂y∂ϕx∂y+∂ϕy∂x,εxx3εyy3γxy3=−c1∂ϕx∂x+∂2w0∂x2∂ϕy∂y+∂2w0∂y2∂ϕx∂y+∂ϕy∂x+2∂2w0∂x∂y,(18)γyz0γzx0=ϕy+∂w0∂yϕx+∂w0∂x,γyz2γzx2=−c2ϕy+∂w0∂yϕx+∂w0∂x,where(19)c1=43h2,c2=3c1.The symmetric cross-ply laminated composite plate is adopted, and we can obtain(20)σxxσyyσyzσzxσxy=Q¯11Q¯12000Q¯21Q¯2200000Q¯4400000Q¯5500000Q¯66εxxεyyγyzγzxγxy,where(21a)Q¯11Q¯12Q¯22Q¯16Q¯26Q¯66=C42C2S2S24C2S2C2S2C4+S4C2S2−4C2S2S42C2S2C44C2S2C3SCS3−C3S−CS3−2CSC2−S2CS3C3S−CS3−C3S2CSC2−S2C2S2−2C2S2C2S2C2−S22Q11Q12Q22Q66,(21b)Q¯44Q¯45Q¯55=C2S2−CSCSS2C2Q44Q55,(22)C=cosθ,S=sinθ,where θ is the angle of laminated layer. The stiffness of each layer is(23)Q11=E11−ν12ν21,Q12=ν12E21−ν12ν21,Q22=E21−ν12ν21,Q44=Q55=Q66=G12,where Ei is the Young modulus, G12 is the shear modulus, and ν12 is Poisson’s ratio.According to Hamilton’s principle,(24)∫0TδK−δU+δWdt=0.The nonlinear governing equations of motion are given as follows:(25a)Nxx,x+Nxy,y=I0u¨0+I1−c1I3φ¨x−c1I3∂w¨0∂x,(25b)Nyy,y+Nxy,x=I0v¨0+I1−c1I3φ¨y−c1I3∂w¨0∂y,(25c)Nyy,y∂w0∂y+Nyy∂2w0∂y2+Nxy,x∂w0∂y+Nxy,y∂w0∂x+Nxx,x∂w0∂x+Nxx∂2w0∂x2+c1Pxx,xx+2Pxy,xy+Pyy,yy+2Nxy∂2w0∂x∂y+Qx,x−c2Rx,x+Qy,y−c2Ry,y+L−γw˙0=c1I3∂u¨0∂x+∂v¨∂y+c1I4−c1I6∂φ¨x∂x+∂φ¨y∂y+I0w¨0−c12I6∂2w¨0∂x2+∂2w¨0∂y2,(25d)Mxx,x+Mxy,y−c1Pxx,x−c1Pxy,y−Qx−c2Rx=I1−c1I3u¨0+I2−2c1I4+c12I6φ¨x−c1I4−c1I6∂w¨0∂x,(25e)Myy,y+Mxy,x−c1Pyy,y−c1Pxy,x−Qy−c2Ry=I1−c1I3v¨0+I2−2c1I4+c12I6φ¨y−c1I4−c1I6∂w¨0∂y,where γ is the damping coefficient.The internal force is expressed as follows:(26)NαβMαβPαβ=∫−h/2h/2σαβ1zz3dz,QαRα=∫−h/2h/2σαz1z2dz.Substituting the internal force of equation (26) into equations (25a)–(25e), we can obtain the expression of the governing equation of motion by the displacements, which are given in Appendix.The boundary conditions of the cantilever plate are obtained at the same time as(27)x=0:w=v=u=φy=φx=0,(28a)x=a:Nxx=Nxy=Mxx=Mxy−c1Pxy=Q¯x=0,(28b)y=0:Nxy=Nyy=Myy=M−xyc1Pxy=Q¯y=0,(28c)y=b:Nxy=Nyy=Myy=Mxy−c1Pxy=Q¯y=0,(28d)Q¯x=Qx+∂Mxy∂y−c2Rx+c1∂Pxx∂x+∂Pxy∂y,(28e)Q¯y=Qy+∂Mxy∂x−c2Ry+c1∂Pyy∂y+∂Pxy∂x.The dimensionless equations can be obtained by introducing the following parameters:(29)w¯=w0h,t¯=tπ2Eabρ1/2,A¯ij=Aijab1/2Eh2,B¯ij=ab1/2Eh3Bij,D¯ij=ab1/2Eh4Dij,E¯ij=ab1/2Eh5Eij,F¯ij=ab1/2Eh6Fij,H¯ij=ab1/2Eh8Hij,I¯i=1abi+1/2ρIi,γ¯=ab2π2h4ρE1/2γ,x¯=xa,y¯=by,Ω¯i=1π2Eabρ1/2Ωi,F¯=b2Eh3F.


## 4. Frequency Analysis

High-dimensional nonlinear dynamic systems contain several types of the internal resonant cases which can lead to different forms of the nonlinear vibrations. When the system exists a special internal resonant relationship between two linear modes, the large amplitude nonlinear responses may suddenly happen which is dangerous in engineering and should be avoided.In order to study the internal resonance relationship between two bending modes, the finite element model of the cantilever laminated composite plate is established. The ply stacking sequence is0∘/90∘S and the material parameters [51] are shown in Table 2.Table 2
Single-layer material setting.
Material
E
11
E
22
S
12
υ
T300/QY8911
135 GPa
88 GPa
45 GPa
0.3
The single layer thickness ishii=1,2,…,6 and they are specified as h1=1mm,h2=2mm,…,h6=6mm.The natural frequencies of bending vibration of laminated composite cantilever plates with different span-chord ratio and different layer thickness are calculated, and the results are shown in Figures5–8. When the span-chord ratio is 2 : 1 and single thickness is 6 mm, the transverse vibration modes of the cantilever laminated composite plate with corresponding frequencies are shown in Figure 9.Figure 5
Natural frequencies with different thicknesses under the span-chord ratio 1 : 1.Figure 6
Natural frequencies with different thicknesses under the span-chord ratio 2 : 1.Figure 7
Natural frequencies with different thicknesses under the span-chord ratio 3 : 1.Figure 8
Natural frequencies with different thicknesses under the span-chord ratio 4 : 1.Figure 9
The transverse modes of the cantilever laminated composite plate when the span-chord ratio is 2 : 1 and single thickness is 6 mm.In frequency analysis, the first mode natural frequency is most important. Table3 is to illustrate the 1st mode natural frequency for 6 (thickness values) × 4 (span-chord ratio values) from Figures 5–8.Table 3
The 1st mode natural frequency for four span-chord ratio values corresponding to the thickness.
Span-chord ratio
Frequency (Hz)
h
1
h
2
h
3
h
4
h
5
h
6
1 : 1
15.15
30.31
45.46
60.61
75.77
90.92
2 : 1
3.83
7.65
11.48
15.31
19.13
22.96
3 : 1
1.69
3.39
5.08
6.78
8.47
10.17
4 : 1
0.95
1.90
2.85
3.80
4.75
5.70Based on the results of numerical simulations, the first six orders natural frequencies of the laminated composite cantilever plate are obtained, as shown in Figures5–8. It is obviously observed that there is a proportional relation between the two bending modes, such as relation 1 : 1 in area c of Figure 7, relation 1 : 2 in Figure 8, and relation 1 : 3 in Figures 5–7. We select 1 : 2 internal resonance relationship between two bending modes and the nonlinear vibrations of the laminated composite cantilever plate are investigated in the following analysis.


## 5. Perturbation Analysis

The discrete equation is derived by the Galerkin method, and the discrete function adopts the following expression:(30)w0=wt1X1xY1y+wt2X2xY2y,where(31a)Xix=sinλiax−sinhλiax+αicoshλiax−cosλiax,(31b)Yjy=sinβmby+sinhβmby−αmcoshβmby+cosβmby,where(32a)cosλiacoshλia−1=0,cosβmbcoshβmb−1=0,(32b)αi=coshλi−cosλisinhλi+sinλi,αm=−coshβm−cosβmsinhβm−sinβm.Similarly, the aerodynamic force is discretized by using the modal function:(33)L=l1X1xY1y+l2X2xY2y,where l1 and l2 represent the amplitudes of the aerodynamic forces corresponding to the two transvers vibration modes and they contain the perturbed items.Substituting equations (30)–(33) into equation (A.3), the governing differential equations of transverse motion of the system are derived as follows:(34a)w¨1+γ11w˙1+ω12w1+γ12f1cosΩ1tw1+γ13w2−γ14w13−γ15w12w2−γ16w1w22−γ17w23=γ18l1,(34b)w¨2+γ21w˙2+ω22w2+γ22f1cosΩ1tw2+γ23w1−γ24w23−γ25w22w1−γ26w2w12−γ27w13=γ28l2.Next, the averaged equations of the system are obtained by using the multiple scale method, and the internal resonant relationship 1 : 2 is considered:(35)ω12=14Ω2+εσ1,ω22=Ω2+εσ2,Ω1=Ω2=1,where ω1 and ω2 are two frequencies of the laminated composite plate and σ1 and σ2 are the tuning parameters.The scale transformations are given as follows:(36)γ11⟶εγ11,γ12⟶εγ12,γ13⟶εγ13,γ14⟶εγ14,γ15⟶εγ15,γ16⟶εγ16,γ17⟶εγ17,γ18⟶εγ18,γ21⟶εγ21,γ22⟶εγ22,γ23⟶εγ23,γ24⟶εγ24,γ25⟶εγ25,γ26⟶εγ26,γ27⟶εγ27,γ28⟶εγ28.The scale transformation can be obtained by introducing equation (36) into equations (34a) and (34b):(37a)w¨1+εγ11w˙1+ω12w1+εγ12f1cosΩ1tw1+εγ13w2−εγ14w13−εγ15w12w2−εγ16w1w22−εγ17w23=εγ18l1,(37b)w¨2+εγ21w˙2+ω22w+2εγ22f1cosΩ1tw2+εγ23w1−εγ24w23−εγ25w22w1−εγ26w2w12−εγ27w13=εγ28l2.The method of multiple scales [52] is used to obtain the averaged equations in following form:(38)wx,t,ε=w0x,T0,T1+εw1x,T0,T1,where T0=t and T1=εt.The derivatives with respect tot become(39a)ddt=∂∂T0dT0dt+∂∂T1dT1dt+⋯=D0+εD1+⋯,(39b)d2dt2=D0+εD1+⋯2=D02+2εD0D1+⋯,where Dn=∂/∂Tn,n=0,1,2,….Substituting equations (38)–(39b) into equations (37a) and (37b) and eliminating the secular terms, we have the averaged equations as follows:(40a)D1A1=−12γ11A1+iσ1A+112iγ12f1A¯1−3iA12A¯1γ13−2iγ16A1A2A¯2,(40b)D1A2=−12γ21A2+12iσ2A2−32iγ23A22A¯2−iγ26A1A2A¯1−14iγ28L2.In order to obtain the averaged equations in the polar coordinate system, we expressA1 and A2 in the following form:(41)A1T1=12a1T1ejϕT1,A2T2=12a2T2ejϕT2.Substituting equation (41) into equations (40a) and (40b) and separating the real and imaginary parts, the four-dimensional averaged equations in the polar coordinate system are obtained:(42a)a˙1=−14γ11a1sinϕ1−2σ1a1cosϕ1−γ12a1f1cosϕ1+32γ14a13cosϕ1+γ16a1a22sinϕ1,(42b)a1ϕ˙1=−14γ11a1cosϕ1+2σ1a1sinϕ1+γ12a1f1sinϕ1−32γ14a13sinϕ1+γ16a1a22cosϕ1,(42c)a˙2=−14γ21a2sinϕ2−σ2a2cosϕ2+34γ24a23cosϕ2+12γ26a2a12sinϕ2+γ28l2,(42d)a2ϕ˙2=−14γ21a2cosϕ2+σ2a2sinϕ2+34γ24a23sinϕ2+12γ26a2a12cosϕ2.In order to obtain the averaged equations in the Cartesian coordinate system, we rewriteA1 and A2 in the following forms:(43)A1T1=x1T1+ix2T1,A2T1=x3T1+ix4T1.According to the same way as the above, the averaged equations in the Cartesian coordinate system are obtained as follows:(44a)x˙1=−12γ11x1−σ1−12γ12f1x2+3γ14x2x12+x22+2γ16x2x32+x42,(44b)x˙2=−12γ11x2+σ1+12γ12f1x1−3γ14x1x12+x22−2γ16x1x32+x42,(44c)x˙3=−12γ21x3−12σ2x4+32γ24x4x32+x42+γ26x4x12+x22,(44d)x˙4=−12γ21x4+12σ2x3−32γ24x3x32+x42−γ26x3x12+x22−14γ28l2.


## 6. Numerical Simulation

In order to study the nonlinear vibration characteristics of the laminated composite cantilever plate with different modal modes, based on equations (44a) and (44d), the frequency-response curves are used to reveal the characteristics of this system. Let a˙1=a˙2=ϕ˙1=ϕ˙2=0, ϕ1=π/4, and ϕ2=3π/4.The frequency-response functions of the system are given as follows:(45a)0=γ11a1−2σ1a1−γ12a1f1+32γ14a13+γ16a1a22,(45b)0=γ21a2+σ2a2−34γ24a23+12γ26a2a12+2γ28l2.From equations (45a)–(45b), we can find that the amplitude a1 and amplitude a2 are coupled. The case of weak coupled form is considered here. Through introducing the proportional relation a1/a2=ε, the frequency-responses between the amplitudes and the tuning parameters are analyzed.According to the geometries and the material properties of the nonlinear system, the basic parameters are chosen asγ11=0.2, γ12=−6, γ14=5, γ16=−5, γ21=0.6, γ24=9, γ26=9, and γ28=−2. The relationship between the amplitude and the tuning parameter in different excitation conditions can be obtained.Figure10 gives the relationship between the amplitude a1 and tuning parameter σ1 in different internal force amplitudes f1, and Figure 11 gives the relationship between the amplitude a2 and the tuning parameter σ2 in different aerodynamic amplitudes l2.Figure 10
The relationship between amplitudea1 and tuning parameter σ1 in different internal force amplitudes f1.Figure 11
The relationship between amplitudea2 and tuning parameter σ2 in different aerodynamic amplitudes l2.The stiffness hardening phenomenon of the system can be seen in the relationship between the tuning parameters and the amplitude. With the increase of external excitation amplitude, the stiffness hardening phenomenon is gradually strengthened. The typical jump phenomenon of the nonlinear oscillations also happened in the system. The jump phenomenon appeared in the frequency-response curves at pointA∗ and point B∗ with the increase of the tuning parameters in Figure 11. The frequency-response curves have the wider resonance interval and the larger oscillation amplitudes under the stronger external excitation amplitude l2.Next, the geometries and the material properties of the nonlinear system are fixed, and the effects of different decoupling parametersε are calculated and the relationship between the amplitude and the tuning parameter can be found in Figures 12 and 13.Figure 12
The relationship between amplitudea1 and tuning parameter σ1 in different decoupling parameters ε.Figure 13
The relationship between amplitudea2 and tuning parameter σ2 in different decoupling parameters ε.In the following investigation, the bifurcation and chaotic motions of the laminated composite cantilever plate based on the averaged equation in the case of one to two internal resonances by using the fourth-order Runge–Kutta algorithm are analyzed. We choose the force excitationl2 as the controlling parameter to study the complicated nonlinear dynamics of the system.The tuning parameters and the initial conditions are chosen asσ1=0.1, σ2=1.1, γ11=0.2, γ12=5.2, γ14=4.2, γ16=−2.5, γ21=0.529, γ24=7, γ26=5, γ28=0.44, f1=0.9, x10=0.44, x20=0.798, x30=1.35, and x40=−0.5.Figure14 presents the two-dimensional bifurcation diagram to show the nonlinear oscillations of the laminated composite plate by varying the force excitation l2. The periodic and chaotic motions of the laminated composite cantilever plate system for x3 alternately occur with the increase of the parameter l2 in the interval of 20,50.Figure 14
The bifurcation diagram of the laminated composite cantilever plate forx3 via the forcing excitation l2.In Figures15–18, (a) and (c) show the waveforms on the planes t,x1 and t,x3, (b) and (d) depict the phase portraits on the planes x1,x2 and x3,x4, and (e) represents the three-dimensional phase portrait in the space x1,x2,x3.Figure 15
The period motion of the laminated composite cantilever plate is obtained whenl2=21.5.
(a)
(b)
(c)
(d)
(e)Figure 16
The multiperiod motion of the laminated composite cantilever plate is obtained whenl2=22.7.
(a)
(b)
(c)
(d)
(e)Figure 17
The quasi-period motion of the laminated composite cantilever plate is obtained whenl2=25.3.
(a)
(b)
(c)
(d)
(e)Figure 18
The chaotic motion of the laminated composite cantilever plate is obtained whenl2=26.5.
(a)
(b)
(c)
(d)
(e)From Figures15–18, the periodic and chaotic motions of the laminated composite cantilever plate occur with the increase of the amplitude of the parameter l2 can be seen. It is noticed that both the power spectrum and the Poincaré map can distinguish the periodic motions and the chaotic motions. When the l2 is in the range of 21,27, the periodic motions and chaotic motions alternately occur, the system experiences periodic motion, multiperiodic motion, quasi-periodic motion, and chaotic motion. It is obvious that the nonlinear vibration characteristics of the system are complicated under the influence of external excitation.


## 7. Conclusions

The nonlinear dynamics of the laminated composite cantilever plate under subsonic excitation force have been investigated in this paper. The aerodynamic force of a three-dimensional flat wing was calculated. Unlike a two-dimensional airfoil or an infinite length wing, the aerodynamic force of a three-dimensional flat wing is calculated. Some results are obtained:(1)
According to the ideal incompressible fluid flow condition and the Kutta–Joukowski lift theorem, the subsonic aerodynamic lift on the three-dimensional finite length flat wing is calculated by using the vortex lattice (VL) method.(2)
The finite length flat wing is modeled as a laminated composite cantilever plate based on Reddy’s third-order shear deformation plate theory and the von Karman geometry nonlinearity is introduced. The nonlinear partial differential governing equations of motion for the laminated composite cantilever plate subjected to the subsonic aerodynamic force are established via Hamilton’s principle.(3)
Through comparing the natural frequencies of the linear system with different materials and geometry parameters, the relationship of 1 : 2 internal resonance is found and considered for the nonlinear vibration analysis.(4)
Corresponding to several selected parameters, the frequency-response curves are obtained. The hardening-spring-type behaviors and jump phenomena are exhibited with the variation of the tuning parameters.(5)
The influence of the force excitation on the bifurcations and chaotic behaviors of the laminated composite cantilever plate have been investigated numerically. The periodic motion, multiperiodic motion, quasi-periodic motion, and chaotic motion of the system occurred with the increasing of the external excitation.



---

*Source: 4601672-2020-04-22.xml*