# Free Vibration Analysis of Laminated Composite Double-Plate Structure System with Elastic Constraints Based on Improved Fourier Series Method

Authors: Ying Zhang; Dongyan Shi; Dongze He; Dong Shao
Journal: Shock and Vibration
 (2021)

---

## Abstract

An analytical model of laminated composite double-plate system (LCDPS) is established, which efficiently analyzes the common 3D plate structure in engineering applications. The proposed model combines the first-order shear deformation theory (FSDT) and the classical delamination theory, and then the LCDPS’s vibration characteristics are investigated. In the process of analysis, the improved Fourier series method (IFSM) is used to describe the displacement admissible function of the LCDPS, which can remove the potential discontinuities at the boundaries. Five sets of artificial springs are introduced to simulate the elastic boundary constraints, and the restraints of the Winkler elastic layer can be adjustable. The improved Fourier series is substituted into the governing equations and boundary conditions; then, applying the Rayleigh–Ritz method, we take all the series expansion coefficients as the generalized coordinates. After that, a set of standard linear algebraic equations was obtained. On this basis, the natural frequency and mode shapes of the LCDPS can be obtained by solving the standard eigenvalue problem. By the discussion of numerical examples and the comparison with those of the reports in the literature, the convergence and the reliability of the present approach are validated. Finally, the parametric investigations of the free vibration with complex boundary conditions are carried out, including the influence of boundary conditions, lamination scheme, plate geometric parameters, and elastic coefficient between two plates.


---

## Body

## 1. Introduction

Plate structure is one of the most important structures in the field of engineering applications. In particular, laminated composite plates have been widely applied in engineering owing to their excellent properties such as higher strength and hardness and lighter weight. The analysis of vibration characteristics for plates has been widely studied for decades. Large numbers of designers, engineers, and scholars have been attracted to conduct in-depth research on the vibration features of plates [1, 2]. However, most researches are based on the classical thin plate theory, studying the vibration characteristics of single-layer plate structures with classical boundary conditions [3–7]. Li et al. [8, 9] studied the modal characteristics of a plate with general edges supports by employing the Rayleigh-Ritz method. Shuyu [10] analyzed the flexural vibration characteristics of a thin rectangular plate under free boundary constraints. Liew et al. [11] researched the frequency of rectangular plates based on the energy principle also. On the basis of the boundary-continuous displacement Fourier approach, Hatami et al. [12] calculated the free vibration features of moving laminated composite plates. Carrera et al. [13] studied the 2D and 3D closed-form solutions of simply supported plates and then they extended this method to the 3D sandwich beams [14]. Kabir et al. [15] carried out boundary-continuous displacement Fourier methods to analyze the frequency characteristics of thin anisotropic plates. Su et al. [16] explored the vibration problem of laminated composite rectangular plates. Rongong et al. [17] used plasma technology to produce constrained layer damping (CLD) coatings on metal substrates. The process involves depositing relatively thick hard ceramic layers on soft polymer damping materials while maintaining the integrity of the two layers. Based on the Reuss model and the Hashin-Shtrikman equation, Yu et al. [18] discussed the damping efficiency of the coating system. Catania and Strozzi [19] proposed a simple model based on the assumption of multilayer sawtooth beams to study the performance of different solutions in free and forced bending vibration responses.A major advancement of the single-layered system would be extended to complex plate systems such as double-beam system [20], double-plate system [21–24], and double-shell system [25]. Kukla [22] used Green’s functions to analyze the free vibration of two elastic connected rectangular plates. Oniszczuk [23] analytically investigated two simply supported rectangular plates connected by a Winkler layer based on the Kirchhoff-Love plate theory. Hedrih [24] used Bernoulli’s particular integral and Lagrange’s method to simulate the transversal vibrations of double-plate systems. The laminated composite double-plate system (LCDPS) is composed of two parallel plates and an elastic layer between them as the Winkler elastic layer [26]. Compared with the single-layered plate structure, double-layered plate structures have better shock absorption and higher impact-damage tolerance and noise-reducing performance, which can be widely used in the maritime field, aeronautics and astronautics, auto part, and petrochemical and other industries [27], such as the application in the bottom plate and making the bulkhead of the hull structure, the skin of the aircraft structure, the inner decorative plate, and the double-layer sound insulation plates [28]. Many shock absorbers also use double-layer plate structures [29]. Although there have been some studies on the vibration characteristics of double-plate systems, the boundary conditions of the plates are limited to the classical boundary conditions (fixed or simply supported). These two-plate systems are made of ordinary isotropic materials. Therefore, it is necessary to extend the research to the composite laminate double-plate system to provide a relevant basis for engineering design.Lim et al. [30] analyzed doubly connected superelliptical laminated composite plates with arbitrary geometric boundary conditions. Oniszczuk [31] studied free transverse vibrations of a simply supported rectangular double-plate system interconnected by a homogeneous Winkler elastic layer first and then analyzed forced transverse vibrations of the same structure based on the Kirchhoff–Love plate theory [32]. Mahi et al. [33] researched bending and free vibration characters of isotropic, functionally various plates by using a new hyperbolic shear deformation theory. Sometimes double-layer plates’ structures are known as sandwich plate structures. Fu et al. [34] presented sound transmission through stiffened double laminated composite sandwich plates. Liu and Jeffers [35] adopted a layer-wise displacement theory to build functionally graded sandwich plates and laminated composite models. Yazdani et al. [36] analyzed the vibration of double-bonded sandwich cylindrical shells under multiphysical loadings. Park and Baz [37] presented the equations of motion for the composite plates with a single active constrained damping layer based on the classical laminated plate theory and investigated transient vibration of the composite plates. Similarly, Cao et al. [38] explored vibrations and the damping behavior of thin constrained composite plates with double piezoelectric layers based on Fourier transformation and classical laminated plate theory. Nastos and Saravanos [39] developed simulation of transient dynamic response with thick composite and sandwich panels. Naik and Sayyad [40] presented a 1D thermal analysis (cylindrical bending) of laminated composite and sandwich plates by using a new fifth-order shear and normal deformation theory. In addition, the research on structural vibration characteristics has been extended to many new composite materials in recent years. For example, Gurses et al. [41] analyzed the free vibration of nanosized annular sector plate. Akgöz and Civalek [42] investigated static bending response of single-walled carbon nanotubes (SWCNTs) embedded in an elastic medium. Jalaei and Civalek [43] presented the dynamic instability of viscoelastic porous functionally graded (FG) nanobeam, which is embedded on visco-Pasternak foundation subjected to an axially oscillating loading as well as magnetic field. Arefi et al. [44] studied size-dependent deflection analysis of composite microplates, which were reinforced by FG graphene nanoplatelets (GNPs) under transverse impact, and Karami et al. [45] calculated the size-dependent free vibration characteristics of the same nanocomposite plates. Therefore, the composite structure has a wide range of application fields and important research significance in engineering practice. It is necessary to carry out theoretical research and practical applications on such structures.According to the above literature review, it is not difficult to find that the current research mainly focuses on the plate structure with classical boundaries. Although there were some methods for vibration analysis of the sandwich plate, there is a lack in research on the laminated composite double-plate structure. Therefore, it is still necessary to develop a practical method to study the vibration behavior of laminated composite double-plate system with general boundary conditions. In the present study, the linear and torsional springs will be introduced on the boundary of the LCDPS, and arbitrarily general boundary conditions can be achieved via setting different spring stiffness parameters. The displacement function of the LCDPS is expressed by the improved Fourier series expansion, and the discontinuity of displacement function at the boundary of the plate is avoided by adding auxiliary functions. Based on the Rayleigh-Ritz method, the natural frequency and mode shape of the LCDPS can be obtained by solving the standard eigenvalue problem. Then the correctness and convergence of the present method are verified by comparison with the experimental results, relevant literature, and FEM results.


## 2. Analysis Model of Vibration Characteristics of the Laminated Composite Double-Plate Structure

### 2.1. Description of the Model

Figure1 shows the structure of a laminated composite double-plate system (LCDPS) with general elastic boundary conditions. There are two composite laminate parallel plates in the structure with the same length (a) and width (b) connected by a uniformly distributed Winkler-type elastic layer. In LCDPS, the boundary weight parameter is introduced to simulate the constraint of boundary conditions, which are equivalent to the effect of load on the boundary. By using the FSDT, there are three kinds of boundary forces that need to be simulated: tensile force, shear force, and bending moment. To meet the required boundary conditions, five groups of spring elastic weight parameters are introduced, which can be divided into two categories: linear springs, ku, kv, and kw, and torsional springs, Kφx and Kφy. The elastic force of a linear spring can be equivalent to shear force and tension force, and the bending moment is described by the coil spring. In order to describe the double-plate structure clearly, a coordinate system is established. Drawing x-axis, y-axis, and z-axis along with the length, width, and thickness directions, suppose that the plate is composed of multiple composite layers, which are closely combined together.Figure 1
Model of LCDPS (laminated composite double-plate system) with general boundary conditions.


### 2.2. Kinematic Relations and Stress Resultants

The mechanical properties of each plate in the double-plate structure are the same. According to the FSDT, the displacement and rotation components of any point of LCDPS can be represented as follows [46–49]:(1)Ux,y,z,t=ux,y,t+zφxx,y,t,Vx,y,z,t=vx,y,t+zφyx,y,t,Wx,y,z,t=wx,y,t,where u, v, and w express the middle surface displacements and t is the time variable. φx and φy are normal rotations, with respect to y-axis and x-axis in the transverse direction. Based on the linear elastic theory, the relationship between strain and displacement of the k-th layer of the laminated plate can be represented as [50](2)εxx=εxx0+zχxx,εyy=εyy0+zχyy,γxy=γxy0+zχxy,γxz=γxz0,γyz=γyz0,where εxx0,εyy0,γxy0 represent the normal and shear strain and χxx,χyy,χxy represent curvature and twist changes component on the middle surface of a laminated plate. The transverse shear strains γxz0 and γyz0 are constant through the thickness. These quantities are expressed according to the generalized displacements as(3)εxx0=∂u∂x,εyy0=∂v∂y,γxy0=∂v∂x+∂u∂y,χxx=∂φx∂x,χyy=∂φy∂x,χxy=∂φx∂y+∂φy∂x,γxz0=φx+∂w∂x,γyz0=φy+∂w∂x,and we have that zzk<z<zk+1 expresses the variation range of plane thickness; zk and zk+1 are the coordinate values on the bottom and top of the k-th layer in the thickness direction. According to the generalized Hooke’s Law, normal stresses σxx and σyy, shear stresses τyz, τxz, and τxy in the x-, y-, and z-coordinate system can be obtained as [47, 49, 51](4)σxxσyyτxy=Q¯11kQ¯12kQ¯16kQ¯12kQ¯22kQ¯26kQ¯16kQ¯26kQ¯66kεxxεyyγxy,τxzτyz=Q¯55kQ¯45kQ¯45kQ¯44kγxzγyz,where Q¯ijki,j=1,2,4,5,6 represents the stiffness coefficient of laminated plates in the O-xyz coordinate system, which can transfer the stress and strain of all layers to the coordinate system O-xyz. Q¯ijk can be determined in the following way [46, 49, 50]:(5a)Q¯11k=Q11km4+2Q12k+2Q66km2n2+Q22kn4,(5b)Q¯12k=Q11k+Q22k−4Q66km2n2+Q12km4+n4,(5c)Q¯22k=Q11kn4+2Q12k+2Q66km2n2+Q22km4,(5d)Q¯16k=Q11k−Q12k−2Q66km3n+Q12k−Q22k+2Q66kmn3,(5e)Q¯26k=Q11k−Q12k−2Q66kmn3+Q12k−Q22k+2Q66km3n,(5f)Q¯66k=Q11k+Q22k−2Q12k−2Q66km2n2+Q66km4+n4,(5g)Q¯44k=Q44km2+Q55kn2,Q¯45k=Q55k−Q44kmn,(5h)Q¯55k=Q55km2+Q44kn2,m=cosαfiberk,n=sinαfiberk,where the fiber laying angle of k-th layer fiber is αfiberk, which represents the included angle between the laying direction of the k-th layer material and the x-axis. The elastic parameters of materials Qijki,j=1,2,4,5,6 can be defined as [52](6)Q11k=E1k1−μ12kμ21k,Q22k=E2k1−μ12kμ21k,Q12k=μ21kQ11k,Q44k=G23k,Q55k=G13k,Q66k=G12k,where E1k and E2k represent the moduli of the k-th layer, respectively, and μ12k and μ21k represent Poisson’s ratios. G13k, G23k, and G12k express the shear moduli. By letting E1k=E2k and G12k=G13k=G23k=E1k/21+1μ12, the isotropy single-layer plate will be simulated. Based on the classical laminated theory and the stress-strain relation of the plane, we can obtain the constitutive model of the laminated composite plate [46, 47].(7)NM=ABBDεχ,QxQy=κA55κA45κA45κA44γxz0γyz0,where N=NxxNyyNxyT, Q=Qx,QyT, and M=MxxMyyMxyT, respectively, represent the internal force, shear force, and moment components of the structure. ε=εxx0εyy0γxy0T is the strain vector and χ=χxxχyyχxyT represents the curvature vector of the middle plane. κ is the shear correction factor. Aiji,j=1,2,6 is the extensional stiffness; Aiji,j=4,5 is the torsional stiffness. Bij and Diji,j=1,2,6 denote the coupling stiffness and bending stiffness, individually. Aij,Bij,Dij can be gained by the following formula:(8)Aij,Bij,Dij=∑zkNk∫zkzk+1Q¯ijk1,z,z2dz.Through the derivation of the above formulas, the kinematic relationships of the plate model and the stress results are given.


### 2.3. Energy Expressions

From the structural feature, we can infer that it is difficult to theoretically solve the boundary conditions and structural differential equations of the laminated composite double-plate structure system (LCDPS) and other more complex structures. Therefore, most researchers study the vibration characteristics of the complex plate structures by employing the numerical analysis method, such as the Rayleigh-Ritz method, variational method [53], and the Galerkin method [54]. In this section, we will derive expressions for various energies stored in the LCDPS and construct the energy function of the structure (Lagrange energy function: L) to prepare the Rayleigh-Ritz energy method. On the basis of the previous text, the strain energy (E) for the double-plate structure can be formulated as(9)Ep=12∫Vσxxpεxxp+σyypεyyp+τyzpγyzp+τxzpγxzp+τxypγxypdV,which can be further deduced as(10)Ep=12∫sNxpεxx0,p+Nypεyy0,p+Nxypεxy0,p+Mxpχxx0,p+Mypχyy0,p+Mxypχxyp+Qxpγxz0,p+Qypγyz0,pdS,where the symbolic variable p represents the number of plates in the laminated composite double-plate structure system. If no special instructions are given later, p=U represents the upper plate on the system (Plate 1 in Figure 1), and p=L denotes the lower plate on the system (Plate 2 in Figure 1). Then the strain energy expression can be written in three components: Ep=ESp+EBp+EBSp, in which ESp, EBp, and EBSp represent the stretching energy, the bending energy, and the bending-stretching coupling energy, separately. By substituting formulas (3) and (7) into formula (10), they can be expressed as follows:(11)Esp=12∫0La∫0LbA11εxx0,p2+2A12εxx0,pεyy0,p+2A16εxx0,pγxy0,p+A22εyy0,p2+2A26εyy0,pγxy0,p+A66γxy0,p2+κA44γxz0,p2+κA55γyz0,p2+2κA45γxz0,pγyz0,pdydx,=12∫0La∫0LbA11∂up∂x2+2A12∂up∂x∂vp∂y+2A16∂up∂x∂up∂y+∂vp∂x+A22∂vp∂y2+2A26∂vp∂y∂up∂y+∂vp∂x+A66∂up∂y+∂vp∂x2+κA44∂wp∂x+φxp2+κA55∂wp∂y+φyp2+2κA45∂wp∂x+φxp∂wp∂y+φypdydx,(12)EBp=12∫0La∫0LbD11χxxp2+2D12χxxpχyyp+D22χxxpχyyp+2D16χxxpχxyp+2D26χyypχxyp+D66χxyp2dydx,=12∫0La∫0LbD11∂φxp∂x2+2D12∂φxp∂x∂φyp∂y+2D16∂φxp∂x∂φxp∂y+∂φyp∂x+D22∂φyp∂y2+2D26∂φyp∂y∂φyp∂y+∂φxp∂x+D66∂φyp∂y+∂φxp∂x2dydx,(13)EBSp=12∫0La∫0LbB11εxx0,pχxxp+B12εxx0,pχyyp+B22εyy0,pχyyp+B16εxx0,pχxyp+B26εyy0,pχyyp+B16χxxp+B26χyyp+B26χxypγxy0,pdydx,=12∫0La∫0LbB11∂up∂x∂φxp∂x+B12∂up∂x∂φyp∂y+B16∂up∂x∂φxp∂y+∂φyp∂x+B22∂vp∂y2+B26∂vp∂y∂φyp∂y+∂φxp∂x+B16∂up∂x+B26∂vp∂y+B66∂φyp∂y+∂φxp∂x∂vp∂x+∂up∂ydydx.The kinetic energy function of the double-plate structure produced by vibration can be expressed as follows:(14)Tp=b2∫0La∫0Lbρpz∂Up∂t2+∂Vp∂t2+∂Wp∂t2dydx,=b2∫0La∫0LbI0p∂up∂t2+2I1p∂up∂t∂ϕxp∂t+I2p∂ϕyp∂t2+I0p∂vp∂t+2I1p∂vp∂t∂ϕyp∂t+I2p∂ϕyp∂t2+I0p∂wp∂tdydx,where the inertia terms can be obtained by the following formula:(15)I0p,I1p,I2p=∑k=1Nkp∫zkpzk+1pρpz1,z,z2dz.In the present model, the upper plate and the lower plate are connected by the middle Winkler-type elastic layer, which consists of a series of face springs, and the stiffness of the middle Winkler layer is represented byKc. When the system vibrates, the following formula denotes energy contained in the middle elastic layer:(16)Ec=12∫0La∫0LbKcwU−wLdydx.The major purpose is to establish a general solution model for the vibration characteristics of LCDPS under general boundary constraints, which are achieved by the penalty parameter method. By introducing five sets of penalty factors, the boundary conditions are added to the total energy function in the form of additional potential energy. Therefore, any specific elastic boundary constraints can be replaced by the proper penalty parameter. The resulting additional potential energy is defined as(17)Ebp=12∫0bku0x,pupxp,yp2+kv0x,pvpxp,yp2+kw0x,pwpxp,yp2+Kφx0x,pφxpxp,yp2+Kψy0x,pφypxp,yp2xp=0ku1x,pupxp,yp2+kv1x,pvpxp,yp2+kw1x,pwpxp,yp2+Kφx1x,pφxpxp,yp2+Kφy1x,pφypxp,yp2xp=adyp+12∫0aku0y,pupxp,yp2+kv0y,pvpxp,yp2+kw0y,pwpxp,yp2+Kφx0y,pφxpxp,yp2+Kφy0y,pφypxp,yp2yp=0ku1y,pupxp,yp2+kv1y,pvpxp,yp2+kw1y,pwpxp,yp2+Kφx1y,pφxpxp,yp2+Kφy1y,pφypxp,yp2yp=bdxp.We have thatkuΩχ,kvΩχ,kwΩχ,KφxΩχandKφxΩχΩ=0,1,χ=x,y represent the five sets of penalty factors of different boundaries at x=0, x=a, y=0, and y=b, which are stiffness coefficients of springs. It can be seen from the above formula that the general boundary conditions can be obtained as long as the penalty factors are given suitable values. It shows the superiority of this method in dealing with boundary conditions, which can deal with boundary conditions as easily as geometric parameters.


### 2.4. Displacement Admissible Functions and Solution Process

Based on the Hamilton principle [46, 47, 55], the governing equations of LCDPS on O-xyz coordinate system can be obtained as(18)∂Nx∂x+∂Nxy∂y=I0∂2u∂t2+I1∂2ϕx∂t2,∂Ny∂y+∂Nxy∂x=I0∂2v∂t2+I1∂2ϕy∂t2,∂Qx∂x+∂Qy∂y=I0∂2w∂t2,∂Mx∂x+∂Mxy∂y−Qx+mx=I1∂2u∂t2+I2∂2ϕx∂t2,∂My∂y+∂Mxy∂x−Qy+my=I1∂2v∂t2+I2∂2ϕy∂t2.It is obvious that there are the second derivatives in different governing equations. In order to avoid a discontinuity at boundary conditions, the admissible function of displacement is expressed by trigonometric Fourier series [56]. In addition, the supplementary function is introduced to ensure the existence of the derivative and the continuity of functions. The modified Fourier series technique is concise and efficient. This method does not need to solve partial differential equations like other numerical methods [57, 58], thus avoiding numerical errors. The whole solving process is based on the energy integration of the structure, and the Fourier series used in this paper can obtain the accurate solution to the integral calculation because of its own integral property [59, 60]. The displacement admissible function can be written as(19)upx,y,t=e−jωt∑m=0M∑n=0NAmnpcosλmxcosλny+∑l=12∑n=0Nalnpplxcosλny∑l=12∑m=0Mblmpplycosλmx,vpx,y,t=e−jωt∑m=0M∑n=0NBmnpcosλmxcosλny+∑l=12∑n=0Nclnpplxcosλny∑l=12∑m=0Mdlmpplycosλmx,wpx,y,t=e−jωt∑m=0M∑n=0NCmnpcosλmxcosλny+∑l=12∑n=0Nelnpplxcosλny∑l=12∑m=0Mflmpplycosλmx,ϕxpx,y,t=e−jωt∑m=0M∑n=0NDmnpcosλmxcosλny+∑l=12∑n=0Nglnpplxcosλny∑l=12∑m=0Mhlmpplycosλmx,ϕypx,y,t=e−jωt∑m=0M∑n=0NEmnpcosλmxcosλny+∑l=12∑n=0Nilnpplxcosλny∑l=12∑m=0Mjlmpplycosλmx.where λm=mπ/La, λn=nπ/Lb, and M and N are the truncation number.Amnp, Bmnp, Cmnp, Dmnp, and Emnp express coefficients expansion of two-dimensional standard cosine function of Fourier series. Accordingly, alnp, blmp, clnp, dlmp, elnp, flmp, glnp, hlmp, ilnp, and jlmp represent the corresponding coefficients of the supplementary terms.(20)p1x=xxa−12,p2x=x2axa−1,p1y=yyb−12,p2y=y2byb−1.It is easy to verify that(21)p10=p1a=p1′a=0,p1′0=1,p20=p2a=pa′0=0,p2′a=1,p10=p1b=p1′a=0,p1′0=1,p20=p2b=p2′0=0,p2′b=1.As explained in [8, 9, 47, 54, 55], the design of modified Fourier series should satisfy the second-order partial differential governing equations of plate. It is easy to see from equation (21) that the function complementary p can make sure that the first derivative of the admission displacement functions is continuous in the whole solution domain and the second derivative of each point exists. So displacement function and the derivative term have a continuous derivative across the edge, and the potential discontinuities with the displacement function and its derivatives at the end points can be removed.Substituting equations (2) and (4) into equation (18), we can get matrix equation as(22)L11UL12U0L14UL15U00000L12UL22U0L14UL15U0000000L33UL34UL35U00LkL00L14UL24UL34UL44UL45U00000L15UL25UL35UL45UL55U0000000000L11LL12L0L14LL15L00000L12LL22L0L24LL25L00LkU0000L33LL34LL35L00000L14LL24LL34LL44LL45L00000L15LL25LL35LL45LL55L0−ω2M11U00M14U0000000M22U00M25U0000000M33U0000000M14U00M44U0000000M25U00M55U0000000000M11L00M14L0000000M22L00M25L0000000M33L0000000M14L00M44L0000000M25L00M55Lupvpwpφxpφyp=00000,where L and MLij=Lji,Mij=Mji indicate the coefficients of the linear operator. They can be represented as(23)L11p=A11p∂2∂x2+2A16p∂2∂x∂y+A66p∂2∂y2,L12p=A16p∂2∂x2+A12p+A66p∂2∂x∂y+A26p∂2∂y2,L13p=0,L23p=0,L14p=B11p∂2∂x2+2B16p∂2∂x∂y+B66p∂2∂y2,L15p=B16p∂2∂x2+B12p+B66p∂2∂x∂y+B26p∂2∂y2,L22p=A66p∂2∂x2+2A26p∂2∂x∂y+A22p∂2∂y2,L24p=B16p∂2∂x2+B12p+B66p∂2∂x∂y+B26p∂2∂y2,L25p=B66p∂2∂x2+2B26p∂2∂x∂y+B22p∂2∂y2,L33p=−A55p∂2∂x2−2A45p∂2∂x∂y−A44p∂2∂y2,L34p=−A55p∂∂x−2A45p∂∂y,L35p=−A45p∂∂x−2A44p∂∂y,L44p=D11p∂2∂x2+2D16p∂∂x∂y+D66p∂2∂y2−A55p,L45p=D16p∂2∂x2+D12p+D66p∂2∂x∂y+D26p∂2∂y2−A45p,L55p=D66p∂2∂x2+2D26p∂2∂x∂y+D22p∂2∂y2−A44p,M11p=M22p=M33p=−I0p∂2∂t,M14p=M25p=−I1p∂2∂t,M44p=M55p=−I2p∂2∂t.One can note that there are second-order derivatives of all the displacements and rotation components of the plate, as we have posted before. After deriving the expressions of dynamic energy and potential energy, the Lagrange function of the double-plate structure can be written as(24)L=Tp+Wp−Ep−Espp−Ec,(25)∂L∂Ξmnp=0,andΞp=Ap,Bp,Cp,Dp,Ep,m=0,1,…M;n=0,1,…N,∂L∂Ψlnp=0,andΨp=ap,cp,ep,gp,ip,l=1,2;n=0,1,…N,∂L∂ϒlmp=0,andϒp=bp,dp,fp,hp,jp,l=1,2;m=0,1,…M.By substituting equations (11)–(15) and (19) into equation (24) and then carrying out the Rayleigh-Ritz energy method, we can obtain a set of linear algebraic equations to determine the unknown Fourier coefficients for eigenfunctions. The number of linear equations is 5∗2∗, M+1, ∗, N+1 + 10∗2∗, M+N+2, which can be turned into the matrix as(26)K−ω2MΠ=0.We have thatΠ=AUBUCUDUEUALBLCLDLELT are the unknown Fourier coefficients; K stands for the stiffness matrix and M represents the mass matrix, separately.Please refer to Appendix A for the detailed expression of the above matrices. Obviously, by solving the eigenproblem of the standard matrix, the natural frequency and eigenvector can be easily obtained. Once the eigenvector of the given frequencyΠ is determined, the displacement of the double-plate structure system can be determined by substituting the coefficients into the equation.


## 2.1. Description of the Model

Figure1 shows the structure of a laminated composite double-plate system (LCDPS) with general elastic boundary conditions. There are two composite laminate parallel plates in the structure with the same length (a) and width (b) connected by a uniformly distributed Winkler-type elastic layer. In LCDPS, the boundary weight parameter is introduced to simulate the constraint of boundary conditions, which are equivalent to the effect of load on the boundary. By using the FSDT, there are three kinds of boundary forces that need to be simulated: tensile force, shear force, and bending moment. To meet the required boundary conditions, five groups of spring elastic weight parameters are introduced, which can be divided into two categories: linear springs, ku, kv, and kw, and torsional springs, Kφx and Kφy. The elastic force of a linear spring can be equivalent to shear force and tension force, and the bending moment is described by the coil spring. In order to describe the double-plate structure clearly, a coordinate system is established. Drawing x-axis, y-axis, and z-axis along with the length, width, and thickness directions, suppose that the plate is composed of multiple composite layers, which are closely combined together.Figure 1
Model of LCDPS (laminated composite double-plate system) with general boundary conditions.


## 2.2. Kinematic Relations and Stress Resultants

The mechanical properties of each plate in the double-plate structure are the same. According to the FSDT, the displacement and rotation components of any point of LCDPS can be represented as follows [46–49]:(1)Ux,y,z,t=ux,y,t+zφxx,y,t,Vx,y,z,t=vx,y,t+zφyx,y,t,Wx,y,z,t=wx,y,t,where u, v, and w express the middle surface displacements and t is the time variable. φx and φy are normal rotations, with respect to y-axis and x-axis in the transverse direction. Based on the linear elastic theory, the relationship between strain and displacement of the k-th layer of the laminated plate can be represented as [50](2)εxx=εxx0+zχxx,εyy=εyy0+zχyy,γxy=γxy0+zχxy,γxz=γxz0,γyz=γyz0,where εxx0,εyy0,γxy0 represent the normal and shear strain and χxx,χyy,χxy represent curvature and twist changes component on the middle surface of a laminated plate. The transverse shear strains γxz0 and γyz0 are constant through the thickness. These quantities are expressed according to the generalized displacements as(3)εxx0=∂u∂x,εyy0=∂v∂y,γxy0=∂v∂x+∂u∂y,χxx=∂φx∂x,χyy=∂φy∂x,χxy=∂φx∂y+∂φy∂x,γxz0=φx+∂w∂x,γyz0=φy+∂w∂x,and we have that zzk<z<zk+1 expresses the variation range of plane thickness; zk and zk+1 are the coordinate values on the bottom and top of the k-th layer in the thickness direction. According to the generalized Hooke’s Law, normal stresses σxx and σyy, shear stresses τyz, τxz, and τxy in the x-, y-, and z-coordinate system can be obtained as [47, 49, 51](4)σxxσyyτxy=Q¯11kQ¯12kQ¯16kQ¯12kQ¯22kQ¯26kQ¯16kQ¯26kQ¯66kεxxεyyγxy,τxzτyz=Q¯55kQ¯45kQ¯45kQ¯44kγxzγyz,where Q¯ijki,j=1,2,4,5,6 represents the stiffness coefficient of laminated plates in the O-xyz coordinate system, which can transfer the stress and strain of all layers to the coordinate system O-xyz. Q¯ijk can be determined in the following way [46, 49, 50]:(5a)Q¯11k=Q11km4+2Q12k+2Q66km2n2+Q22kn4,(5b)Q¯12k=Q11k+Q22k−4Q66km2n2+Q12km4+n4,(5c)Q¯22k=Q11kn4+2Q12k+2Q66km2n2+Q22km4,(5d)Q¯16k=Q11k−Q12k−2Q66km3n+Q12k−Q22k+2Q66kmn3,(5e)Q¯26k=Q11k−Q12k−2Q66kmn3+Q12k−Q22k+2Q66km3n,(5f)Q¯66k=Q11k+Q22k−2Q12k−2Q66km2n2+Q66km4+n4,(5g)Q¯44k=Q44km2+Q55kn2,Q¯45k=Q55k−Q44kmn,(5h)Q¯55k=Q55km2+Q44kn2,m=cosαfiberk,n=sinαfiberk,where the fiber laying angle of k-th layer fiber is αfiberk, which represents the included angle between the laying direction of the k-th layer material and the x-axis. The elastic parameters of materials Qijki,j=1,2,4,5,6 can be defined as [52](6)Q11k=E1k1−μ12kμ21k,Q22k=E2k1−μ12kμ21k,Q12k=μ21kQ11k,Q44k=G23k,Q55k=G13k,Q66k=G12k,where E1k and E2k represent the moduli of the k-th layer, respectively, and μ12k and μ21k represent Poisson’s ratios. G13k, G23k, and G12k express the shear moduli. By letting E1k=E2k and G12k=G13k=G23k=E1k/21+1μ12, the isotropy single-layer plate will be simulated. Based on the classical laminated theory and the stress-strain relation of the plane, we can obtain the constitutive model of the laminated composite plate [46, 47].(7)NM=ABBDεχ,QxQy=κA55κA45κA45κA44γxz0γyz0,where N=NxxNyyNxyT, Q=Qx,QyT, and M=MxxMyyMxyT, respectively, represent the internal force, shear force, and moment components of the structure. ε=εxx0εyy0γxy0T is the strain vector and χ=χxxχyyχxyT represents the curvature vector of the middle plane. κ is the shear correction factor. Aiji,j=1,2,6 is the extensional stiffness; Aiji,j=4,5 is the torsional stiffness. Bij and Diji,j=1,2,6 denote the coupling stiffness and bending stiffness, individually. Aij,Bij,Dij can be gained by the following formula:(8)Aij,Bij,Dij=∑zkNk∫zkzk+1Q¯ijk1,z,z2dz.Through the derivation of the above formulas, the kinematic relationships of the plate model and the stress results are given.


## 2.3. Energy Expressions

From the structural feature, we can infer that it is difficult to theoretically solve the boundary conditions and structural differential equations of the laminated composite double-plate structure system (LCDPS) and other more complex structures. Therefore, most researchers study the vibration characteristics of the complex plate structures by employing the numerical analysis method, such as the Rayleigh-Ritz method, variational method [53], and the Galerkin method [54]. In this section, we will derive expressions for various energies stored in the LCDPS and construct the energy function of the structure (Lagrange energy function: L) to prepare the Rayleigh-Ritz energy method. On the basis of the previous text, the strain energy (E) for the double-plate structure can be formulated as(9)Ep=12∫Vσxxpεxxp+σyypεyyp+τyzpγyzp+τxzpγxzp+τxypγxypdV,which can be further deduced as(10)Ep=12∫sNxpεxx0,p+Nypεyy0,p+Nxypεxy0,p+Mxpχxx0,p+Mypχyy0,p+Mxypχxyp+Qxpγxz0,p+Qypγyz0,pdS,where the symbolic variable p represents the number of plates in the laminated composite double-plate structure system. If no special instructions are given later, p=U represents the upper plate on the system (Plate 1 in Figure 1), and p=L denotes the lower plate on the system (Plate 2 in Figure 1). Then the strain energy expression can be written in three components: Ep=ESp+EBp+EBSp, in which ESp, EBp, and EBSp represent the stretching energy, the bending energy, and the bending-stretching coupling energy, separately. By substituting formulas (3) and (7) into formula (10), they can be expressed as follows:(11)Esp=12∫0La∫0LbA11εxx0,p2+2A12εxx0,pεyy0,p+2A16εxx0,pγxy0,p+A22εyy0,p2+2A26εyy0,pγxy0,p+A66γxy0,p2+κA44γxz0,p2+κA55γyz0,p2+2κA45γxz0,pγyz0,pdydx,=12∫0La∫0LbA11∂up∂x2+2A12∂up∂x∂vp∂y+2A16∂up∂x∂up∂y+∂vp∂x+A22∂vp∂y2+2A26∂vp∂y∂up∂y+∂vp∂x+A66∂up∂y+∂vp∂x2+κA44∂wp∂x+φxp2+κA55∂wp∂y+φyp2+2κA45∂wp∂x+φxp∂wp∂y+φypdydx,(12)EBp=12∫0La∫0LbD11χxxp2+2D12χxxpχyyp+D22χxxpχyyp+2D16χxxpχxyp+2D26χyypχxyp+D66χxyp2dydx,=12∫0La∫0LbD11∂φxp∂x2+2D12∂φxp∂x∂φyp∂y+2D16∂φxp∂x∂φxp∂y+∂φyp∂x+D22∂φyp∂y2+2D26∂φyp∂y∂φyp∂y+∂φxp∂x+D66∂φyp∂y+∂φxp∂x2dydx,(13)EBSp=12∫0La∫0LbB11εxx0,pχxxp+B12εxx0,pχyyp+B22εyy0,pχyyp+B16εxx0,pχxyp+B26εyy0,pχyyp+B16χxxp+B26χyyp+B26χxypγxy0,pdydx,=12∫0La∫0LbB11∂up∂x∂φxp∂x+B12∂up∂x∂φyp∂y+B16∂up∂x∂φxp∂y+∂φyp∂x+B22∂vp∂y2+B26∂vp∂y∂φyp∂y+∂φxp∂x+B16∂up∂x+B26∂vp∂y+B66∂φyp∂y+∂φxp∂x∂vp∂x+∂up∂ydydx.The kinetic energy function of the double-plate structure produced by vibration can be expressed as follows:(14)Tp=b2∫0La∫0Lbρpz∂Up∂t2+∂Vp∂t2+∂Wp∂t2dydx,=b2∫0La∫0LbI0p∂up∂t2+2I1p∂up∂t∂ϕxp∂t+I2p∂ϕyp∂t2+I0p∂vp∂t+2I1p∂vp∂t∂ϕyp∂t+I2p∂ϕyp∂t2+I0p∂wp∂tdydx,where the inertia terms can be obtained by the following formula:(15)I0p,I1p,I2p=∑k=1Nkp∫zkpzk+1pρpz1,z,z2dz.In the present model, the upper plate and the lower plate are connected by the middle Winkler-type elastic layer, which consists of a series of face springs, and the stiffness of the middle Winkler layer is represented byKc. When the system vibrates, the following formula denotes energy contained in the middle elastic layer:(16)Ec=12∫0La∫0LbKcwU−wLdydx.The major purpose is to establish a general solution model for the vibration characteristics of LCDPS under general boundary constraints, which are achieved by the penalty parameter method. By introducing five sets of penalty factors, the boundary conditions are added to the total energy function in the form of additional potential energy. Therefore, any specific elastic boundary constraints can be replaced by the proper penalty parameter. The resulting additional potential energy is defined as(17)Ebp=12∫0bku0x,pupxp,yp2+kv0x,pvpxp,yp2+kw0x,pwpxp,yp2+Kφx0x,pφxpxp,yp2+Kψy0x,pφypxp,yp2xp=0ku1x,pupxp,yp2+kv1x,pvpxp,yp2+kw1x,pwpxp,yp2+Kφx1x,pφxpxp,yp2+Kφy1x,pφypxp,yp2xp=adyp+12∫0aku0y,pupxp,yp2+kv0y,pvpxp,yp2+kw0y,pwpxp,yp2+Kφx0y,pφxpxp,yp2+Kφy0y,pφypxp,yp2yp=0ku1y,pupxp,yp2+kv1y,pvpxp,yp2+kw1y,pwpxp,yp2+Kφx1y,pφxpxp,yp2+Kφy1y,pφypxp,yp2yp=bdxp.We have thatkuΩχ,kvΩχ,kwΩχ,KφxΩχandKφxΩχΩ=0,1,χ=x,y represent the five sets of penalty factors of different boundaries at x=0, x=a, y=0, and y=b, which are stiffness coefficients of springs. It can be seen from the above formula that the general boundary conditions can be obtained as long as the penalty factors are given suitable values. It shows the superiority of this method in dealing with boundary conditions, which can deal with boundary conditions as easily as geometric parameters.


## 2.4. Displacement Admissible Functions and Solution Process

Based on the Hamilton principle [46, 47, 55], the governing equations of LCDPS on O-xyz coordinate system can be obtained as(18)∂Nx∂x+∂Nxy∂y=I0∂2u∂t2+I1∂2ϕx∂t2,∂Ny∂y+∂Nxy∂x=I0∂2v∂t2+I1∂2ϕy∂t2,∂Qx∂x+∂Qy∂y=I0∂2w∂t2,∂Mx∂x+∂Mxy∂y−Qx+mx=I1∂2u∂t2+I2∂2ϕx∂t2,∂My∂y+∂Mxy∂x−Qy+my=I1∂2v∂t2+I2∂2ϕy∂t2.It is obvious that there are the second derivatives in different governing equations. In order to avoid a discontinuity at boundary conditions, the admissible function of displacement is expressed by trigonometric Fourier series [56]. In addition, the supplementary function is introduced to ensure the existence of the derivative and the continuity of functions. The modified Fourier series technique is concise and efficient. This method does not need to solve partial differential equations like other numerical methods [57, 58], thus avoiding numerical errors. The whole solving process is based on the energy integration of the structure, and the Fourier series used in this paper can obtain the accurate solution to the integral calculation because of its own integral property [59, 60]. The displacement admissible function can be written as(19)upx,y,t=e−jωt∑m=0M∑n=0NAmnpcosλmxcosλny+∑l=12∑n=0Nalnpplxcosλny∑l=12∑m=0Mblmpplycosλmx,vpx,y,t=e−jωt∑m=0M∑n=0NBmnpcosλmxcosλny+∑l=12∑n=0Nclnpplxcosλny∑l=12∑m=0Mdlmpplycosλmx,wpx,y,t=e−jωt∑m=0M∑n=0NCmnpcosλmxcosλny+∑l=12∑n=0Nelnpplxcosλny∑l=12∑m=0Mflmpplycosλmx,ϕxpx,y,t=e−jωt∑m=0M∑n=0NDmnpcosλmxcosλny+∑l=12∑n=0Nglnpplxcosλny∑l=12∑m=0Mhlmpplycosλmx,ϕypx,y,t=e−jωt∑m=0M∑n=0NEmnpcosλmxcosλny+∑l=12∑n=0Nilnpplxcosλny∑l=12∑m=0Mjlmpplycosλmx.where λm=mπ/La, λn=nπ/Lb, and M and N are the truncation number.Amnp, Bmnp, Cmnp, Dmnp, and Emnp express coefficients expansion of two-dimensional standard cosine function of Fourier series. Accordingly, alnp, blmp, clnp, dlmp, elnp, flmp, glnp, hlmp, ilnp, and jlmp represent the corresponding coefficients of the supplementary terms.(20)p1x=xxa−12,p2x=x2axa−1,p1y=yyb−12,p2y=y2byb−1.It is easy to verify that(21)p10=p1a=p1′a=0,p1′0=1,p20=p2a=pa′0=0,p2′a=1,p10=p1b=p1′a=0,p1′0=1,p20=p2b=p2′0=0,p2′b=1.As explained in [8, 9, 47, 54, 55], the design of modified Fourier series should satisfy the second-order partial differential governing equations of plate. It is easy to see from equation (21) that the function complementary p can make sure that the first derivative of the admission displacement functions is continuous in the whole solution domain and the second derivative of each point exists. So displacement function and the derivative term have a continuous derivative across the edge, and the potential discontinuities with the displacement function and its derivatives at the end points can be removed.Substituting equations (2) and (4) into equation (18), we can get matrix equation as(22)L11UL12U0L14UL15U00000L12UL22U0L14UL15U0000000L33UL34UL35U00LkL00L14UL24UL34UL44UL45U00000L15UL25UL35UL45UL55U0000000000L11LL12L0L14LL15L00000L12LL22L0L24LL25L00LkU0000L33LL34LL35L00000L14LL24LL34LL44LL45L00000L15LL25LL35LL45LL55L0−ω2M11U00M14U0000000M22U00M25U0000000M33U0000000M14U00M44U0000000M25U00M55U0000000000M11L00M14L0000000M22L00M25L0000000M33L0000000M14L00M44L0000000M25L00M55Lupvpwpφxpφyp=00000,where L and MLij=Lji,Mij=Mji indicate the coefficients of the linear operator. They can be represented as(23)L11p=A11p∂2∂x2+2A16p∂2∂x∂y+A66p∂2∂y2,L12p=A16p∂2∂x2+A12p+A66p∂2∂x∂y+A26p∂2∂y2,L13p=0,L23p=0,L14p=B11p∂2∂x2+2B16p∂2∂x∂y+B66p∂2∂y2,L15p=B16p∂2∂x2+B12p+B66p∂2∂x∂y+B26p∂2∂y2,L22p=A66p∂2∂x2+2A26p∂2∂x∂y+A22p∂2∂y2,L24p=B16p∂2∂x2+B12p+B66p∂2∂x∂y+B26p∂2∂y2,L25p=B66p∂2∂x2+2B26p∂2∂x∂y+B22p∂2∂y2,L33p=−A55p∂2∂x2−2A45p∂2∂x∂y−A44p∂2∂y2,L34p=−A55p∂∂x−2A45p∂∂y,L35p=−A45p∂∂x−2A44p∂∂y,L44p=D11p∂2∂x2+2D16p∂∂x∂y+D66p∂2∂y2−A55p,L45p=D16p∂2∂x2+D12p+D66p∂2∂x∂y+D26p∂2∂y2−A45p,L55p=D66p∂2∂x2+2D26p∂2∂x∂y+D22p∂2∂y2−A44p,M11p=M22p=M33p=−I0p∂2∂t,M14p=M25p=−I1p∂2∂t,M44p=M55p=−I2p∂2∂t.One can note that there are second-order derivatives of all the displacements and rotation components of the plate, as we have posted before. After deriving the expressions of dynamic energy and potential energy, the Lagrange function of the double-plate structure can be written as(24)L=Tp+Wp−Ep−Espp−Ec,(25)∂L∂Ξmnp=0,andΞp=Ap,Bp,Cp,Dp,Ep,m=0,1,…M;n=0,1,…N,∂L∂Ψlnp=0,andΨp=ap,cp,ep,gp,ip,l=1,2;n=0,1,…N,∂L∂ϒlmp=0,andϒp=bp,dp,fp,hp,jp,l=1,2;m=0,1,…M.By substituting equations (11)–(15) and (19) into equation (24) and then carrying out the Rayleigh-Ritz energy method, we can obtain a set of linear algebraic equations to determine the unknown Fourier coefficients for eigenfunctions. The number of linear equations is 5∗2∗, M+1, ∗, N+1 + 10∗2∗, M+N+2, which can be turned into the matrix as(26)K−ω2MΠ=0.We have thatΠ=AUBUCUDUEUALBLCLDLELT are the unknown Fourier coefficients; K stands for the stiffness matrix and M represents the mass matrix, separately.Please refer to Appendix A for the detailed expression of the above matrices. Obviously, by solving the eigenproblem of the standard matrix, the natural frequency and eigenvector can be easily obtained. Once the eigenvector of the given frequencyΠ is determined, the displacement of the double-plate structure system can be determined by substituting the coefficients into the equation.


## 3. Numerical Calculation and Analysis

In this part, the fast convergence and accuracy of this method are verified by a series of LCDPS examples under different boundary conditions. There are two main parts in this section. First of all, in order to prove that the methods described can be applied to LCDPS stably and reliably, the convergence characteristics and calculation accuracy of the present method are studied. Then, the vibration characteristics of LCDPS with changing parametric are presented. Some new calculation results would enrich the conclusions in the field. The solution and data analysis show that the model is feasible and effective. Besides, the symbolic representations of boundary conditions for LCDPS should be noted. C, S, F, and E represent clamped boundary, simply supported boundary, free boundary, and elastic boundary, respectively. Set an example of FFFF-CCCC boundary conditions, which show that the four boundary conditions of the upper plate (Plate 1) in anticlockwise order from the borderx=0, y=0, x=a, and y=b are all free constraints. CCCC represents that four boundaries of the lower plate are all clamped boundaries. Moreover, the dimensionless natural frequency parameters, if not specified, are expressed as Ω=ωa2ρu/E2uhu2.


### 3.1. Convergence and Validation Study of the LCDPS

In this section, the material parameters of the upper plate in the double-plate structure system are the same as those of the lower plate. A kind of carbon fiber composite material was taken as the sample, which has the advantages of light weight, high strength, and low cost. The material characteristics parameters are defined asG12=G13=5GPa, G23=6GPa, E1=150GPa, E2=10GPa, μ12=0.25, ρ=1500kg/m3, μ12E1=μ21E2, and αfiber=30°−30°30°−30°. In engineering applications, different laminates will be used according to different structural forms and design requirements. This kind of typical material is often used in numerical simulation, such as [61]. The geometric characteristic parameters of the plate are set as b/a=0.5 and h/a=0.1. The elasticity coefficient of the intermediate Winkler layer is set to Kc=1014. In the calculation, the truncation indexes M and N in the extended repressions are finite. Therefore, it is necessary to explore the reasonable truncation index to meet the need for accuracy. Table 1 records the first five frequency parameters of the LCDPS under various classical boundary conditions when the displacement truncation index is different.Table 1
Convergence of the first five frequency parameters Ω for an [30°/−30°/30°/−30°] LCDPS with CCCC-CCCC and FFFF-FFFF conditions.
BCΩM and N2∗23∗34∗48∗812∗1214∗1416∗1618∗1820∗20CCCC-CCCC132.87732.59432.52532.50332.50232.50232.50132.50132.501246.03645.71945.66145.63245.62945.62945.62945.62945.629366.06864.64462.60362.50962.50462.50362.50362.50362.503467.64765.13164.61964.57164.56964.56964.56964.56964.569576.25574.99274.95274.88774.88374.88374.88274.88274.882FFFF-FFFF113.19113.08413.06713.05913.05813.05813.05813.05813.058217.48816.81216.80216.74916.74716.74716.74616.74616.746329.18429.11029.05929.03729.03529.03529.03429.03429.034431.11731.01330.90230.87030.86730.86730.86730.86730.867533.88833.39533.24233.17133.16633.16633.16533.16533.165It can be seen from Table1 that, for CCCC-CCCC boundary conditions, when the truncation values of displacement admissible functions M and N are 8∗8 and 20∗20, respectively, the difference of frequency results is very small. The maximum difference is only 0.01%. Then, the vibration frequency of FFFF-FFFF boundary conditions was presented. The maximum difference of calculation results when the truncation values of displacement admissible functions M and N are 8∗8 and 20∗20 is 0.02%. It can also be got from Table 1 that when the truncation number of displacement expansions is greater than 16∗16, the dimensionless natural frequency remains unchanged within the set of significant number. Thus, from Table 1, we get that the frequency parameters of the LCDPS with different borderline conditions always can converge fast to a certain solution and not be sensitive to increase of the cutoff values. However, when the truncation number is larger and larger, the calculation cost is also higher and higher. Thus, in the following calculation, in order to ensure the accuracy of calculation results, the cutoff value of displacement function is taken as M=N=14.In this structural model, the constraints on the boundary are simulated by weight parameters. Based on the FSDT, boundary parameters can be changed by setting five sets of spring stiffness values. When the boundary condition is fixed constraint, the spring stiffness should be infinite in theory. But, in the actual calculation, it is difficult to realize. Therefore, it is necessary to determine the spring stiffness through numerical calculation, so that the boundaries of plates are clamped approximately. The convergent curves of boundary weight parameters are given in Figure2. The boundary condition of LCDPS is defined as two edges of x=0 and x=a are clamped boundaries (C-C, C-C), and two ends of y=0 and y=b are elastic boundaries (E-E, E-E). In the calculation, there is only one set of spring stiffness values on the elastic boundary varying between 101 and 1018, while the stiffness values of other groups of springs remain unchanged (1014) as an infinite constraint. For instance, the stiffness of linear springs ku increased from 100 to 1018, and the other stiffness values of springs (ky, kw, Kφx, and Kφy) are uniformly set as 1014. The geometric parameters and material parameters of the LCDPS are consistent with Table 1, except the different laying schemes. The fiber laying scheme in this example is αfiber=0°90°0°90°.Figure 2
Variation of dimensionless natural frequency versus spring parameters for an LCDPS with different mode order. (a) 1st mode; (b) 2nd mode; (c) 3rd mode; (d) 4th mode.
(a)(b)(c)(d)As we can see from Figure2, due to changes in the boundary spring stiffness coefficient, the frequency parameters of LCDPS change accordingly. It can be found that when the spring stiffness value is less than 107, the dimensionless natural frequency parameter remains unchanged. This is mainly because the spring stiffness value is too small to affect the vibration of the structure. We should also pay attention from the calculation results that the change of vibration frequency is very sharp when the spring stiffness is in a certain range of 107 to 1012. The stiffness of the boundary spring is large enough to affect the structure, so the resonance frequency increases rapidly with the increase of the stiffness parameter. When the spring stiffness value is larger than 107, with the increase of stiffness parameters, the natural dimensionless frequency parameter of the structure does not increase. This shows that when the boundary spring stiffness value is large enough, the boundary of the structure will be restricted and cannot be displaced. Thus, the approximate spring stiffness value of the clamped restrained boundary can be deduced, and the reasonable selection range is suggested. In the subsequent calculation, the spring stiffness coefficient can be selected as 1014 to simulate the clamped boundary condition. It should also be noted that different types of boundary springs represent different boundary constraints, and their effects on the vibration characteristics of LCDPS are quite different. It is found that the stiffness of the spring with boundary constraint kΩw, that is, the transverse shear force, has the most important effect on the vibration frequency of the LCDPS. Based on the above analysis, the classical boundary conditions and corresponding boundary spring stiffness values of the LCDPS model in this paper are given in Table 2.Table 2
Weight parameters of LCDPS under different conditions.
BCkΩukΩvkΩwKΩxKΩyF00000S10141014101401014SD01014101401014C10141014101410141014E110810810810141014E21014101410141081014E3108108108108108On the basis of convergence verification in the previous section, this section further verifies the accuracy and effectiveness of the established model by comparing the calculation data with the results in the literature. Table3 shows the calculation results of this paper and compares them with the solutions in literature with the same boundary conditions of CCCC-CCCC and SSSS-SSSS.Table 3
Comparison of the natural frequency parameters Ω for an [45°/60°/45°/60°] LCDPS with various conditionsΩ=ωa2/π2ρu/E2uhu2.
CCCC-CCCCSSSS-SSSSKc=0ΩRef [60]Ref [57]Ref [58]Kc=0ΩRef [60]Ref [57]Ref [58]ΩPresentΩPresent12.284912.28552.28572.278511.880111.88031.87881.846622.2849————21.8801————33.736123.73633.73923.738333.374023.37633.37763.377443.7361————43.3740————53.978733.97923.98133.958353.692233.69233.69213.642563.9787————63.6922————75.177345.17775.17995.183674.965244.96654.96804.966185.1773————84.9652————Kc=1e14————Kc=1e14————12.284912.28552.28572.278511.880111.88031.87881.846623.736123.73633.73923.738323.374023.37633.37763.377433.978733.97923.98133.958333.692233.69233.69213.642545.177345.17775.17995.183644.965244.96654.96804.966155.695255.69665.70195.680855.479855.48015.48345.434865.840265.84165.84545.806665.592865.59855.60005.518776.712276.71296.71666.735876.541676.54396.54736.562887.040087.04017.04497.004686.909386.91096.91356.8375For LCDPS, the upper plate and the lower plate are the same, and the carbon fiber composite material parameters and geometric parameters were of the following values:E1=400GPa, E1/E2=40, G12=G13=5/6, G23=6GPa, μ12=0.25, ρ=1600kg/m3, αfiber=45°−45°45°−45°45°, a=b=1m, and hU=hL=a/10. The stiffness of the middle elastic layer Kc=0 represents free connection and Kc=1014 represents rigidity connection, respectively. The dimensionless natural frequency is Ω=ωa2/π2ρu/E2uhu2. The results of Ye et al. [60], Karami et al. [57], and Liew et al. [58] were selected as comparative data. From the calculation results in Table 3, when the rigidity of the middle layer Kc=0, the phenomenon of “repeated frequency” appears obviously. It can be seen that the 1st order and the 2nd order have the same frequencies and the 3rd order and the 4th order are in the same situation. The reason for this phenomenon is that when the stiffness of the middle elastic layer is zero, the system degenerates into two independent and parallel composite laminates. As a result, the frequencies of the two plates will alternate. The “repetition frequency” of the LCDPS can be eliminated in such a situation that the stiffness coefficient of the Winkler layer Kc=1014. Then the rigid connected double-layer plate can be regarded as a single plate. At the same time, it is found that the natural frequency obtained at present is in accordance with the literature results of the study. The maximum error occurs when the boundary is SSSS-SSSS and Kc=1014. The difference between calculation results in this paper and the calculation results in [58] is only 1.325%, which appeared in the 7th-order frequency, and other differences between them are minor. As mentioned above, the validity of the analysis model in this paper is verified.As an effective numerical simulation method, the experimental verification is necessary. In this part, we will select a fully clamped rectangular plate as experimental objects to obtain its natural frequency. As depicted in Figure3, the method of single point vibration pick-up is used to excite the rectangular plate with a force hammer. The alloy steel plate is taken as the subject of research, which has the advantages of high strength, hardness, and toughness. The geometric parameters of this plate are a=0.6m, b=0.35m, and h=0.01m. The material properties are given as E=216GPa, μ=0.3, and ρ=7800kg/m3. For the isotropic plate with single-layered material, the setting of material parameters is given in Section 2.2. In Figure 4, the vibration modes and corresponding natural frequencies of the plate obtained from experimentation and the solution by the present method are compared. As can be seen from Figure 4, both modal and frequency results are in good agreement.Figure 3
The experimental setting for a rectangular plate. (a) Experimental plate. (b) Single vibration excitation.
(a)(b)Figure 4
Comparison between the present solutions and experimental results. (a) Present solutions. (b) Experimental results.
(a)(b)


### 3.2. Parametric Analysis of the LCDPS

In the last section, the convergence and efficiency of this method are obtained. In this part, a parametric approach to free vibration characteristics of LCDPS is carried out. Tables4 and 5 show the frequency parameters of LCDPS under different laying schemes and boundary conditions. The stiffness of the middle elastic layer between two plates is set as Kc=1014. In the calculation, the material parameters of the double plates are the same and consistent with Table 3. Geometric parameters are defined as a=b=1m and hU=hL=a/50, and the laying schemes of composite laying are as follows: αfiber=0°90°90°0°, αfiber=30°60°30°60°, and αfiber=45°−45°45°−45°. The model constructed in this paper can easily simulate different boundary conditions by setting different spring stiffness coefficients, such as CCCC-CCCC, FCFC-FCFC, and CCCC-FFFF. There is no need to change the model or rewrite the program. Such an approach shows great flexibility and adaptability. It can be seen from the calculation results in Table 4 that the vibration frequency of LCDPS changes with the stiffness of boundary springs. Under CCCC-CCCC boundary condition, the system has the maximum frequency parameter. When the boundary condition is FCFC-FCFC, the vibration frequency parameters of the system are obviously reduced. Different boundary conditions of top and bottom plates can also be set, such as CCCC-FFFF, and, with the increase of boundary degrees of freedom, the frequency decreases further. This phenomenon is due to the change of plate stiffness caused by different boundary constraints. The structural stiffness of the LCDPS in the clamped restrained boundary is the largest and decreases with boundary constraints relaxing gradually. In addition, we can observe from Table 4 that different composite layer design schemes will cause marked difference in the vibration characteristics of the structure. This phenomenon is mainly due to the fact that the layer scheme has a significant impact on the structural stiffness of the LCDPS.Table 4
Frequency parameters Ω for an LCDPS with different layer scheme and classical conditions.
LayerΩBCCCCC-CCCCCSCS-CSCSCFCF-CFCFSSSS-SSSSSSSS-FFFFCCCC-FFFFSSSS-CFCFSSSS-CSCSSSSS-CCCC[0°/90°/90°/0°]139.62539.62436.05223.3690.0090.00434.78135.26535.265256.46456.46236.39446.4480.0650.02751.28953.44453.444390.21190.21040.60766.9420.0790.09783.93585.55985.559494.95894.95857.93579.29122.63034.78185.33788.31188.3115104.017104.01692.35484.25343.94751.28994.23095.48495.4846126.654126.65393.150107.30966.65083.934116.436119.713119.714[30°/60°/30°/60°]131.38930.39318.35124.1840.0000.00527.61328.98329.554251.61350.05219.01040.1940.0300.19945.97547.96348.738371.37969.53230.38459.2130.3470.35664.38566.82967.805471.83870.07243.95761.20722.29428.42064.78967.26168.312594.37192.18950.15180.68138.48147.17485.97888.89890.0336105.633103.18754.91689.83856.57165.72596.58399.712100.827[45°/−45°/45°/−45°]135.38435.31017.85930.0450.0000.00332.75534.09334.129269.49569.27425.91056.5100.0420.24664.83265.97266.080369.49569.32546.97166.2210.0930.39864.83568.39168.4714108.122107.56748.71496.8560.39828.198101.761104.911105.1835116.321115.96859.67999.04228.10554.528109.399110.602110.7866117.213116.95279.567114.60454.50762.114110.069115.873115.992Table 5
Frequency parameters Ω for an LCDPS with different middle spring parameters, layer scheme, and classical conditions.
LayerΩKc=0Kc=109Kc=1014CCCC-CCCCCSDCSD-CSDCSDCFCF-CFCFCCCC-CCCCCSDCSD-CSDCSDCFCF-CFCFCCCC-CCCCCSDCSD-CSDCSD|CFCF-CFCF[45°/−45°]1114.26414.25211.93714.26414.25211.93714.26414.25211.937214.26414.25211.93719.56019.52713.12019.56019.52713.120319.56019.52713.12027.27427.21917.04727.27427.21917.047419.56019.52713.12034.24734.22022.83334.24734.22022.833527.27427.21917.04737.05536.97830.46337.05536.97830.463[45°/−45°]2120.74920.69417.00220.74920.69417.00220.74920.69417.002220.74920.69417.00228.26928.11518.62428.26928.11518.624328.26928.11518.62438.82438.58624.59738.82438.58624.597428.26928.11518.62445.01944.99432.86247.47547.37632.862538.82438.58624.59747.47547.37642.76251.76851.47442.762[45°/−45°]3121.60421.54017.68921.60421.54017.68921.60421.54017.689221.60421.54017.68929.37529.19719.35829.37529.19719.358329.37529.19719.35840.24939.97825.57440.24939.97825.574429.37529.19719.35845.42445.39434.13149.09248.98234.131540.24939.97825.57449.09248.98243.69953.53953.21044.189[45°/−45°]4121.88721.82017.91821.88721.82017.91821.88721.82017.918221.88721.82017.91829.73929.55319.60229.73829.55319.602329.73829.55319.60240.71740.43425.89740.71740.43425.897429.73829.55319.60245.56145.52934.54949.62149.50734.549540.71740.43425.89749.62149.50743.79354.11853.77744.658As mentioned above, with the solution of the vibration characteristic equation, the eigenvector corresponding to the natural frequency can be obtained as well. Then substitute the eigenvectors into the displacement admissible function of the system. From the first four modes simulated in Figure5 corresponding to Table 4, we can comprehend the natural vibration characteristics directly and vividly. It is shown in the figure that the mode shapes are different with varying boundary conditions, but the vibration modes of the top and bottom plate structures of the LCDPS are the same.Figure 5
The lowest four mode shapes for an LCDPS with several conditions.
(a)(b)(c)(d)(e)In the next section, Tables5 and 6 show the first five natural frequency parameters of the LCDPS under the classical and elastic boundary conditions with a changing number of layers of the layering scheme. The geometric parameters and physical properties are exactly the same as in Table 4. Set the stiffness of the middle elastic layer to Kc=0, Kc=109, and Kc=1014 respectively. From Table 5, we can see the trends from the calculation results that the vibration behaviors of its structure changed somewhat with the increase of the fiber layers. This effect will be more significant with lower layers number. As fiber layer counts increase, the properties of composite materials are close to those of isotropic materials. Table 6 presents the vibration characteristics of the LCDPS with elastic boundary conditions. By changing the stiffness of the Winkler layer, different frequency values can be obtained. While the stiffness coefficient of the middle layer has little effect on the lower-order frequency and the influence becomes larger with the increase of the order, we can infer that the stiffness of the middle layer maybe has an obvious effect on the higher-order frequency of the structure that will be referred to later.Table 6
Frequency parameters Ω for an LCDPS with different middle spring parameters, layer scheme, and elastic conditions.
LayerΩKc=0Kc=109Kc=1014E1E1E2E2-CCCCE1E1E3E3-CCCCE2E2E3E3-CCCCE1E1E2E2-CCCCE1E1E3E3-CCCCE2E2E3E3-CCCCE1E1E2E2-CCCCE1E1E3E3-CCCCE2E2E3E3-CCCC[45°/−45°]1111.95310.57011.87813.68413.22413.65614.22314.21814.178214.26414.19914.26418.65214.92818.60919.52414.92819.473316.41914.26416.28625.77615.30225.70927.24215.30327.184419.56014.92819.56030.08217.88030.04234.14919.50834.048522.52215.29822.35834.16320.36734.07237.02820.36836.960[45°/−45°]2115.05012.45114.90119.15214.92819.08620.63414.92820.505220.74914.92820.74925.86715.30425.74328.18015.30428.045321.09915.30320.79234.53618.05534.34638.75820.36938.612428.26916.76828.26936.83220.36936.71547.23920.61846.997529.26120.36828.88442.64123.96242.42349.83728.13249.817[45°/−45°]3115.38212.62015.22319.82214.92919.74921.47514.92921.331221.60414.92921.30026.70815.30426.56929.27615.30429.127321.63315.30421.60435.45118.62335.24440.17820.36940.019429.37517.01629.37537.54420.36937.40648.83621.45748.573530.11120.36929.69343.51324.61743.27549.84829.22249.839[45°/−45°]4115.08812.72814.96018.64215.77718.58719.85815.77719.755219.95015.77719.95025.28016.21025.18427.20716.21027.098320.96916.21020.72234.09317.70233.94437.49919.84637.381427.28017.08727.28037.09821.55337.02045.74321.55345.544528.85719.95028.55842.82723.66942.65850.13327.16950.003In addition, Figures6 and 7 show the influence of the fiber laying angle on the LCDPS under classical or elastic boundary conditions. Maintaining the original material parameters as in Table 6, the fiber angle laying scheme of the composite plates in the structure is defined as αfiberp=90°ααfiberp90°, in which ααfiberpp=UandL represent the laying angles of fiber in the middle surface of top and bottom plate, respectively. Changing the fiber angle of the lower plate from 10° to 180°, the fiber angle has a different influence on the vibration characteristics of LCDPS with different stiffness coefficients of intermediate Winkler elastic layers. It can be concluded that, at any boundary condition, in both the lower-order frequency range and the higher-order frequency range, the change of natural frequency parameter is symmetrical about the frequency of the system, whose fiber angle is 90°.Figure 6
Variation of dimensionless natural frequency versus fiber angle for an LCDPS with different conditions. (a) CCCC-CCCC, 1st mode; (b) SSSS-SSSS, 2nd mode; (c) CCCC-CCCC, 3rd mode; (d) SSSS-SSSS, 4th mode.
(a)(b)(c)(d)Figure 7
Variation of dimensionless natural frequency versus fiber angle for an LCDPS with different boundary conditions. (a) FCFC-FCFC, 5th mode; (b) FFCC-FFCC, 6th mode; (c) FCCC-FCCC, 7th mode; (d) FFFF-FFFF, 8th mode.
(a)(b)(c)(d)In the following examples, the influences of the geometric parameters of the LCDPS on the vibration characteristics are demonstrated. Keep material parameters of the top and bottom plates always be the same carbon fiber composite material, and the material parameters are assumed as follows:E1=740GPa, E1/E2=40, G12=G13=111GPa, G23=E1/8, μ12=0.25, ρ=1600kg/m3, and αfiber=30∘0∘−30∘0∘30∘. The geometrical parameters of LCDPS are defined as follows: a=b=1m, hU=a/50, and hL/hU=1∼20. Firstly, the effect of the thickness of composite plates on the vibration frequency of LCDPS is investigated. From Figure 8, we can see that in a certain thickness range the vibration frequency of LCDPS increases rapidly with the thickness. There is a unified trend with any boundary condition; when the thickness coefficient hL/hU increases to a certain value, the vibration frequency of the LCDPS will not change obviously. The structural vibration characteristics tend to be stable. Secondly, the effects of plate aspect ratio on the vibration behavior are studied. The material characteristics of the plates are the same and consistent with Figure 7. There are four layers for one laminated plate, and the geometric parameters are as follows: b/a=1∼10, hU=0.05, hL/hU=1∼20, and αfiber=0°90°0°90°. The result demonstrates in Figure 9 that the increase of plate length ratio causes the raise of vibration frequency of the system. It shows that the change of geometric parameters affects the vibration characteristic stiffness matrix of the LCDPS in a certain range. Taking advantage of this feature, the vibration characteristics can be adjusted by changing different geometry and material parameters. But, over the sensitive range, the effect of adjustment by geometrical parameters will be limited.Figure 8
Variation of dimensionless natural frequency versus thickness for an LCDPS with different conditions. (a) CCCC-CCCC; (b) SDSDSDSD-SDSDSDSD; (c) E1E1E1E1-CCCC.
(a)(b)(c)Figure 9
Variation of dimensionless natural frequency versus length-width ratio for an LCDPS with different conditions. (a) CCCC-CCCC; (b) SSSS-SSSS; (c) E1E1E1E1-CCCC.
(a)(b)(c)As we have mentioned earlier, Figure10 shows the regularities of an LCDPS on frequency parameter with the stiffness of the middle Winkler layer changing. The geometric and materials parameters are consistent with Figure 7. From the first 10 modal frequencies given, the change rule of the frequency parameters of the system can be observed intuitively. Firstly, with the increasing stiffness coefficient of the middle elastic layer, the frequency parameters from the first order to the tenth order have a common increasing trend. Secondly, the elastic stiffness coefficient of the middle Winkler elastic layer has a varied sensitive area, which is similar to the boundary condition. In this region, the vibration characteristics of LCDPS are very sensitive to the stiffness coefficient of the Winkler layer, and the influence of stiffness parameters on higher-order frequency is more apparent. This view improves the point mentioned earlier. Outside this region, the influence is weak. Therefore, the adjustment of the elastic parameters of the middle Winkler layer is another way to obtain the ideal vibration characteristic besides certain parameters such as boundary conditions, thickness coefficient, and aspect length-width ratio.Figure 10
Variation of frequency parameter Ω versus middle spring stiffness for an LCDPS with different boundary conditions. (a) CCCC-CCCC; (b) CCCC-SSSS; (c) CFCF-CFCF; (d) CCCC-E1E1E1E1; (e) CCCC-E2E2E2E2; (f) CCCC-E3E3E3E3.
(a)(b)(c)(d)(e)(f)


## 3.1. Convergence and Validation Study of the LCDPS

In this section, the material parameters of the upper plate in the double-plate structure system are the same as those of the lower plate. A kind of carbon fiber composite material was taken as the sample, which has the advantages of light weight, high strength, and low cost. The material characteristics parameters are defined asG12=G13=5GPa, G23=6GPa, E1=150GPa, E2=10GPa, μ12=0.25, ρ=1500kg/m3, μ12E1=μ21E2, and αfiber=30°−30°30°−30°. In engineering applications, different laminates will be used according to different structural forms and design requirements. This kind of typical material is often used in numerical simulation, such as [61]. The geometric characteristic parameters of the plate are set as b/a=0.5 and h/a=0.1. The elasticity coefficient of the intermediate Winkler layer is set to Kc=1014. In the calculation, the truncation indexes M and N in the extended repressions are finite. Therefore, it is necessary to explore the reasonable truncation index to meet the need for accuracy. Table 1 records the first five frequency parameters of the LCDPS under various classical boundary conditions when the displacement truncation index is different.Table 1
Convergence of the first five frequency parameters Ω for an [30°/−30°/30°/−30°] LCDPS with CCCC-CCCC and FFFF-FFFF conditions.
BCΩM and N2∗23∗34∗48∗812∗1214∗1416∗1618∗1820∗20CCCC-CCCC132.87732.59432.52532.50332.50232.50232.50132.50132.501246.03645.71945.66145.63245.62945.62945.62945.62945.629366.06864.64462.60362.50962.50462.50362.50362.50362.503467.64765.13164.61964.57164.56964.56964.56964.56964.569576.25574.99274.95274.88774.88374.88374.88274.88274.882FFFF-FFFF113.19113.08413.06713.05913.05813.05813.05813.05813.058217.48816.81216.80216.74916.74716.74716.74616.74616.746329.18429.11029.05929.03729.03529.03529.03429.03429.034431.11731.01330.90230.87030.86730.86730.86730.86730.867533.88833.39533.24233.17133.16633.16633.16533.16533.165It can be seen from Table1 that, for CCCC-CCCC boundary conditions, when the truncation values of displacement admissible functions M and N are 8∗8 and 20∗20, respectively, the difference of frequency results is very small. The maximum difference is only 0.01%. Then, the vibration frequency of FFFF-FFFF boundary conditions was presented. The maximum difference of calculation results when the truncation values of displacement admissible functions M and N are 8∗8 and 20∗20 is 0.02%. It can also be got from Table 1 that when the truncation number of displacement expansions is greater than 16∗16, the dimensionless natural frequency remains unchanged within the set of significant number. Thus, from Table 1, we get that the frequency parameters of the LCDPS with different borderline conditions always can converge fast to a certain solution and not be sensitive to increase of the cutoff values. However, when the truncation number is larger and larger, the calculation cost is also higher and higher. Thus, in the following calculation, in order to ensure the accuracy of calculation results, the cutoff value of displacement function is taken as M=N=14.In this structural model, the constraints on the boundary are simulated by weight parameters. Based on the FSDT, boundary parameters can be changed by setting five sets of spring stiffness values. When the boundary condition is fixed constraint, the spring stiffness should be infinite in theory. But, in the actual calculation, it is difficult to realize. Therefore, it is necessary to determine the spring stiffness through numerical calculation, so that the boundaries of plates are clamped approximately. The convergent curves of boundary weight parameters are given in Figure2. The boundary condition of LCDPS is defined as two edges of x=0 and x=a are clamped boundaries (C-C, C-C), and two ends of y=0 and y=b are elastic boundaries (E-E, E-E). In the calculation, there is only one set of spring stiffness values on the elastic boundary varying between 101 and 1018, while the stiffness values of other groups of springs remain unchanged (1014) as an infinite constraint. For instance, the stiffness of linear springs ku increased from 100 to 1018, and the other stiffness values of springs (ky, kw, Kφx, and Kφy) are uniformly set as 1014. The geometric parameters and material parameters of the LCDPS are consistent with Table 1, except the different laying schemes. The fiber laying scheme in this example is αfiber=0°90°0°90°.Figure 2
Variation of dimensionless natural frequency versus spring parameters for an LCDPS with different mode order. (a) 1st mode; (b) 2nd mode; (c) 3rd mode; (d) 4th mode.
(a)(b)(c)(d)As we can see from Figure2, due to changes in the boundary spring stiffness coefficient, the frequency parameters of LCDPS change accordingly. It can be found that when the spring stiffness value is less than 107, the dimensionless natural frequency parameter remains unchanged. This is mainly because the spring stiffness value is too small to affect the vibration of the structure. We should also pay attention from the calculation results that the change of vibration frequency is very sharp when the spring stiffness is in a certain range of 107 to 1012. The stiffness of the boundary spring is large enough to affect the structure, so the resonance frequency increases rapidly with the increase of the stiffness parameter. When the spring stiffness value is larger than 107, with the increase of stiffness parameters, the natural dimensionless frequency parameter of the structure does not increase. This shows that when the boundary spring stiffness value is large enough, the boundary of the structure will be restricted and cannot be displaced. Thus, the approximate spring stiffness value of the clamped restrained boundary can be deduced, and the reasonable selection range is suggested. In the subsequent calculation, the spring stiffness coefficient can be selected as 1014 to simulate the clamped boundary condition. It should also be noted that different types of boundary springs represent different boundary constraints, and their effects on the vibration characteristics of LCDPS are quite different. It is found that the stiffness of the spring with boundary constraint kΩw, that is, the transverse shear force, has the most important effect on the vibration frequency of the LCDPS. Based on the above analysis, the classical boundary conditions and corresponding boundary spring stiffness values of the LCDPS model in this paper are given in Table 2.Table 2
Weight parameters of LCDPS under different conditions.
BCkΩukΩvkΩwKΩxKΩyF00000S10141014101401014SD01014101401014C10141014101410141014E110810810810141014E21014101410141081014E3108108108108108On the basis of convergence verification in the previous section, this section further verifies the accuracy and effectiveness of the established model by comparing the calculation data with the results in the literature. Table3 shows the calculation results of this paper and compares them with the solutions in literature with the same boundary conditions of CCCC-CCCC and SSSS-SSSS.Table 3
Comparison of the natural frequency parameters Ω for an [45°/60°/45°/60°] LCDPS with various conditionsΩ=ωa2/π2ρu/E2uhu2.
CCCC-CCCCSSSS-SSSSKc=0ΩRef [60]Ref [57]Ref [58]Kc=0ΩRef [60]Ref [57]Ref [58]ΩPresentΩPresent12.284912.28552.28572.278511.880111.88031.87881.846622.2849————21.8801————33.736123.73633.73923.738333.374023.37633.37763.377443.7361————43.3740————53.978733.97923.98133.958353.692233.69233.69213.642563.9787————63.6922————75.177345.17775.17995.183674.965244.96654.96804.966185.1773————84.9652————Kc=1e14————Kc=1e14————12.284912.28552.28572.278511.880111.88031.87881.846623.736123.73633.73923.738323.374023.37633.37763.377433.978733.97923.98133.958333.692233.69233.69213.642545.177345.17775.17995.183644.965244.96654.96804.966155.695255.69665.70195.680855.479855.48015.48345.434865.840265.84165.84545.806665.592865.59855.60005.518776.712276.71296.71666.735876.541676.54396.54736.562887.040087.04017.04497.004686.909386.91096.91356.8375For LCDPS, the upper plate and the lower plate are the same, and the carbon fiber composite material parameters and geometric parameters were of the following values:E1=400GPa, E1/E2=40, G12=G13=5/6, G23=6GPa, μ12=0.25, ρ=1600kg/m3, αfiber=45°−45°45°−45°45°, a=b=1m, and hU=hL=a/10. The stiffness of the middle elastic layer Kc=0 represents free connection and Kc=1014 represents rigidity connection, respectively. The dimensionless natural frequency is Ω=ωa2/π2ρu/E2uhu2. The results of Ye et al. [60], Karami et al. [57], and Liew et al. [58] were selected as comparative data. From the calculation results in Table 3, when the rigidity of the middle layer Kc=0, the phenomenon of “repeated frequency” appears obviously. It can be seen that the 1st order and the 2nd order have the same frequencies and the 3rd order and the 4th order are in the same situation. The reason for this phenomenon is that when the stiffness of the middle elastic layer is zero, the system degenerates into two independent and parallel composite laminates. As a result, the frequencies of the two plates will alternate. The “repetition frequency” of the LCDPS can be eliminated in such a situation that the stiffness coefficient of the Winkler layer Kc=1014. Then the rigid connected double-layer plate can be regarded as a single plate. At the same time, it is found that the natural frequency obtained at present is in accordance with the literature results of the study. The maximum error occurs when the boundary is SSSS-SSSS and Kc=1014. The difference between calculation results in this paper and the calculation results in [58] is only 1.325%, which appeared in the 7th-order frequency, and other differences between them are minor. As mentioned above, the validity of the analysis model in this paper is verified.As an effective numerical simulation method, the experimental verification is necessary. In this part, we will select a fully clamped rectangular plate as experimental objects to obtain its natural frequency. As depicted in Figure3, the method of single point vibration pick-up is used to excite the rectangular plate with a force hammer. The alloy steel plate is taken as the subject of research, which has the advantages of high strength, hardness, and toughness. The geometric parameters of this plate are a=0.6m, b=0.35m, and h=0.01m. The material properties are given as E=216GPa, μ=0.3, and ρ=7800kg/m3. For the isotropic plate with single-layered material, the setting of material parameters is given in Section 2.2. In Figure 4, the vibration modes and corresponding natural frequencies of the plate obtained from experimentation and the solution by the present method are compared. As can be seen from Figure 4, both modal and frequency results are in good agreement.Figure 3
The experimental setting for a rectangular plate. (a) Experimental plate. (b) Single vibration excitation.
(a)(b)Figure 4
Comparison between the present solutions and experimental results. (a) Present solutions. (b) Experimental results.
(a)(b)


## 3.2. Parametric Analysis of the LCDPS

In the last section, the convergence and efficiency of this method are obtained. In this part, a parametric approach to free vibration characteristics of LCDPS is carried out. Tables4 and 5 show the frequency parameters of LCDPS under different laying schemes and boundary conditions. The stiffness of the middle elastic layer between two plates is set as Kc=1014. In the calculation, the material parameters of the double plates are the same and consistent with Table 3. Geometric parameters are defined as a=b=1m and hU=hL=a/50, and the laying schemes of composite laying are as follows: αfiber=0°90°90°0°, αfiber=30°60°30°60°, and αfiber=45°−45°45°−45°. The model constructed in this paper can easily simulate different boundary conditions by setting different spring stiffness coefficients, such as CCCC-CCCC, FCFC-FCFC, and CCCC-FFFF. There is no need to change the model or rewrite the program. Such an approach shows great flexibility and adaptability. It can be seen from the calculation results in Table 4 that the vibration frequency of LCDPS changes with the stiffness of boundary springs. Under CCCC-CCCC boundary condition, the system has the maximum frequency parameter. When the boundary condition is FCFC-FCFC, the vibration frequency parameters of the system are obviously reduced. Different boundary conditions of top and bottom plates can also be set, such as CCCC-FFFF, and, with the increase of boundary degrees of freedom, the frequency decreases further. This phenomenon is due to the change of plate stiffness caused by different boundary constraints. The structural stiffness of the LCDPS in the clamped restrained boundary is the largest and decreases with boundary constraints relaxing gradually. In addition, we can observe from Table 4 that different composite layer design schemes will cause marked difference in the vibration characteristics of the structure. This phenomenon is mainly due to the fact that the layer scheme has a significant impact on the structural stiffness of the LCDPS.Table 4
Frequency parameters Ω for an LCDPS with different layer scheme and classical conditions.
LayerΩBCCCCC-CCCCCSCS-CSCSCFCF-CFCFSSSS-SSSSSSSS-FFFFCCCC-FFFFSSSS-CFCFSSSS-CSCSSSSS-CCCC[0°/90°/90°/0°]139.62539.62436.05223.3690.0090.00434.78135.26535.265256.46456.46236.39446.4480.0650.02751.28953.44453.444390.21190.21040.60766.9420.0790.09783.93585.55985.559494.95894.95857.93579.29122.63034.78185.33788.31188.3115104.017104.01692.35484.25343.94751.28994.23095.48495.4846126.654126.65393.150107.30966.65083.934116.436119.713119.714[30°/60°/30°/60°]131.38930.39318.35124.1840.0000.00527.61328.98329.554251.61350.05219.01040.1940.0300.19945.97547.96348.738371.37969.53230.38459.2130.3470.35664.38566.82967.805471.83870.07243.95761.20722.29428.42064.78967.26168.312594.37192.18950.15180.68138.48147.17485.97888.89890.0336105.633103.18754.91689.83856.57165.72596.58399.712100.827[45°/−45°/45°/−45°]135.38435.31017.85930.0450.0000.00332.75534.09334.129269.49569.27425.91056.5100.0420.24664.83265.97266.080369.49569.32546.97166.2210.0930.39864.83568.39168.4714108.122107.56748.71496.8560.39828.198101.761104.911105.1835116.321115.96859.67999.04228.10554.528109.399110.602110.7866117.213116.95279.567114.60454.50762.114110.069115.873115.992Table 5
Frequency parameters Ω for an LCDPS with different middle spring parameters, layer scheme, and classical conditions.
LayerΩKc=0Kc=109Kc=1014CCCC-CCCCCSDCSD-CSDCSDCFCF-CFCFCCCC-CCCCCSDCSD-CSDCSDCFCF-CFCFCCCC-CCCCCSDCSD-CSDCSD|CFCF-CFCF[45°/−45°]1114.26414.25211.93714.26414.25211.93714.26414.25211.937214.26414.25211.93719.56019.52713.12019.56019.52713.120319.56019.52713.12027.27427.21917.04727.27427.21917.047419.56019.52713.12034.24734.22022.83334.24734.22022.833527.27427.21917.04737.05536.97830.46337.05536.97830.463[45°/−45°]2120.74920.69417.00220.74920.69417.00220.74920.69417.002220.74920.69417.00228.26928.11518.62428.26928.11518.624328.26928.11518.62438.82438.58624.59738.82438.58624.597428.26928.11518.62445.01944.99432.86247.47547.37632.862538.82438.58624.59747.47547.37642.76251.76851.47442.762[45°/−45°]3121.60421.54017.68921.60421.54017.68921.60421.54017.689221.60421.54017.68929.37529.19719.35829.37529.19719.358329.37529.19719.35840.24939.97825.57440.24939.97825.574429.37529.19719.35845.42445.39434.13149.09248.98234.131540.24939.97825.57449.09248.98243.69953.53953.21044.189[45°/−45°]4121.88721.82017.91821.88721.82017.91821.88721.82017.918221.88721.82017.91829.73929.55319.60229.73829.55319.602329.73829.55319.60240.71740.43425.89740.71740.43425.897429.73829.55319.60245.56145.52934.54949.62149.50734.549540.71740.43425.89749.62149.50743.79354.11853.77744.658As mentioned above, with the solution of the vibration characteristic equation, the eigenvector corresponding to the natural frequency can be obtained as well. Then substitute the eigenvectors into the displacement admissible function of the system. From the first four modes simulated in Figure5 corresponding to Table 4, we can comprehend the natural vibration characteristics directly and vividly. It is shown in the figure that the mode shapes are different with varying boundary conditions, but the vibration modes of the top and bottom plate structures of the LCDPS are the same.Figure 5
The lowest four mode shapes for an LCDPS with several conditions.
(a)(b)(c)(d)(e)In the next section, Tables5 and 6 show the first five natural frequency parameters of the LCDPS under the classical and elastic boundary conditions with a changing number of layers of the layering scheme. The geometric parameters and physical properties are exactly the same as in Table 4. Set the stiffness of the middle elastic layer to Kc=0, Kc=109, and Kc=1014 respectively. From Table 5, we can see the trends from the calculation results that the vibration behaviors of its structure changed somewhat with the increase of the fiber layers. This effect will be more significant with lower layers number. As fiber layer counts increase, the properties of composite materials are close to those of isotropic materials. Table 6 presents the vibration characteristics of the LCDPS with elastic boundary conditions. By changing the stiffness of the Winkler layer, different frequency values can be obtained. While the stiffness coefficient of the middle layer has little effect on the lower-order frequency and the influence becomes larger with the increase of the order, we can infer that the stiffness of the middle layer maybe has an obvious effect on the higher-order frequency of the structure that will be referred to later.Table 6
Frequency parameters Ω for an LCDPS with different middle spring parameters, layer scheme, and elastic conditions.
LayerΩKc=0Kc=109Kc=1014E1E1E2E2-CCCCE1E1E3E3-CCCCE2E2E3E3-CCCCE1E1E2E2-CCCCE1E1E3E3-CCCCE2E2E3E3-CCCCE1E1E2E2-CCCCE1E1E3E3-CCCCE2E2E3E3-CCCC[45°/−45°]1111.95310.57011.87813.68413.22413.65614.22314.21814.178214.26414.19914.26418.65214.92818.60919.52414.92819.473316.41914.26416.28625.77615.30225.70927.24215.30327.184419.56014.92819.56030.08217.88030.04234.14919.50834.048522.52215.29822.35834.16320.36734.07237.02820.36836.960[45°/−45°]2115.05012.45114.90119.15214.92819.08620.63414.92820.505220.74914.92820.74925.86715.30425.74328.18015.30428.045321.09915.30320.79234.53618.05534.34638.75820.36938.612428.26916.76828.26936.83220.36936.71547.23920.61846.997529.26120.36828.88442.64123.96242.42349.83728.13249.817[45°/−45°]3115.38212.62015.22319.82214.92919.74921.47514.92921.331221.60414.92921.30026.70815.30426.56929.27615.30429.127321.63315.30421.60435.45118.62335.24440.17820.36940.019429.37517.01629.37537.54420.36937.40648.83621.45748.573530.11120.36929.69343.51324.61743.27549.84829.22249.839[45°/−45°]4115.08812.72814.96018.64215.77718.58719.85815.77719.755219.95015.77719.95025.28016.21025.18427.20716.21027.098320.96916.21020.72234.09317.70233.94437.49919.84637.381427.28017.08727.28037.09821.55337.02045.74321.55345.544528.85719.95028.55842.82723.66942.65850.13327.16950.003In addition, Figures6 and 7 show the influence of the fiber laying angle on the LCDPS under classical or elastic boundary conditions. Maintaining the original material parameters as in Table 6, the fiber angle laying scheme of the composite plates in the structure is defined as αfiberp=90°ααfiberp90°, in which ααfiberpp=UandL represent the laying angles of fiber in the middle surface of top and bottom plate, respectively. Changing the fiber angle of the lower plate from 10° to 180°, the fiber angle has a different influence on the vibration characteristics of LCDPS with different stiffness coefficients of intermediate Winkler elastic layers. It can be concluded that, at any boundary condition, in both the lower-order frequency range and the higher-order frequency range, the change of natural frequency parameter is symmetrical about the frequency of the system, whose fiber angle is 90°.Figure 6
Variation of dimensionless natural frequency versus fiber angle for an LCDPS with different conditions. (a) CCCC-CCCC, 1st mode; (b) SSSS-SSSS, 2nd mode; (c) CCCC-CCCC, 3rd mode; (d) SSSS-SSSS, 4th mode.
(a)(b)(c)(d)Figure 7
Variation of dimensionless natural frequency versus fiber angle for an LCDPS with different boundary conditions. (a) FCFC-FCFC, 5th mode; (b) FFCC-FFCC, 6th mode; (c) FCCC-FCCC, 7th mode; (d) FFFF-FFFF, 8th mode.
(a)(b)(c)(d)In the following examples, the influences of the geometric parameters of the LCDPS on the vibration characteristics are demonstrated. Keep material parameters of the top and bottom plates always be the same carbon fiber composite material, and the material parameters are assumed as follows:E1=740GPa, E1/E2=40, G12=G13=111GPa, G23=E1/8, μ12=0.25, ρ=1600kg/m3, and αfiber=30∘0∘−30∘0∘30∘. The geometrical parameters of LCDPS are defined as follows: a=b=1m, hU=a/50, and hL/hU=1∼20. Firstly, the effect of the thickness of composite plates on the vibration frequency of LCDPS is investigated. From Figure 8, we can see that in a certain thickness range the vibration frequency of LCDPS increases rapidly with the thickness. There is a unified trend with any boundary condition; when the thickness coefficient hL/hU increases to a certain value, the vibration frequency of the LCDPS will not change obviously. The structural vibration characteristics tend to be stable. Secondly, the effects of plate aspect ratio on the vibration behavior are studied. The material characteristics of the plates are the same and consistent with Figure 7. There are four layers for one laminated plate, and the geometric parameters are as follows: b/a=1∼10, hU=0.05, hL/hU=1∼20, and αfiber=0°90°0°90°. The result demonstrates in Figure 9 that the increase of plate length ratio causes the raise of vibration frequency of the system. It shows that the change of geometric parameters affects the vibration characteristic stiffness matrix of the LCDPS in a certain range. Taking advantage of this feature, the vibration characteristics can be adjusted by changing different geometry and material parameters. But, over the sensitive range, the effect of adjustment by geometrical parameters will be limited.Figure 8
Variation of dimensionless natural frequency versus thickness for an LCDPS with different conditions. (a) CCCC-CCCC; (b) SDSDSDSD-SDSDSDSD; (c) E1E1E1E1-CCCC.
(a)(b)(c)Figure 9
Variation of dimensionless natural frequency versus length-width ratio for an LCDPS with different conditions. (a) CCCC-CCCC; (b) SSSS-SSSS; (c) E1E1E1E1-CCCC.
(a)(b)(c)As we have mentioned earlier, Figure10 shows the regularities of an LCDPS on frequency parameter with the stiffness of the middle Winkler layer changing. The geometric and materials parameters are consistent with Figure 7. From the first 10 modal frequencies given, the change rule of the frequency parameters of the system can be observed intuitively. Firstly, with the increasing stiffness coefficient of the middle elastic layer, the frequency parameters from the first order to the tenth order have a common increasing trend. Secondly, the elastic stiffness coefficient of the middle Winkler elastic layer has a varied sensitive area, which is similar to the boundary condition. In this region, the vibration characteristics of LCDPS are very sensitive to the stiffness coefficient of the Winkler layer, and the influence of stiffness parameters on higher-order frequency is more apparent. This view improves the point mentioned earlier. Outside this region, the influence is weak. Therefore, the adjustment of the elastic parameters of the middle Winkler layer is another way to obtain the ideal vibration characteristic besides certain parameters such as boundary conditions, thickness coefficient, and aspect length-width ratio.Figure 10
Variation of frequency parameter Ω versus middle spring stiffness for an LCDPS with different boundary conditions. (a) CCCC-CCCC; (b) CCCC-SSSS; (c) CFCF-CFCF; (d) CCCC-E1E1E1E1; (e) CCCC-E2E2E2E2; (f) CCCC-E3E3E3E3.
(a)(b)(c)(d)(e)(f)


## 4. Conclusions

In this paper, on the basis of the first-order shear deformation theory (FSDT), an effective free vibration analysis model of the laminated composite double-plate system (LCDPS) with general boundary conditions is established. By adding auxiliary function to the traditional Fourier series, we successfully overcome the jump phenomenon of displacement function on the boundary. Based on the Rayleigh-Ritz principle, the vibration problem of LCDPS is transformed into a standard eigenvalue problem. The following conclusions are given:(a)
The consistent model of the laminated composite double-plate system is established. By introducing the uniformly distributed artificial spring technology, the proposed model can simply change the spring stiffness parameters of the middle elastic layer and on plates boundary, which can meet the complex engineering practical needs flexibly. The examples show that the method outperforms other existing ways.(b)
The results show that the method has good convergence and high efficiency. It can be concluded that the present method provides good flexibility in various applications compared with the FEM results and the solutions from references. Some new vibration results of the LCDBS are given, which can provide reference for future researchers in this field.(c)
There is a sensitive area of numerical influence of the stiffness of the middle elastic layer on the vibration features of the system. In the sensitive region, changes in the relevant parameters have a significant impact on the vibration characteristics of the system. The effects of different material parameters and geometric parameters on the vibration characteristics of the LCDPS are studied. These characteristics can be used to improve the design of some engineering structures, such as double-layer vibration isolation system and sound insulation structure.



---

*Source: 8811747-2021-03-03.xml*