Three-Dimensional Vibration Analysis of Laminated Composite Rectangular Plate with Cutouts

The main purpose of this paper is to establish an analysis model of laminated composite rectangular plate with/without cutouts on the basis of three-dimensional elasticity theory and provide the exact three-dimensional solution. In the present work, the effect of the cutout is considered by subtracting the energies of cutouts from the total energy of the entire plate. The standard three-dimensional trigonometric cosine Fourier series and auxiliary Fourier series are chosen as admissible functions, and the Hamilton’s principle and Rayleigh-Ritz procedure are used to obtain the exact solution. In order to verify the effectiveness and reliability of the proposed method, some numerical results are obtained, and the results are compared with the ones available in the literature or finite element analysis. Finally, the effects of some key parameters which will affect the vibration characteristics are analyzed, and the non-dimensional frequency parameters are obtained, which can serve as benchmark data for the future research.


Introduction
A composite material is a heterogeneous material formed by combining two or more constituent materials with different properties integrated together to achieve enhanced structural performance. Owing to the high specific strength and stiffness, high fatigue resistance, light weight, and good thermal stability compared with conventional heterogeneous materials, applications of composite structures are employed in diverse fields that include vehicles, aircraft, turbines, architecture, submarines, and so on. [1,2]. With the continuous development of composite material manufacturing technology, composite materials have developed into one of the four major material systems in parallel with metal materials, polymer materials, inorganic non-metal materials, and laminated composite plates have become the basic structural components of most engineering structures. Therefore, it is necessary for engineers to have a thorough understanding of the vibration characteristics of laminated composite plates to ensure the reliability of structural performance.
In recent decades, the prediction of vibration characteristics of composite structure has attracted great interest from engineers and experts, and significant research advances have been achieved. Furthermore, a large amount of analytical theories have been established, such as the classical laminated plate theory (CLPT), first-order shear deformation theory (FSDT), higher-order shear deformation theory (HSDT), zigzag theory (ZZT), and three-dimensional elasticity theory (3DT).
The classical laminated plate theory is the simplest theory of all plate theories, Love [3] proposed this theory which is based on the Kirchhoff's hypothesis in the late 19th century. According to the theory, the displacement field is expressed as three displacement components along the x, y, z coordinate directions of a point on the midplane, thus the three-dimensional structures are reduced to two-dimensional structures. The simple form of displacement fields and reduction of variables result in simple equations and high calculation efficiency, it is applicable to the thin multi-layer composite structures. By neglecting the influence of the transverse shear deformation and rotational inertia, the vibration characteristics of the rather thick plates obtained by this theory are inaccurate, which will lead to smaller displacement and stress and higher natural frequency of the plates. Even for thin laminated composite plates, when the Young's modulus of the material is small, the transverse shear deformation cannot be neglected, and the resultant deviation of the theory increases sharply. CLPT remains a popular approach to solving the vibration problems of thin structures, and many research groups have proposed some theories in order to achieve better results based on CLPT. Gorman and Ding [4,5] exploited the superposition method and Galerkin method to obtain the exact analytical solution of laminated composite plates under free boundary conditions. Nallim et al. [6][7][8] used the Ritz method in conjunction with natural coordinates to express the geometry of general plates by using a set of beam characteristic orthogonal polynomials. Vibration characteristics of general quadrilateral plates, trapezoidal plates, skew plates, and rhomboidal plates with angle-ply layers have been studied. Secgin et al. [9] proposed a discrete singular convolution (DSC) approach by using a grid discretization based on distribution theory and wavelets. Cosentino and Weaver [10] developed a mixed theory by means of Reissner's variational approach based on Castigliano's principle of least work in conjunction with a Lagrange multiplier method to assess the effect of transverse shear stresses of thick composite laminates and sandwich plates.
In order to obtain more accurate vibration characteristics of laminated composite structures, Reissner [11] and Mindlin [12] proposed the first-order shear deformation theory by including the transverse shear deformation, and assumed that the shear force and strain are constant in the thickness coordinate direction, which is not consistent with the parabola distribution of shear force and strain. Therefore, the shear correction factor has been introduced into this theory to correct the shear strain and stress [13][14][15], which is related to the material, geometry, loads, boundary conditions, and laminations. The results obtained by the first-order shear deformation theory largely depend on the shear correction factor which is generally selected as 5/6 or π 2 /12. Subsequently, the higher-order shear deformation theory was proposed by Reddy [16,17] to overcome the shortcomings of CLPT or FSDT. Liew et al. [18] adopted the first-order shear deformation theory in the moving least squares differential quadrature procedure by using the moving least squares shape functions and their partial derivatives to obtain the weighting coefficients to predict the free vibration behavior of moderately thick symmetrically laminated composite plates. Karami et al. [19] applied the differential quadrature method (DQM) for free vibration analysis of moderately thick composite plates with edges elastically restrained against translation and rotation. The governing differential equations with their boundary conditions are transformed into algebraic equations by using differential quadrature rules in order to establish the eigenvalue equations, and the results showed that the DQM can yield an accurate solution using few grid points. Ngo-Cong et al. [20] employed one-dimensional integrated RBF networks instead of conventional differentiated RBF networks to approximate the field variables, and the rectangular or non-rectangular plates are discretized by means of Cartesian grids. Carrera [21] reduced the third-order model from five displacement variables to three by imposing homogeneous stress conditions with correspondence to the plate top-surface, and then modified it to apply to non-homogeneous stress conditions. Closed form solution has been obtained for both stresses and displacements in the case of harmonic loadings and simply supported boundary conditions. Chen et al. [22] employed the p-Ritz method by using the polynomials as the admissible trial displacement and rotation functions to study the vibration of laminated plates. Aydogdu [23] chose different types of shape functions to determine the distribution of the transverse shear strains and stresses along the thickness according to 3-D results to present a new higher order shear deformation theory. The new shear model was used to analyze bending, free vibration, and bulking of laminated composite plates. Ferreira et al. [24] applied a multi-quadrics radial basis function method and a third-order shear deformation theory to solve the analysis of isotropic and symmetric laminated composite thick beams and plates. Liu et al. [25] used Materials 2020, 13, 3113 3 of 23 the radial basis functions with polynomial reproduction to present the problem domain by a set of scattered nodes in its support domain. The static deflection, free vibration, and bulking analysis of laminated composite plates were presented.
Generally, the solution of higher-order shear deformation theory is more accurate than that obtained by using the first-order shear deformation theory due to the shear strain being distributed as a cubic function in the thickness direction; however, they introduce rather sophisticated formulations and boundary terms that are not easily applicable or yet understood [26]. All of the above methods belong to the equivalent two-dimensional methods, and the two-dimensional methods cannot satisfy the piecewise continuous displacement requirement and insufficient considerations for the transverse stress fields between layers. Murakami [27] developed the zigzag theory by including a zigzag shaped function to approximate the thickness variation of in-plane displacement. Sciuva [28,29] proposed a piecewise linear zigzag theory and then a non-linear third-order zigzag theory which accounts for continuous inter-laminar transverse shearing stresses at the interfaces between any two adjacent layers. Kapuria et al. [30][31][32][33] applied a new zigzag theory to analyze the static, dynamic, bulking, and thermal problems of laminated composite beams or plates. Pandit et al. [34] proposed an improved higher order zigzag theory for the static analysis of laminated sandwich plates with soft compressible core. The variation of in-plane displacements is assumed to be cubic for both the face sheets and the core, the transverse displacement is assumed to vary quadratically within the core while it remains constant through the faces. Khandelwal et al. [35] developed an improved FE plate model based on refined higher order shear deformation theory (RHSDT) and a least square error method (LSE) to accurately predict the deflections and stresses of composite and sandwich laminates. The C 1 zigzag theory was proposed based on the nine node C 0 element which satisfies the inter-laminar shear stress continuity conditions at the layer interfaces and zero transverse shear stress conditions at the top and bottom of the plate. The development of zigzag theory can be found in [36]. Generally, there are many variables in the zigzag theory and the displacement field functions are relatively complicated.
The three-dimensional elasticity theory does not rely on any hypothesis, each sublayer of the structure is regarded as a three-dimensional entity which is described as an independent function, and it satisfies the continuity of interlaminar displacement and stress requirement. Liew et al. [37] presented a continuum three-dimensional Ritz formulation by selecting sets of orthogonally generated polynomial functions as shape functions. The frequency parameters and mode shapes of thick plates have been solved with five practical groups of boundary conditions. Senthil et al. [38] presented a three-dimensional exact analytical solution and benchmark results for the free and forced vibrations of simply supported functionally graded rectangular plates by using suitable displacement functions in order to reduce the governing partial differential equations to a set of coupled ordinary differential equations in the thickness direction, which are solved by the power series method. Zhou et al. [39][40][41][42][43] used Chebyshev polynomials and static beam functions as admissible functions in conjunction with Rayleigh-Ritz method to study the three-dimensional vibrations of rectangular, skew, elliptical, and circular annular plates, as well as solid and hollow circular cylinders with different boundary conditions. Qu et al. [44] employed a multilevel partitioning hierarchy, viz., multilayered parallelepiped, individual layer and layer segment, which are based on the exact three-dimensional elasticity theory to analyze the free and transverse vibrations of multilayered laminated composite and sandwich beams, plates, and solids with various boundary conditions. The displacement components of each layer segment are approximated as the product of orthogonal polynomials and/or trigonometric functions.
The structure with cutouts is common in the practical engineering, the cutouts are used to reduce the weight or provide operational and maintenance access. The existence of cutouts changes the dynamic characteristics of the structure and damages the service life of the structures. It is necessary to study the effect of the cutouts on the vibration response of the structures. Liew et al. [45] employed a domain decomposition method by using a basic L-shaped element which is divided into appropriate sub-domains that are dependent upon the location of the cutouts as the basic element to solve the free vibration of plates with central cutouts. Kumar et al. [46] developed a finite element formulation based on higher-order shear deformation theory and Hamilton's principle to study the vibration response of thick square composite plates with a central rectangular cutout. The effect of material orthotropy, boundary conditions, side-to-thickness ratio, delamination size, and location around the cutout is investigated. Sakiyama et al. [47] proposed an approximate method and transformed the problem into an equivalent square plate with non-uniform thickness by considering the cutout as an extremely thin part of the equivalent plate. The Green function is used to obtain the discrete solution and characteristic equation of the plate. Laura et al. [48][49][50][51] have made some achievements in the vibration characteristics of plates with cutouts by applying the Rayleigh-Ritz variational method. Shufrin et al. [52] presented a new semi-analytical variational extended the Kantorovich method to model rectangular plates with variable thickness and cutouts. The plate thickness and deflections are represented as a finite sum of multiplications of one-dimensional functions. Kwak et al. [53] developed an independent coordinated coupling method (ICCM) in which the energies of the plate domain and the cutout domain are derived independently and the two independent coordinates are coupled by imposing kinematic relations. The vibration analysis of a rectangular plate with a rectangular cutout or a circular cutout is investigated. Huang [54,55] then applied this method to solve the problem of orthotropic composite laminated plate and functionally graded carbon nanotube-reinforced plate. The existing literature surveys show that most of the studies on the vibration of rectangular plates with cutouts are based on the traditional two-dimensional theories, and that three-dimensional exact solutions are very rare.
In this paper, a unified analysis model of the vibration characteristics of the laminated composite plates with/without cutouts is established, and the three-dimensional exact solution is provided with different boundary conditions. According to the energy principle, the Lagrangian energy functions of the plate and the cutout are obtained, and then a standard three-dimensional trigonometric cosine function and auxiliary Fourier series are selected as the admissible functions. The basic idea of solving the free vibration of a plate with a cutout is to subtract the energy of the cutout part from the total energy of the plate. The governing equations are solved by using the Hamilton's principle and Rayleigh-Ritz method. In order to verify the reliability of the method, some numerical examples are carried out to compare with those available in the literature and the finite element analysis. Furthermore, the effect of some key parameters which may affect the vibration characteristics is conducted, including the position of the cutouts, the size of the cutouts, the laminations, the boundary conditions, and the layer fiber direction angle. The non-dimensional frequency parameters are shown which can serve as a benchmark solution for future research.

Description of the Model
Consider a multilayered laminated composite rectangular plate with a rectangular cutout, as shown in Figure 1. The global coordinate system (x, y, z) which located in the middle surface of the structure is considered to present the dimensions. The length, width, and thickness of the plate are assumed as a, b, and h, the dimensions of the rectangular cutout in the x and y directions are c and d, the distances of the cutout in the x and y directions with respect to the coordinated point O are x a and y b , the coordinates of the cutout center are x h and y h , whereas x h = x a + c/2 and y h = y b + d/2, respectively. The middle surface displacements of the plate in the x, y, and z directions are denoted by u, v, and w, respectively. The number of the layers of the laminated composite plate is assumed to be N, and the kth layer has a thickness of h k , where the h k = z k+1 − z k , and the z k+1 and z k are the distances from the top surface and the bottom surface of the layer to the referenced middle surface. The principal coordinates of the laminated composite plate in the kth layer are denoted as 1, 2, and 3, the angle between the material axis and x coordinate is denoted by ϑ k . The artificial springs are introduced in order to simulate different kinds of boundary conditions, three groups of springs k u , k v , and k w , which are distributed and kw, which are distributed uniformly along the edges. The stiffness of the springs can take any value from zero to infinity to model the classical boundary conditions or elastic boundary conditions.

Stress-Strain Relations and Stress Resultants
Based on the three-dimensional elasticity theory and small deformation theory, the straindisplacement relations of the kth layer for the plate are of the form where the u k , v k , and w k are the middle surface displacement components of an arbitrary point of the plate in the x, y, and z directions, respectively. k  are the normal and shear strain components in the x, y, and z coordinate system. According to the generalized Hooke's law, the corresponding stress-strain relations in the kth layer of the plate are obtained as

Stress-Strain Relations and Stress Resultants
Based on the three-dimensional elasticity theory and small deformation theory, the straindisplacement relations of the kth layer for the plate are of the form where the u k , v k , and w k are the middle surface displacement components of an arbitrary point of the plate in the x, y, and z directions, respectively. ε k x , ε k y , ε k z , γ k yz , γ k xz , and γ k xy are the normal and shear strain components in the x, y, and z coordinate system.
According to the generalized Hooke's law, the corresponding stress-strain relations in the kth layer of the plate are obtained as where the σ k x , σ k y , σ k z , τ k yz , τ k xz , and τ k xy are the normal and shear stress components of the kth layer of the plate. The constants C k ij (i, j = 1, 2, · · · , 6) are elastic stiffness coefficients of the layer, which can be expressed as where C k ij (i, j = 1, 2, · · · , 6) are material constants in the principal coordinate system (1, 2, and 3), which are expressed as where E k 1 ,E k 2 , and E k 3 are moduli of elasticity of the kth layer in the principle directions, G k 12 , G k 13 , and G k 23 are shear moduli and µ k ij (i, j = 1, 2, 3, i j) are Poisson's ratios. According to the symmetry of coefficients in the stress-strain relationship, the constants satisfy the following An isotropic model can be obtained by letting: , where E and µ are elastic modulus and Poisson's ratio of the isotropic material, respectively.

Energy Functions
As shown in Figure 1, according to the three-dimensional elasticity theory, the energy functions of the plate and the cutouts are established independently in the global coordinate system, which can be expressed as where T p is the kinetic energy of the entire plate, T h is the kinetic energy of the cutout part, U p is the strain energy of the entire plate, U h is the strain energy of the cutout part, U sp is the potential energy stored in the boundary springs, and W e is the work done by the external force.
The basic idea to solve the vibration of a plate with cutouts is to subtract the energies of the cutouts part from the total energy of the plate without cutouts. In that sense, the total energy of a plate with a cutout is express as Equation (38), a plate without cutouts can be seen as a special case of a plate with cutouts, the total energy of a plate without cutouts is expressed as Equation (39)

Admissible Displacement Functions
The selection of the admissible displacement functions has a significant effect on the convergence and efficiency of the results. According to [56], the author expressed the displacement field functions as Fourier series with auxiliary polynomials, the integration process is complicated due to the non-orthogonality between Fourier series and polynomials. With the increase of the cutout number, the solution process will become more and more complex, therefore, the selection of the admissible displacement functions is significant to reduce the complexity of the calculation. In this paper, the admissible displacement functions are expressed in the form of complete trigonometric Fourier series expansions due to the orthogonality and completeness, as well as the excellent stability in numerical calculations when using the Rayleigh-Ritz method. The admissible functions are expressed as the component segment standard cosine Fourier series functions, and the potential derivative discontinuities of the displacements at the boundary are transferred to the added auxiliary trigonometric terms. The admissible displacement functions are expressed as three variables separated along the x, y, and z directions.
where λ m = mπ/a, λ n = nπ/b, λ q = qπ/h, and A mnq , B mnq , and C mnq are the unknown Fourier coefficients of the three-dimensional Fourier series expansions for the admissible displacement It is easy to verify that ξ 1b , ξ 2b , ξ 1c and ξ 2c have the same expression as ξ 1a and ξ 2a . The six supplementary functions are introduced to solve the potential discontinuities of displacement functions of three-dimensional plate structure at the boundary conditions. The displacement functions can satisfy arbitrary boundary conditions and will significantly improve the convergence of the solution compared with the polynomials and reduce the calculation time due to the orthogonality of the trigonometric Fourier series.

Solution Procedure
In order to obtain the exact solution of the governing equations and boundary conditions of the plate with a cutout, the Hamilton's principle is applied. Substituting the admissible displacement functions Equation (40) to Equation (43) into the Lagrangian energy functions Equation (32) to Equation (39), and minimizing the function L with respect to the unknown Fourier coefficients of the admissible displacement functions, one can easily obtain The Equations of motion for the plate with a cutout can be yielded and are summed up in the matrix form as where K is the stiffness matrix of the plate and M is the mass matrix. The stiffness matrix and the mass matrix both are symmetric and they can be expressed as X is the column vector which contains the unknown Fourier coefficients. The natural frequencies and mode shapes of laminated composite rectangular plate with cutouts can be obtained by solving the eigenvalues of Equation (47). Due to the large number of formula terms of K and M, the detailed calculation expressions are not given here.

Results and Discussion
In this part, numerical examples pertaining to the free vibration analysis of multilayered laminated composite rectangular or square plates with/without cutouts are presented in order to demonstrate the validity and efficiency of the described method. The results are compared with the solutions available in the existing literature and the finite element analysis. The effect of eccentricity of cutout positions and cutout sizes are studied for free vibration behavior with different boundary conditions. Numerical results are obtained with different layered laminations of symmetric and antisymmetric cross-ply and angle-ply.
Four groups of springs are introduced to simulate different kinds of boundary conditions, including free (F), simply supported (S), clamped supported (C), and elastic supported (E) boundary conditions, the stiffness of the springs are given in Table 1. Unless other stated, all layers of the plate have equal thickness and the same material parameters. The geometrical dimensions and material properties of the plate are given as follows The non-dimensional frequency parameters are expressed as Ω = ωb 2 ρ/E 2 h 2 .

Validity and Convergence of The Method
In this section, some numerical examples are carried out to verify the validity and the convergence of the proposed method in terms of thickness, boundary conditions, and cutout sizes. The accuracy of the method is confirmed by a comparison study with the data published in the literature of rectangular plate with/without considering cutouts, the data are given in Tables 2 and 3. The symbol '-' in Tables 2  and 3 indicate that the frequencies were not considered in the reference work.  First, for the case of plate without cutouts, the first six natural frequencies of a three-layered, cross-ply (0 • /90 • /0 • ) laminated square plate without cutouts under simply supported boundary conditions are presented in Table 2. Three kinds of thickness-length ratios, i.e., h/a = 0.01, 0.1, 1 corresponding to thin, moderately thick, and thick plates are performed, whereas h/a = 1 represents a solid cube, respectively. From Table 2, it is obvious that the proposed method has fast convergence behavior. The maximum error of the 9 × 9 × 9 truncated configuration compared with [44] is 0.086%, which proves the validity and accuracy of this method. In the following examples, the Fourier series is truncated to M × N × Q = 9 × 9 × 9. In Table 3, the first six natural frequencies of a four-layered, angle-ply (45 • /−45 • /45 • /−45 • ) laminated square plate without cutouts are presented, and the FFFF, SFSF, FFCF boundary conditions with different thickness-length ratios are performed in the comparison. The geometrical dimensions and material properties are the same as for Table 2. It can be seen from the table that the maximum error with the first-four natural frequencies provided by Qu et al. [44] is 0.155%, and the present results are in good agreement with those results of Ref. [44].
Second, in order to verify the reliability of solving the vibration characteristics of a plate with cutouts, a two-layered, cross-ply (0 • /90 • ) laminated square plate with central square cutouts or rectangular cutouts under simply supported boundary condition has been considered, and the first six natural frequency parameters are performed in Table 4. The results of the proposed method are compared with those obtained by Sheikh et al. [57]. The geometrical dimensions and material properties of the plate are given as follows: b = 1 m, b/a = 1, h/b = 0.01; E 1 = 25E 2 , E 2 = E 3 = 2 GPa, G 23 /E 2 = 0.2, G 12 = G 13 = 0.5E 2 , µ 12 = µ 12 = µ 12 = 0.25, ρ = 1500 kg/m 3 . From the table, it can be seen that the difference between the data of the proposed method and Sheikh et al. [57] is less than 0.556%, which shows that the present method is able to accurately predict the frequency parameters of a rectangular plate with cutouts. The small differences are related to different laminated plate theory and solution program methods.

Effect of Eccentricity of Cutout Position
In this section, the effect of eccentricity of cutout positions is studied for antisymmetric cross-ply (0 • /90 • /0 • /90 • ) laminated rectangular plate with CCCC (four edges clamped) boundary condition. The geometry dimensions are: h/b = 0.2, b = 1 m, a/b = 2, c/a = 0.2, d/b = 0.2. The coordinates of the cutout center ratio, thus eccentricity ratio, x h /a and y h /b are varied from 0.2 to 0.8 along the x direction and y direction, separately. Tables 5 and 6 show some of the first 20 non-dimensional frequency parameters compared with the solution of finite element method. The finite element method results are obtained from 3D model by using the software ANSYS Workbench, and the element size is selected as 5 × 5 mm. The present solution agrees well with the results of finite element method, and the data obtained by finite element method are slightly larger than the data obtained by the present method because of the finite element method needs to refine mesh to obtain enough accuracy. The data are symmetric about the geometric centers due to the symmetry of the geometry.
In Figure 2, variation of the lowest four order frequency parameters against the cutout position ratio are shown for the eccentricity in both x and y directions. From Figure 2a, we can see that the frequency firstly increases and then decreases for the first and second order frequency parameters when the eccentricity in the x direction goes from 0.2 to 0.5, while the third and fourth orders have reverse behavior, respectively. From Figure 2b, there is little variation in the second and third orders with the eccentricity ratio varied in the y direction. However, for the fourth order, the frequency parameter has gone through the process of decrease after going up for y h /b in the range of 0.2-0.5. For the first order, the fundamental frequency increases and then declines both in the x direction and y direction, and reaches its wave crest at 0.5. It can conclude that, for a rectangular plate with different aspect ratios, the cutout positions have different effects on the vibration characteristics. Nevertheless, for a square plate with a square cutout, it is conceivable that the frequency parameter variation is consistent in the x direction and y direction.

Effect of Cutout Size
Three examples are selected to demonstrate the effect of the cutout size with different boundary conditions, the cutout size ratio is varied from 0 to 0.8. In Table 7, some of the first 20 non-dimensional frequency parameters for the symmetric cross-ply (0°/90°/0°/90°) lamination with different cutout sizes under CCCC boundary conditions are examined. From Table 7, it can be seen that increasing of the cutout size of the plate increases the fundamental frequency parameters. The increasing trend is due to the introduction of the cutout; however, the loss of mass gets more significance over the loss of stiffness and hence the frequency parameters increase. Table 8 shows the data of antisymmetric cross-ply (0°/90°/90°/0°) laminated square plate with different cutout sizes under CSCS boundary condition, and Table 9 for angle-ply (45°/−45°/45°/−45°) laminated square plate with different cutout sizes under CFCF boundary condition. It is evident that from Table 8 and Table 9, the fundamental frequency parameters increase and then decrease when the cutout size ratio is varied from 0 to 0.8. The typical mode shapes corresponding to the first four orders of vibration for the plate compared with those obtained from the finite element method are given in Figures 3-5.

Effect of Cutout Size
Three examples are selected to demonstrate the effect of the cutout size with different boundary conditions, the cutout size ratio is varied from 0 to 0.8. In Table 7, some of the first 20 non-dimensional frequency parameters for the symmetric cross-ply (0 • /90 • /0 • /90 • ) lamination with different cutout sizes under CCCC boundary conditions are examined. From Table 7, it can be seen that increasing of the cutout size of the plate increases the fundamental frequency parameters. The increasing trend is due to the introduction of the cutout; however, the loss of mass gets more significance over the loss of stiffness and hence the frequency parameters increase. Table 8 shows the data of antisymmetric cross-ply (0 • /90 • /90 • /0 • ) laminated square plate with different cutout sizes under CSCS boundary condition, and Table 9 for angle-ply (45 • /−45 • /45 • /−45 • ) laminated square plate with different cutout sizes under CFCF boundary condition. It is evident that from Tables 8 and 9, the fundamental frequency parameters increase and then decrease when the cutout size ratio is varied from 0 to 0.8. The typical mode shapes corresponding to the first four orders of vibration for the plate compared with those obtained from the finite element method are given in Figures 3-5.

Effect of Lamination
To further study the effect of the number of layers on the frequency characteristic, the fundamental frequency parameters of two-layered, three-layered, four-layered, six-layered, and eight-layered, cross-ply and angle-ply plates with different cutout sizes are presented, the data are given in Appendix A1. Figure 6 shows the fundamental frequency parameters of cross-ply and angle-ply laminations with even numbers of layers under CCCC and SSSS boundary conditions. From Figure 6, it can be noticed that the frequency parameters increase with the increase of the number of the layers for both cross-ply and angle-ply. From Figure 6a,b, it is obvious that for CCCC boundary condition, the frequency parameters increase with the increase of cutout size ratio, and when the cutout size ratio exceeds 0.6, the frequency parameters increase sharply. However, it is somewhat complicated for SSSS boundary condition. From Figure 6c, for cross-ply lamination, the trend of frequency parameters declines first and then rises, when the cutout size ratio reaches a certain value which is respectively equal to 0.7, the frequency parameter decreases sharply. For angle-ply lamination, there is only a little bit of an upward trend when the cutout size ratio is less than 0.3, while after that, the frequency parameter trend declines with the increase of cutout sizes.

Effect of Lamination
To further study the effect of the number of layers on the frequency characteristic, the fundamental frequency parameters of two-layered, three-layered, four-layered, six-layered, and eight-layered, cross-ply and angle-ply plates with different cutout sizes are presented, the data are given in Table A1. Figure 6 shows the fundamental frequency parameters of cross-ply and angle-ply laminations with even numbers of layers under CCCC and SSSS boundary conditions. From Figure 6, it can be noticed that the frequency parameters increase with the increase of the number of the layers for both cross-ply and angle-ply. From Figure 6a,b, it is obvious that for CCCC boundary condition, the frequency parameters increase with the increase of cutout size ratio, and when the cutout size ratio exceeds 0.6, the frequency parameters increase sharply. However, it is somewhat complicated for SSSS boundary condition. From Figure 6c, for cross-ply lamination, the trend of frequency parameters declines first and then rises, when the cutout size ratio reaches a certain value which is respectively equal to 0.7, the frequency parameter decreases sharply. For angle-ply lamination, there is only a little bit of an upward trend when the cutout size ratio is less than 0.3, while after that, the frequency parameter trend declines with the increase of cutout sizes.

Effect of Boundary Condition
Next, to research the effect of boundary conditions, a four-layered, symmetric, and antisymmetric composite plate under six groups of boundary conditions, i.e., CCCC, CSCC, CSSC, CSCS, CSSS, and SSSS are considered. The fundamental frequency parameters are calculated for various cutout sizes and the data are detailed in Table 10.

Effect of Boundary Condition
Next, to research the effect of boundary conditions, a four-layered, symmetric, and antisymmetric composite plate under six groups of boundary conditions, i.e., CCCC, CSCC, CSSC, CSCS, CSSS, and SSSS are considered. The fundamental frequency parameters are calculated for various cutout sizes and the data are detailed in Table 10. It can be seen that the fundamental frequency parameters of antisymmetric structures are slightly larger than those of symmetric structures and the maximum frequency parameters are under the CCCC boundary conditions. Comparing the boundary conditions of CCCC and CSCC, for the same cutout size ratios, the boundary constrain is released from clamp supported to simply supported, the stiffness becomes smaller and the mass remains unchanged, thus the frequency parameter decreases both for symmetric and antisymmetric structures, cross-ply and angle-ply laminations. For different boundary conditions, with the increase of the cutout sizes, both of the mass matrix and the stiffness matrix get smaller, the frequency parameters tend to rise or drop. When the frequency parameters increase, mainly because the loss mass is more significant than the loss of stiffness, whereas the reason is reverse for the frequency parameters decreasing.

Effect of Layer Fiber Direction Angle
The effect of layer fiber direction angle on the frequency parameters of laminated composite plate with cutouts is studied in the following section. The detailed data are given in Table A2. The first six frequency parameters of a three-layered angle-ply (0 • /ϑ • /0 • ) square plate with a center cutout under CFCF boundary condition are investigated, where ϑ is the fiber angle in the range of 0 • -90 • . In Figure 7, the curves of frequency parameters changing with the fiber orientations under various cutout sizes are plotted. The fiber direction angle is varied from 0 • to 180 • with a step of 15 • . All of the curves are symmetric about ϑ = 90 • , and the trend of frequency parameters is not changing in the same way with the increase of cutout size ratio. The fundamental frequency parameters increase and then decrease with the increase of the cutout size ratio. The trend can be regrouped in three groups according to the performance: Group I, a smaller cutout size ratio, where 0 ≤ c/a ≤ 0.2; Group II, 0.2 < c/a ≤ 0.5; Group III, a larger cutout size ratio, where 0.5 < c/a ≤ 0.8. The trend in each group is generally consistent. It can be observed that the increment of the fiber direction angle from 0 • to 90 • results in the increase of the lowest two frequency parameters of the plate. For Group III, the frequency parameters climb up with the fiber direction angle changed from 0 • to 90 • . When the cutout size ratio is small, the trend of the lower-order modal is simpler, but for higher-order modal, the trend is more complex. This is mainly because the bending stiffness plays a dominant role for the lower-order modal, while for the higher-order modal, the bending and shear deformation are coupled with each other, which make the vibration behavior more complex.

Conclusions
In this paper, a unified three-dimensional solution is provided to study the free vibration characteristics of laminated composite rectangular plate with/without cutouts by using the energy principles in conjunction with Rayleigh-Ritz method. The energy equations of the plate and the cutouts which are based on the three-dimensional elasticity theory are established independently. The governing equation was obtained by subtracting the energies of cutouts from the total energies of the entire plate. The displacement field is represented as a standard three-dimensional trigonometric cosine Fourier series and modified auxiliary Fourier series which can satisfy the boundary conditions and reduce the computational complexity. The three-dimensional solution is obtained by using the Hamilton's principle to solve the unknown Fourier coefficients. The advantage of this method is that the unified model can simulate the different boundary conditions by changing the values of spring stiffness and there is no need to establish various models separately. For the free vibration characteristics of laminated composite rectangular plate with/without cutouts under various boundary conditions, the accuracy of this method is better than the results of the literature and finite element analysis. Furthermore, some numerical examples are carried out to show the performance of the proposed method and we can arrive at the following conclusions: First, the frequency parameters of laminated composite square plate with/without cutouts under various boundary conditions are calculated and compared with those available in the literature. It can be seen from the comparison that the proposed method is accurate in predicting the vibration characteristics of plates with cutouts.
Second, the influence of some key parameters, such as the cutout positions, the cutout sizes, the boundary conditions, the laminations, and the fiber orientations on the frequency parameters are discussed. The cutout positions have a significant effect on the frequency parameters, and the maximum fundamental frequency occurs when the cutout is located in the center of plate. From the variations of the cutout size ratio, it can be concluded that with the increase of cutout size ratio under various boundary conditions, the loss of mass and stiffness play a different role in determining the frequency parameters. When the total thickness is constant, the frequency parameter grows as the number of the layers increases. The fundamental frequency parameters climb up as the fiber direction angle varies from 0 • -90 • .

Conflicts of Interest:
The authors declare that there are no conflict of interest regarding the publication of this paper.