Vibration Characteristics Analysis of Moderately Thick Laminated Composite Plates with Arbitrary Boundary Conditions

In this study, an improved Fourier series method is presented for the vibration modeling and analysis of moderately thick laminated composite plates with arbitrary boundary conditions, in which the vibration displacements are sought as the linear combination of a double Fourier cosine series and auxiliary series functions. The vibration model was established using the Hamilton energy principle. To study the vibration characteristics of laminated composite plates more comprehensively, firstly, the accuracy of the current results were validated via comparison with previous results and finite element method data. A parametric study was conducted on the effects of several key parameters, such as the h/b ratio, orientation and number of layers. In this section, both solutions are applicable to various combinations of boundary constraints, including classical boundary conditions and elastic-supported boundary conditions. Secondly, in order to identify the action position of vibration and the transmission of vibration energy, the response analysis of laminated plates was studied, and the power flow field for laminated plates was analyzed. Finally, a modal test was introduced to further verify the accuracy of the method in this paper.


Introduction
It is generally known that laminated composite plates, as basic structural components, are commonly applied in mechanical, aerospace, astronautic and civil engineering, among other fields. This is mainly because they are superior to conventional materials when high strength-to-weight and stiffness-to-weight ratios are required. Therefore, it is attractive for engineers and designers to research the material characteristics of the laminated composite plates.
According to Yu et al. [14], the vibration characteristics of anisotropic rectangular plates under the mixed boundary conditions of solid support and simple support are analyzed based on the superposition principle. Ge et al. [15] studied the dynamic response of symmetrically orthogonal composite laminates on elastic foundations under in-plane preloading and transverse impact loads based on the modal superposition method. In this paper, the equation of motion for displacement description is deduced by using the higher-order shear deformation theory. Huang et al. [16] studied the free vibration characteristics of orthotropic rectangular plates subjected to internal compression on a two-parameter elastic foundation by the method of separated variables. Kshirsagar et al. [17] analyzed the free vibration and buckling of rectangular plates under arbitrary classical boundary conditions based on the superposition principle of infinite series with infinite truncation. Liu et al. [18] presented the exact closed solution of a rectangular plate by using the method of separating variables. The boundary conditions of plate structure in this paper are simply supported on one side, and have arbitrary classical boundary conditions on the other. Chung et al. [19] studied the vibration characteristics of rectangular plates under elastic boundary conditions by the Rayleigh-Ritz method, where the Timoshenko beam function is used as the permissible deflection function. Liew et al. [20,21] presented the vibration characteristics of anisotropic plates and symmetrical composite laminates under mixed boundary conditions based on the subdomain and orthogonal polynomial methods. On this basis, they use the Reissner-Mindlin theory and p-Ritz method to analyze the vibration characteristics of composite laminates and calculate the vibration frequencies of composite laminates under different boundary conditions, length-width ratios and width-thickness ratios. Cheung et al. [22] studied the free vibration characteristics of symmetrically laminated plates under point-supported boundary conditions by using a new displacement admissible function and the Rayleigh-Ritz method. The allowable function is composed of a static beam function, which is different from the existing admissible function. The functions they set can satisfy both geometric boundary conditions and point-supported boundary conditions. Matsunaga et al. [23] analyzed the natural frequency and buckling stress of laminates by considering the influence of shear deformation thickness variation and moment of inertia, and expanded the vibration displacement function in power series. Mbakogu et al. [24] studied the bending problem of orthotropic rectangular plates under uniformly distributed loads by the Galerkin method. Dalaei et al. [25] solved the vibration problem of cantilever plates with anisotropic materials for the first time by using the extended Kantorovich method; they also obtained a closed-form high-precision solution in their paper. Bercin et al. [26] studied the bending and vibration of fully clamped plates by the Kantorovich method. Chen et al. [27][28][29][30] established a series of composite laminate models and studied the vibration, stability and large deformation of composite laminate plates. Luccioni et al. [31] presented the free vibration characteristics and stability of composite laminates by the finite element method, combining the classical laminate theory with the first-order shear deformation theory. Rao et al. [32] presented several finite element models to study the static, stability, impact and nonlinearity of laminated plates and shells. Shafiee et al. [33] analyzed the vibration characteristics of composite coupled plates by the finite element method. Huang et al. [34] studied the free vibration characteristics of rectangular plates with variable thickness based on the discrete method, and the characteristic equation of free vibration was obtained by the Green function. Liew et al. [35] presented the free vibration of symmetrical laminated plates by the moving least squares differential integral method. Song et al. [36][37][38] studied the free vibration of laminated plates by the local Radial Point Interpolation method; however, it is difficult to deal with the displacement boundary using a meshless method. Zhang et al. [39,40] presented an improved Fourier series method for the free vibration analysis of a moderately thick laminated composite rectangular plate with non-uniform boundary conditions. They also established a unified analysis model for vibration characteristics of composite laminated annular sector plates, circular sector plates, annular plates and circular plates with various elastic boundary conditions. Qin et al. [41] studied the free vibration analysis of composite laminated plates based on the Jacobi-Ritz method.
Romanelli et al. [42] presented the dynamic response of a simply supported rectangular composite plate under local pressure based on the series method. Shen et al. [43,44] discussed the dynamic response of a cantilever plate on a Pastemak elastic foundation under temperature and transverse dynamic loads by the Rayleigh-Ritz method, considering the influence of first-order shear deformation, foundation stiffness coefficient and temperature change. Khan et al. [45] presented the dynamic response of a composite cantilever plate under a uniform load by the variational method. Niyogi et al. [46] analyzed the free and forced vibration responses of composite laminates based on the nine-node element method, considering the influence of first-order transverse shear deformation and moment of inertia. Biswas et al. [47] studied the transient dynamic response behavior of multi-layered hybrid composite plates by a modified higher-order refined zigzag theory (HRZT). Zhang et al. [48] presented the free and forced vibration behaviors of thin, three-dimensionally coupled plate structures based on the dynamic stiffness method (DSM). Nefske et al. [49] presented a power flow finite element method, which is applied to the calculation of beam structures. In subsequent studies, power flow finite element method has been developed and gradually applied to the power flow analysis of coupled structures. Hambric et al. [50] calculated the structural sound intensity of a cantilever plate with stiffened members by the finite element method, and the bending, torsion and axial power flow were obtained. Li et al. [51] obtained the structural sound intensity vectors in three cases by finite element harmonic response analysis through calculating the surface admittance of a thin plate by the structural sound intensity method. Cieslik et al. [52] presented a method to introduce the bending moment and external force into the theoretical expression of structural vibration and sound intensity by using the complex mode theory. The vibration power flow of simply supported stiffened plates is calculated by the finite element method, and the energy flow distribution of stiffened plates is analyzed in this paper. Xing et al. [53] studied the power flow of structures under fluid-structure interaction based on the substructure method. The energy flow density vector was used to represent the energy flow transmission path between substructures, and the energy flow transmission was visually shown by graphics. Wang et al. [54] studied the structural sound intensity characteristics of composite laminates under dynamic concentrated forces. An example of structural sound intensity was analyzed by using finite element software; in this paper, the results show that orthotropic laminates have different characteristics from isotropic laminates, and that the structural sound intensity characteristics of orthotropic laminates are influenced by the boundary conditions, number of layers and stacking sequence. Zhang et al. [55][56][57][58][59] presented the free and forced vibration behaviors of the laminated plate-cavity coupling system by means of the improved Fourier series method.
There have been many studies on thin plates with arbitrary boundary conditions. However, these studies have limitations for moderately thick composite plates with special supported boundary conditions, for instance, an arbitrary elastic boundary. In addition, because the Mindlin theory considers the shear deformation along the thickness direction compared with the equivalent single-layer theory used in thin plates when studying the vibration of plates, the results obtained by the Mindlin theory are more accurate. On the other hand, most of the recent research merely discusses the vibration analysis, while the research on energy transfer is insufficient. Nevertheless, in many engineering applications, such as in the design of vehicles, it is necessary to comprehend the energy transfer path in order to reduce noise and vibration. Yan [60] points out that, compared with the traditional path analysis method, power flow analysis not only describes the speed of energy transfer or transformation, but also can be used to characterize the flow of vibration energy in the system. An et al. [61] studied the transmission path in hydraulic pipelines by using vibration power flow. Inoue et al. [62][63][64] studied the relationship between the structural noise response and the total power flow transmitted to the receiving structure, and they presented that there is a certain similarity between the response spectrum of structural noise and the total power spectrum transmitted to the receiving structure in a wide frequency range. Lee et al. [65,66] used the method of measuring vibration power flow to analyze the transmission paths of occupant interior structural noise caused by vehicle transmission system; they identified the main transmission paths according to the power flow, and modified the body structure with modal analysis to reduce the structural noise. Wang [67] studied the transmission characteristic of energy flows of micro-vibrations in spacecraft structures. Therefore, the study of power flow analysis of structural component is essential.
In this paper, a mechanics model is presented, based on the Mindlin theory, in order to research its vibration characteristics with arbitrary boundary conditions. Three kinds of restraining springs (translational, rotational and torsional), attached to each edge, are introduced to establish the general structure model of laminated composite plates. In addition, a modified Fourier series method is introduced to describe the vibration displacement function of the structure to study the influences of the structural style and the boundary conditions on the vibration characteristics of laminated composite plates. Moreover, the vibration model of the laminated composite plates with arbitrary boundary conditions is established, which is based on the Hamilton energy variation principle. This paper also focuses on the effect on vibration frequency while changing various key parameters (such as the ratio between thickness and width, the number of layers, as well as the laying angle between two layers). In the following section, several numerical examples of the free vibration of a laminated composite plate with arbitrary boundary conditions are presented, which can serve as references for future engineering. The current results were checked against previous results and achieved good agreement. Moreover, to analyze the vibration characteristics of laminated composite plates more comprehensively, the present work discusses the vibration response characteristic of the laminate, in terms of energy, through harmonic response and power flow analysis. In this section, we discuss the periodic response of a continuous periodic load in a structural system; in addition, we found the distribution character of energy transfer through changing the position of actions and boundary conditions. Finally, we introduce a test to further verify the accuracy of the presented method, in which the experimental data coincide with the theoretical value.

Theoretical Formulations
The mechanical properties of the composite laminate structure were related to the mechanical properties and thickness of the single-layer board, as well as the fiber laying direction and sequence and the number of layers. The laminated plate shown in Figure 1 consisted of five layers. The angle of the first layer was 0 • , and the other four layers were expressed as α, 90 • , −α and 0 • , respectively. The principal axis coordinate system of each layer in laminated composite plates was O 12 , and the global coordinate system was O xy . The ply-angle of the fiber can be expressed by θ, and the positive direction was identified when the x-axis rolled towards the 1 axis in an anticlockwise direction. The material of the plate structure in this chapter was a fiber-reinforced composite material, which was composed of matrix and fiber. The matrix was mainly used to support and protect the fiber, while the fiber was enhanced and mainly supported. In addition, the coupling vibration of bending and stretching of composite laminate plates is not considered in this chapter. introduced to describe the vibration displacement function of the structure to study the influences of the structural style and the boundary conditions on the vibration characteristics of laminated composite plates. Moreover, the vibration model of the laminated composite plates with arbitrary boundary conditions is established, which is based on the Hamilton energy variation principle. This paper also focuses on the effect on vibration frequency while changing various key parameters (such as the ratio between thickness and width, the number of layers, as well as the laying angle between two layers). In the following section, several numerical examples of the free vibration of a laminated composite plate with arbitrary boundary conditions are presented, which can serve as references for future engineering. The current results were checked against previous results and achieved good agreement. Moreover, to analyze the vibration characteristics of laminated composite plates more comprehensively, the present work discusses the vibration response characteristic of the laminate, in terms of energy, through harmonic response and power flow analysis. In this section, we discuss the periodic response of a continuous periodic load in a structural system; in addition, we found the distribution character of energy transfer through changing the position of actions and boundary conditions. Finally, we introduce a test to further verify the accuracy of the presented method, in which the experimental data coincide with the theoretical value.

Theoretical Formulations
The mechanical properties of the composite laminate structure were related to the mechanical properties and thickness of the single-layer board, as well as the fiber laying direction and sequence and the number of layers. The laminated plate shown in Figure 1 consisted of five layers. The angle of the first layer was 0°, and the other four layers were expressed as α, 90°, -α and 0°, respectively. The principal axis coordinate system of each layer in laminated composite plates was O12, and the global coordinate system was Oxy. The ply-angle of the fiber can be expressed by θ, and the positive direction was identified when the x-axis rolled towards the 1 axis in an anticlockwise direction. The material of the plate structure in this chapter was a fiber-reinforced composite material, which was composed of matrix and fiber. The matrix was mainly used to support and protect the fiber, while the fiber was enhanced and mainly supported. In addition, the coupling vibration of bending and stretching of composite laminate plates is not considered in this chapter. For the plate model based on the assumptions of the first-order shear deformation theory, the deformation of the plate was continuous. Thus, on the basis of not considering the in-plane vibration of plates, the displacement field can be expressed as ( , , ) ( , ) x u x y z z x y y v x y z z x y For the plate model based on the assumptions of the first-order shear deformation theory, the deformation of the plate was continuous. Thus, on the basis of not considering the in-plane vibration of plates, the displacement field can be expressed as v(x, y, z) = zψ y (x, y), w(x, y, z) = w(x, y) where u, v and w are the displacement in the x-, yand zdirections, respectively, and ψ x and ψ y are the rotation in the xand ydirections, respectively. Therefore, the stress-strain relationship in coordinates, based on the small deflection theory, can be written as follows: The corresponding stresses are obtained in terms of the generalized Hooke's law: where the normal stresses are σ x and σ y in the x and y directions, respectively, and the shear stresses are τ yz , τ xz and τ xy in the x, y and z coordinate system, respectively. The lamina stiffness coefficients are Q k ij (i, j = 1, 2, 4, 5, 6), which can be expressed as follows: where m = cos θ k ,n = sin θ k , in which the included angle is simply represented as θ between the principal direction of the layer and the x-axis. The lamina elastic coefficients Q k ij in Equation (6) can be obtained from the material properties of the kth orthotropic lamina layer: where E 1 and E 2 denote the longitudinal modulus and the transverse modulus, µ 12 is the major Poisson's ratio and the other Poisson's ratio µ 21 can be obtained by µ 12 E 2 = µ 21 E 1 . G 12 , G 13 and G 23 are the shear moduli. Through taking the unit body of unit length from the plate in the x/y direction, we can synthesize the stress component on the cross section of the unit body into an internal force per unit width. Therefore, the bending moment M x and M y , the torque M xy and M yx and the shear force Q x and Q y can be expressed as where κ is the shear correction factor, and σ k x , σ k y , τ k xy , τ k xz and τ k yz are the corresponding stresses of the kth layer. Consider a structural model of laminated composite plates that is elastically restrained along its four edges, as shown in Figure 2. The laminated composite plates may also be loaded with three kinds of restraining springs (translational, rotational and torsional) at arbitrary locations, and the boundary conditions can be easily obtained by changing the stiffness of the springs.
According to Hamilton's principle, the vibration equation can be easily established. For a holonomic mechanical system with n generalized degrees of freedom, the generalized coordinates are q s (s = 1, 2, . . . n), the Lagrange function is L(q s , . q s , t) = T − U − W, where T, U and W are the kinetic energy, potential energy and the external work, respectively. The Hamilton action is defined as follows: where κ is the shear correction factor, and Consider a structural model of laminated composite plates that is elastically restrained along its four edges, as shown in Figure 2. The laminated composite plates may also be loaded with three kinds of restraining springs (translational, rotational and torsional) at arbitrary locations, and the boundary conditions can be easily obtained by changing the stiffness of the springs.
According to Hamilton's principle, the vibration equation can be easily established. For a holonomic mechanical system with n generalized degrees of freedom, the generalized coordinates are qs (s = 1,2,…n), the Lagrange function is ( , , ) For the laminated composite plate in Figure 2, the Hamilton equation of the plate can be expressed as follows: where V, T and Wext, denote the total potential, the total kinetic energy and the external potential energy, respectively. The total potential includes the bending potential energy of the plate and the elastic potential energy of springs. The total potential energy can be written as follows: plate spring where Vspring is the elastic potential energy of springs, Vplate is the bending potential energy of the plate and these terms can be explicitly expressed as follows.
The elastic potential energy of springs on the four edges can be written as follows: For the laminated composite plate in Figure 2, the Hamilton equation of the plate can be expressed as follows: where V, T and W ext , denote the total potential, the total kinetic energy and the external potential energy, respectively. The total potential includes the bending potential energy of the plate and the elastic potential energy of springs. The total potential energy can be written as follows: where V spring is the elastic potential energy of springs, V plate is the bending potential energy of the plate and these terms can be explicitly expressed as follows.
The elastic potential energy of springs on the four edges can be written as follows: Materials 2019, 12, 2829 The stress relationship can be evaluated according to the stress-strain relationship. On the basis of the elastic mechanic theory, the relationship between elastic potential energy and stress-strain is expressed as follows: In this equation, Θ is given by and ε 1 is given by The stiffness can be defined as Therefore, Equation (5) can be expressed as follows: The total kinetic energy of the laminated composite plate is given by The external potential energy of the laminated composite plate can be expressed as where w is the flexural displacement; a and b are the plate dimension in x direction and the plate dimension in y direction, respectively; ρ and h are the mass density and the thickness of the plate, respectively, and f (x, y) = F δ (x − x 0 ) (y − y 0 ), where δ is the delta function, F is the amplitude of the stresses and x 0 and y 0 are the position coordinates of force. k x0 and k xa (k y0 and k yb ) are the transverse spring constants, K x0 and K xa (K y0 and K yb ) are the rotational spring constants and K yx0 and K yxa (K xy0 and K xyb ) are the torsional spring constants at x = 0 and x = a (y = 0 and y = b), respectively.
The displacement function is expanded as a single Fourier series plus an auxiliary polynomial function. The auxiliary function is used to overcome the discontinuities of the resilient boundary, the transverse displacement function w(x, y).
where p(x,y) can be expressed as where the auxiliary polynomial based on trigonometric function can be written as follows: where λ m = mπ/a, λ n = nπ/b; a is length and b is width and A mn , d 1 1m ,d 1 2m , f 1 1n and f 1 2n are the expansion coefficients.
Therefore, function w(x, y) can be rewritten as follows: A mn cos(λ m x) cos(λ n y) Similarly, the rotational displacement function Ψ x (x, y) and Ψ y (x, y) can be easily obtained as follows: where B mn , d 2 1m ,d 2 2m , f 2 1n , f 2 2n ,C mn , d 3 1m ,d 3 2m , f 3 1n and f 3 2n are the expansion coefficients. Substitution of Equations (17)- (19) into Equation (11) leads to a series of equations, and these equations can be rewritten in matrix form as where K is the stiffness matric, M is the mass matric, A is the unknown coefficient vector and F is the force vector. The detailed expressions of K and M are shown in the Appendix A.
In order to determine the modal characteristics of the laminated composite plate, one needs to solve the characteristic equation by setting the external force vector F in Equation (20) to zero. Obviously, the natural frequencies and eigenvectors of the laminated composite plate can now be easily obtained by solving a standard matrix Eigen problem. Since the components of each eigenvector are actually the expansion coefficients of the Fourier series, the corresponding mode shape can be directly determined by substituting the eigenvector in Equation (17) into Equation (19).
On the other hand, when the external force vector F is not zero, it can be used to study the response analysis of laminates. The response obtained is a harmonic response when the external excitation force is a simple harmonic force.
The lateral displacement field and the corner of the bending vibration of the plate structure at the excitation frequency by substituting the coefficient vector A in Equation (17) into Equation (19).
To solve the structural vibration and noise radiation problems, vibration power flow leads to a more comprehensive evaluation of the vibration response characteristics of the structure. It can intuitively identify the position of the vibration and the transmission path of the vibration energy, providing a better evaluation index and basis for the vibration control of the structure.
The power flow intensity I(x, y) at any point in the Mindlin plate structure can be expressed as follows: where I x (x,y) and I y (x,y) are the components of I(x, y) along the x and y axes, respectively.

Modal Characteristics of Laminated Composite Plates
Consider a laminated composite plate with arbitrary elastic boundary support, with the dimensions of a × b, which is illustrated in Figure 2. If three kinds of restraining springs (translational, rotational and torsional) are assumed to be distributed uniformly along each edge, then all the classical boundary conditions, as well as their combinations, can be easily obtained by simply setting the spring coefficients to zero or infinity. For the sake of convenience, C denotes the clamped cases, E denotes the elastic bearing cases, F denotes the free-boundary condition and S represents the simply supported cases. Taking the edge of x = 0 as an example, the spring stiffness is shown in Table 1 when the boundary condition is F, S, C or E, respectively. (Note that the elastic boundary can be chosen arbitrarily, and the exact value of E can be transformable in different examples.) First, the relationship between the vibration solution of a rectangular plate based on the energy principle and the analytical solution based on the governing equation and the boundary condition equation was studied. The laminate was a single-layer board and the material was an isotropic material, the same as in [18].
In Table 2, the first five non-dimensional frequency parameters are shown (Ω = (ωb 2 /π 2 )(ρh/D) 1/2 ). It can be seen from the table that the deviations of the results obtained by the present method and the analytical method were all within 0.022. The results obtained were almost identical, which verified that the energy method resulted in an exact solution when the permissible displacement function satisfied both the displacement boundary condition and the force boundary condition. In the following calculations, the laminated composite was an orthogonal symmetrical laminate, which was composed of three layers of unit plates of equal thickness-that is, the layup direction of the unit plate was 0 • /90 • /0 • . The physical parameters are specified as E x /E y = 40, G 12 = 3/5E y , G 23 = 1/2E y , To verify the convergence of the solution, Table 3 compares the first seven non-dimensional modal frequency parameters of the laminated composite plate under different cutoff values. For illustrative purposes, one assumes that a SSFF laminated composite plate system shown in Table 3 is defined by the dataset: with E representing the larger value between E x and E y . At the same time, a comparison of the solution in [21] is also presented. It is seen that the maximum deviation of the first seven non-dimensional modal frequency parameters obtained in both cases (M = N = 4 and M = N = 16) was 1.4%; on the other hand, when the cutoff value was more than 10, the results were almost invariant. Thus, it was evident that just a few terms can lead to excellent prediction and the current solution shows remarkable convergence. To further verify the accuracy of the method from a modal perspective, the first five non-dimensional frequency parameters are shown for the plates of various thickness-width ratios and boundary conditions in Tables 4 and 5. Meanwhile, the results calculated by using the separate variable method in [21], the Rayleigh-Ritz method in [18] and the results calculated by FEA are also presented in Tables 4 and 5. It is seen that the present method has led to excellent agreement with the classical solution.    To further verify the accuracy of the method from a modal perspective, the first five nondimensional frequency parameters are shown for the plates of various thickness-width ratios and boundary conditions in Tables 4 and 5. Meanwhile, the results calculated by using the separate variable method in [21], the Rayleigh-Ritz method in [18] and the results calculated by FEA are also presented in Tables 4 and 5. It is seen that the present method has led to excellent agreement with the classical solution. Figure 3 shows the mode shapes of the CCCC and FFFF moderately thick laminated composite plate. The mode shapes obtained by the present method coincide with those obtained by the finite element method. The material properties were as follows: a/b = 1, h/b = 0.2; the layup direction of the unit plate was 0°/90°/0°.  In order to show the consistency between the results of the present method and finite element method clearly, Figure 4 shows the frequency parameter curves of the SFCF laminated composite plates. The layup direction of the unit plate was 45°/−45°/45°; a/b = 1; h/b = 0.1. Figure 5 shows the frequency parameter curves of the EFCS laminated composite plates. The spring stiffness coefficient on the elastic edge were as follows: kx0 = D0, Kx0 = 0, Kyx0 = 0.  In order to show the consistency between the results of the present method and finite element method clearly, Figure 4 shows the frequency parameter curves of the SFCF laminated composite plates. The layup direction of the unit plate was 45 • /−45 • /45 • ; a/b = 1; h/b = 0.1. Figure 5 shows the frequency parameter curves of the EFCS laminated composite plates. The spring stiffness coefficient on the elastic edge were as follows: k x0 = D 0 , K x0 = 0, K yx0 = 0.  In order to show the consistency between the results of the present method and finite element method clearly, Figure 4 shows the frequency parameter curves of the SFCF laminated composite plates. The layup direction of the unit plate was 45°/−45°/45°; a/b = 1; h/b = 0.1. Figure 5 shows the frequency parameter curves of the EFCS laminated composite plates. The spring stiffness coefficient on the elastic edge were as follows: kx0 = D0, Kx0 = 0, Kyx0 = 0.  In Figures 4 and 5, the frequency curves obtained by this method coincided with those obtained by the finite element method at medium and low frequencies. At high frequencies, the results obtained by the finite element method were larger than presented results. This was because the finite element method needs finer meshes in order to get more accurate results in the high frequency band, with the workload increasing sharply. In Figures 4 and 5, the frequency curves obtained by this method coincided with those obtained by the finite element method at medium and low frequencies. At high frequencies, the results obtained by the finite element method were larger than presented results. This was because the finite element method needs finer meshes in order to get more accurate results in the high frequency band, with the workload increasing sharply.
To research the influence of different ratios between thickness and width, the first six non-dimensional frequency parameters for the SSSS and CCCC laminated composite plates are shown in Figures 6 and 7, where the number of piles in the figures were three layers and five layers, respectively. As can be seen in the figures, compared with low-order frequency parameters, high-order frequency parameters showed a remarkable decrease with an increase in the ratio between thickness and width; moreover, the change was more significant when h/b was less than 0.10. In Figures 4 and 5, the frequency curves obtained by this method coincided with those obtained by the finite element method at medium and low frequencies. At high frequencies, the results obtained by the finite element method were larger than presented results. This was because the finite element method needs finer meshes in order to get more accurate results in the high frequency band, with the workload increasing sharply.
To research the influence of different ratios between thickness and width, the first six nondimensional frequency parameters for the SSSS and CCCC laminated composite plates are shown in Figures 6 and 7, where the number of piles in the figures were three layers and five layers, respectively. As can be seen in the figures, compared with low-order frequency parameters, highorder frequency parameters showed a remarkable decrease with an increase in the ratio between thickness and width; moreover, the change was more significant when h/b was less than 0.10.   In Figures 4 and 5, the frequency curves obtained by this method coincided with those obtained by the finite element method at medium and low frequencies. At high frequencies, the results obtained by the finite element method were larger than presented results. This was because the finite element method needs finer meshes in order to get more accurate results in the high frequency band, with the workload increasing sharply.
To research the influence of different ratios between thickness and width, the first six nondimensional frequency parameters for the SSSS and CCCC laminated composite plates are shown in Figures 6 and 7, where the number of piles in the figures were three layers and five layers, respectively. As can be seen in the figures, compared with low-order frequency parameters, highorder frequency parameters showed a remarkable decrease with an increase in the ratio between thickness and width; moreover, the change was more significant when h/b was less than 0.10.  To further verify the effect of the number of layers on the frequency, Figures 8 and 9 show the first six non-dimensional frequency parameters for the SSSC and CCSS laminated composite plates of various layers, where h/b = 0.05 and the layup directions of the unit plates in Figures 8 and 9 were 0 • /45 • /0 • and 0 • /90 • /0 • , respectively (when the number of layers was three), and the layup directions in Figures 8 and 9, respectively were 0 • /45 • /0 • /45 • /0 • and 0 • /90 • /0 • /90 • /0 • (when the number of layers is five), etc. In the following figures, the data were collected when the number of layers was odd (3, 5, 7, 9, etc.).
It can be seen in Figures 8 and 9 that the frequency parameters increased with the increase in the number of layers, so long as the number of layers was less than 13, in which case, it stayed constant as the number of layers changed. It should be mentioned that the effect was more noticeable for higher frequencies.
It can be seen in Figures 8 and 9 that the frequency parameters increased with the increase in the number of layers, so long as the number of layers was less than 13, in which case, it stayed constant as the number of layers changed. It should be mentioned that the effect was more noticeable for higher frequencies.  In order to study the influence of the number of layers laid on the vibration frequency, Figure  10 shows the variation of the dimensionless frequency of the first order with the number of layers and the aspect ratio of the normal symmetric angled plywood structure under different boundary conditions. Among them, the laying angle of the laminate was 45°, and h/b = 0.1.
In Figure 10c  is five), etc. In the following figures, the data were collected when the number of layers was odd (3, 5, 7, 9, etc.).
It can be seen in Figures 8 and 9 that the frequency parameters increased with the increase in the number of layers, so long as the number of layers was less than 13, in which case, it stayed constant as the number of layers changed. It should be mentioned that the effect was more noticeable for higher frequencies. In order to study the influence of the number of layers laid on the vibration frequency, Figure  10 shows the variation of the dimensionless frequency of the first order with the number of layers and the aspect ratio of the normal symmetric angled plywood structure under different boundary conditions. Among them, the laying angle of the laminate was 45°, and h/b = 0.1.
In Figure 10c  In order to study the influence of the number of layers laid on the vibration frequency, Figure 10 shows the variation of the dimensionless frequency of the first order with the number of layers and the aspect ratio of the normal symmetric angled plywood structure under different boundary conditions. Among them, the laying angle of the laminate was 45 • , and h/b = 0.1.

(c) EEEE_1
(d) EEEE_2 The laying angle was another important design parameter of the laminate. Figure 11 shows the variation of the first-order vibration frequency with the laying angle of the composite laminate under different layers.

Vibration Response Analysis of Laminated Composite Plates
In this section, the following example is focused on the vibration harmonic response of the plate. In order to avoid the numerical instability under the excitation of the plate structure resonance frequency in the numerical calculation, the structural damping factor η is introduced in the analysis, and the elastic modulus of the plate structure is correspondingly the complex elastic modulus, E ' = E(1 + jη). To avoid the instability of the response at the resonant frequency, the damping factor η (loss factor value) of the plate structure herein is η = 0.01.

Vibration Response Analysis of Laminated Composite Plates
In this section, the following example is focused on the vibration harmonic response of the plate. In order to avoid the numerical instability under the excitation of the plate structure resonance frequency in the numerical calculation, the structural damping factor η is introduced in the analysis, and the elastic modulus of the plate structure is correspondingly the complex elastic modulus, E ' = E(1 + jη).
To avoid the instability of the response at the resonant frequency, the damping factor η (loss factor value) of the plate structure herein is η = 0.01.
The material in this section is an orthotropic plate, and the material properties of the plate were as follows: a = 1 m, b = 1 m, h = 0.1 m, E 1 = 128.8 GPa, E 2 = 8.3 GPa, G 12 = 4.1 GPa, G 23 = G 12 , G 13 = G 12 , v 2 = 0.355. Figure 12 shows the displacement harmonic response curve at different points when the point force acted on the center of the cantilever plate. This part of the response curve was based on the frequency domain analysis, as the comparison map obtained by the frequency domain analysis was more intuitive, and the full text always analyzed the frequency domain part of the research, which was also to maintain the consistency of the whole article.  Figure 12 shows the displacement harmonic response curve at different points when the point force acted on the center of the cantilever plate. This part of the response curve was based on the frequency domain analysis, as the comparison map obtained by the frequency domain analysis was more intuitive, and the full text always analyzed the frequency domain part of the research, which was also to maintain the consistency of the whole article.
The solid support edge was the edge on x = 0, and the force amplitude was F = 1 N. Figure 12a is a curve where the response point was at the center of the plate (0.5 m, 0.5 m), and Figure 12b is a curve at which the response point was at (0.8 m, 0.8 m). It can also be seen from the figure that the frequency corresponding to the peak value of the response was the modal frequency of the plate structure, and the response curves obtained by the finite element method were also given in the two figures. The mesh size was 0.02 m × 0.02 m. The solid support edge was the edge on x = 0, and the force amplitude was F = 1 N. Figure 12a is a curve where the response point was at the center of the plate (0.5 m, 0.5 m), and Figure 12b is a curve at which the response point was at (0.8 m, 0.8 m). It can also be seen from the figure that the frequency corresponding to the peak value of the response was the modal frequency of the plate structure, and the response curves obtained by the finite element method were also given in the two figures. The mesh size was 0.02 m × 0.02 m.
In order to study the effect of the excitation on the response of the rectangular plates with different positions, Figure 13 shows the displacement response and velocity response of the composite laminate under the action of point force. The boundary condition of the laminate was CFFC. The amplitude of the point force was 1 N, the action direction was along the z-axis, the action position was (0.5 m, 0.5 m) and the action time was 0.5 s. The displacement response in the figure was the same as the position of the velocity response, with both points being (0.5 m, 0.5 m). The calculation curve obtained by the finite element method is also shown in the figure. It can be seen that the results obtained by the two methods agreed well, thus verifying the feasibility of the method for solving transient problems. The material in this section is an orthotropic plate, and the material properties of the plate were as follows: a = 1 m, b = 1 m, h = 0.1 m, E1 = 128.8 GPa, E2 = 8.3 GPa, G12 = 4.1 GPa, G23 = G12, G13 = G12, v2 = 0.355. Figure 12 shows the displacement harmonic response curve at different points when the point force acted on the center of the cantilever plate. This part of the response curve was based on the frequency domain analysis, as the comparison map obtained by the frequency domain analysis was more intuitive, and the full text always analyzed the frequency domain part of the research, which was also to maintain the consistency of the whole article.
The solid support edge was the edge on x = 0, and the force amplitude was F = 1 N. Figure 12a is a curve where the response point was at the center of the plate (0.5 m, 0.5 m), and Figure 12b is a curve at which the response point was at (0.8 m, 0.8 m). It can also be seen from the figure that the frequency corresponding to the peak value of the response was the modal frequency of the plate structure, and the response curves obtained by the finite element method were also given in the two figures. The mesh size was 0.02 m × 0.02 m. In order to study the effect of the excitation on the response of the rectangular plates with different positions, Figure 13 shows the displacement response and velocity response of the composite laminate under the action of point force. The boundary condition of the laminate was CFFC. The amplitude of the point force was 1 N, the action direction was along the z-axis, the action position was (0.5 m, 0.5 m) and the action time was 0.5 s. The displacement response in the figure was the same as the position of the velocity response, with both points being (0.5 m, 0.5 m). The calculation curve obtained by the finite element method is also shown in the figure. It can be seen

Vibration Power Flow Analysis of Laminated Composite Plates
In some engineering applications, studying the energy transfer path in structural systems is helpful for solving the problems of structural vibration and noise radiation. In recent years, as a new energy transfer analysis method, power flow studies the transmission and control of vibration in structures from the perspective of energy. Because the power flow takes into account the inherent information of both force and motion of the structure, it can evaluate the vibration response characteristics of the structure more comprehensively than the previous single description of force or displacement. In addition, energy dissipation and collection can be characterized through the study of power flow fields, having practical significance for engineering.
In order to study the influence factors and characteristics of energy transfer, the dynamical characteristics of the structures are described as the power flow field. The power flow intensity at any point can be obtained by extracting the force and displacement from Equations (37)- (39), as shown in that the results obtained by the two methods agreed well, thus verifying the feasibility of the method for solving transient problems.

Vibration Power Flow Analysis of Laminated Composite Plates
In some engineering applications, studying the energy transfer path in structural systems is helpful for solving the problems of structural vibration and noise radiation. In recent years, as a new energy transfer analysis method, power flow studies the transmission and control of vibration in structures from the perspective of energy. Because the power flow takes into account the inherent information of both force and motion of the structure, it can evaluate the vibration response characteristics of the structure more comprehensively than the previous single description of force or displacement. In addition, energy dissipation and collection can be characterized through the study of power flow fields, having practical significance for engineering.
In order to study the influence factors and characteristics of energy transfer, the dynamical characteristics of the structures are described as the power flow field. The power flow intensity at any point can be obtained by extracting the force and displacement from Equations (37) The material in this part was normal symmetric orthogonally laid laminate, and the material  The excitation frequency in Figure 16 was not the first four natural frequencies of the CSSE rectangular plate, but the excitation frequency from Figure 14. It can be seen from the figure that the distribution of the power flow field had strong dependence on the position of the excitation force and the boundary condition of the rectangular plate, but there was no fixed distribution between them.
In order to study the energy transfer characteristics of plate structures subjected to multiple point forces, Figure 17 shows the distribution of vibration power flow in the whole rectangular plate structure under two-point excitation. The boundary conditions of the plate structure were CSSS. The excitation frequencies of the two point forces were the same, and the direction of action was along the z-axis, with the amplitude of 1 N. The positions of action for the two points were (0.2 m, 0.2 m) and (0.7 m, 0.7 m), respectively. Figure 17a-d shows the power flow vector graphs at first to fourth natural frequencies. It can be seen from the figure that the vibration energy was transmitted from the excitation to the surrounding plate, but the excitation position was not always the source of vibration energy output. Under some excitation frequencies, the excitation position may also have the flow of vibration energy, which further verifies that the energy transfer was strongly dependent on the excitation frequency.  The excitation frequency in Figure 16 was not the first four natural frequencies of the CSSE rectangular plate, but the excitation frequency from Figure 14. It can be seen from the figure that the distribution of the power flow field had strong dependence on the position of the excitation force and the boundary condition of the rectangular plate, but there was no fixed distribution between them.
In order to study the energy transfer characteristics of plate structures subjected to multiple point forces, Figure 17 shows the distribution of vibration power flow in the whole rectangular plate structure under two-point excitation. The boundary conditions of the plate structure were CSSS. The excitation frequencies of the two point forces were the same, and the direction of action was along the z-axis, with the amplitude of 1 N. The positions of action for the two points were (0.2 m, 0.2 m) and (0.7 m, 0.7 m), respectively. Figure 17a-d shows the power flow vector graphs at first to fourth natural frequencies. It can be seen from the figure that the vibration energy was transmitted from the excitation to the surrounding plate, but the excitation position was not always the source of vibration energy output. Under some excitation frequencies, the excitation position may also have the flow of vibration energy, which further verifies that the energy transfer was strongly dependent on the excitation frequency.
the z-axis, with the amplitude of 1 N. The positions of action for the two points were (0.2 m, 0.2 m) and (0.7 m, 0.7 m), respectively. Figure 17a-d shows the power flow vector graphs at first to fourth natural frequencies. It can be seen from the figure that the vibration energy was transmitted from the excitation to the surrounding plate, but the excitation position was not always the source of vibration energy output. Under some excitation frequencies, the excitation position may also have the flow of vibration energy, which further verifies that the energy transfer was strongly dependent on the excitation frequency. In Figures 14-17, it can be seen that the energy transfer in the structure had an important relationship with the frequency, position, number of excitation forces and boundary conditions of the plate structure. As the excitation frequency increased, the flow path of the vibration energy in the structure became more and more complicated. When the boundary condition of one side became soft, the vibration energy rapidly flowed to this side. The vibration energy did not always transfer energy from the excitation position to the boundary according to the shortest path principle. There was also a certain form of energy circulation inside the plate structure. Through the description of the power flow, the transmission of the vibration energy could be clearly seen, thereby providing a basis for the effective control of the vibration of the structure.

Modal Test of Laminated Plates
Aiming at strengthening our understanding of the vibration behaviors of the laminated composite plate, we verified the accuracy of the present method by experimental means. Figure 18 shows the experiment to study the model characteristics of the plate. A hammer was used for excitation, and the experiment was carried out by single-point pock-up and point-by-point excitation. In Figures 14-17, it can be seen that the energy transfer in the structure had an important relationship with the frequency, position, number of excitation forces and boundary conditions of the plate structure. As the excitation frequency increased, the flow path of the vibration energy in the structure became more and more complicated. When the boundary condition of one side became soft, the vibration energy rapidly flowed to this side. The vibration energy did not always transfer energy from the excitation position to the boundary according to the shortest path principle. There was also a certain form of energy circulation inside the plate structure. Through the description of the power flow, the transmission of the vibration energy could be clearly seen, thereby providing a basis for the effective control of the vibration of the structure.

Modal Test of Laminated Plates
Aiming at strengthening our understanding of the vibration behaviors of the laminated composite plate, we verified the accuracy of the present method by experimental means. Figure 18 shows the experiment to study the model characteristics of the plate. A hammer was used for excitation, and the experiment was carried out by single-point pock-up and point-by-point excitation. In this experiment, the boundary condition of the laminated composite plate was CFFF, and the clamped boundary condition was simulated by the pressure plate and the Φ6 bolt fixed on the foundation. The laminated plate was divided into a uniform mesh, the number of grids was 9 × 9 and the excitation position of the hammer was at each grid point. According to the experimental requirements, the sampling frequency was selected to be 5000 Hz and the analysis frequency was 1950 Hz, where the sampling mode was transient and the triggering mode was signal triggering. The Donghua DH5922 data acquisition analyzer and dynamic signal acquisition system were used to analyze the excitation signal of the hammer and the response signal of the acceleration sensor. Finally, a modal analysis was performed to obtain the natural frequency in Table 6 and the mode shape of the rectangular plate in Figure 19. Table 6. Natural frequency of the laminate plate. In this experiment, the boundary condition of the laminated composite plate was CFFF, and the clamped boundary condition was simulated by the pressure plate and the Φ6 bolt fixed on the foundation. The laminated plate was divided into a uniform mesh, the number of grids was 9 × 9 and the excitation position of the hammer was at each grid point. According to the experimental requirements, the sampling frequency was selected to be 5000 Hz and the analysis frequency was 1950 Hz, where the sampling mode was transient and the triggering mode was signal triggering. The Donghua DH5922 data acquisition analyzer and dynamic signal acquisition system were used to analyze the excitation signal of the hammer and the response signal of the acceleration sensor. Finally, a modal analysis was performed to obtain the natural frequency in Table 6 and the mode shape of the rectangular plate in Figure 19.  Table 6 shows the first eight natural frequency values obtained from the experimental composite laminate structure. The theoretical results are also given in the table. By comparing the experimental values and the theoretical values, we see that the results were in good agreement, and the maximum deviation was less than 4%, which verified the accuracy of the proposed modeling method. The results of the two methods showed a certain deviation, mainly because during the experiment the clamped boundary conditions were simulated by using multiple bolt joints, which had certain inconsistencies with the complete solidification in the theoretical calculation process. This may cause the stiffness to be less than the clamped end stiffness, resulting in the frequency of the experiment being higher than the theoretical value. Furthermore, a certain error was also generated, as the acceleration sensor introduced additional mass when it was connected to the plate structure. Figure 19 shows the first four natural modes of the composite laminate obtained by the experimental method and theoretical method.

Summary
In this paper, on the basis of the improved Fourier series method, vibration modeling of a laminated composite plate was proposed to study vibration characteristics of laminates. According to the Mindlin plate theory, the mechanical models of the laminates were established and a number of functions were deduced, such as the potential energy function, as well as the function of kinetic energy. Simultaneously, according to the Hamilton principle, the vibration equations of the plates with arbitrary boundary conditions were derived by expressing the displacement as a superposition of a Fourier cosine series and four auxiliary polynomials.
Through the numerical analysis, the effectiveness of the present method was verified by a comparison between the FEA results and the results of the method proposed in relevant references and mode shapes. The natural frequency was analyzed in terms of the different boundary conditions, and the present method had the better accuracy. The effects of boundary conditions, the ratio between thickness and width, the number of layers of laminates and the laying angle on the vibration

Summary
In this paper, on the basis of the improved Fourier series method, vibration modeling of a laminated composite plate was proposed to study vibration characteristics of laminates. According to the Mindlin plate theory, the mechanical models of the laminates were established and a number of functions were deduced, such as the potential energy function, as well as the function of kinetic energy. Simultaneously, according to the Hamilton principle, the vibration equations of the plates with arbitrary boundary conditions were derived by expressing the displacement as a superposition of a Fourier cosine series and four auxiliary polynomials.
Through the numerical analysis, the effectiveness of the present method was verified by a comparison between the FEA results and the results of the method proposed in relevant references and mode shapes. The natural frequency was analyzed in terms of the different boundary conditions, and the present method had the better accuracy. The effects of boundary conditions, the ratio between thickness and width, the number of layers of laminates and the laying angle on the vibration characteristics of composite laminates were studied.
It can be seen that the vibration frequency of laminate increased as the spring stiffness on the boundary increased. As the number of layers increased, the frequencies, except for the first-order frequency, showed a noticeable increase, and the frequency eventually became stable when the number of layers increased to a certain value. It can also be concluded that the effect of frequency regarding the ratio between thickness and width is significant; with the increase in h/b, the frequency of the plate vibration decreased.
Furthermore, in order to study the influencing factors and characteristics of the energy transfer, the harmonic response analysis and power flow analysis of the composited plate vibration showed that the vibration energy does not always transfer energy to the boundary according to the shortest path principle, and that there is also a certain form of energy cycle inside the plate structure. The distribution field and frequency of the energy flow had a significant influence. When the excitation force, excitation position and boundary conditions remain unchanged, the frequency change may also cause a complete change in the energy transfer.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The detailed expressions of the submatrices K ij and M ij in Equations (34)- (35) are provided as follows.