Finite Element Analyses of the Modiﬁed Strain Gradient Theory Based Kirchhoff Microplates

: In this contribution, the variational problem for the Kirchhoff plate based on the modiﬁed strain gradient theory (MSGT) is derived, and the Euler-Lagrange equations governing the equation of motion are obtained. The Galerkin-type weak form, upon which the ﬁnite element method is constructed, is derived from the variational problem. The shape functions which satisfy the governing homogeneous partial differential equation are derived as extensions of Adini-Clough-Melosh (ACM) and Bogner-Fox-Schmit (BFS) plate element formulations by introducing additional curvature degrees of freedom (DOF) on each node. Based on the proposed set of shape functions, 20-, 24-, 28- and 32-DOF modiﬁed strain gradient theory-based higher-order Kirchhoff microplate element are proposed. The performance of the elements are demonstrated in terms of various tests and representative boundary value problems. Length scale parameters for gold are also proposed based on experiments reported in literature.

The mechanical behavior of MEMS and NEMS devices diverge from that of macrodevices as their sizes get smaller. This phenomenon is known as size effect or scale effect, resulting in stiffening in mechanical response as the structural length-scale tends to the material length-scale describing the microstructure. The theoretical modelling of size effects in terms of nonlocal continuum models dates back to the early works of [11][12][13][14][15][16][17][18][19][20][21][22]. The concepts of couple stresses, micropolar elasticity, micromorphic theory, nonlocal elasticity originate from these pioneering works. Within this context, the term micromorphic theory is coined for the first time by Eringen and Suhubi [23,24]. A more specific form, the Strain Gradient Theories (SGT) originate from the works of Mindlin et al. [13][14][15]20], Toupin [22] and Koiter [19]. These are in general called the higher-order theories in which the internal potential of a solid body depends only on strains but also higher-order strains as strain gradients. Fleck and Hutchinson [25][26][27] extended SGT, pioneering the framework of modern higherorder theories. The Modified Strain Gradient Theory (MSGT) of Lam et al. [28] decreased the number of length-scale parameters from five to three in SGT using new equilibrium conditions. Another variation with one length-scale parameter is also used as the Modified Couple Stress Theory (MCST) of Yang et al. [29]. Although many studies exist regarding application and implementation of higher-order theories to microbeams, microplates did not receive much attention so far [30][31][32]. Reddy [33] proposed a nonlocal plate formulation extending the nonlocal continuum theory of Eringen [34]. Phadikar and Pradhan developed finite element formulations of Euler-Bernoulli beams and Kirchhoff plates for nonlocal elasticiy using one nonlocal parameter [35]. Wang et al. [36] developed a model for SGT for Kirchhoff microplates and solved various bending problems for specific boundary conditions [37]. SGT has also been implemented in Kirchhoff plates via isogeometric analyses successfully [38][39][40]. Arefi and Zenkour analyzed magneto-electro-thermo-mechanical bending and free vibration behavior of a sandwich microbeam [41,42] and microplate using SGT [43]. They also analyzed a three-layered microbeam using SGT and three-unknown shear and normal deformation theory [44]. Sohby and Zenkour [45] analyzed hygrothermal behavior of exponentially graded microplates on elastic foundations with MCST, whereas nonlocal SGT is also utilized [46]. Movassagh and Mahmoodi [47] derived the formulation for application of Kirchhoff plate theory to MSGT and further developed a Kirchhoff plate model in MSGT using extended Kantarovich method. Sahmani and Ansari analyzed FGM's based on shear deformation plate theory and MSGT [48]. Li et al. studied the bending behavior of bi-layered Kirchhoff microplates based on SGT [49].
The above-mentioned contributions provided great insight towards the essense of size-dependent theories, and are oriented towards analytical modelling applicable for either ideal geometries or ideal boundary conditions that would be applicable to extremely limited conditions. However, stable and convergent finite element implementations allow the analysis and design of MEMS switches of various geometries and boundary conditions. Recently, there is considerable effort regarding the development of finite element formulations for size-dependent theories using Kirchhoff plates. Babu and Patel [50,51] developed a quadrilateral Kirchhoff plate element for second-order strain gradient theory proposed by Ru and Aifantis [52,53]. Beheshti [54] used the same formulation to develop four types of rectangular elements with different nodal degree of freedoms. Bacciocchi et al. used nonlocal strain gradient theory to formulate a finite element method and used it for laminated nanoplates in bending [55] and in hygro-thermal environment [56].
The main objective of this contribution is to propose novel finite element formulations for Kirchhoff microplates in a specific branch of strain gradient theories, i.e., MSGT proposed by Lam et al. [28], which contains three length scale parameters for modelling size effects. The key features and contributions of this study are; (i) to introduce a variational formulation for microplates based on the Kirchhoff kinematics and the modified strain gradient theory, (ii) to formulate 20-and 24-degree of freedom (DOF) rectangular micro-plate element formulations from serendipity-type higher-order basis functions enriched by sinh(x, y) and cosh(x, y) terms as extensions of the Adini-Clough-Melosh [57,58] element and Bogner-Fox-Schmit [59] element respectively, (iii) to further formulate 28-and 32-degree of freedom (DOF) rectangular micro-plate element formulations from 20-and 24-degree of freedom (DOF) rectangular microplate element formulations respectively, (iv) to assess the convergence of the derived element formulations and elaborate on the continuity requirements, (v) investigate the performance of the element through realistic MEMS switch geometries, and (vi) to contribute to the analysis of MEMS devices with various boundary conditions, particularly in the domain of RF-MEMS made of gold by proposing length scale parameters based on the proposed formulations.
Particularly, we propose a set of shape functions that satisfies the exact solution to the homogeneous PDE governing the MSGT based higher-order Kirchhoff plate theory. To this end, higher-order Kirchhoff plate finite elements for bending are developed based on Modified Strain Gradient Theory (MSGT). Then the weak form, from which the finite element formulations are constructed, is derived from variational principle or the principle of minimum potential energy. Using these elements, length scale parameters of gold are identified based on experiments reported in literature [60]. As a result, considerable size effect for gold microplate structures is revealed.
The proposed set of shape functions for the 20-DOF element naturally depends on a twenty-term basis function consisting of twelve monomials applicable to conventional rectangular Kirchhoff finite elements, and eight additional terms based on hyperbolic sine and hyperbolic cosine terms for the higher-order portion of the problem. The said twelve monomials belong to the simplest classical rectangular element which has twelve degrees of freedom, as known as Adini-Clough-Melosh (ACM) element [57,58]. This element does not pass the patch test but is convergent when used as rectangular finite elements [61]. For 24-DOF element, sixteen monomials applicable to a different classical element are used. This element has sixteen degrees of freedom and is known as Bogner-Fox-Schmit (BFS) element [59]. The additional eight terms are the same hyperbolic terms discussed above for 20-DOF element. The remaning 28-and 32-DOF elements are extensions of the 20-and 24-DOF elements -each with additional eight DOFs.
In order to assess the performance of the proposed higher-order finite element formulations, realistic NEMS switch geometries from literature are analyzed and the results are compared to those obtained from classical Kirchhoff plate theory. Yet, the dependence of the results on the aspect ratio is also investigated for the quadrilateral element formulation for the classical and higher-order Kirchhoff plates. It must be emphasized that the higherorder finite elements developed within this study are applicable to rectangular elements, as they fail to pass the patch test. However, this does not reduce the implementation of these elements in the community, given most of the MEMS and NEMS devices can be modelled by rectangular elements [6,[62][63][64]. As a novel aspect, it is shown that 20-DOF rectangular elements can be used to model these types of MEMS structures. Chamfers at the edges and corners can also be found in some MEMS and NEMS devices [2,[65][66][67]. Circular shapes are even possible in several applications such as micropumps and pressure sensors [68][69][70]. Again, as a novel aspect, it is shown that 24-DOF rectangular elements can be used in rectangular-triangular meshes, complementing with the readily available Bell elements [71]. With this method, any planar finite element domain of any shape can be modelled with the said rectangular and triangular elements [72].
This manuscript is organized as follows: In Section 2, the variational principle governing the higher-order Kirchhoff plate theory-based on the modified strain gradient theory is described and the respective Euler-Lagrange equations are derived. Section 3 is devoted to the novel finite element formulation of microplates based on the modified strain gradient theory of nonlocal elasticity. In Section 4, the performance of the proposed element is assessed through representative benchmark numerical examples. The paper ends with concluding remarks in Section 5.

Theory: Local and Nonlocal
where is the strain tensor. Herein, λ and µ are the first Lamé constant and the shear modulus respectively [32]. The nonlocal elasticity is implemented by the addition of higher-order strains into the free energy function in the sense of Mindlin [15] ψ =ψ(ε, η) where η = ∇∇u (2) where is the second gradient of the displacement field where u is the displacement field. For isotropic materials, the free energy function can be written in an additive format and where are the local and nonlocal parts of the free energy function, respectively, see references [20,26]. In the most general case, five additional material parameters a i exist. From (3) 1 , the first variation of the free energy can be derived as where are the second order and third order stress tensors as work conjugates of the second order strain tensor ε and the third order strain tensor η, respectively.

The Modified Strain Gradient Theory
Lam et al. [28] proposed a more specific form of the ansatz (2) and (5) where σ, p, τ, m are the local stress tensor, pressure gradient vector, the double stress tensor, and the couple stress tensor, as the work conjugates of the strain tensor ε, dilatation gradient ∇ = ∇ tr ε, deviatoric stretch gradient tensor η 1 , and curvature tensor χ such that T + 2(1 ⊗ tr(∇ε)) 13 T + (1 ⊗ ∇ ) 23 T + 2(1 ⊗ tr(∇ε)) 23 T , and Therein, ∇ p is the permutational gradient defined as ∇ p ε = ε jk,i + ε ki,j + ε ij,k . The higher-order strain metrics are the deviatoric stretch gradient tensor η 1 in addition to the rotation gradient tensor or the curvature tensor χ. ∇ is the dilatation gradient vector. The total internal energy for the modified strain gradient theory of Lam et. al [28] then takes the following form The Euler-Lagrange equations of the minimization principle then reads Herein, the corresponding stress measures as work conjugates of the strain measures ε, ∇ and χ respectively are where l 0 , l 1 , l 2 are the three length scale parameters of the modified strain gradient theory associated with the higher-order strains ∇ , η 1 , and χ, respectively.

Classical Kirchhoff Plate Theory
Kirchhoff plate theory is the extension of Euler-Bernoulli beam theory to two-dimensional space. Accordingly, rotations are assumed to be small and plane sections are assumed to remain plane and perpendicular to the neutral axis. Consequently, the displacement field simplifies to where the plate axes and geometry are as defined as in Figure 1a and w = w(x, y). Therein, L is the length, W is the width, and h is the thickness of the plate. The out-of-plane rotations and the curvatures can be described as The small-strain linear isotropic material response leads to the constitutive relation between the moment and the curvature Herein, D, E, ν are the isotropic plate rigidity, modulus of elasticity and Poisson's ratio, respectively. The displacement vector u consists of vertical displacements w and rotations θ and surface traction vector t consists of shear forces V and moments M as given in Figure 1b. Due to kinematic equations, shearing, twisting and membrane effects are neglected in the strain energy equation. Therefore strain energy is governed by bending action solely. For a Kirchhoff plate, the internal and the external potentials read The first variation of the total potential Π has to vanish at the equilibrium state. Hence The Euler-Lagrange equation of the minimization principle can then be derived as (18) where ∇ 2 ∇ 2 (•) is the biharmonic operator. The Dirichlet (essential) and Neumann (natural) boundary conditions are respectively,

Modified Strain Gradient Theory for Kirchhoff Plates
With the displacement field given in Equation (12) and kinematic relations given in Equation (13) valid, an additional set of higher-order out-of-plane curvatures are defined as The constitutive relations between the moment & curvature and higher-order forces & higher-order curvatures are, respectively Derivation of the above equations are given in Appendix A. For the MSGT based higher-order Kirchhoff plate, the internal potential reads The external potential reads Q is therefore the work conjugate of the curvature κ, similar to the work conjugate couples V-w and M-θ. The first variation of the internal potential leads to Similarly, the first variation of the external work potential can be derived as The principle of minimum potential energy also requires the variation of the total potential to vanish at equilibrium. The Euler-Lagrange equation of the minimization principle can be derived as follows The Dirichlet and Neumann boundary conditions are respectively,

The Microplate Finite Element Formulation
The finite element formulation for the classical plate theory and the modified strain gradient theory will be obtained using the variational formulations presented in Section 2. The proposed finite element formulation is derived by first proposing additional higherorder DOF's. Then, homogeneous solutions that satisfy Equation (18) are proposed and the shape functions are checked. From the relevant shape functions, local element stiffmess matrices are derived.

12-Dof Adini-Clough-Melosh Element for Classical Kirchhoff Plates
Let us consider a classical Kirchhoff plate element domain B e as depicted in Figure 1a. The generalized displacements at the element nodes are prescribed as 1. w = w 1 , θ x = θ x1 and θ y = θ y1 at x = 0 and y = 0 , 2. w = w 2 , θ x = θ x2 and θ y = θ y2 at x = L and y = 0 , 3. w = w 3 , θ x = θ x3 and θ y = θ y3 at x = L and y = W , 4. w = w 4 , θ x = θ x4 and θ y = θ y4 at x = 0 and y = W , see also Figure 2b. Accordingly, the vector field containing the nodal displacement vector d and the nodal force vector f read The displacement field w(x, y) within the element domain B e is interpolated as where is the row vector including the set of interpolation/shape functions. n DOF = 3 is the number of degrees of freedom (DOFs) per node and n nodes = 4 is the number of nodes per element, with N j i , i denoting the relevant DOF of the node and j denoting the relevant node. The number of DOFs per element is hence 12. The homogenous solution of the partial differential Equation (18) is w(x, y) = a 1 + a 2 x + a 3 y + a 4 x 2 + a 5 xy + a 6 y 2 + a 7 x 3 + a 8 x 2 y + a 9 xy 2 + a 10 y 3 + a 11 x 3 y + a 12 xy 3 .
The 12 DOF plate element with such shape functions is known as the ACM quadrilateral [73]. The first three of the twelve shape functions are given in Figure 3, where the remaining nine are symmetric with respect to two principal centroidal axesx andȳ around the geometric center of the element. The analytic expressions for the ACM shape functions are given in Appendix B, see also [73]. The ACM element does not satisfy the C 1 continuity requirement, and therefore it is a non-conforming element [61]. The displacement w(x, y), rotation {θ x , θ y }, and the curvature {κ x , κ y , κ xy } fields (13) can be approximated as Consequently, the variation of the field variables (34) are Incorporation of the discrete counterpart of the curvature in Equation (34) and their variations (35) into (17), we obtain where where are the element stiffnes matrix and the element nodal force vectors. For the sake of convenience the loading terms on the boundary are omitted for f e . Therein, the operator ∇ C is defined as and A refers to the standard assembly of element contributions at the local element nodes where n denotes the total number of elements. In Equation (47) the strain-displacement matrix is ∇ C N as given in Appendix D.
The global stiffness matrix, generalized nodal displacement vector and the generalized force vector assembled from local force vectors read respectively. No variation exists δd = 0 at essential boundary ∂B u where the displacements are prescribed d =d. The equilibrium is satisfied for arbitrary variation of the displacement δd leading to the set of linear algebraic equations

16-Dof Bogner-Fox-Schmit Element for Classical Kirchhoff Plates
For a classical Kirchhoff plate element the generalized displacements at the element nodes are prescribed as see also Figure 4a. Similarly, the generalized nodal force resultants are prescribed as The displacement field w(x, y) within the element domain is interpolated as where N is the row vector including the set of interpolation/shape functions, n DOF = 4 is the number of degrees of freedom (DOFs) per node and n nodes = 4 is the number of nodes per element. For N j i , i denotes the relevant DOF of the node and j denotes the relevant node. The number of DOFs per element is hence 16. The homogenous solution of the partial differential biharmonic plate equation is w(x, y) = a 1 + a 2 x + a 3 y + a 4 x 2 + a 5 xy + a 6 y 2 + a 7 x 3 + a 8 x 2 y + a 9 xy 2 + a 10 y 3 + a 11 x 3 y + a 12 x 2 y 2 + a 13 xy 3 + a 14 x 2 y 3 + a 15 x 3 y 2 + a 16 x 3 y 3 .
Incorporation of these into the total functional, we obtain where are the element stiffness matrix and the element nodal force vectors. The global stiffness matrix, generalized nodal displacement vector and the generalized force vector are assembled from local force vectors in a similar manner as described in Section 3.1.

New 20-Dof Finite Element Formulation for Msgt-Based Kirchhoff Microplates
Let us now consider a MSGT-based Kirchhoff plate element domain B e . The generalized displacements at the element nodes are prescribed as see also Figure 5a. Similarly, the generalized nodal force resultants are prescribed as see Figure 5b for relevant nomenclature. The element nodal displacement vector d and the element nodal force vector f of the higher-order Kirchhoff plate are, respectively, where the nodal bi-axial curvatures and higher-order moments as their work conjugates enter the formulation as additional nodal field variables. The displacement field w(x, y) within the element domain B e is interpolated as is the row vector including the set of interpolation/shape functions with number of DOFs per node is n DOF = 5 and number of nodes is n nodes = 4. We propose a homogenous solution of the partial differential Equation (27) in the form w(x, y)= a 1 + a 2 x + a 3 y + a 4 x 2 + a 5 xy + a 6 y 2 + a 7 x 3 + a 8 x 2 y + a 9 xy 2 + a 10 y 3 + a 11 x 3 y + a 12 xy 3 + a 13 sinh(Ax) + a 14 cosh(Ax) + a 15 sinh(By) + a 16 cosh(By) + a 17 sinh(Ax)y + a 18 cosh(Ax)y The use of hyperbolic sine and hyperbolic cosine terms are in fact motivated by the nature of the plate equation with fourth and sixth order terms (27). The choice of the additional degrees of freedom is also based on the study of Kahrobaiyan et al. [74]. This solution naturally extends the MSGT based Euler-Bernoulli beam solution to Kirchhoff plate solution, see Reference [32].
We propose 20 shape functions derived from the basis function that satisfies the homogeneous solution (65). The first five of these shape functions (i.e., those for the first node) are shown in Figure 6. The remaining fifteen are symmetric with respect to two centroidal principal axesx andȳ around the geometric center of the element. The analytic expressions for the shape functions N  The higher-order derivatives appearing in Equation (A12) in Appendix A are interpolated in the element domain B e as follows Consequently, the variation fields for the relevant terms in Equation (25) read Incorporation of the discrete counterpart of the curvature in Equation (34), higherorder curvature (53) and their variations (35) and (54) along with (25) and (26) into minimization principle, we obtain where are the element stiffnes matrix and the element nodal force vector (Figure 5b), respectively. Therein, the higher-order in-plane Laplacian operator ∇ H is There are two strain-displacement matrices in Equation (56) that can be expressed as ∇ C N and ∇ H N as given in Appendix D.
By changing of variables (x, y) with (ξ 1 , ξ 2 ), the stiffness matrix can be written as Therein, the Jacobian or the transformation between the natural parameter space ξ and physical space x is Making use of two-point Gaussian quadrature, the continuous integral (58) for the element stiffness matrix can be replaced with the discrete representation where ξ 1 (i), ξ 2 (j) and Ω(ξ 1 (i), ξ 2 (j)) are the Gaussian quadrature points and the weight factors, respectively. The 28-DOF finite element formulation variant derived from ACM and 20-DOF finite elements is given in Appendix E.

New 24-Dof Finite Element Formulation for Msgt-Based Kirchhoff Microplates
We develop the 24-DOF element from the 16-DOF BFS element in exactly the same manner the 20-DOF element is developed from the 12-DOF ACM element in Section 3.3. Herein the generalized displacements at the element nodes are prescribed as , κ yy = κ yy1 and κ xy = κ xy1 , at at node 1 , 2. w = w 2 , θ x = θ x2 , θ y = θ y2 , κ xx = κ xx2 , κ yy = κ yy2 and κ xy = κ xy2 , at at node 2 , 3. w = w 3 , θ x = θ x3 , θ y = θ y3 , κ xx = κ xx3 , κ yy = κ yy3 and κ xy = κ xy3 , at at node 2 , 4. w = w 4 , θ x = θ x4 , θ y = θ y4 , κ xx = κ xx4 , κ yy = κ yy4 and κ xy = κ xy4 , at at node 4 , see Figure 7a. Similarly, the generalized nodal force resultants are prescribed as see Figure 7b. The element nodal displacement vector d and the element nodal force vector f are, respectively, The displacement field w(x, y) within the element domain is then interpolated as with number of DOFs per node is n DOF = 6 and number of nodes is n nodes = 4. We propose a homogenous solution of the partial differential Equation (27) in the form w(x, y)= a 1 + a 2 x + a 3 y + a 4 x 2 + a 5 xy + a 6 y 2 + a 7 x 3 + a 8 x 2 y + a 9 xy 2 + a 10 y 3 + a 11 x 3 y + a 12 x 2 y 2 With similar hyperbolic sine and cosine terms, the equation above satisfies the homogeneous solution for Equation (27). Hence we also propose twenty-four shape functions as given in Equation (64). The first six of these shape functions (i.e., those for the first node) are shown in Figure 8. The remaining eighteen are symmetric with respect to two centroidal principal axesx andȳ around the geometric center of the element. The shape functions N j 1 , N j 2 , N j 3 , and N j 6 recover their classical counterparts in BFS element for j'th node whereas N j 4 and N j 5 vanish as length scale parameters tend to zero. By changing of variables and the use of two-point Gaussian quadrature, the discrete representation for the element stiffness matrix k e similar to Section 3.3 is similar to the formulation in Section 3.3. The 32-DOF finite element formulation variant derived from BFS and 24-DOF finite elements is given in Appendix F.

Conformity
Since the energy Equation (A12) in Appendix A involves terms with third derivatives of displacement w, C 2 -continuity is required for conformity. Herein, 20-DOF element is demonstrated and 24-, 28-, and 32-DOF versions behave similarly. Let us consider a rectangular element as depicted in Figure 9. For an element boundary AB, the C 2 -continuity would require displacements, rotations and curvatures to be uniquely defined in terms of the nodal degrees of freedoms, i.e., w, θ x , θ y , κ x and κ y at points A and B respectively (ten nodal variables). Similarly for an element boundary DA given in the same figure, these need to be defined by the same nodal variables at D and A. For the boundary along the element edge AD, where y is constant and zero with 18 unkown constants that can not be defined by invoking 10 nodal displacement boundary conditions. Hence, just like its conventional counterpart [61], it is not possible to specify a polynomial set for the shape functions that ensure compatibility. The proposed elements are non-conforming elements whereas the convergence and numerical perfomance of the elements are investigated in Section 4. Cirak et al. [75,76] have overcome this problem in the context of classical Kirchfoff plates and shells by introducing a concept of triangular elements with subdivision surfaces where nonlocal interpolations are imposed on the subdivison elements, fully conforming to C 1 continuity requirements. The current contribution, however, shall be considered as a first step towards incorporating MSGT into Kirchhoff plates in the context of finite element method. In the subsequent part of the work, we also develop several elements with more degrees of freedom, particularly from the conforming Bogner-Fox-Schmit (BFS) element [59], in addition to adding more higher-order terms.

Representative Numerical Examples
In this section, the performances of the proposed 20-, 24-, 28-, and 32-DOF MSGT-based Kirchhoff plate elements are assessed under various boundary and loading conditions for different geometries. To this end, the finite element method outlined in Section 3.3 is implemented into our in-house finite element code that is developed via Matlab. First, the performance comparison is made between all the proposed elements. Then, we find the length scale parameters of gold based on the experiments of Espinosa et al. [60].
Four benchmark examples concerned with the assessment of the performance of the proposed element formulation with respect to mesh irregularity due to element size, mesh convergence, and response to various displacement boundary conditions are given. We then study the convergence of the numerical solutions to the approximate analytical solutions for a square and rectangular plate subjected to evenly distributed surface loads per the study of Movassagh and Moahmoodi [47]. Finally, several representative boundary value problems involving realistic MEMS switch geometries are investigated, and the results are compared with those obtained from the classical Kirchhoff plate theory. For all the problems, the number of elements and mesh size are chosen based on the studies given in Section 4.3.3.

Comparison
A basic comparison case study is performed by using a fixed-fixed plate with dimensions 20 µm × 5 µm × 1 µm [32]. A concentrated midpoint load of F z = 1 mN is applied as shown in Figure 10. Therein, square elements are used with a mesh of 32 × 8 elements, i.e., with a mesh density of 1.6 elements/µm at the edges. The resulting deflection field is shown in Figure 11. The deflection profiles of the principal centroidal axes of the plates are depicted in Figure 12a,b respectively. The difference between the models constructed with 20-, 24-, 28-, and 32-DOF elements is insignificant -i.e., in the order of magnitude of 10 −5 µm and smaller than 0.2% of the total tip deflection-, and the deflection profile looks almost the same in Figure 12a. In Figure 12b the difference can be seen, along with a slight saddle effect that is expected due to Poisson effect [37].
Along with these results and similar results found for the forthcoming sections for the analysis of length scale parameters and assessment of element performance, only the results for the 20-DOF element are included for Sections 4.2-4.4.

Length Scale Parameters for Gold Microplates
Experiments in the study of Espinosa et al. [60] are modelled by microplates. The specimen dimensions and force-displacement values are given in Table 1, and the boundary conditions are shown in Figure 13. Force-displacement behavior is assumed as linear elastic and elastic parameters E and ν for Specimens 1 and 2 are assumed to be the same. The methodology is based on our previous study [32], i.e., for quantification of elastic modulus E and length scale parameter l, an error parameter Err is defined as the L 2 -norm of the residual vector, for the quantification of the best fit at which Err is minimum.
Herein w sim 1 and w sim 2 are the midpoint deflections predicted by higher order theories, w exp 1 and w exp 2 are the actual midpoint deflections from experiments for specimens 1 and 2. Err is evaluated for different values of E and l in order to find the minimum error. It is seen that the error function is minimum along the dark blue region in Figure 14. For bulk elastic modulus of gold i.e., E = 80 GPa, l 0 is found as 3.60 µm. Table 1. Experiments taken from [60]. W, t, L denote the width, thickness and length of the specimen respectively. F z and w ma are the vertical force applied at the midpoint and the actual displacement of the midpoint respectively. (*gage dimension).   Throughout the following numerical examples demonstrating element performance, these material parameters specific to gold are used. Those of epoxy are taken from the pioneering work of Lam et al. [28]. These are summarized in Table 2.

Assessment of Element Performance
The proposed MSGT-based Kirchhoff microplate elements can be used for rectangular elements similar to their classical counterparts. The results obtained with the developed MSGT-based Kirchhoff microplate are exactly the same with those with the ACM element, with vanishing length scale parameters. Although the ACM element passes the patch test, the usage area should be confined to rectangular meshes [61,77]. The proposed formulation shows similar performance to ACM element under distorted element geometries. In this regard, we confine ourselves to the investigation of the convergence behaviour upon mesh refinement and the element performance to irregular rectangular meshes. To this end, the convergence of the displacement field for a square block subjected to a point load is investigated for various boundary conditions, Additionally, the sensitivity of the displacement field to mesh irregularity under prescribed displacement/rotation field is investigated.

Microplate Response to Point Load
A 6 µm × 6 µm × 1 µm microplate is subjected to a concentrated load of 1 mN applied at the centroid. The plate is double-cantilevered as depicted in Figure 15a. The problem is investigated with several meshes as shown in Figure 15 in order to assess its sensitivity to mesh irregularity. The corresponding deflection profiles are depicted in Figure 16a-d. In order to assess the difference in displacement fields, the midline deflections for the relevant nodes at x = 0 and y = 0 for the regular mesh and the irregular mesh given in Figure 16b are given in Figure 17a. The points correspond to the nodal displacement values along x = 0 and y = 0. The problem is also solved with ACM element and similar midline deflections are indicated in Figure 17b. Therein, curves in between nodal values are interpolated with the element shape functions. It is seen that the deflections of the nodes for the regular and irregular meshes complement each other. The difference between centroidal deflections is smaller than 1% as seen in Figure 16, with the displacement fields aligned to a reasonably acceptable extent. In fact the average difference between the deformation results at the nodes is less than 1%, meaning it performs slightly better than ACM element does for midline deflections in both directions. Although it is suggested that the element size variation such as aspect ratio change to be minimum as a best practice, varying aspect ratios do not yield erroneous results at least up to some extent, regarding displacement results.

Microplate Response to Displacement and Rotation
In order to check the integrity of the formulation and consistency of the numerical implementation in xand y-directions, the 6 µm × 6 µm × 1 µm microplate is subjected to a unit displacement and rotation at two perpendicular edges, respectively, see Figure 18 (first column). Therein, the left edge is fixed along y-axis and the right edge is displaced 1 µm in z-direction, see Figure 18a. Then, the same plate is subjected to a unit rotation (1 rad) about y-axis, see Figure 18b. The same procedure is repeated for the perpendicular direction in Figure 18c,d, respectively. The simulation is first carried out with 4 × 4 regular mesh (second column) and for an irregular mesh (third column). Although the proposed element is shown to be non-conforming in Section 3.5, the displacement fields obtained from the regular and irregular meshes are nearly identical. The average difference between the deformation results at the nodes is less than 1%, hence the performance is deemed to be satisfactory. Figure 18. Geometry and boundary conditions (left) and corresponding displacement fields for the aspect ratio test considering response to prescribed displacements, for regular mesh (middle) and irregular mesh (right). The thickness of the plates is 1 µm for (a-d).

Mesh-Refinement and Convergence
The 6 µm × 6 µm × 1 µm is now subjected to a concentrated midpoint load of 1 mN with both classical ACM elements and the proposed higher-order microplate element. The the boundary conditions are specified as (i) double-cantilevered (two opposite edges clamped, two opposite edges free) and (ii) wholly-cantilevered (all four sides clamped), see Figure 19 (first column). The displacement profile for these boundary conditions on a 4 × 4 mesh are given in Figure 19 (second column). The midpoint deflections versus element per edge results are depicted in Figure 20a-d, for the proposed element formulation and the ACM element, respectively. The proposed element formulation is converging slightly faster than the classical counterpart. The convergence behaviour of the proposed element upon mesh-refinement is satisfactory for rectangular meshes given that the differ-ences between the maximum displacement for consequent mesh refinements are less than 0.1%.

Square Microplate Subjected to Different Boundary Conditions
A 20 µm × 20 µm × 1 µm microplate subjected to a point load is analyzed under various boundary conditions with the classical ACM plate element and the proposed MSGT-based microplate element formulations, respectively. Four different boundary conditions are considered, see Figure 21: (a) CFCF, (b) CFFF, (c) CCFC and (d) CCFC. The boundary conditions are abbreviated by "C" for clamped ends and by "F" for free ends. In the example (a) a centroidal, (b) midpoint of the free edge, (c) centroidal and (d) midpoint of the free edge, are subjected to point load, respectively. The domain is discretized with 20 × 20 elements. The results obtained from the classical ACM plate finite element and the porposed finite element formulation are visualized in Figure 21. In order to study the difference between the deformation patterns obtained from the classical and MSGT-based Kirchhoff theory, the contour plots are normalized with respect to the peak values. The numerical results show that, not only the maximum deflections but also the deformation patterns change significantly by considering the size effect in terms of the modified strain gradient theory.
It is seen in Figure 22 that the MSGT approaches asymptotically to classical theory as the thickness, which is related to structural length scale, increases. With the assumed length-width-thickness (aspect) ratio of 20:20:1 for gold microstructures, the error of using classical theory with macroscopic material parameters vary with boundary conditions. It still seems to be more than 10% if thickness is reduced below ca. 40 µm at best case. This error increases with decreasing plate thickness -it becomes more than 25% if the relevant thickness is taken smaller than ca. 25 µm at best case. It is seen that for the particular cases with 20 µm × 20 µm × 1 µm microplate, the size effect yields deflections in the range of 0.3-0.5% of the deflections that are found with classical theory.

Benchmark Example: Rectangular Microplates Subjected to Evenly Distributed Load
Three benchmark examples, which are analytically solved by Movassagh and Mahmoodi [47] using extended Kantarovich method (EKM), are numerically reproduced with the proposed higher-order microplate elements. The boundary conditions and the dimensions in terms of the length scale parameters l = l 0 = l 1 = l 2 are depicted in Figure 23a-c. A distributed load of 1 kN/m 2 is applied in all three cases. Therein, the material parameters specific to epoxy are used [28], see Table 2. The normalized midpoint deflections (w/l) versus number of elements per length results are depicted in Figure 24a-c. The numerical results recover the analytical solution of Movassagh and Mahmoodi upon mesh refinement for all three cases. The difference between the results for consequent three refinements are less than 0.05%.    [47] replicated numerically with the proposed higher order microplate elements. The dimensions are in terms of the length scale parameters l = l 0 = l 1 = l 2 , with thickness h = l in all cases. A distributed load of 1 kN/m 2 is applied for all cases as in the said study. Therein, the material parameters specific to epoxy are used [28], see Table 2.

Analysis of Realistic Mems Switches with 20-Dof Elements
In this example, we consider three real MEMS switch structures from Patel and Rebeiz [66] and Stefanini et al. [65], respectively. The chamfers of the plates and tethers are ignored in this section. The geometry and the boundary conditions are depicted in Figure 25. Stefanini et al. [65] discusses an actuation electrode and the corresponding MEMS structure (see Figure 25a) to transfer the majority of the electrostatic force to the contact force, i.e., where F c is the contact force, F r is the release force, and F e is the electrostatic force. Therein, F c = 34.7 µN, F r = 15.5 µN, and hence the electrostatic force is found as F e = 77.2 µN. This force is equally distributed to the nodes that are electrostatically actuated, see Figure 25a.
The deflection of the plate should be equal to the clearance of 0.3 µm for contact condition. The electrostatic loads are applied to the structures and the deflected shapes, which are obtained from classical and the MSGT-based Kirchhoff plate theory, are shown in Figure 26.  Figure 25a,b, can now be designed and analyzed more effectively making use of the new MSGT plate elements. The MEMS community traditionally use higher elasticity parameters such as Young's modulus µ and shear modulus µ. This choice, for uniform thickness and under pure bending deformations leads to satisfactory results in line with the modified couple stress theory (MCST). This is mainly due to the fact that, MCST, when applied to Kirchhoff plate theory, leads to the same differential equation as the classical counterpart, where the nonlocal effects are merely reflected to the material parameters. In order to assess the difference between two theories, we depict the normalized tip deflections corresponding to each switch structure in Figure 28. The normalized tip deflections of the classical and the MCST-based Kirchhoff theory will lead to equivalent result. The idea here is to show, how the deflection pattern changes as we switch to the MSGT-based Kircchhoff plate theory from the classical counterpart or the MCST-based Kirchhoff plate theory. As seen from, Figure 28, for cantilever geometries, where the bending governs the deformation, normalized results overlap. However, for the second geometry, where highly complex local and nonlocal deformations exist due to the complex geometry and boundary conditions, the normalized deflection patterns are quite dissimilar, revealing the necessity for the MSGT-based Kirchhoff plate theory.

Analysis of Realistic Mems Switches with 24-Dof Elements
We consider the two real MEMS switch structures from Stefanini et al. [65] and Patel and Rebeiz [66] that had been investigated in the previous section. In this section, chamfers are added, making the geometric models more realistic considering their original shapes in respective studies. The geometry and the boundary conditions are depicted in Figure 29.
We solve the problems by making use of the rectangular mesh (figures in the middle in Figure 25) as well as with a mesh consisting of 18-DOF Bell triangular and 24-DOF rectangular elements (figures on the right in Figure 25). For the triangular-rectangular mesh, very few 18-DOF Bell elements are used as seen in the relevant figure. Herein, in determining the higher-order stiffness matrix for the 18-DOF Bell element, Equation (58) is used that incorporate the higher-order d 1 -d 6 terms and hence the length scale parameters l 0 , l 1 , and l 2 for the readily available classical Bell triangular element with the explicitly defined shape functions [78]. The length scale parameters are chosen as found in Section 4.2. The mix of rectangular-triangular elements yield the deflection fields as given in Figure 30. The deflection profile of the end tips are also given in Figure 31 for the rectangular mesh in the first study and the mix of rectangular-triangular elements. The average difference of all data points when compared with the purely rectangular mesh is ca. 0.2% for the case in Stefanini et al. [65] and ca. 0.02% for the case in Patel and Rebeiz [66]. It can then be concluded that the use of triangular Bell elements especially in the truncated parts of the microplates does not result in a considerable difference in the analyses when used with the same length scale parameters.

Conclusions
In this paper, we presented novel MSGT-based element formulations for Kirchhoff plate theory-based on the variational principle dictating the minimum potential energy. Finite element discretization is implemented based on this framework and relevant codes are developed. Herein, 2 × 2 Gaussian quadrature is implemented. Despite being nonconforming elements and only applicable to rectangular meshes similar to the original classical plate elements, it is demonstrated that the new elements yield expected results in various numerical examples. It is explicitly demonstrated that the convergence and aspect ratio test performances of the proposed elements are equivalent the classical elements. The numerical problems investigated clearly show that using the classical plate theory in predicting microplate behavior results in significant errors, especially if the thickness of the structures are smaller than 40 µm.
In the future, discrete techniques, isogeometric methods, or similar will be utilized to come up with a conforming higher-order element applicable to general quadrilaterals, satisfying C 2 -continuity requirement. However, it is also worth mentioning that many, if not most, of the MEMS-NEMS plate structures can be modelled with only rectangular elements, similar to the ones modelled and analyzed in this study.
All simulations are carried out on a standard Laptop with Intel I7 processor having 8 × 2.4 GHz cores and 8 GB Ram, without any parallelization, requiring several minutes of computation time for the most demanding simulation. Computational timing benchmarks for the analyses mentioned under Section 4.1 reveal that analysis with MSGT with the newly developed elements takes up to approximately 2 min, 14 × longer than classical analyses which takes less than 10 s, ceteris paribus. This is in fact why 2 × 2 Gaussian quadrature method is adopted, which decreases the computational duration ratio up to 6 at most, when compared to the classical theory, again ceteris paribus, yielding a computational duration of less than 1 min. This computational time ratio of approximately 6 is found in other analyses within this study as well. It is noted however, most MEMS-NEMS structures can be modelled with very small number of nodes and elements than macrostructures, hence the increase in the computational time can be tolerated with the state of the art CPUs and parallelization techniques.

Conflicts of Interest:
The authors declare no conflict of interest.
The indicial expression for deviatoric stretch gradient terms (η 1 ) is: and where δ ij is the Kronecker's delta. After several steps, the terms of the nonlocal deviatoric stretch gradient tensor are derived as with The curvature or the rotation gradient terms in indicial notation are e imn u n,mj + e jmn u n,mi , where (A9) e ijk is the permutation symbol. After some manipulations, the components of the curvature tensor read Taking the variations of above expressions and multiplying with their stress conjugates as appearing in Equation (A3), the following are acquired: where B is 2D domain of the mid-plane of the plate bounded by a piecewise smooth curve Γ. Then Inserting the stress terms, stress resultant terms read Now, using the divergence theorem for Equation (A12) and minimum potential energy principle with we obtain the Euler-Lagrange equation governing the minimization principle M xx,xx + M xy,xy + M yy,yy − Q xxx,xxx − Q xxy,xxy − Q xyy,xyy − Q yyy,yyy = q(x, y) . (A28) Insertion of Equations (A20)-(A26) into (A28) leads to the Equation (27).

Appendix B. Shape Functions for Acm Element Based on Classical Kirchhoff Plate Theory
The first three elements of ACM shape functions N 1 corresponding to the first node for classical theory, which are shown in Figure 3, are as given as:

Appendix C. Shape Functions of the Proposed Element for Msgt-Based Kirchhoff Microplates
The first five elements of the shape function N for MSGT, which are shown in Figure 6, are as given below. where with the terms in the denominators defined as (A32)

Appendix D. Strain-Displacement Matrices for Msgt-Based Kirchhoff Microplates
The expression for the elements of the strain-displacement matrices ∇ C N (3 × 20 in size) and ∇ H N (4 × 20 in size) are given in this section. Null elements are also indicated. ∇ C N converges to the classical strain-displacement matrix, and ∇ H N, i.e., the higher order strain-displacement matrix converges to zero, when l 0 = l 1 = l 2 = 0.
The strain-displacement matrix similar to the classical and nonlocal counterpart are with Therein, the superscript j corresponds to the node number. Herein the number of DOFs per node is n DOF = 7. We propose a homogenous solution of the partial differential Equation (27) in the form w(x, y)= a 1 + a 2 x + a 3 y + a 4 x 2 + a 5 xy + a 6 y 2 + a 7 x 3 + a 8 x 2 y + a 9 xy 2 + a 10 y 3 + a 11 x 3 y + a 12 xy 3 +a 13 sinh(Ax) + a 14 cosh(Ax) + a 15 sinh(By) + a 16 cosh(By) + a 17 sinh(Ax)y + a 18 cosh(Ax)y +a 19 sinh(By)x + a 20 cosh(By)x + a 21 sinh(Ax)y 2 + a 22 cosh(Ax)y 2 + a 23 sinh(By)x 2 +a 24 cosh(By)x 2 + a 25 sinh(Ax)y 3 + a 26 cosh(Ax)y 3 + a 27 sinh(By)x 3 + a 28 cosh(By)x 3 . (A39) The additional shape functions to those in 20-DOF element are as depicted in Figure A2. This can also be interpreted as the additional shape functions to the first five (N 1 1 to N 1 5 ) given in Figure 8. The element stiffness matrix can be found in a similar fashion as given in Section 3.3.

Appendix F. New 32-Dof Finite Element Formulation for Msgt-Based Kirchhoff-Love Microplates
The two additional generalized displacements and generalized nodal force resultants to those of the 24-DOF element are the same as those given in Equations (A35) and (A36), see Figure A3a,b.
The element nodal displacement vector d and the element nodal force vector f are, The additional shape functions to those in 24-DOF element are as depicted in Figure A4. This can also be interpreted as the additional shape functions to the six shape functions (N 1 1 to N 1 6 ) given in Figure 8. The element stiffness matrix can be found in a similar fashion as given in Section 3.3.