Meshless Computational Strategy for Higher Order Strain Gradient Plate Models

: The present research focuses on the use of a meshless method for the solution of nanoplates by considering strain gradient thin plate theory. Unlike the most common ﬁnite element method, meshless methods do not rely on a domain decomposition. In the present approach approximating functions at collocation nodes are obtained by using radial basis functions which depend on shape parameters. The selection of such parameters can strongly inﬂuences the accuracy of the numerical technique. Therefore the authors are presenting some numerical benchmarks which involve the solution of nanoplates by employing an optimization approach for the evaluation of the undetermined shape parameters. Stability is discussed as well as numerical reliability against solutions taken for the existing literature.


Introduction
In the broad context of meshless methods Meshless Local Petrov-Galerkin (MLPG) [1] demonstrated to have strong capabilities to solve problems where weakly-singular traction and displacement boundary integral equations are involved. Moreover, other problems in this context have been analyzed and solved [2,3]. The main idea of meshless methods is to go beyond the current limitations of finite element models for analyzing problems in mechanics [4]. Unlike classical FEM shape functions are developed for a scattered set of collocation nodes [4,5] and their generation defines different meshless approaches. The one considered in the present work is named Radial Point Interpolation Method (RPIM) [6][7][8][9] which has the fundamental property of Kronecker delta function [10][11][12][13]. This approach makes the RPIM extremely easy to be implemented and used for the solution of any problem in structural mechanics.
It has been demonstrated to be extremely relevant in industrial applications for nanoengineering that nanoelectromechanical (NEMS) systems have been widely considered in the recent literature with several gradient elasticity problems [14,15]. It has been demonstrated in recent literature that nano-structures are used for micro-sized systems and devices such as biosensors, nano-actuators and nano-electro-mechanical systems [16].
Among the vast literature of nonlocal theories stress-driven models have been recently presented for nanobeams presenting analytical solutions. For instance closed-form solution of Bernoulli and Timoshenko type for Erigen-like formulations and stress field and a bi-exponential averaging kernel functions characterized by a scale parameter is presented in [17][18][19]. Nonlocal strain gradient Bernoulli like beam models by considering special bi-exponential averaging kernels and functionally graded materials has been presented in [20]. Nonlocal beam formulations have been presented within a thermodynamic framework, variational formulation within its analytical solution has been provided in [21]. Nanobeams can be subjected to axial loads which leads to buckling, such effects must be considered for a proper nanoengineering design, thus stress-driven buckling of nanobeams can be found in [22,23]. In the context of dynamic problems vibrations in nonlocal integral elasticity has been recently considered for beams and plates in [24,25].
Most of the nonlocal theories relies on homogenization approaches which aim at simplifying the problem by considering less modelling parameters in composite materials [26,27]. It is remarked that size effects and microstructures paved the way in presenting innovative and multiscale approches in solid mechanics [28,29].
It has been demonstrated by several researchers that higher-order elasticity is becoming of paramount importance for solid mechanics as mentioned in [30][31][32]. In most researchers the term nonlocality has been brought by Eringen [33] and by Eringen and Edelen [34] where it can be found that the constitutive relations have to be modified to take into account their dependency on the mechanical properties of the entire body and not only of the properties in the neighborhood of the material point. In this regard, in the present work nonlocal effects have be meaning introduced by Altan and Aifantis [35,36], which considered a simplified nonlocal model where all nonlocalities are concentrated in a gradient model similar to the one proposed by Mindlin [37]. Such approach, known as strain gradient theory, has been utilized also by others in other contexts in solid mechanics [38][39][40]. Strain gradient theory [41] has been demonstrated to be a constrained version of other higher-order unconstrained versions available in the literature such as couple stress [42][43][44] as well as micropolar theories [45][46][47].
In the framework of plate theories thin plate model is very common within the area of structural mechanics of investigating thin-walled structures [48,49]. Such approach can be easily extended in order to consider composite structures such as laminated [50,51] or sandwich ones [52].
To this day, the application of Mesh Free Methods for the analysis of structures modeled according to the strain gradient theory appears to be limited. Static analysis as well as buckling and free vibration of strain gradient beams and plates have been studied by means of different Mesh Free techniques [53][54][55][56]. However, none of the works mentioned makes use of the Mesh Free RPIM to perform the analysis.
In the present work the advantages of meshless methods are considered in the framework of strain gradient composite thin plate theory to solve such numerical problem.

Theoretical Background
The present work considers the problem of laminated thin plates of rectangular plane form of size a × b, where a, b indicate lengths with respect to x, y, respectively. The plate thickness is indicated with h and the plate represents the middle plane of the actual plate with z axis pointed normal to the x, y plate. Each ply which constitutes the stacking sequence is indicated with h k and the total thickness is computed as h = ∑ N L k=1 h k , where N L denotes the number of plies [57].
The present displacement field follows the classical thin plate theory and Cartesian displacements can be represented as [58] where u is the vector collecting the three middle plane displacements of the material point and D (s) is the derivative operation defined in [59]. The strain vector ε is given by in which ε (m) and ε (b) denote respectively the membrane and bending strains, respectively. They can be evaluated as follows where the meaning of the differential operators D (m) , D (b) is reported in [58,59]. Linear constitutive law is considered within the strain gradient theory as suggested in [40] allows to relate the membrane stresses in the k-th layer σ (k) to the corresponding strain components ε as shown below [60,61] in which the nonlocal parameter includes the micro/macro-scale interaction effects.
The dependency of the stresses on the strain distribution within the medium is emphasized by the presence of the Laplacian in Cartesian coordinate system: ∇ 2 = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 . On the other hand,Q (k) represents the plane stress-reduced stiffness coefficients matrix of the k-th layer. The termsQ (k) ij of this matrix depend on the orthotropic properties of the layer (Young's moduli E 1 , E 2 , Poisson's ratio ν 12 and shear modulus G 12 ), as well as by an arbitrary orientation θ (k) as indicated in the book [62]. the S N and S M are the stress resultants that can be defined as follows The differential operators D yy are defined in [58]. The constitutive operators A, B, D represent instead the membrane, membrane-bending coupling and bending stiffness matrices of the laminated composite plates [62].
In the current paper, isotropic and composite schemes are considered, therefore for some configurations B = 0, thus, membrane and bending behaviors are coupled. The variational form of the present equilibrium is represented by [58,63] The present meshless technique is applied to such variational statements of the problem.

Mesh Free Methods
The present numerical solution is obtained using a Mesh Free approach. As explained in detail in [4], with respect to traditional FEM, the two methods diverge at the stage of geometry discretization. In conventional FEM, the continuum is discretized by means of a mesh which also introduces predefined relationship between the nodes. Shape functions and system equations are then written for each element in the mesh.
In a Mesh Free method the domain is not discretized. It is represented by a set of nodes scattered in the domain and on its boundaries. Since no elements exist in this description, the shape functions and the system equations are written for the nodes. This characteristics allows to introduce a huge degree of flexibility, with respect to FEM, in term of adding or removing nodes from the domain.
The absence of the mesh implies the the relationship between the nodes representing the domain are not defined a priori. Therefore, a new concept has to be introduced to explain how the shape functions are constructed.
In Mesh Free methods, the field variable is interpolated using only the nodes falling within a local domain, within the problem domain, called support domain. The support domain can have different shapes, the most common being circular or rectangular, the latter being the choice for this work. The numerical procedure for the shape functions calculation is as follows. The support domain is first centered in a point of interest. The nodes falling inside it are identified and the shape functions are constructed using only those nodes. Once all the necessary calculations are performed, the support domain is moved to be centered in the next point of interest and the procedure starts anew. The so-called point of interest in which this local domain is centered can be either a node or, as chosen for this work, a Gauss point.
It has been explained how in mesh free methods the domain is represented by a set of scattered nodes rather than discretized by means of a mesh. Hence, it may seem contradictory to talk about Gauss points in such context. As a matter of fact, the Gauss points come from the necessity of numerically evaluate the integrals appearing in the weak form of the equation of motion, as shown in Equation (6).
For the sole purpose of integrals evaluation, in fact, a so-called background mesh needs to be introduced. The background mesh and nodal distribution are two independent objects, as shown in Figure 1, and its sole purpose is the integrals evaluation. As stated, the only nodes used to construct the shape functions are the ones falling within the support domain of the i-th point of interest. It follows that the size of the support domain is an important aspect to ensure that the shape functions are computed properly. The size of the local domain is evaluated as: where d c = ∆x 2 + ∆y 2 is the average nodal spacing and α s is a nondimensional parameter. This α s parameter has to be tuned to ensure that the appropriate amount of nodes is contained within the support domain.

Radial Point Interpolation Method (RPIM)
In the present approach solution is obtained in scattered points located in the given plate rectangular domain. The approximating polynomials involved possess the Kronecker delta property which allows a straightforward implementation of boundary conditions. In the following implementation, in plane displacements u and v, transverse displacement w as well as rotations w x , w y are included in the numerical implementation even though more primary variables are involved in the present formulation [58,59].
Since both isotropic and laminated plates are studied, the degrees of freedom for the different cases are taken as shown in Table 1. In the context of the strain gradient theory, the essential boundary conditions for the plate are shown in Table 2 where the four edges are identified by the values of the physical coordinates x and y.

Material
Degrees of Freedom Any other higher-order derivative is carried out by numerical derivation of the aforementioned parameters. Since both the deflection and its first derivatives are considered unknown, the Hermite-RPIM formulation is here presented. Let's consider a domain enclosing n arbitrarily scattered nodes. The approximation of the generic displacement w(x, y) can be expressed as: where R, R ,x and R ,y are the vectors including the radial basis functions (RBF) and their derivatives. The correspondent coefficients are indicated using vectors a, a x and a y . For the following numerical applications the well-known multi-quadrics (MQ) RBF is used in its general form where C = α C d c . Both q and α C are shape parameters that have to be tuned while d c is the average nodal spacing. The vectors of coefficients in Equation (8) can be obtained by enforcing the field function and its derivatives to be satisfied at all the n nodes falling within the support domain of the point of interest (x, y). The support domain is a local domain, typically circular or rectangular, centered in a point of interest which can either be a node or an integration point.
This leads to 3n linear equations where W is a vector containing the function values of the degrees of freedom considered in the collocation nodes and listed in Table 1. Thus, the independent parameter can be carried out as An example of what the shape functions look like as computed with this method is given in Figure 2. A squared domain, represented by 3 × 3 regularly distributed nodes, is considered and all of the nodes are used to construct the shape functions for this domain.    By following the same computational strategy of conventional finite element method [58] the algebraic form of the variational statement can be carried out because shape functions are evaluated at the collocation nodes. The needed integration is performed by following well-known Gauss integration rules. Solution of the present static problem is provided by Gauss elimination algorithm. The following section is dedicated to numerical results, stability and accuracy.

Applications
In the following sections, the numerical results of the analysis are shown and commented. Discussion about the degrees of freedom taken as variables in the problem and numerical convergence analysis are introduced, as well as plots of both convergence and deformed configurations which are also reported for clarity of representation.

Isotropic Plates
This section shows the results of the numerical analysis of isotropic Kirchhoff nanoplates modeled according to the second-order strain gradient theory and analyzed by means of a mesh free RPIM. In this case, there is no coupling between the in plane and out of plane behavior. Hence, the only unknown variables considered in the numerical implementation the are transverse displacement w and the corresponding rotations w x , w y . The numerical codes are developed in MATLAB.
The results in terms of mid transverse displacement are presented in the non-dimensional form as follows:w = 1000wD q z a 4 (12) where w is the central plate deflection, q z is the magnitude of the transverse external load and D is the bending rigidity D = Eh 3 /12(1 − ν 2 ). Different nodal densities are also taken into account. Nanoplates represented by 3 × 3, 5 × 5, 7 × 7 and 11 × 11 equally spaced node grids are studied to analyze the convergence of the method. The non-local parameter also varies according to the analysis and results presented in the available literature [64].
The support domain used for the mesh free implementation has rectangular shape and is centered in the Gauss points. Its dimensions are considered in the classical way [4] d s = α s d c where d c is the average nodal spacing d c = ∆x 2 + ∆y 2 and α s is a dimensionless parameter which, in this work, varies from 2 to 3. Note that, in this work, 2 × 2 Gauss integration points are used in each cell of the background integration mesh.
Results listed in Tables 3-6 compare the available analytical solutions [64] with the present ones in terms of percentage error: where w e is the exact solution taken from the aforementioned references.  The provided comparison is performed for different boundary conditions, number of collocation nodes and non-local parameter values.
As mentioned, in Mesh Free methods the background mesh and nodal distribution which represents the domain are two independent objects. However, in this work, the analysis is performed considering the collocation nodes to be coincident with nodes of the background mesh. This choice influences the convergence analysis in terms of mesh refinement. A small number of nodes implies a coarser background mesh while high number of nodes implies a finer mesh.
The α C and q coefficients characterizing the MQ radial basis functions, as well as the non-dimensional α s support domain parameter, vary as the number of nodes changes. On the contrary, they have been kept constant as the boundary conditions change. This choice, although not very convenient, is made to show how the accuracy of the results changes based only on the number of nodes used to represent the domain and on how refined the background mesh is. The parameters are changed to provide the most suitable set for the specific nodal distribution, so that the analysis is performed under some of the best achievable conditions in each case. It is true that an optimal combination of shape parameters may be found for all the different nodal distributions considered, but its application on the different nodal densities will introduce new sources of error. Therefore, in this work, we only focus on the influence that the nodal density and the integration have on the accuracy of the results, while performing the analysis with the most convenient set of parameters in each case.
In addition, a visualization of the convergence trend is shown in Figure 5 as well as a visualization of the deformed configuration of the nanoplates here analyzed shown in Figure 6.
Observing the results, it should be remarked that while in the case of the non-local parameter = 0 the results converge monotonically in most of the cases, this trend is no longer followed as the value assumed by changes. To explain this behavior, it is necessary to acknowledge the relative influence that both the nodal coordinates and the nodal density have on the numerical results.
In a case as the one under analysis, of a plate with fixed geometry, it comes natural to understand that representing the problem domain with a small or large number of nodes influences the accuracy of the results. In mesh free methods, this numerical phenomenon occurs because the nodal spacing directly enters the calculation of derivatives of any order. As explained in previous sections, the strain gradient theory requires the computation of high order derivatives of the RBF. This implies that the spacial coordinates of the nodes are multiplied several times, once per each derivative order, for each other. For small nodal densities, meaning less nodes, these coordinates have larger values which implies that larger numbers appear in the shape functions. However, as the nodal density increases, the coordinates of the nodes get smaller and smaller. When these numbers enter the shape functions, they yields larger errors. The balancing of these two aspects, nodal density and magnitude of the nodal coordinates, may cause the alteration in the trend exhibited by the convergence.
Overall, results of the analysis are quite accurate, never reaching 2% error value for the highest nodal density. As expected, the errors increase as the value of the non-local parameter increases, for the reason mentioned above.

Composite Plates
A similar analysis is also performed on squared cross-ply laminates.
Simply-supported (SSSS) laminates with lamination schemes 0, (0/90), (0/90) 2 and (0/90) 4 , subjected to a sinusoidal load are analyzed. The material property used are taken as E 1 /E 2 = 25, G 12 = G 13 = 0.5E 2 , G 23 = 0.2E 2 . Moreover, the thickness is given by h = a/100. As in the previous section, the results are compared in terms of non-dimensional mid-deflection:w where q s is the magnitude of the sinusoidal load, taken as 1. The analysis is performed using 11 × 11 regularly distributed nodes to represent the domain and 2 × 2 Gauss points to perform the numerical integration. The results of the analysis, compared in terms of percentage error as shown in Equation (13), are shown in Table 7. In this case, since the nodal density is kept constant, there is no need to change the values of the non-dimensional parameters.
Even in this case, the values of the percentage error remain contained, never even reaching 4%.

Conclusions
In this work, strain gradient nanoplates, both isotropic and laminated, have been analyzed by means of the Radial Point Interpolation Method. The aim was to apply a RPIM formulation to thin plates modeled via strain gradient theory. The paper provides detailed theoretical notions on mesh free methods in general as well as accurate explanation on the RPIM theoretical and numerical implementation. Moreover, theoretical notions are here reported in both explicit and matrix formulation.
This work proves the validity of the RPIM for problems with higher order of derivatives involved. Isotropic second order strain gradient Kirchhoff nanoplates with various boundary conditions are first analyzed. More specifically, SSSS, CCCC, SCSC, SFSF and SCSF boundary conditions are accounted for. Numerical convergence with the analytical results achieved in recent literature was studied. Results variation as the non-local parameter , introduced by the strain gradient theory, varies has also been studied. In a similar way, cross-ply composite plates subjected to a sinusoidal load are also analyzed. In this case, 11 × 11 nodes are used to represent the domain and different lamination sequences are considered. The paper provides a detailed explanation of the RPIM method theoretical and numerical implementation as well as theoretical notions in both explicit and matrix form. This work proves the validity of the RPIM for problems with higher order of derivatives involved.