High frequency modes meshfree analysis of Reissner-Mindlin plates

Finite element method (FEM) is well used for modeling plate structures. Meshfree methods, on the other hand, applied to the analysis of plate structures lag a little behind, but their great advantages and potential benefits of no meshing prompt continued studies into practical developments and applications. In this work, we present new numerical results of high frequency modes for plates using a meshfree shearlocking-free method. The present formulation is based on ReissnereMindlin plate theory and the recently developed moving Kriging interpolation (MK). High frequencies of plates are numerically explored through numerical examples for both thick and thin plates with different boundaries. We first present formulations and then provide verification of the approach. High frequency modes are compared with existing reference solutions and showing that the developed method can be used at very high frequencies, e.g. 500th mode, without any numerical instability. © 2016 The Authors. Publishing services by Elsevier B.V. on behalf of Vietnam National University, Hanoi. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Eigenvalue analysis of plate structures is an important research area to designers and researchers because of their wide applications in engineering such as aerospace, marine, ship building, and civil. Many different theories accounting for such plate structures have been developed, see e.g., [1e5]. One of the most successful theories is based on the Kirchhoff hypothesis for thin plates neglecting the transverse shear strains [1,5], but this strong assumption causes the main reason for the inaccuracy of the solutions, especially at high modes. In order to accommodate the transverse shear strain effect, a theory, which is based on the ReissnereMindlin's assumption, has been introduced as a remarkable candidate and commonly used for thick plate analysis [2e5].
Analytical solutions to free vibration of thick plates are certainly available and extended to analyze the vibration of functionally graded material plates [46e48] but unfortunately they are limited to structures which consist of simple geometries and boundary conditions and often, the exact solutions are very difficult to obtain. Thus, approximate solutions of eigenfrequency plates problems at high modes derived from numerical approaches are often chosen. The development of numerical approaches, in particular, for plates has led the invention of some important computational methods such as Ritz method [6], isogeometric analysis [7], finite strip method [8], spline finite strip method [9e11], finite element method (FEM) [12e16], discrete singular convolution (DSC) method [17,18], and DSC-Ritz method [19,20]. The FEM is well-advanced and is one of the most popular techniques for practice, but till has some inherent disadvantages, e.g., mesh distortion. In order to avoid such disadvantages, meshfree or meshless methods have been developed, and some superior advantages over the classical numerical ones have illustrated, see e.g., [21e25]. Unlike the conventional approaches, the entire domain of interest is discretized by a set of scattered nodes in meshfree methods irrespective of any connectivity.
Plate structures with high frequency modes have been analyzed using numerical methods, for instance, by FEM [26]; DSC method [17,19]; DSC-Ritz method [20]. The hierarchical FEM by Beslin et al. [27] was to reduce the well-known numerical instability of the conventional p-version FEM [28], due to computer's round-off errors. For more information related to this issue, readers can refer to an elegant review done by Langley et al. [26]. This work is devoted to the numerical investigation of high frequency modes of plates. A meshfree method is thus adopted. We numerically demonstrate the applicability and performance of our meshfree moving Kriging interpolation method (MK) [29] to high frequency mode analysis of ReissnereMindlin plates without numerical instability. The meshfree MK [29] has recently been extended to other problems such as two-dimensional plane problems [30,31], shell structures [32], static deflections of thin plates [33], piezoelectric structures [34] and dynamic analysis of structures [35]. Another important shear-locking issue, which occurs when using thick plate theories to analyze for thin plates, is taken into account in the present formulation. To this end, a special technique proposed in [36], using the approximation functions for the rotational degrees of freedom (DOF) as the derivatives of the approximation function for the translational DOF, is incorporated into the present formulation to eliminate the shear-locking effect.
Most recent meshfree methods still have the same problem in dealing with the essential boundary conditions, although many efforts have been devoted to overcoming that subject and some particular techniques have been proposed to eliminate this difficulty in several ways, such as the Lagrange multipliers [22], penalty methods [21,37], coupling with the traditional FEM [38e42], and transformation method [43,44]. In other words, the MK is a wellknown geostatistical technique for spatial interpolation in geology and mining. The basic idea of the MK interpolation is that any unknown nodes can be interpolated from known scatter nodes in a sub-domain and move over any sub-domain [29]. The procedure is similar to the moving least-square (MLS) method [22,45], but the formulation employs the stochastic process instead of least-square process. The MK is smooth and continuous over the global domain and one of the superior advantages of the present method over the traditional ones. The Kronecker delta property is satisfied automatically. Hence, the essential boundary conditions are exactly imposed without any requirement of special treatment techniques as the conventional FEM.
Because the MSL approximation is not the interpolation function, this is a major drawback of the standard EFG method. Hence, the present work describes a means using the MK interpolation technique to high vibration modes analysis of plates. As far as the present authors' knowledge goes, no such task has been studied when this work is being reported. The paper is structured as follows. A meshless formulation for free vibration of Reiss-nereMindlin plates is presented in the next section, showing a brief description of governing equations and their weak form. Approximation of displacements is then presented in Section 3 and the corresponding discrete equation systems are given in Section 4. Numerical examples are presented and discussed in Section 5 dealing with natural frequencies of the square and circular plates at high modes. We shall end with a conclusion.

Formulation of ReissnereMindlin plates for high frequency variation analysis
In this section, formulation of ReissnereMindlin plates for the analysis high frequency modes is briefly presented [29]. A FSDT plate as depicted in Fig. 1 with two-dimensional mid-surface U3< 2 , boundary G ¼ vU, the thickness t and the transverse coordinate z is considered. The displacements u and v can be expressed as [43] u ¼ zb x (x) and v ¼ zb y (x), with x ¼ {x,y} Τ and independent  ; The assumption for displacement of three independent field variables u2H 1 0 ðUÞ Â ðH 1 0 ðUÞÞ 2 for ReissnereMindlin plates is The linear elastic material is assumed with Young's modulus E and Poisson's ratio n, strong form for free vibration of plates is given by [12,13] V$D b kðbÞ þ ltg þ t 3 12 where D t ¼ Et 3 /12(1 À n 2 ) is the flexural rigidity. The bending k and transverse shear g strains are expressed as where L b and L s are differential operator matrices and are explicitly given by The bending k and transverse shear g strains in Eqs. (6) and (7) can be rewritten as Fig. 1. Geometric notation of a FSDT plate [29].
High frequency modes of ReissnereMindlin plate are derived from the principle of virtual work under the assumptions of the FSDT plate theory [12,13,43]: find the natural frequencies u2< þ and 0s(w,b)2S such that aðb; hÞ þ ltðVw þ b; where Q is a set of the essential boundary conditions and the L 2 inner-products is [29] aðb; hÞ ¼ In meshfree methods implementation, the bounded domain U is discretized into a set of scattered n nodes, and each node is covered by a sub-domain U x associated with an appropriate influence domain such that U x 4U. The meshfree solution of high modes for ReissnereMindlin plate is to find natural frequencies u h 2< þ and where the meshfree approximation spaces, S h and S h 0 , are expressed as with P 1 (U x ) being the set of polynomials for each variable within the sub-domain U x 4U.
Dynamic equation by a minimization form of energy principle of virtual displacements incorporating the FSDT plate theory is [43] where du is the variation of displacement field u, € u is the secondorder derivatives of displacement over time or acceleration, B m is the matrix consisting of the mass density r and the thickness t B m ¼ r while D s is the tensor of shear modulus as 3. Meshfree approximation of field variables and treatment of shear-locking In this section, the MK meshfree approximation for field variables (i.e., deflection and rotations) for ReissnereMindlin plates and a technique for treatment of shear-locking effect are briefly presented [29]. Field variables of plates are the deflection w(x) and the two rotation components b x (x) and b y (x) at all nodes. The approximation is utilized through parameters of nodal values expressed in a group of nodes within a compact domain of support. This means that these values can be interpolated based on all nodal values x i (i2 [1,n]), where n is the total number of the nodes in U x so that U x 4U. Thus, the meshfree approximation The superscript h in Eq. (21) is omitted without loss of generality, i.e., where u I ¼ (w I ,b xI ,b yI ) Τ is the vector of nodal variables at node I whereas f I ,f xI and f yI are the MK shape functions defined by In this work formulations using the first-order derivatives of shape functions presented in [36] to eliminate the shear-locking is taken The matrices A and B are determined via where I is an unit matrix and p(x) in Eq. (20) is the polynomial with m basis functions For n coupling nodes, the n Â m matrix P is expressed as (20) is also given by where R(x i ,x j ) is the correlation function between any pair of the n nodes x i and x j . The correlation matrix R½Rðx i ; x j Þ nÂn is given by A Gaussian function with a correlation parameter q is employed where r ij ¼ x i À x j and q > 0 is a correlation parameter.
The quadratic basic function p T ðxÞ ¼ Â 1 x y x 2 y 2 xy Ã is taken throughout the study.
The firstand second-order derivatives of the shape function can be computed as The influence domain radius is determined by with d c being a characteristic length relative to the nodal spacing close to the interest point while a standing for a scaling factor. The MK shape functions f I (x j ) at node x I for interpolation node x j possess the Kronecker delta function property The order continuity of the MK interpolation is mostly dependent on the continuity of semivariogram. Since the Gaussian function Eq. (31) used in interpolation has high continuity, leading to that the MK interpolation also has high continuity. Other properties of the MK shape functions such as consistency can also be found in Refs. [29e31].

Meshfree discrete equations for high frequency analysis
Based on the preceding section on the variational form in Eq. (17), the bending strain and transverse shear strain for plates are By inserting Eqs. (22) and (35) into Eq. (17), discrete system of equations for vibration problems is obtained as where the global stiffness matrix K, which consists bending K b and transverse shear K s forms  (40) and the global mass matrix M A general solution of such a homogeneous equation is where i is the imaginary unit, b t indicates time and u is the eigenvector. Substituting Eq. (42) into Eq. (37), natural frequencies u is obtained solving the following eigenvalue equation For numerical integration, a background cell with 16 Gaussian points is used [29e31].

Numerical results of high frequency modes and discussion
High frequency modes results of FSDT plates with various boundary conditions derived from the proposed meshless are analyzed here. The boundaries of the plates, for convenience, are denoted as (F) completely free, (S) simply supported and (C) fully clamped edges. Throughout the paper, if not specified otherwise, Fig. 3. The rate convergence study with the SSSS square plates (t/a ¼ 0.1) for the first six modes using the proposed meshfree method.    5. Influence of the correlation parameter q on the natural dimensionless frequencies of the square plate (t/a ¼ 0.1) at low modes. This result is similar to that presented in [29]. Fig. 6. Influence of the correlation parameter q on the natural dimensionless frequencies of the square plate (t/a ¼ 0.1) at high modes. the following parameters are used: the Young's modulus E ¼ 200 Â 10 9 N/m 2 , the Poisson's ratio n ¼ 0.3 and the mass density r ¼ 8000kg/m 3 , the shear correction factor m ¼ 5/6 and the dimensionless frequency coefficient 6 ¼ ðu 2 a 4 rt=D t Þ 1=4 .

Convergence study
A square plate with a ¼ b ¼ 10m is considered. Since analytical solutions of this plate are available at low frequency modes, a convergence study of the method at low frequencies is explored. The dimensionless frequencies of a square plate accounting for CCCC [29] and SSSS boundaries are computed for different sets of regular distributed nodes, e.g., 7 Â 7, 9 Â 9, 11 Â 11, 13 Â 13 and 15 Â 15. The first six modes results of non-dimensional frequencies compared with exact solutions [45] are reported in Table 1. The frequency convergence u num /u exact (u num : meshfree solutions, u exact : analytical solutions) of square plates for the first six modes is also depicted in Fig. 2. Here Dh is the average spacing of scattered nodes in the domain. Compared with theoretical solutions, the frequencies obtained by the present method are in good agreement. Sufficient accuracy can be found for both the considered boundaries with a regular density of 13 Â 13 nodes, especially even for a course set of 9 Â 9 nodes the solution of the CCCC plate matches well with the exact one. Thus, we decide to use a pattern of 13 Â 13 nodes for all implementations unless specified.
Further convergence study is made to again verify the convergence rate of this meshfree method. The SSSS boundary associated with three regularly distributed nodes 7 Â 7(49), 9 Â 9(81) and 13 Â 13(169) is used. The first six modes are considered and their relative error plotted in a logelog plot is depicted in Fig. 3, showing a good convergence.

Shear-locking examination
Square plates under SSSS and CCCC boundaries are considered. The same parameters as above are used, except the thickness-span aspect ratio t/a ¼ 0.005 (thin plate). Table 2 presents the results of the first six modes calculated by the proposed method in comparison with the analytical solutions [45]. In Table 2, results obtained by using the elimination technique of the shear-locking are   named as "Shear-MK". The percentage errors in normalized frequencies estimated over the exact solutions are visualized in Fig. 4. As expected, the free of shear-locking is achieved when the Shear-MK is employed and large errors are found for the standard MK.

Effects of the correlation and scaling parameters
The correlation parameter q has some effects on the solutions,    Fig. 11. Influence of the thickness-span ratios on the dimensional frequencies 6 1 for the SSSS square plate. parameter varies from 0.1 to 50 for low frequencies and this range is wider for high frequencies. We examine low frequencies because of exact solutions, and thus it is easy to validate the results. The SSSS boundary is used here. Fig. 5 represents the percentage errors in non-dimensional natural frequencies at low modes estimated over the exact solutions [45], it can be seen that acceptable solution are gained if 1 q 10 is taken. Fig. 6 depicts dimensionless natural frequencies at high modes for each value of the correlation parameter. We found that 1 q < 10 can be selected for free vibration analysis of plates at high modes. We now decide to use q ¼ 5 for all implementations if not specified, otherwise.
Similarly, the scaling factor altering the high modes is analyzed, a correlation parameter of q ¼ 5 is used, and several scaling factors from 2.5 to 6 are considered for low modes and other higher values are for high modes. The results calculated for low and high modes are represented in Figs. 7 and 8, respectively. According to our own numerical experiments, we found that a range of 2.8 a 4 would be possible to be used for analyzing both low and high modes, and we thus decided to use a ¼ 3 for all implementations if not specified, otherwise.

Comparison study
A comparison of high frequencies of a square plate (a/b ¼ 1) among the present method and other existing reference solutions is explored.
The dimensionless natural frequencies the SSSS boundary and the thickness-span ratio t/a ¼ 0.1 are used. Table 3 and Fig. 9 show the frequency results at high modes up to 1500th obtained from the present MK meshfree method, the DSC-Ritz method with both the Shannon and the de la Vallee Poussin kernels [20] and the KirchhoffeMindlin relationship [20]. It can be seen that the frequencies calculated by the proposed method match well with the DSC-Ritz method for both given kernels. However, the results obtained from the Kirch-hoffeMindlin relationship also match very well with the DSC-Ritz and the present approach only at the modes below 112th and beyond that mode 112th the solutions of the KirchhoffeMindlin failed. The absence of shear deformation modes may cause such inaccuracy. As shown in Fig. 9 at the same modes after 112th, the KirchhoffeMindlin relationship offers higher frequencies than other methods, implying that less accuracy can be found for the KirchhoffeMindlin relationship when high frequency modes of thick plates are considered.

Effect of the length-to-width and the thickness-span ratios
The influence of length-to-width ratio for thick plates (t/a ¼ 0.2) on high frequencies is analyzed. This is because the natural frequencies may have significant variation when varying this aspect ratio.
The non-dimensional frequency coefficient  Fig. 10(aed), respectively. The high frequencies behave the same situation for all the considered boundaries, i.e., the frequencies increase with increasing the aspect ratios a/b. Unlike the length-to-width ratios, it can be observed that when the thickness-span ratio increases, the corresponding frequencies decrease.

Effect of the boundary
The influence of the different boundaries on the high modes is studied. A thick square plate (t/a ¼ 0.1) with different boundaries CCCC, SSSS, SCSC, CCCF, SFSF, CFFF and CFCF is studied. The nondimensional natural frequency coefficient is estimated by Fig. 12 represents the dimensionless frequencies calculated by the present method up to 450th modes and     Fig. 21. Influence of the thickness-span ratio on the dimensionless frequencies of a CCCC circular plate at high modes. Fig. 13 shows the calculated results for the lower modes. Highest frequency modes are found for the CCCC plate, whereas the modes for other SCSC, SSSS, CCCF, CFCF, SFSF and CFFF plates gradually decrease.

Mode shape analysis
In free vibration analysis of plates, mode shapes are often considered to view how vibratory structures look like especially at high modes. In this section, a thick square plate associated with SSSS boundary is taken with the thickness-span aspect ratio t/ a ¼ 0.1. Four different sets of six mode shapes are picked up typically from the low to the high frequencies that are plotted in a series of Fig. 14 (for modes 1st to 6th), Fig. 15 (for modes 90th to 95th), Fig. 16 (for modes 200th to 205th) and Fig. 17 (for modes 495th to 500th), respectively. It is easy to see that the wavelengths are decreased from the low to the high modes.

Circular plate
A circular plate shown in Fig. 18 is also considered to illustrate the applicability of the proposed method to arbitrary geometries at high modes. The SSSS and CCCC boundaries are taken into account. The problem parameters are taken the same as used in the rectangular above. The radius of the circular plate is indicated by R parameter.
A non-dimensional frequency coefficient is also employed. Figs. 19 and 20 present the influence of the radius of the plates on the dimensionless frequencies at high (up to 450th) and low (20th) modes undergone by both  CCCC and SSSS boundaries. In these results, the thickness-span ratio t/(2R) ¼ 0.1 is employed and the radius R is varied i.e. R ¼ 3, 5, 7 and 9, respectively, while other problem parameters are unchanged. It is found for both boundaries at high and low modes that the frequencies are increased when increasing the radius. Additionally, the influence of the thickness-span aspect ratio on the frequencies at high modes is considered and its results undergone by the CCCC boundary are given in Fig. 21. Various thickness-span ratios are taken the same as above such as t/a ¼ 0.01, 0.03, 0.06, 0.09, 0.1, 0.15 and 0.2. Likewise the results accounted for the rectangular plate above, it again confirms that the frequencies are decreased once the thickness-span ratio increases. Furthermore, several mode shapes shown in Figs. 22e25, from low to high frequencies are also provided in order to get a better observation. The wavelengths are decreased from the low to the high modes as expected.

Conclusions
In this paper, we present new numerical results of high frequency modes of ReissnereMindlin plates using an effective meshless method eliminating the shear-locking. The accuracy of the proposed formulation is demonstrated through numerical examples and the obtained results are analyzed and discussed in detail. The achieved results are compared with existing reference solutions and very good agreements are obtained. The influences of various aspect ratios and different boundaries dealt with both relatively thick and thin plates are considered. The developed method is efficient, robust, stable, accurate and free from the shearlocking effect. It has a good convergence and allows predicting at high frequency modes of the FSDT plates without numerical instabilities or numerical round-off errors. The details for the computational time (CPU-time) of the proposed formulation  compared over the conventional methods, e.g., the MLS-based EFG method, can be found in our previous work [33]. The application of the method to other complex problems is also possible.