# A Gradient Stable Node-Based Smoothed Discrete Shear Gap Method for Analysis of Reissner–Mindlin Plates

Authors: Yadong Xu; Guangsong Chen; Jinsong Tang
Journal: Shock and Vibration
 (2020)

---

## Abstract

In this paper, a gradient stable node-based smoothed discrete shear gap method (GS-DSG) using 3-node triangular elements is presented for Reissner–Mindlin plates in elastic-static, free vibration, and buckling analyses fields. By applying the smoothed Galerkin weak form, the discretized system equations are obtained. In order to carry out the smoothing operation and numerical integration, the smoothing domain associated with each node is defined. The modified smoothed strain with gradient information is derived from the Hu–Washizu three-field variational principle, resulting in the stabilization terms in the system equations. The stabilized discrete shear gap method is also applied to avoid transverse shear-locking problem. Several numerical examples are provided to illustrate the accuracy and effectiveness. The results demonstrate that the presented method is free of shear locking and can overcome the temporal instability issues, simultaneously obtaining excellent solutions.


---

## Body

## 1. Introduction

Thin-walled structures (shells) render a majority of engineering structures, and as one special case of shells, the plate has been widely used in mechanical, civil, marine, aerospace, and other engineering science fields. The analyses of plate structures in elastic-static, free vibration, and buckling fields stand for the key three aspects in their engineering application. There exist two first-order plate theories, namely, the Kirchhoff plate theory and the Reissner–Mindlin one. Kirchhoff theory is usually applied to thin structures with negligible shear strain, and the C1-continuous shape function is required. In view of its simplicity and efficiency, the lower-order Reissner–Mindlin plate which considered the shear effects is appealing in practical and only requires the C0-continuity shape function for both translational and rotational displacement fields. However, the shear-locking phenomenon of Reissner–Mindlin plate elements emerges when the thickness reaches the thin limit, and this is due to spurious transverse shear strains/stresses in bending. A development residing on node-based kinematics that aims at alleviation of the shear-locking effect and local improvement of stress recovery has been recently presented by Valvano et al. [1] Besides, many researchers proposed large amounts of effective elements to address this difficulty, such as the assumed natural discrete shear gap (DSG) method [2–5], strain (ANS) methods [6, 7], and also the methods in [8–16]. All the methods show excellent performance in reducing the shear-locking deficiency and increasing the solution accuracy. The DSG method, similar to the ANS method, owns a property of “glue” because of no additional collocation points and fits the combination with other novel element techniques. Marinkovic et al. [17]. applies DSG for the plate part of a flat shell element together with a strain smoothing technique implemented so as to make it independent from node numbering, and it was made available to users through ABAQUS implementation of the element.Up to now, the finite element method (FEM) still holds its place as the most widely used numerical tool to simulate different behaviors of plates [18, 19] and other structures due to its robustness, reliability, and effectiveness. Unlike traditional FEM in which element connectivity should be established to form the discretized equations, another form of numerical method development called the meshfree or meshless method [20–44] has attracted much attention. Regardless of the element connectivity among the nodes and mesh, only a set of nodes scattered in the problem domain are required. Although the meshfree indeed can overcome some drawbacks of FEM, it still cannot overcome all the deficiencies of FEM. Some key limitations are the difficulties in essential boundary condition implementation, high computational cost, and overly complex trial function construction processes. In an effort to make use of both advantages of FEM and meshless methods, Liu et al. have extended the concept of smoothing domains to formulate a family of smoothed finite element methods (SFEM) [45–47] by using the strain smoothing technique [48]. Researchers further proposed the edge-based smoothed finite element method (ES-FEM) [49] and node-based smoothed finite element method (NS-FEM) [50] based on the concept of SFEM.NS-FEM can be regarded as a modified model of FEM. It has very attractive properties and can be easily applied to tetrahedral or triangular elements without any modification of formulas and procedures. NS-FEM wins the favor recently for its prominent inherent properties [51], such as its insensitivity to element distortion and its immunity to volumetric locking. Moreover, the computation efficiency of NS-FEM has been studied in previous works using bandwidth solvers for linear elastic-statics [52, 53]. It is, however, found that the NS-FEM behaves “overly soft” resulting from correction to the “overly stiff” behavior of the compatible FEM. Such an “overly soft” behavior leads to the so-called temporal instability [54]. In addition, spatial instability, another kind of instability, is also a common problem in node integration. The spatial instability can be successfully eliminated by smoothing operation. Temporal instability can be reflected in the modal frequency analysis of structures, which often leads to spurious nonzero energy modes in free vibration analyses and is still a problem to be solved. In [55], Beissel and Belytschko pointed out that by adding a stabilization term that contains the square of the residual of the equilibrium equation to the potential energy functional, the problem of nodal integration which suffers from spurious single modes due to underintegration of the weak form can be solved. Chai et al. [56] also proposed a stable NS-FEM to cure the “overly soft” of NS-FEM for analysis of underwater acoustic scattering problems. To overcome temporal instability of nodal integration in metal-forming simulations, Bonet and Kulasegaram [57] presented a least-square stabilization procedure based on these previous works, and Zhang and Liu [53] further developed a stabilization procedure for NS-FEM and then provided a recommended range for the stabilization parameter. By expanding the Taylor series of the function of the displacement field [58], it can be used to reduce the instability in direct node integration. However, since the high-order derivatives appear in underlying formulations, the computational cost will increase. Other forms of stabilization consisting of the Taylor expansion and displacement smoothing have been proposed [59], wherein the nodal integration technique is directly applied to obtain stable solutions. Puso et al. [60] developed a nodal integration technique by adding integration points, of which the effectiveness has been proved for both small and large deformation problems. Feng et al. [61] proposed a stable nodal integration method with strain gradient for dynamic analyses of solid structures based on NS-FEM. The proposed method can achieve appropriate system stiffness in train energy between FEM and NS-FEM solutions and indeed provide temporally stable results. There still exist a variety of gradient term constructions available for different cases [37, 40, 62–67].In this work, a gradient stable node-based smoothed discrete shear gap method (GS-DSG) using 3-node triangular elements is formulated for elastic-static, free vibration, and buckling analyses of Reissner–Mindlin plates. In order to overcome the temporal instability problem encountered in the nodal integration process, the smoothed Galerkin weak form is applied by using the strain smoothing technique with gradient information, which is derived from the Hu–Washizu three-field variation principle. The stabilized discrete shear gap method is also incorporated into the presented method to avoid the transverse shear locking and improve the accuracy of the present formulation. The numerical examples presented herein demonstrate that the present method is both free of shear locking and temporal instability. It also achieves high accuracy compared with the exact solutions and other existing methods in the literature.The outline of this paper is as follows.Section 2 describes the weak form of the governing equations and the formulation of the 3-node plate element. In Section 3, the gradient stable integration formulation for Reissner–Mindlin plates with the stabilized discrete shear gap technique is introduced. Section 4 demonstrates the effectiveness of the presented method through numerical examples. Finally, the paper closes with concluding remarks.


## 2. Theoretical Formulations

### 2.1. Basic Equations for Reissner–Mindlin Plates

Based on the assumption of the first-order shear-deformation plate theory, the displacements in the Cartesian coordinate system can be expressed as follows:(1)ux,y,z=u0x,y+zθyx,y,(2)vx,y,z=v0x,y−zθxx,y,(3)wx,y,z=w0x,y,where u0, v0, and w0 are the displacements of the plate midplane in the x, y, and z directions and θx and θy denote the rotations with respect to x and y axes, respectively, as shown in Figure 1.Figure 1
Positive directions of the deflection and two rotations.The relevant strain vectorε can be written in terms of the midplane deformations of equations (1)–(3), which gives(4)ε=εxxεyyγxyγxzγyz=εm0+εb0+0εs,where εm, εb, and εs are the membrane strain, the bending strain (curvature), and the shear strain, respectively:(5)εm=∂u0∂x∂v0∂y∂u0∂y+∂v0∂x,εb=z∂θy∂x−∂θx∂y∂θy∂y−∂θx∂x,εs=∂w0∂x+θy∂w0∂y−θx.By applying the principle of virtual work, the weak form can be stated as follows:(6)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ−∫ΩδuTf˜dΩ−∫ΓδuTt˜dΓ=0,where Dm is the membrane stiffness constitutive coefficients, Db represents the bending stiffness constitutive coefficients, and Ds denotes the transverse shear stiffness constitutive coefficients defined as(7)Dm=∫−t/2t/2D0dz=tD0,(8)Db=∫−t/2t/2z2D0dz=t312D0,(9)Ds=∫−t/2t/2κG1001dz=κtG1001,in which G is the shear modulus, κ=5/6 is the shear correction factor, and the matrix D0 contains the constitutive coefficients:(10)D0=E1−v21v0v10001−v2,where E and v are Young’s modulus and Poisson’s ratio, respectively.Based on the assumption of the first-order shear-deformation plate theory, the weak form for the free vibration analysis of the Reissner–Mindlin plate can be derived from the dynamic form of energy principle, i.e.,(11)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ+∫ΩδdTmd¨dΩ=0,where δd is the variation of the displacement field d=u0v0w0θxθyT and m is the inertia matrix containing the mass density ρ and thickness t:(12)m=ρt00000t00000t00000t31200000t312.For the buckling analysis, when the plate is subjected to in-plane prebuckling stressesσ0, the corresponding weak form can be reformulated as(13)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ+t∫Ωδ∇uTσ0000σ0000σ0∇udΩ+t312∫Ωδ∇θTσ000σ0∇θdΩ=0,where(14)σ0=σxx0τxy0τxy0σyy0,∇u=∇u0∇v0∇w0,∇θ=∇θx∇θy,∇=∂∂x∂∂y.Equation (13) can be rewritten in a compact form as(15)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ+∫ΩδεgTτεgdΩ=0,in which(16)τ=tσ000000tσ000000tσ000000t3σ01200000t3σ012,εgT=u0,xu0,y0000000000v0,xv0,y0000000000w0,xw0,y0000000000θx,xθx,y0000000000θy,xθy,y,where “u0,x” represents the derivative of u0 with respect to x.


### 2.2. Discrete Formulation for Reissner–Mindlin Plates

The bounded domainΩ is discretized into Ne triangular elements such that Ω=∪e=1NeΩe and Ωi∩Ωj=∅, i≠j. For any point in a 3-node triangular element, by using the nodal displacements at the nodes of the element using the shape functions, the generalized displacement field d in the element is interpolated. The same shape functions are used for both displacements and rotations as(17)d=u0,v0,w0,θx,θyT=∑i=1NnNidi,where Nn is the total number of nodes, di=u0i,v0i,w0i,θxi,θyiT is the generalized nodal displacement at node i, and Nix is a diagonal matrix of shape functions given by(18)Nix=diagNix,Nix,Nix,Nix,Nix.Substituting equation (17) into (5), the membrane strain εm, the bending strain εb, and the shear strain εs can be written as(19)εm=∑i=1NnBmidi,(20)εb=∑i=1NnBbidi,(21)εs=∑i=1NnBsidi,in which(22)Bmi=Ni,x00000Ni,y000Ni,yNi,x000,(23)Bbi=0000Ni,x000−Ni,y0000−Ni,xNi,y,(24)Bsi=00Ni,x0Ni00Ni,y−Ni0,where the subscript i=1,2,3.From equations (14)–(16), the geometrical strain εg can be written as(25)εg=∑i=1NnBgidi,in which(26)BgiT=Ni,xNi,y0000000000Ni,xNi,y0000000000Ni,xNi,y0000000000Ni,xNi,y0000000000Ni,xNi,y.Substituting equations (22)–(24) into (6), a set of discretized algebraic system equations of Reissner–Mindlin plates for static analysis can be obtained in the following matrix form:(27)Kd^−f^=0,where d^ denotes the vector of global nodal displacement at all of the nodes and f^ is the force vector (including forces and torques) defined as(28)f^=∫ΩNxfdΩ+∫ΓNxtdΓ,where f and t denote the distributed load and prescribed boundary load, respectively.In equation (27), the global stiffness matrix K can be expressed as(29)K=∑e=1NeKme+Kbe+Kse,in which(30)Kme=∫ΩeBmTDmBmdΩ,(31)Kbe=∫ΩeBbTDbBbdΩ,(32)Kse=∫ΩeBsTDsBsdΩ.For free vibration, we have(33)K−ω2Md^=0,where ω denotes the natural frequency and M is the global inertia matrix:(34)M=∑k=1Ne∫ΩeNTmNdΩ.For the buckling analysis, we have(35)K−λcrKgd^=0,where(36)Kg=∑k=1Ne∫ΩeBgTτBgdΩ,which is the geometrical stiffness matrix, and λcr denotes the critical buckling load. In addition, it is noted that the summation in equations (29), (34), and (36) means an assembly process.


## 2.1. Basic Equations for Reissner–Mindlin Plates

Based on the assumption of the first-order shear-deformation plate theory, the displacements in the Cartesian coordinate system can be expressed as follows:(1)ux,y,z=u0x,y+zθyx,y,(2)vx,y,z=v0x,y−zθxx,y,(3)wx,y,z=w0x,y,where u0, v0, and w0 are the displacements of the plate midplane in the x, y, and z directions and θx and θy denote the rotations with respect to x and y axes, respectively, as shown in Figure 1.Figure 1
Positive directions of the deflection and two rotations.The relevant strain vectorε can be written in terms of the midplane deformations of equations (1)–(3), which gives(4)ε=εxxεyyγxyγxzγyz=εm0+εb0+0εs,where εm, εb, and εs are the membrane strain, the bending strain (curvature), and the shear strain, respectively:(5)εm=∂u0∂x∂v0∂y∂u0∂y+∂v0∂x,εb=z∂θy∂x−∂θx∂y∂θy∂y−∂θx∂x,εs=∂w0∂x+θy∂w0∂y−θx.By applying the principle of virtual work, the weak form can be stated as follows:(6)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ−∫ΩδuTf˜dΩ−∫ΓδuTt˜dΓ=0,where Dm is the membrane stiffness constitutive coefficients, Db represents the bending stiffness constitutive coefficients, and Ds denotes the transverse shear stiffness constitutive coefficients defined as(7)Dm=∫−t/2t/2D0dz=tD0,(8)Db=∫−t/2t/2z2D0dz=t312D0,(9)Ds=∫−t/2t/2κG1001dz=κtG1001,in which G is the shear modulus, κ=5/6 is the shear correction factor, and the matrix D0 contains the constitutive coefficients:(10)D0=E1−v21v0v10001−v2,where E and v are Young’s modulus and Poisson’s ratio, respectively.Based on the assumption of the first-order shear-deformation plate theory, the weak form for the free vibration analysis of the Reissner–Mindlin plate can be derived from the dynamic form of energy principle, i.e.,(11)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ+∫ΩδdTmd¨dΩ=0,where δd is the variation of the displacement field d=u0v0w0θxθyT and m is the inertia matrix containing the mass density ρ and thickness t:(12)m=ρt00000t00000t00000t31200000t312.For the buckling analysis, when the plate is subjected to in-plane prebuckling stressesσ0, the corresponding weak form can be reformulated as(13)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ+t∫Ωδ∇uTσ0000σ0000σ0∇udΩ+t312∫Ωδ∇θTσ000σ0∇θdΩ=0,where(14)σ0=σxx0τxy0τxy0σyy0,∇u=∇u0∇v0∇w0,∇θ=∇θx∇θy,∇=∂∂x∂∂y.Equation (13) can be rewritten in a compact form as(15)∫ΩδεmTDmεmdΩ+∫ΩδεbTDbεbdΩ+∫ΩδεsTDsεsdΩ+∫ΩδεgTτεgdΩ=0,in which(16)τ=tσ000000tσ000000tσ000000t3σ01200000t3σ012,εgT=u0,xu0,y0000000000v0,xv0,y0000000000w0,xw0,y0000000000θx,xθx,y0000000000θy,xθy,y,where “u0,x” represents the derivative of u0 with respect to x.


## 2.2. Discrete Formulation for Reissner–Mindlin Plates

The bounded domainΩ is discretized into Ne triangular elements such that Ω=∪e=1NeΩe and Ωi∩Ωj=∅, i≠j. For any point in a 3-node triangular element, by using the nodal displacements at the nodes of the element using the shape functions, the generalized displacement field d in the element is interpolated. The same shape functions are used for both displacements and rotations as(17)d=u0,v0,w0,θx,θyT=∑i=1NnNidi,where Nn is the total number of nodes, di=u0i,v0i,w0i,θxi,θyiT is the generalized nodal displacement at node i, and Nix is a diagonal matrix of shape functions given by(18)Nix=diagNix,Nix,Nix,Nix,Nix.Substituting equation (17) into (5), the membrane strain εm, the bending strain εb, and the shear strain εs can be written as(19)εm=∑i=1NnBmidi,(20)εb=∑i=1NnBbidi,(21)εs=∑i=1NnBsidi,in which(22)Bmi=Ni,x00000Ni,y000Ni,yNi,x000,(23)Bbi=0000Ni,x000−Ni,y0000−Ni,xNi,y,(24)Bsi=00Ni,x0Ni00Ni,y−Ni0,where the subscript i=1,2,3.From equations (14)–(16), the geometrical strain εg can be written as(25)εg=∑i=1NnBgidi,in which(26)BgiT=Ni,xNi,y0000000000Ni,xNi,y0000000000Ni,xNi,y0000000000Ni,xNi,y0000000000Ni,xNi,y.Substituting equations (22)–(24) into (6), a set of discretized algebraic system equations of Reissner–Mindlin plates for static analysis can be obtained in the following matrix form:(27)Kd^−f^=0,where d^ denotes the vector of global nodal displacement at all of the nodes and f^ is the force vector (including forces and torques) defined as(28)f^=∫ΩNxfdΩ+∫ΓNxtdΓ,where f and t denote the distributed load and prescribed boundary load, respectively.In equation (27), the global stiffness matrix K can be expressed as(29)K=∑e=1NeKme+Kbe+Kse,in which(30)Kme=∫ΩeBmTDmBmdΩ,(31)Kbe=∫ΩeBbTDbBbdΩ,(32)Kse=∫ΩeBsTDsBsdΩ.For free vibration, we have(33)K−ω2Md^=0,where ω denotes the natural frequency and M is the global inertia matrix:(34)M=∑k=1Ne∫ΩeNTmNdΩ.For the buckling analysis, we have(35)K−λcrKgd^=0,where(36)Kg=∑k=1Ne∫ΩeBgTτBgdΩ,which is the geometrical stiffness matrix, and λcr denotes the critical buckling load. In addition, it is noted that the summation in equations (29), (34), and (36) means an assembly process.


## 3. The Formulation of Gradient Stable Node-Based Smoothed Integration and Discrete Shear Gap Technique

According to the introduction inSection 2.2, the structural stiffness matrix is composed of three parts, namely, Km, Kb, and Ks. On the one hand, for Km and Kb, we take the method of gradient stable node-based smoothed integration to avoid temporal instability and spatial instability; for Ks, on the other hand, the discrete shear gap technique is employed to avoid the shear-locking problem. The formula for the smoothed integral is derived from the Hu–Washizu three-field variational principle as shown in Section 3.1, and its discrete form is given in Section 3.2. The discrete shear technology is introduced in Section 3.3.


### 3.1. Gradient Stable Smoothed Derivative Correction

Based on the Hu–Washizu three-field variational principle, Duan [68] gives the corrected nodal derivative using more rigorous mathematics, and the quadratically consistent nodal integration is proposed for second-order meshfree Galerkin methods. However, the proposal is more complex and requires two-order Gauss integration of the boundary integral. In this paper, a simplified scheme is provided by using another correct method.Assume thatu and σ^ are the displacement and assumed Cauchy stress, respectively, ε˜ is the interpolated strain or smoothed strain, and ε is the actual strain. The Hu–Washizu three-field weak form for the elastic-static problem can be written as(37)δΠ∗u,ε˜,σ˜=∫Ωδε˜TDε˜dΩ−∫ΓtδuTtdΓ−∫ΩδuTfdΩ+δ∫Ωσ^Tε−ε˜dΩ.Clearly, if the interpolated strainε˜ can be somehow constructed from the displacement u and meet the following orthogonality condition:(38)∫Ωσ^Tε−ε˜dΩ=0.Then, a form containing only independent variables can be obtained as simple as the classical one:(39)δΠ∗∗u,ε˜,σ˜=∫Ωδε˜TDε˜dΩ−∫ΓtδuTt¯dΓ−∫ΩδuTf¯dΩ.Deriving by reference [68], in order to meet the orthogonality condition as equation (38), let the following equation be satisfied for each subdomain ΩL and expressed as a finite element form:(40)∫ΩLσ^T∂NI∂xi−∂N˜I∂xidΩ=0,where NI is the shape function corresponding to strain ε in element I and N˜I is the shape function corresponding to smoothed strain ε˜ in element I. Since stress and strain are related to the first partial derivative of displacement, which is a polynomial combination of coordinates, stress can also be regarded as a polynomial of position. Equation (40) can be equivalent to the following:(41)∫ΩLqx∂NI∂xi−∂N˜I∂xidΩ=0,where qx is the space base which is one order lower than the space for the displacement u.Different from reference [68] which chose the quadratic base for displacement, the first base is selected for displacement as the shape function of the first-order triangular element, and qx=1. Then, the following equations can be obtained:(42)∫ΩL∂N˜I∂xidΩ=∫ΓLNInidΓ,(43)∫ΩL∂N˜2I∂xi∂xjdΩ=∫ΓLNI,injdΓ.In the nodal integration scheme, nodeL is the only point in ΩL, and there exists only one unknown N˜I,xxL. Hence, equations (42) and (43) cannot be satisfied at the same time. To this end, we introduce the derivatives of the function N˜I,x by means of Taylor’s expansion such that N˜I,xxxL and N˜I,xyxL are introduced and can serve as the other two unknowns. Taylor’s expansion for N˜I,x can be formulated as(44)N˜I,xx=N˜I,xxL+x−xLN˜I,xxxL+y−yLN˜I,xyxL+H.O.T,(45)N˜I,xxx=N˜I,xxxL+H.O.T,(46)N˜I,xyx=N˜I,xyxL+H.O.T,where H.O.T means higher-order terms.Substitution of equations (56)–(58) into (42) and (43) leads to(47)ALILxILy0AL000ALN˜I,xxLN˜I,xxxLN˜I,xyxL=∫ΓLNInxdΓ∫ΓLNI,xnxdΓ∫ΓLNI,xnydΓ,where AL=∫ΩLdΩ is the area of ΩL and(48)ILx=∫ΩLx−xLdΩ,ILy=∫ΩLy−yLdΩ.By solving equation (47), the corrected nodal derivative N˜I,xxL and its derivatives, i.e., N˜I,xxxL and N˜I,xyxL, are obtained. Following the same derivation, Taylor’s expansion for N˜I,y is(49)N˜I,yx=N˜I,yxL+x−xLN˜I,yxxL+y−yLN˜I,yyxL+H.O.T,(50)N˜I,yxx=N˜I,yxxL+H.O.T,(51)N˜I,yyx=N˜I,yyxL+H.O.T.The equation fory-derivatives can be written as(52)ALILxILy0AL000ALN˜I,yxLN˜I,yxxLN˜I,yyxL=∫ΓLNInydΓ∫ΓLNI,ynxdΓ∫ΓLNI,ynydΓ.To simplify the calculation, the equivalent circle domain can be assumed for the subdomainΩL. The area AL, first moments ILx and ILy, and second moments of inertia ILxx, ILxy, and ILyy of each nodal domain can be expressed as(53)AL=∫ΩLdΩ,ILx=ILy=ILxy=0,ILxx=ILyy=AL24π.It should be noted that the concept of equivalent circles is only introduced to simplify the calculation of these integrals, and the actual smooth region is still a polygon composed of elements and nodes.


### 3.2. Discrete Form of Gradient Stabilized Nodal Integration

In this part, the nodal integration formulation will be introduced. We can discretize the problem domainΩ with Ne triangular elements as in standard FEM, but the integral required in this work is now based on the node and utilizes strain smoothing operations. In the process of such node integration, the middle edge points are connected with the center points of the surrounding triangular elements in order to form the smoothed domain of each node sequentially, as shown in Figure 2, such that Ω=∪k=1NnΩk and Ωi∩Ωj=ϕ for i≠j, in which Nn is the total number of nodes of the problem domain.Figure 2
Triangular elements and smoothing cells associated with nodes.From equations (44)–(46) and equations (49)–(51), the smoothed strain ε˜mx and ε˜bx in ΩL surrounding node L can be expressed as(54)ε˜mx=ε˜mL+ε˜mxLx−xL+ε˜myLy−yL,(55)ε˜bx=ε˜bL+ε˜bxLx−xL+ε˜byLy−yL,where(56)ε˜mL=B˜mLd^,ε˜mxL=B˜mxLd^,ε˜myL=B˜myLd^,(57)ε˜bL=B˜bLd^,ε˜bxL=B˜bxLd^,ε˜byL=B˜byLd^,in which(58)B˜mL=B˜m1L⋯B˜mIL⋯B˜mNnL,B˜bL=B˜b1L⋯B˜bIL⋯B˜bNnL,(59)B˜mxL=B˜mx1L⋯B˜mxIL⋯B˜mxNnL,B˜bxL=B˜bx1L⋯B˜bxIL⋯B˜bxNnL,(60)B˜myL=B˜my1L⋯B˜myIL⋯B˜myNnL,B˜byL=B˜by1L⋯B˜byIL⋯B˜byNnL,(61)B˜mIL=N˜I,x00000N˜I,y000N˜I,yN˜I,x000,B˜bIL=0000N˜I,x000−N˜I,y0000−N˜I,xN˜I,y,(62)B˜mxIL=N˜I,xx00000N˜I,yx000N˜I,yxN˜I,xx000,B˜bxIL=0000N˜I,xx000−N˜I,yx0000−N˜I,xxN˜I,yx,(63)B˜myIL=N˜I,xy00000N˜I,yy000N˜I,yyN˜I,xy000,B˜byIL=0000N˜I,xy000−N˜I,yy0000−N˜I,xyN˜I,yy.


### 3.3. Formulation of the Stabilized Discrete Shear Gap Technique

The discrete shear gap method is adopted here to eliminate the shear locking. In each triangular element, the nodes are denoted anticlockwise asi, j, and k, respectively. The shear strain can be given as(64)γxz=∑mNm,xxΔwxm+∑mNm,xxΔwym,(65)γyz=∑mNm,yxΔwxm+∑mNm,yxΔwym,where Δwxm and Δwym are the discrete shear gaps at the node m given by(66)Δwxi=Δwxk=Δwyi=Δwyj=0,(67)Δwxj=wj−wi−12bθxi+θxj+12aθyi+θyj,(68)Δwyk=wk−wi−12dθxi+θxk+12cθyi+θyk,where(69)a=xj−xi,b=yj−yi,(70)c=xk−xi,d=yk−yi.From equations (64)–(70), the shear strain εs in each element can be rewritten as(71)εs=γxzγyz=Bsd=BsiBsjBskdidjdk,in which(72)Bsi=12Ae00b−d0Ae00c−a−Ae0,(73)Bsj=12Ae00d−bd2ad200−cbc2−ac2,(74)Bsk=12Ae00−bbd2−bc200a−ad2ac2,where Ae is the area of the element.To improve significantly the accuracy of approximate solutions and to stabilize shear force oscillations presenting the triangular element, a stabilization technique [69, 70] needs to be added to the original discrete shear gap element. Therefore, the transverse shear stiffness constitutive coefficients Ds should be corrective as D¯s:(75)D¯s=κt3Gt2+αhe21001,where α is a positive constant [69, 71] and the characteristic length he can be estimated as the diameter of the equation circle domain:(76)he=2Aeπ.


### 3.4. Discrete Formulation for GS-DSG Method

We now seek for a weak form solution of the generalized displacement fieldd that satisfies the following smoothed Galerkin weak form:(77)∫Ωδε˜mTDmε˜mdΩ+∫Ωδε˜bTDbε˜bdΩ+∫ΩεsTDsεsdΩ−∫ΩδdTfdΩ−∫ΓδdTtdΓ=0.Substituting equations (74)–(76) into (77), a set of discretized algebraic system equations can be obtained in the following matrix form:(78)K˜d^−f^=0,where K˜ is the global smoothed stiffness matrix assembled in the form of(79)K˜=∑L=1NnK˜mL+K˜bL+∑e=1NeKse.The summation in equation (79) means an assembly process same as the practice in the FEM, and Nn is the number of the nodes of the whole problem domain Ω. Kse is given from equation (32) by using the stabilized discrete shear gap method to eliminate the shear locking. K¯mL and K¯bL are the stiffness matrices associated with node L given as(80)K˜mL=B˜mLTDmB˜mLAL+B˜mxLTDmB˜mxLILxx+B˜myLTDmB˜myLILyy,(81)K˜bL=B˜bLTDbB˜bLAL+B˜bxLTDbB˜bxLILxx+B˜byLTDbB˜byLILyy.For free vibration, we have(82)K˜−ω2Md^=0.For the buckling analysis, we have(83)K˜−λcrK˜gd^=0,where K˜g denotes the geometrical stiffness matrix assembled in the form of(84)K˜g=∑L=1NnK˜gL.The nodal geometrical stiffness matrixK˜gL in equation (84) can be calculated by(85)K˜gL=B˜gLTτB˜gLAL+B˜gxLTτB˜gxLILxx+B˜gyLTτB˜gyLILyy,with(86)B˜gLT=N˜i,xN˜i,y0000000000N˜i,xN˜i,y0000000000N˜i,xN˜i,y0000000000N˜i,xN˜i,y0000000000N˜i,xN˜i,y,(87)B˜gxLT=N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx,(88)B˜gyLT=N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy.


## 3.1. Gradient Stable Smoothed Derivative Correction

Based on the Hu–Washizu three-field variational principle, Duan [68] gives the corrected nodal derivative using more rigorous mathematics, and the quadratically consistent nodal integration is proposed for second-order meshfree Galerkin methods. However, the proposal is more complex and requires two-order Gauss integration of the boundary integral. In this paper, a simplified scheme is provided by using another correct method.Assume thatu and σ^ are the displacement and assumed Cauchy stress, respectively, ε˜ is the interpolated strain or smoothed strain, and ε is the actual strain. The Hu–Washizu three-field weak form for the elastic-static problem can be written as(37)δΠ∗u,ε˜,σ˜=∫Ωδε˜TDε˜dΩ−∫ΓtδuTtdΓ−∫ΩδuTfdΩ+δ∫Ωσ^Tε−ε˜dΩ.Clearly, if the interpolated strainε˜ can be somehow constructed from the displacement u and meet the following orthogonality condition:(38)∫Ωσ^Tε−ε˜dΩ=0.Then, a form containing only independent variables can be obtained as simple as the classical one:(39)δΠ∗∗u,ε˜,σ˜=∫Ωδε˜TDε˜dΩ−∫ΓtδuTt¯dΓ−∫ΩδuTf¯dΩ.Deriving by reference [68], in order to meet the orthogonality condition as equation (38), let the following equation be satisfied for each subdomain ΩL and expressed as a finite element form:(40)∫ΩLσ^T∂NI∂xi−∂N˜I∂xidΩ=0,where NI is the shape function corresponding to strain ε in element I and N˜I is the shape function corresponding to smoothed strain ε˜ in element I. Since stress and strain are related to the first partial derivative of displacement, which is a polynomial combination of coordinates, stress can also be regarded as a polynomial of position. Equation (40) can be equivalent to the following:(41)∫ΩLqx∂NI∂xi−∂N˜I∂xidΩ=0,where qx is the space base which is one order lower than the space for the displacement u.Different from reference [68] which chose the quadratic base for displacement, the first base is selected for displacement as the shape function of the first-order triangular element, and qx=1. Then, the following equations can be obtained:(42)∫ΩL∂N˜I∂xidΩ=∫ΓLNInidΓ,(43)∫ΩL∂N˜2I∂xi∂xjdΩ=∫ΓLNI,injdΓ.In the nodal integration scheme, nodeL is the only point in ΩL, and there exists only one unknown N˜I,xxL. Hence, equations (42) and (43) cannot be satisfied at the same time. To this end, we introduce the derivatives of the function N˜I,x by means of Taylor’s expansion such that N˜I,xxxL and N˜I,xyxL are introduced and can serve as the other two unknowns. Taylor’s expansion for N˜I,x can be formulated as(44)N˜I,xx=N˜I,xxL+x−xLN˜I,xxxL+y−yLN˜I,xyxL+H.O.T,(45)N˜I,xxx=N˜I,xxxL+H.O.T,(46)N˜I,xyx=N˜I,xyxL+H.O.T,where H.O.T means higher-order terms.Substitution of equations (56)–(58) into (42) and (43) leads to(47)ALILxILy0AL000ALN˜I,xxLN˜I,xxxLN˜I,xyxL=∫ΓLNInxdΓ∫ΓLNI,xnxdΓ∫ΓLNI,xnydΓ,where AL=∫ΩLdΩ is the area of ΩL and(48)ILx=∫ΩLx−xLdΩ,ILy=∫ΩLy−yLdΩ.By solving equation (47), the corrected nodal derivative N˜I,xxL and its derivatives, i.e., N˜I,xxxL and N˜I,xyxL, are obtained. Following the same derivation, Taylor’s expansion for N˜I,y is(49)N˜I,yx=N˜I,yxL+x−xLN˜I,yxxL+y−yLN˜I,yyxL+H.O.T,(50)N˜I,yxx=N˜I,yxxL+H.O.T,(51)N˜I,yyx=N˜I,yyxL+H.O.T.The equation fory-derivatives can be written as(52)ALILxILy0AL000ALN˜I,yxLN˜I,yxxLN˜I,yyxL=∫ΓLNInydΓ∫ΓLNI,ynxdΓ∫ΓLNI,ynydΓ.To simplify the calculation, the equivalent circle domain can be assumed for the subdomainΩL. The area AL, first moments ILx and ILy, and second moments of inertia ILxx, ILxy, and ILyy of each nodal domain can be expressed as(53)AL=∫ΩLdΩ,ILx=ILy=ILxy=0,ILxx=ILyy=AL24π.It should be noted that the concept of equivalent circles is only introduced to simplify the calculation of these integrals, and the actual smooth region is still a polygon composed of elements and nodes.


## 3.2. Discrete Form of Gradient Stabilized Nodal Integration

In this part, the nodal integration formulation will be introduced. We can discretize the problem domainΩ with Ne triangular elements as in standard FEM, but the integral required in this work is now based on the node and utilizes strain smoothing operations. In the process of such node integration, the middle edge points are connected with the center points of the surrounding triangular elements in order to form the smoothed domain of each node sequentially, as shown in Figure 2, such that Ω=∪k=1NnΩk and Ωi∩Ωj=ϕ for i≠j, in which Nn is the total number of nodes of the problem domain.Figure 2
Triangular elements and smoothing cells associated with nodes.From equations (44)–(46) and equations (49)–(51), the smoothed strain ε˜mx and ε˜bx in ΩL surrounding node L can be expressed as(54)ε˜mx=ε˜mL+ε˜mxLx−xL+ε˜myLy−yL,(55)ε˜bx=ε˜bL+ε˜bxLx−xL+ε˜byLy−yL,where(56)ε˜mL=B˜mLd^,ε˜mxL=B˜mxLd^,ε˜myL=B˜myLd^,(57)ε˜bL=B˜bLd^,ε˜bxL=B˜bxLd^,ε˜byL=B˜byLd^,in which(58)B˜mL=B˜m1L⋯B˜mIL⋯B˜mNnL,B˜bL=B˜b1L⋯B˜bIL⋯B˜bNnL,(59)B˜mxL=B˜mx1L⋯B˜mxIL⋯B˜mxNnL,B˜bxL=B˜bx1L⋯B˜bxIL⋯B˜bxNnL,(60)B˜myL=B˜my1L⋯B˜myIL⋯B˜myNnL,B˜byL=B˜by1L⋯B˜byIL⋯B˜byNnL,(61)B˜mIL=N˜I,x00000N˜I,y000N˜I,yN˜I,x000,B˜bIL=0000N˜I,x000−N˜I,y0000−N˜I,xN˜I,y,(62)B˜mxIL=N˜I,xx00000N˜I,yx000N˜I,yxN˜I,xx000,B˜bxIL=0000N˜I,xx000−N˜I,yx0000−N˜I,xxN˜I,yx,(63)B˜myIL=N˜I,xy00000N˜I,yy000N˜I,yyN˜I,xy000,B˜byIL=0000N˜I,xy000−N˜I,yy0000−N˜I,xyN˜I,yy.


## 3.3. Formulation of the Stabilized Discrete Shear Gap Technique

The discrete shear gap method is adopted here to eliminate the shear locking. In each triangular element, the nodes are denoted anticlockwise asi, j, and k, respectively. The shear strain can be given as(64)γxz=∑mNm,xxΔwxm+∑mNm,xxΔwym,(65)γyz=∑mNm,yxΔwxm+∑mNm,yxΔwym,where Δwxm and Δwym are the discrete shear gaps at the node m given by(66)Δwxi=Δwxk=Δwyi=Δwyj=0,(67)Δwxj=wj−wi−12bθxi+θxj+12aθyi+θyj,(68)Δwyk=wk−wi−12dθxi+θxk+12cθyi+θyk,where(69)a=xj−xi,b=yj−yi,(70)c=xk−xi,d=yk−yi.From equations (64)–(70), the shear strain εs in each element can be rewritten as(71)εs=γxzγyz=Bsd=BsiBsjBskdidjdk,in which(72)Bsi=12Ae00b−d0Ae00c−a−Ae0,(73)Bsj=12Ae00d−bd2ad200−cbc2−ac2,(74)Bsk=12Ae00−bbd2−bc200a−ad2ac2,where Ae is the area of the element.To improve significantly the accuracy of approximate solutions and to stabilize shear force oscillations presenting the triangular element, a stabilization technique [69, 70] needs to be added to the original discrete shear gap element. Therefore, the transverse shear stiffness constitutive coefficients Ds should be corrective as D¯s:(75)D¯s=κt3Gt2+αhe21001,where α is a positive constant [69, 71] and the characteristic length he can be estimated as the diameter of the equation circle domain:(76)he=2Aeπ.


## 3.4. Discrete Formulation for GS-DSG Method

We now seek for a weak form solution of the generalized displacement fieldd that satisfies the following smoothed Galerkin weak form:(77)∫Ωδε˜mTDmε˜mdΩ+∫Ωδε˜bTDbε˜bdΩ+∫ΩεsTDsεsdΩ−∫ΩδdTfdΩ−∫ΓδdTtdΓ=0.Substituting equations (74)–(76) into (77), a set of discretized algebraic system equations can be obtained in the following matrix form:(78)K˜d^−f^=0,where K˜ is the global smoothed stiffness matrix assembled in the form of(79)K˜=∑L=1NnK˜mL+K˜bL+∑e=1NeKse.The summation in equation (79) means an assembly process same as the practice in the FEM, and Nn is the number of the nodes of the whole problem domain Ω. Kse is given from equation (32) by using the stabilized discrete shear gap method to eliminate the shear locking. K¯mL and K¯bL are the stiffness matrices associated with node L given as(80)K˜mL=B˜mLTDmB˜mLAL+B˜mxLTDmB˜mxLILxx+B˜myLTDmB˜myLILyy,(81)K˜bL=B˜bLTDbB˜bLAL+B˜bxLTDbB˜bxLILxx+B˜byLTDbB˜byLILyy.For free vibration, we have(82)K˜−ω2Md^=0.For the buckling analysis, we have(83)K˜−λcrK˜gd^=0,where K˜g denotes the geometrical stiffness matrix assembled in the form of(84)K˜g=∑L=1NnK˜gL.The nodal geometrical stiffness matrixK˜gL in equation (84) can be calculated by(85)K˜gL=B˜gLTτB˜gLAL+B˜gxLTτB˜gxLILxx+B˜gyLTτB˜gyLILyy,with(86)B˜gLT=N˜i,xN˜i,y0000000000N˜i,xN˜i,y0000000000N˜i,xN˜i,y0000000000N˜i,xN˜i,y0000000000N˜i,xN˜i,y,(87)B˜gxLT=N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx0000000000N˜i,xxN˜i,yx,(88)B˜gyLT=N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy0000000000N˜i,xyN˜i,yy.


## 4. Numerical Examples

In this section, static, free vibration, and buckling analyses of square, T-shape, elliptical, and rectangular plates are considered. In addition, the present method is compared with other three methods, the FEM-DSG, NS-FEM, and NS-DSG methods. To examine the numerical error precisely, the displacement error norm is defined as(89)eu=∑i=1Nnuiexact−uinum2∑i=1Nnuiexact2,where the superscript “exact” denotes the exact solution (if an exact solution does not exist, “exact” is the reference solution) and “num” represents the displacement vector obtained using numerical methods including the present method.In the following example, material parameters’ Young's modulus is expressed asE, Poisson’s ratio is expressed as v, and mass density is expressed as ρ.


### 4.1. Static Analysis

#### 4.1.1. Square Plate

Consider the model of a simply supported square plate subjected to a uniform loadq=1 as shown in Figure 3. The geometric and material parameters are length L=10 and thickness t=2.0; E=2.1×109 and v=0.3. Due to symmetry, only a quarter of the plate is modeled to reduce the computation cost, and uniform meshes are employed. The center deflection wc is normalized as w^=wcD/qL4, where D=Et3/121−v2 is the bending stiffness. In order to test the performance of the mentioned numerical methods, the numerical results obtained using the present method are compared with other three methods. The result calculated by the ABAQUS software is used for reference, using S4R elements and a large number (37,249) of nodes.Figure 3
A square plate with a simply supported boundary condition.Table1 shows the numerical results of the normalized center deflection. Figure 4 shows the relative error, and the label “mesh density” on the horizontal axis shows the number of cells on each side. Figure 5 shows the convergence status of the displacement error norm eu, where h is the average nodal spacing of the node distribution. From the results, it can be seen that the deflection obtained by NS-FEM and NS-DSG is larger than the reference solution, whereas the deflection solved by the proposed method and FEM-DSG is smaller than the reference solution. Meanwhile, the proposed method is more accurate than the others. By comparing the convergence rate of the methods as shown in Figure 5, the proposed method has a higher convergence rate than the others, as far as the average nodal spacing trail off is concerned.Table 1
Numerical results of normalized central deflectionw^ for the simply supported square plate subjected to uniform load.
MeshFEM-DSGNS-FEMNS-DSGGS-DSG8 × 80.0054420.0056780.0055670.00551916 × 160.0055160.0055850.0055540.00554024 × 240.0055320.0055630.0055500.00554432 × 320.0055380.0055550.0055490.005545Reference0.005545Figure 4
Central deflection.Figure 5
Convergence properties.


#### 4.1.2. T-Shaped Plate

In this section, a T-shaped plate with clamped edges and subjected to two kinds of loads, i.e., uniform and concentrated loads, is analyzed to further examine the efficiency of the present method. The geometric parameters are shown in Figure6. Two thickness are considered, t=0.01 and t=0.5. Material parameters are E=2.1×109 and v=0.3. The uniform load applied to the entire plate is given by q=0.01, and the concentrated load applied to point A is taken as p=1. Due to the symmetry of the plate, only half of the model is studied to reduce the calculation cost.Figure 6
Geometry illustration of a T-plate.Figure7 shows the mesh model, which is discretized using 154 nodes with 243 triangular elements. Numerical results of the present method are compared with other three methods. Since the analytical solution is unavailable for this problem, the result calculated by the ABAQUS software with a very fine mesh (2751 nodes and 5200 elements) was used for reference. The deflections along the line OA are plotted as shown in Figures 8 and 9. From the results, it can be seen that, for the thick plate, the results are almost identical, and the result is close to the reference solution. That is, because the shear-locking phenomenon does not appear in the thick plate, all the four methods can obtain high accuracy and is hard to distinguish which is higher. However, for the thin plate, the difference in the accuracy of the four methods is obvious. The accuracy of FEM-DSG and NS-DGS is lower, and both NS-FEM and the proposed method can achieve high accuracy compared with the reference result.Figure 7
Domain discretization using triangular.Figure 8
Deflections along the line OA with different thickness (uniformly distributed load): (a)t=0.01; (b) t=0.5.
(a)(b)Figure 9
Deflections along the line OA with different thickness (concentrated load): (a)t=0.01; (b) t=0.5.
(a)(b)


### 4.2. Free Vibration Analysis

In this section, numerical examples of free vibration for various plates are given. The nondimensional frequency parameter is normalized byλ=ωa2ρt/D1/2 [35], where ω is the circle frequency value, a is the geometry size given in each problem, ρ is the mass density, t is the thickness, and D=Et3/121−v2 is the bending stiffness.


#### 4.2.1. Square Plates

Square plates of lengtha, width b, and thickness t are considered. The material parameters are E=2.0×1011, v=0.3, and ρ=8000. The plate is modeled with uniform meshes of 4, 8, 16, and 24 elements each side. The boundary conditions are simply supported (S), clamped (C), and free (F). SSSS means that all four sides are simply supported, and the others are similar.First, the SSSS plate corresponding to length-to-width ratiosa/b=1 is considered. The thickness-to-length for the thin plate is t/a=0.005 and for the thick plate is t/a=0.1, respectively. Figures 10(a), 10(c), and 10(d) show the geometry of the plate and its mesh grid, respectively. Table 2 shows the values of the nondimensional frequency parameter λ corresponding to the six frequencies using 4 × 4, 8 × 8, 16 × 16, and 24 × 24 meshes. It is observed that the accuracy of the presented method increases with the decreasing size of the mesh elements, and the results of GS-DSG agree well with the analytical ones. For the same mesh, the presented method is more accurate than FEM-DSG, NS-FEM, and NS-DSG elements for both thin and thick plates. Table 3 shows four methods for the first six modes under the 16×16 mesh. The spurious nonzero energy mode is marked with the blue wireframe. It can be clearly found in the table that NS-FEM has spurious nonzero energy modes, that is, severe time instability. The mode obtained by the NS-DSG method cannot eliminate spurious nonzero energy modes completely even using the discrete shear gap technology. In contrast, the advantage of the GS-DSG method is particularly obvious, and there are no spurious nonzero energy modes, which indicates its stability in the time domain.Figure 10
(a) Supported plate, (b) clamped plate, and (c, d) triangular meshes.
(a)(b)(c)(d)Table 2
A nondimensional frequency parameterλ of an SSSS plate (a/b) = 1).
t/aMethodMode sequence number1234560.005FEM-DSG4.80678.50389.263112.335316.657317.94004.51907.35297.51239.705811.207111.27984.45457.09437.13449.074810.214810.21804.44407.04997.06808.959110.048110.0521NS-FEM3.85234.61804.66965.48035.58055.70684.28734.74644.95285.42435.97026.24744.39314.81824.82385.28625.84636.10544.41564.76594.76605.48646.28916.2921NS-DSG4.19526.24896.37817.40327.51548.29444.31916.74456.75738.43709.31319.37584.38536.92346.92488.70099.73669.74624.40536.96586.96678.77539.83089.8348GS-DSG4.28936.78886.93138.38169.58179.60354.36376.88886.89668.67079.69919.74034.41016.96876.97298.78469.82609.83934.42286.99186.99458.82499.87389.8807Reference [72]4.4437.0257.0258.8869.9359.9350.1FEM-DSG4.67007.97138.581110.326110.601110.60174.39646.99717.13288.96159.753210.18494.32006.75656.78748.42019.41199.41854.30126.70866.72168.30939.27339.2794NS-FEM3.73664.41024.43975.14035.22705.37014.14384.43894.58625.01525.45765.66164.23274.45344.45534.89235.34345.53594.25344.40574.40585.02805.69715.6991NS-DSG4.10265.94186.03666.75216.78437.15404.21756.46116.47287.38217.67957.79824.26756.61586.61587.62057.84467.84514.27596.64396.64587.70187.86037.8604GS-DSG4.18636.46436.54167.74828.50518.64284.25466.58596.58778.08668.94138.96374.28166.64636.64878.18179.09679.09794.28306.65766.65938.20029.13149.1365Reference [72]4.376.746.748.359.229.22Table 3
First 6 modes of the square plate (SSSS) solved by different methods with16×16 meshes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGt/a = 0.005123456t/a = 0.1123456The second problem is that the mesh of CCCC square plate as shown in Figure10(b) is the same as that of the SSSS plate. The nondimensional frequency parameter λ corresponding to the first six frequencies of the CCCC plate is shown in Table 4. The corresponding modes under 16×16 are given in Table 5, and the spurious nonzero energy mode is marked with the blue wireframe. It can be found again that the GS-DSG method is superior to the other three methods.Table 4
A nondimensional frequency parameterλ of a CCCC plate (a/b = 1).
t/aMethodMode sequence number1234560.005FEM-DSG8.674314.607118.478122.121528.389129.17946.38399.46939.773212.410313.889214.23656.07408.74228.809710.808911.947311.98236.02848.63558.655410.565311.662911.6897NS-FEM4.866227.299427.403030.548531.246833.71825.51135.63695.69086.15736.73456.85785.14635.15675.65975.89896.21536.54284.96374.96395.74115.95206.51026.5153NS-DSG6.23997.783314.514316.878718.550419.45848.84888.09458.26089.961310.415710.82715.93568.39818.417810.174311.114411.15765.96608.48108.487910.277011.289011.3201GS-DSG7.139110.478814.736217.063918.806019.53715.99028.46718.617410.622011.286811.64865.96948.48568.501910.311211.305811.33985.98078.51918.524910.336311.372411.4008Reference [73]5.9998.5688.56810.40711.47211.4980.1FEM-DSG7.025510.228411.353911.385011.834012.80155.95768.40118.648910.377411.208411.36915.76097.99688.05259.583810.427810.48035.72747.92717.95149.436910.251710.2995NS-FEM4.61546.34326.78676.85687.08587.64054.81194.97495.37965.47435.91445.93804.55014.55725.11115.51085.62065.72904.44944.44945.25455.66585.75245.7533NS-DSG5.71317.06307.84268.22228.38018.46985.55407.48527.55087.76387.91538.25005.65447.75587.76757.92947.93108.03545.67977.81907.82427.91327.91337.9705GS-DSG6.13058.39878.998610.205910.309510.35115.65057.73777.76839.23859.85429.89665.67877.81307.82259.24779.998510.04215.69027.84337.84799.283110.057510.1036Reference [73]5.717.887.889.3310.1310.18Table 5
First 6 modes of the square plate (CCCC) solved by different methods with16×16 meshes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGt/a = 0.005123456t/a = 0.1123456In this example, we further studied five different boundary conditions: SSSF, SFSF, CCCF, CFCF, and CFSF. In this case, a 16 × 16 regular mesh is adopted for square plates with different boundary conditions. The square of the nondimensional frequency parameterλ corresponding to the first four lowest frequencies is shown in Table 6. As a result, the GS-DSG method is almost superior to the other three methods and is consistent with the exact solution of all frequencies in this problem.Table 6
A nondimensional frequency parameterλ2 of the square plate (t/a) = 0.005) with various boundary conditions.
Plate typeMethodMode sequence number1234SSSFFEM-DSG11.706627.937642.134760.7093NS-FEM11.508423.168823.216826.8704NS-DSG11.396026.612840.293556.4514GS-DSG11.549827.091940.714557.6083Reference [74]11.68527.75641.19759.066SFSFFEM-DSG9.679716.111237.038439.8117NS-FEM9.549115.728823.123023.1655NS-DSG9.575715.298934.826438.2995GS-DSG9.599615.746835.494338.5924Reference [74]9.63116.13536.72638.945CCCFFEM-DSG24.306541.262165.672280.5606NS-FEM23.426325.586525.658430.8248NS-DSG23.485839.034160.859973.7164GS-DSG23.698139.485261.986074.9583Reference [74]24.02040.03963.49376.761CFCFFEM-DSG22.444426.989045.168663.3315NS-FEM21.754824.781124.821725.7548NS-DSG21.776825.856942.446558.8471GS-DSG21.963126.131742.909759.9521Reference [74]22.27226.52943.66464.466CFSFFEM-DSG15.325120.797740.543350.8558NS-FEM14.995820.082623.908923.9517NS-DSG15.026719.896538.104848.1404GS-DSG15.103920.250338.707148.7522Reference [74]15.28520.67339.88249.500


#### 4.2.2. Elliptical Plate

In this section, a simply supported elliptical plate is considered. The geometric parameters of the plate are shown in Figure11 with thickness t=0.5. The material properties are E=2.0×1011, v=0.3, and ρ=8000. Since the analytical solution is unavailable for this problem, the result calculated by the ABAQUS software with a very fine mesh (33345 nodes) is employed as a reference.Figure 11
Geometric illustration.To illustrate the benefits of triangular grids, we use an unstructured mesh layout with 446 nodes, as shown in Figure12. The first 12 natural modes solved by the proposed method are plotted in Table 7. The spurious nonzero energy modes are marked by the blue border, which indicates that the NS-FEM method has the time instability problem. The natural frequencies corresponding to the modes together with other three method solutions are listed in Table 8. Natural frequencies marked by a black border in Table 8 denote the spurious nonzero modes. The relative errors of natural frequencies solved by GS-DSG together with solutions obtained using other three methods are plotted in Figure 13. The following conclusions can be drawn from the results: (1) the results obtained by GS-DSG calculation have no spurious nonzero energy modes and therefore no temporal instability problem; (2) the natural frequencies obtained using the FEM-DSG method are all higher than the reference solutions, with accuracy severely decreasing with increasing frequencies; (3) NS-DSG has been improved in the mode calculation of this problem, but its frequency calculation accuracy is still slightly lower than that of GS-DSG. Especially with increasing frequencies, the accuracy of NS-DSG decreases, whereas the presented method maintains accuracy not only for lower frequencies, but also for higher ones; (4) the presented method provides more accurate natural frequencies than FEM-DSG, NS-FEM, and NS-DSG and can effectively eliminate singular modes.Figure 12
Domain discretization.Table 7
First modes of the elliptical plate solved by different methods.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGReference123456789101112Table 8
First 12 parameterized natural frequenciesω of the elliptical plate with 446 nodes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGReference16.60936.39516.4636.48576.3595211.9388.58711.51311.57211.369319.6119.15118.49618.66218.425423.9849.57922.04222.38822.203529.91610.67127.36727.78227.628633.24410.90329.87330.40230.158743.22411.27338.24139.08939.012844.83311.66139.06439.93439.874954.67213.26845.67447.26947.5001059.02613.96349.71951.11951.4211159.87115.24650.81852.40552.5821269.46515.65556.69058.61858.565Figure 13
Relative errors for the first 12 parameterized natural frequencies of the elliptical plate.


### 4.3. Buckling Analysis

In the following examples, we use the proposed method to study the critical buckling load of rectangular plates with different length-width ratios and different edge loads, as shown in Figure14. For all cases considered here, the nondimensional buckling load factor is defined as K=λcrb2/π2D [35], where the edge width of the plate is expressed as b, the critical buckling load is expressed as λcr, and the bending stiffness is expressed as D. For the material parameters, E=2.0×1011 and v=0.3.Figure 14
Rectangular plates: (a) model; (b) axial compression; (c) biaxial compression; (d) shear in-plane.
(a)(b)(c)(d)


#### 4.3.1. Rectangular Plates Subjected to Different Edge Loading

Firstly, a square plate with thicknesst under a simply supported boundary condition is considered. Uniaxial compression (UC), biaxial compression (BC), and shear in-plane (SP) are studied. The problem domain is discretized with a 16 × 16 uniformly distributed triangular mesh.The critical buckling load factorK solved using the GS-DSG together with the solutions calculated by FEM-DSG, NS-FEM, and NS-DSG is listed in Table 9. The reference solutions are given by Timoshenko [72], Tham [75], and Wang [40]. From the table, it can be clearly seen that (1) the FEM-DSG and NS-FEM are slightly less accurate than the NS-DSG and GS-DSG due to the resulting overly stiff and soft system, respectively; (2) numerical tests demonstrate that GS-DSG can provide a relatively good accuracy of the critical buckling load factor compared with other three methods. Figure 15 shows the buckling modes of simply supported square plates under different edge loading. It is clear that the results of the GS-DSG can provide very stable buckling modes.Table 9
The buckling load factorK of a simply supported rectangular plate with the length-to-width ratio a/b=1, thickness-to-width ratio t/b=0.01, and various loading types.
LoadFEM-DSGNS-FEMNS-DSGGS-DSGWang (α=0.2) [40]Tham [75]Timoshenko [72]UC3.983.913.883.983.934.004.00BC1.991.961.941.991.962.002.00SP9.298.959.029.379.429.409.33Figure 15
The buckling modes of the simply supported square with different types of edge loading: (a) axial compression; (b) biaxial compression; (c) shear in-plane.
(a)(b)(c)


#### 4.3.2. Rectangular Plates with Different Length-to-Width and Thickness-to-Width Ratios

In this section, uniaxial compression rectangular plates with different length-to-width ratios and thickness-to-width ratios are considered. Simply supported boundary conditions are assumed. Four types of length-to-width ratios,a/b=1.0, 1.5, 2.0, and 2.5, and three types of thickness-to-width ratios, t/b = 0.05, 0.1, and 0.2, are investigated. The problem domain is discretized using a uniformly distributed triangular mesh with 16 elements along the y-axis.The critical buckling load factors solved by different schemes are given in Table10. The reference solutions are given by Liew [73], Kitipornchai [74], and Nguyen-Xuan [37]. As in the previous section, the GS-DSG method can obtain very good results compared with other methods. The axial buckling modes of simply supported rectangular plates with thickness-to-width ratios a/b=1.0; 1.5; 2.0; 2.5 are shown in Figure 16. Again, very stable buckling modes can be observed.Table 10
The factorK of axial buckling loads along the x axis of rectangular plates with various length-to-width ratios and various thickness-to-width ratios.
a/bt/bFEM-DSGNS-FEMNS-DSGGS-DSGMeshfree [68]Pb-2 Ritz [75]ES-DSG [37]1.00.053.83343.71673.76533.84143.92933.94443.94120.13.50123.35843.46373.51843.72703.78653.74020.22.63792.56622.65262.68393.14713.26373.12631.50.054.19984.04224.15794.23374.21164.25704.28520.13.80133.63083.78713.84193.89824.02503.98440.22.78792.70252.81762.84903.10323.30483.14612.00.053.88663.75973.85343.90883.86573.94443.98110.13.58673.43773.57173.61313.67973.78653.77110.22.76722.68432.79682.81433.07833.26373.14152.50.054.02733.88014.00494.06023.96004.06454.16910.13.68853.52553.68273.72523.73113.86833.89240.22.79182.70192.81672.84453.03063.24213.1234Figure 16
Axial buckling modes of simply supported rectangular plates with various length-to-width ratios: (a)a/b=1.0; (b) a/b=1.5; (c) a/b=2.0; (d) a/b=2.5.
(a)(b)(c)(d)


## 4.1. Static Analysis

### 4.1.1. Square Plate

Consider the model of a simply supported square plate subjected to a uniform loadq=1 as shown in Figure 3. The geometric and material parameters are length L=10 and thickness t=2.0; E=2.1×109 and v=0.3. Due to symmetry, only a quarter of the plate is modeled to reduce the computation cost, and uniform meshes are employed. The center deflection wc is normalized as w^=wcD/qL4, where D=Et3/121−v2 is the bending stiffness. In order to test the performance of the mentioned numerical methods, the numerical results obtained using the present method are compared with other three methods. The result calculated by the ABAQUS software is used for reference, using S4R elements and a large number (37,249) of nodes.Figure 3
A square plate with a simply supported boundary condition.Table1 shows the numerical results of the normalized center deflection. Figure 4 shows the relative error, and the label “mesh density” on the horizontal axis shows the number of cells on each side. Figure 5 shows the convergence status of the displacement error norm eu, where h is the average nodal spacing of the node distribution. From the results, it can be seen that the deflection obtained by NS-FEM and NS-DSG is larger than the reference solution, whereas the deflection solved by the proposed method and FEM-DSG is smaller than the reference solution. Meanwhile, the proposed method is more accurate than the others. By comparing the convergence rate of the methods as shown in Figure 5, the proposed method has a higher convergence rate than the others, as far as the average nodal spacing trail off is concerned.Table 1
Numerical results of normalized central deflectionw^ for the simply supported square plate subjected to uniform load.
MeshFEM-DSGNS-FEMNS-DSGGS-DSG8 × 80.0054420.0056780.0055670.00551916 × 160.0055160.0055850.0055540.00554024 × 240.0055320.0055630.0055500.00554432 × 320.0055380.0055550.0055490.005545Reference0.005545Figure 4
Central deflection.Figure 5
Convergence properties.


### 4.1.2. T-Shaped Plate

In this section, a T-shaped plate with clamped edges and subjected to two kinds of loads, i.e., uniform and concentrated loads, is analyzed to further examine the efficiency of the present method. The geometric parameters are shown in Figure6. Two thickness are considered, t=0.01 and t=0.5. Material parameters are E=2.1×109 and v=0.3. The uniform load applied to the entire plate is given by q=0.01, and the concentrated load applied to point A is taken as p=1. Due to the symmetry of the plate, only half of the model is studied to reduce the calculation cost.Figure 6
Geometry illustration of a T-plate.Figure7 shows the mesh model, which is discretized using 154 nodes with 243 triangular elements. Numerical results of the present method are compared with other three methods. Since the analytical solution is unavailable for this problem, the result calculated by the ABAQUS software with a very fine mesh (2751 nodes and 5200 elements) was used for reference. The deflections along the line OA are plotted as shown in Figures 8 and 9. From the results, it can be seen that, for the thick plate, the results are almost identical, and the result is close to the reference solution. That is, because the shear-locking phenomenon does not appear in the thick plate, all the four methods can obtain high accuracy and is hard to distinguish which is higher. However, for the thin plate, the difference in the accuracy of the four methods is obvious. The accuracy of FEM-DSG and NS-DGS is lower, and both NS-FEM and the proposed method can achieve high accuracy compared with the reference result.Figure 7
Domain discretization using triangular.Figure 8
Deflections along the line OA with different thickness (uniformly distributed load): (a)t=0.01; (b) t=0.5.
(a)(b)Figure 9
Deflections along the line OA with different thickness (concentrated load): (a)t=0.01; (b) t=0.5.
(a)(b)


## 4.1.1. Square Plate

Consider the model of a simply supported square plate subjected to a uniform loadq=1 as shown in Figure 3. The geometric and material parameters are length L=10 and thickness t=2.0; E=2.1×109 and v=0.3. Due to symmetry, only a quarter of the plate is modeled to reduce the computation cost, and uniform meshes are employed. The center deflection wc is normalized as w^=wcD/qL4, where D=Et3/121−v2 is the bending stiffness. In order to test the performance of the mentioned numerical methods, the numerical results obtained using the present method are compared with other three methods. The result calculated by the ABAQUS software is used for reference, using S4R elements and a large number (37,249) of nodes.Figure 3
A square plate with a simply supported boundary condition.Table1 shows the numerical results of the normalized center deflection. Figure 4 shows the relative error, and the label “mesh density” on the horizontal axis shows the number of cells on each side. Figure 5 shows the convergence status of the displacement error norm eu, where h is the average nodal spacing of the node distribution. From the results, it can be seen that the deflection obtained by NS-FEM and NS-DSG is larger than the reference solution, whereas the deflection solved by the proposed method and FEM-DSG is smaller than the reference solution. Meanwhile, the proposed method is more accurate than the others. By comparing the convergence rate of the methods as shown in Figure 5, the proposed method has a higher convergence rate than the others, as far as the average nodal spacing trail off is concerned.Table 1
Numerical results of normalized central deflectionw^ for the simply supported square plate subjected to uniform load.
MeshFEM-DSGNS-FEMNS-DSGGS-DSG8 × 80.0054420.0056780.0055670.00551916 × 160.0055160.0055850.0055540.00554024 × 240.0055320.0055630.0055500.00554432 × 320.0055380.0055550.0055490.005545Reference0.005545Figure 4
Central deflection.Figure 5
Convergence properties.


## 4.1.2. T-Shaped Plate

In this section, a T-shaped plate with clamped edges and subjected to two kinds of loads, i.e., uniform and concentrated loads, is analyzed to further examine the efficiency of the present method. The geometric parameters are shown in Figure6. Two thickness are considered, t=0.01 and t=0.5. Material parameters are E=2.1×109 and v=0.3. The uniform load applied to the entire plate is given by q=0.01, and the concentrated load applied to point A is taken as p=1. Due to the symmetry of the plate, only half of the model is studied to reduce the calculation cost.Figure 6
Geometry illustration of a T-plate.Figure7 shows the mesh model, which is discretized using 154 nodes with 243 triangular elements. Numerical results of the present method are compared with other three methods. Since the analytical solution is unavailable for this problem, the result calculated by the ABAQUS software with a very fine mesh (2751 nodes and 5200 elements) was used for reference. The deflections along the line OA are plotted as shown in Figures 8 and 9. From the results, it can be seen that, for the thick plate, the results are almost identical, and the result is close to the reference solution. That is, because the shear-locking phenomenon does not appear in the thick plate, all the four methods can obtain high accuracy and is hard to distinguish which is higher. However, for the thin plate, the difference in the accuracy of the four methods is obvious. The accuracy of FEM-DSG and NS-DGS is lower, and both NS-FEM and the proposed method can achieve high accuracy compared with the reference result.Figure 7
Domain discretization using triangular.Figure 8
Deflections along the line OA with different thickness (uniformly distributed load): (a)t=0.01; (b) t=0.5.
(a)(b)Figure 9
Deflections along the line OA with different thickness (concentrated load): (a)t=0.01; (b) t=0.5.
(a)(b)


## 4.2. Free Vibration Analysis

In this section, numerical examples of free vibration for various plates are given. The nondimensional frequency parameter is normalized byλ=ωa2ρt/D1/2 [35], where ω is the circle frequency value, a is the geometry size given in each problem, ρ is the mass density, t is the thickness, and D=Et3/121−v2 is the bending stiffness.


### 4.2.1. Square Plates

Square plates of lengtha, width b, and thickness t are considered. The material parameters are E=2.0×1011, v=0.3, and ρ=8000. The plate is modeled with uniform meshes of 4, 8, 16, and 24 elements each side. The boundary conditions are simply supported (S), clamped (C), and free (F). SSSS means that all four sides are simply supported, and the others are similar.First, the SSSS plate corresponding to length-to-width ratiosa/b=1 is considered. The thickness-to-length for the thin plate is t/a=0.005 and for the thick plate is t/a=0.1, respectively. Figures 10(a), 10(c), and 10(d) show the geometry of the plate and its mesh grid, respectively. Table 2 shows the values of the nondimensional frequency parameter λ corresponding to the six frequencies using 4 × 4, 8 × 8, 16 × 16, and 24 × 24 meshes. It is observed that the accuracy of the presented method increases with the decreasing size of the mesh elements, and the results of GS-DSG agree well with the analytical ones. For the same mesh, the presented method is more accurate than FEM-DSG, NS-FEM, and NS-DSG elements for both thin and thick plates. Table 3 shows four methods for the first six modes under the 16×16 mesh. The spurious nonzero energy mode is marked with the blue wireframe. It can be clearly found in the table that NS-FEM has spurious nonzero energy modes, that is, severe time instability. The mode obtained by the NS-DSG method cannot eliminate spurious nonzero energy modes completely even using the discrete shear gap technology. In contrast, the advantage of the GS-DSG method is particularly obvious, and there are no spurious nonzero energy modes, which indicates its stability in the time domain.Figure 10
(a) Supported plate, (b) clamped plate, and (c, d) triangular meshes.
(a)(b)(c)(d)Table 2
A nondimensional frequency parameterλ of an SSSS plate (a/b) = 1).
t/aMethodMode sequence number1234560.005FEM-DSG4.80678.50389.263112.335316.657317.94004.51907.35297.51239.705811.207111.27984.45457.09437.13449.074810.214810.21804.44407.04997.06808.959110.048110.0521NS-FEM3.85234.61804.66965.48035.58055.70684.28734.74644.95285.42435.97026.24744.39314.81824.82385.28625.84636.10544.41564.76594.76605.48646.28916.2921NS-DSG4.19526.24896.37817.40327.51548.29444.31916.74456.75738.43709.31319.37584.38536.92346.92488.70099.73669.74624.40536.96586.96678.77539.83089.8348GS-DSG4.28936.78886.93138.38169.58179.60354.36376.88886.89668.67079.69919.74034.41016.96876.97298.78469.82609.83934.42286.99186.99458.82499.87389.8807Reference [72]4.4437.0257.0258.8869.9359.9350.1FEM-DSG4.67007.97138.581110.326110.601110.60174.39646.99717.13288.96159.753210.18494.32006.75656.78748.42019.41199.41854.30126.70866.72168.30939.27339.2794NS-FEM3.73664.41024.43975.14035.22705.37014.14384.43894.58625.01525.45765.66164.23274.45344.45534.89235.34345.53594.25344.40574.40585.02805.69715.6991NS-DSG4.10265.94186.03666.75216.78437.15404.21756.46116.47287.38217.67957.79824.26756.61586.61587.62057.84467.84514.27596.64396.64587.70187.86037.8604GS-DSG4.18636.46436.54167.74828.50518.64284.25466.58596.58778.08668.94138.96374.28166.64636.64878.18179.09679.09794.28306.65766.65938.20029.13149.1365Reference [72]4.376.746.748.359.229.22Table 3
First 6 modes of the square plate (SSSS) solved by different methods with16×16 meshes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGt/a = 0.005123456t/a = 0.1123456The second problem is that the mesh of CCCC square plate as shown in Figure10(b) is the same as that of the SSSS plate. The nondimensional frequency parameter λ corresponding to the first six frequencies of the CCCC plate is shown in Table 4. The corresponding modes under 16×16 are given in Table 5, and the spurious nonzero energy mode is marked with the blue wireframe. It can be found again that the GS-DSG method is superior to the other three methods.Table 4
A nondimensional frequency parameterλ of a CCCC plate (a/b = 1).
t/aMethodMode sequence number1234560.005FEM-DSG8.674314.607118.478122.121528.389129.17946.38399.46939.773212.410313.889214.23656.07408.74228.809710.808911.947311.98236.02848.63558.655410.565311.662911.6897NS-FEM4.866227.299427.403030.548531.246833.71825.51135.63695.69086.15736.73456.85785.14635.15675.65975.89896.21536.54284.96374.96395.74115.95206.51026.5153NS-DSG6.23997.783314.514316.878718.550419.45848.84888.09458.26089.961310.415710.82715.93568.39818.417810.174311.114411.15765.96608.48108.487910.277011.289011.3201GS-DSG7.139110.478814.736217.063918.806019.53715.99028.46718.617410.622011.286811.64865.96948.48568.501910.311211.305811.33985.98078.51918.524910.336311.372411.4008Reference [73]5.9998.5688.56810.40711.47211.4980.1FEM-DSG7.025510.228411.353911.385011.834012.80155.95768.40118.648910.377411.208411.36915.76097.99688.05259.583810.427810.48035.72747.92717.95149.436910.251710.2995NS-FEM4.61546.34326.78676.85687.08587.64054.81194.97495.37965.47435.91445.93804.55014.55725.11115.51085.62065.72904.44944.44945.25455.66585.75245.7533NS-DSG5.71317.06307.84268.22228.38018.46985.55407.48527.55087.76387.91538.25005.65447.75587.76757.92947.93108.03545.67977.81907.82427.91327.91337.9705GS-DSG6.13058.39878.998610.205910.309510.35115.65057.73777.76839.23859.85429.89665.67877.81307.82259.24779.998510.04215.69027.84337.84799.283110.057510.1036Reference [73]5.717.887.889.3310.1310.18Table 5
First 6 modes of the square plate (CCCC) solved by different methods with16×16 meshes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGt/a = 0.005123456t/a = 0.1123456In this example, we further studied five different boundary conditions: SSSF, SFSF, CCCF, CFCF, and CFSF. In this case, a 16 × 16 regular mesh is adopted for square plates with different boundary conditions. The square of the nondimensional frequency parameterλ corresponding to the first four lowest frequencies is shown in Table 6. As a result, the GS-DSG method is almost superior to the other three methods and is consistent with the exact solution of all frequencies in this problem.Table 6
A nondimensional frequency parameterλ2 of the square plate (t/a) = 0.005) with various boundary conditions.
Plate typeMethodMode sequence number1234SSSFFEM-DSG11.706627.937642.134760.7093NS-FEM11.508423.168823.216826.8704NS-DSG11.396026.612840.293556.4514GS-DSG11.549827.091940.714557.6083Reference [74]11.68527.75641.19759.066SFSFFEM-DSG9.679716.111237.038439.8117NS-FEM9.549115.728823.123023.1655NS-DSG9.575715.298934.826438.2995GS-DSG9.599615.746835.494338.5924Reference [74]9.63116.13536.72638.945CCCFFEM-DSG24.306541.262165.672280.5606NS-FEM23.426325.586525.658430.8248NS-DSG23.485839.034160.859973.7164GS-DSG23.698139.485261.986074.9583Reference [74]24.02040.03963.49376.761CFCFFEM-DSG22.444426.989045.168663.3315NS-FEM21.754824.781124.821725.7548NS-DSG21.776825.856942.446558.8471GS-DSG21.963126.131742.909759.9521Reference [74]22.27226.52943.66464.466CFSFFEM-DSG15.325120.797740.543350.8558NS-FEM14.995820.082623.908923.9517NS-DSG15.026719.896538.104848.1404GS-DSG15.103920.250338.707148.7522Reference [74]15.28520.67339.88249.500


### 4.2.2. Elliptical Plate

In this section, a simply supported elliptical plate is considered. The geometric parameters of the plate are shown in Figure11 with thickness t=0.5. The material properties are E=2.0×1011, v=0.3, and ρ=8000. Since the analytical solution is unavailable for this problem, the result calculated by the ABAQUS software with a very fine mesh (33345 nodes) is employed as a reference.Figure 11
Geometric illustration.To illustrate the benefits of triangular grids, we use an unstructured mesh layout with 446 nodes, as shown in Figure12. The first 12 natural modes solved by the proposed method are plotted in Table 7. The spurious nonzero energy modes are marked by the blue border, which indicates that the NS-FEM method has the time instability problem. The natural frequencies corresponding to the modes together with other three method solutions are listed in Table 8. Natural frequencies marked by a black border in Table 8 denote the spurious nonzero modes. The relative errors of natural frequencies solved by GS-DSG together with solutions obtained using other three methods are plotted in Figure 13. The following conclusions can be drawn from the results: (1) the results obtained by GS-DSG calculation have no spurious nonzero energy modes and therefore no temporal instability problem; (2) the natural frequencies obtained using the FEM-DSG method are all higher than the reference solutions, with accuracy severely decreasing with increasing frequencies; (3) NS-DSG has been improved in the mode calculation of this problem, but its frequency calculation accuracy is still slightly lower than that of GS-DSG. Especially with increasing frequencies, the accuracy of NS-DSG decreases, whereas the presented method maintains accuracy not only for lower frequencies, but also for higher ones; (4) the presented method provides more accurate natural frequencies than FEM-DSG, NS-FEM, and NS-DSG and can effectively eliminate singular modes.Figure 12
Domain discretization.Table 7
First modes of the elliptical plate solved by different methods.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGReference123456789101112Table 8
First 12 parameterized natural frequenciesω of the elliptical plate with 446 nodes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGReference16.60936.39516.4636.48576.3595211.9388.58711.51311.57211.369319.6119.15118.49618.66218.425423.9849.57922.04222.38822.203529.91610.67127.36727.78227.628633.24410.90329.87330.40230.158743.22411.27338.24139.08939.012844.83311.66139.06439.93439.874954.67213.26845.67447.26947.5001059.02613.96349.71951.11951.4211159.87115.24650.81852.40552.5821269.46515.65556.69058.61858.565Figure 13
Relative errors for the first 12 parameterized natural frequencies of the elliptical plate.


## 4.2.1. Square Plates

Square plates of lengtha, width b, and thickness t are considered. The material parameters are E=2.0×1011, v=0.3, and ρ=8000. The plate is modeled with uniform meshes of 4, 8, 16, and 24 elements each side. The boundary conditions are simply supported (S), clamped (C), and free (F). SSSS means that all four sides are simply supported, and the others are similar.First, the SSSS plate corresponding to length-to-width ratiosa/b=1 is considered. The thickness-to-length for the thin plate is t/a=0.005 and for the thick plate is t/a=0.1, respectively. Figures 10(a), 10(c), and 10(d) show the geometry of the plate and its mesh grid, respectively. Table 2 shows the values of the nondimensional frequency parameter λ corresponding to the six frequencies using 4 × 4, 8 × 8, 16 × 16, and 24 × 24 meshes. It is observed that the accuracy of the presented method increases with the decreasing size of the mesh elements, and the results of GS-DSG agree well with the analytical ones. For the same mesh, the presented method is more accurate than FEM-DSG, NS-FEM, and NS-DSG elements for both thin and thick plates. Table 3 shows four methods for the first six modes under the 16×16 mesh. The spurious nonzero energy mode is marked with the blue wireframe. It can be clearly found in the table that NS-FEM has spurious nonzero energy modes, that is, severe time instability. The mode obtained by the NS-DSG method cannot eliminate spurious nonzero energy modes completely even using the discrete shear gap technology. In contrast, the advantage of the GS-DSG method is particularly obvious, and there are no spurious nonzero energy modes, which indicates its stability in the time domain.Figure 10
(a) Supported plate, (b) clamped plate, and (c, d) triangular meshes.
(a)(b)(c)(d)Table 2
A nondimensional frequency parameterλ of an SSSS plate (a/b) = 1).
t/aMethodMode sequence number1234560.005FEM-DSG4.80678.50389.263112.335316.657317.94004.51907.35297.51239.705811.207111.27984.45457.09437.13449.074810.214810.21804.44407.04997.06808.959110.048110.0521NS-FEM3.85234.61804.66965.48035.58055.70684.28734.74644.95285.42435.97026.24744.39314.81824.82385.28625.84636.10544.41564.76594.76605.48646.28916.2921NS-DSG4.19526.24896.37817.40327.51548.29444.31916.74456.75738.43709.31319.37584.38536.92346.92488.70099.73669.74624.40536.96586.96678.77539.83089.8348GS-DSG4.28936.78886.93138.38169.58179.60354.36376.88886.89668.67079.69919.74034.41016.96876.97298.78469.82609.83934.42286.99186.99458.82499.87389.8807Reference [72]4.4437.0257.0258.8869.9359.9350.1FEM-DSG4.67007.97138.581110.326110.601110.60174.39646.99717.13288.96159.753210.18494.32006.75656.78748.42019.41199.41854.30126.70866.72168.30939.27339.2794NS-FEM3.73664.41024.43975.14035.22705.37014.14384.43894.58625.01525.45765.66164.23274.45344.45534.89235.34345.53594.25344.40574.40585.02805.69715.6991NS-DSG4.10265.94186.03666.75216.78437.15404.21756.46116.47287.38217.67957.79824.26756.61586.61587.62057.84467.84514.27596.64396.64587.70187.86037.8604GS-DSG4.18636.46436.54167.74828.50518.64284.25466.58596.58778.08668.94138.96374.28166.64636.64878.18179.09679.09794.28306.65766.65938.20029.13149.1365Reference [72]4.376.746.748.359.229.22Table 3
First 6 modes of the square plate (SSSS) solved by different methods with16×16 meshes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGt/a = 0.005123456t/a = 0.1123456The second problem is that the mesh of CCCC square plate as shown in Figure10(b) is the same as that of the SSSS plate. The nondimensional frequency parameter λ corresponding to the first six frequencies of the CCCC plate is shown in Table 4. The corresponding modes under 16×16 are given in Table 5, and the spurious nonzero energy mode is marked with the blue wireframe. It can be found again that the GS-DSG method is superior to the other three methods.Table 4
A nondimensional frequency parameterλ of a CCCC plate (a/b = 1).
t/aMethodMode sequence number1234560.005FEM-DSG8.674314.607118.478122.121528.389129.17946.38399.46939.773212.410313.889214.23656.07408.74228.809710.808911.947311.98236.02848.63558.655410.565311.662911.6897NS-FEM4.866227.299427.403030.548531.246833.71825.51135.63695.69086.15736.73456.85785.14635.15675.65975.89896.21536.54284.96374.96395.74115.95206.51026.5153NS-DSG6.23997.783314.514316.878718.550419.45848.84888.09458.26089.961310.415710.82715.93568.39818.417810.174311.114411.15765.96608.48108.487910.277011.289011.3201GS-DSG7.139110.478814.736217.063918.806019.53715.99028.46718.617410.622011.286811.64865.96948.48568.501910.311211.305811.33985.98078.51918.524910.336311.372411.4008Reference [73]5.9998.5688.56810.40711.47211.4980.1FEM-DSG7.025510.228411.353911.385011.834012.80155.95768.40118.648910.377411.208411.36915.76097.99688.05259.583810.427810.48035.72747.92717.95149.436910.251710.2995NS-FEM4.61546.34326.78676.85687.08587.64054.81194.97495.37965.47435.91445.93804.55014.55725.11115.51085.62065.72904.44944.44945.25455.66585.75245.7533NS-DSG5.71317.06307.84268.22228.38018.46985.55407.48527.55087.76387.91538.25005.65447.75587.76757.92947.93108.03545.67977.81907.82427.91327.91337.9705GS-DSG6.13058.39878.998610.205910.309510.35115.65057.73777.76839.23859.85429.89665.67877.81307.82259.24779.998510.04215.69027.84337.84799.283110.057510.1036Reference [73]5.717.887.889.3310.1310.18Table 5
First 6 modes of the square plate (CCCC) solved by different methods with16×16 meshes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGt/a = 0.005123456t/a = 0.1123456In this example, we further studied five different boundary conditions: SSSF, SFSF, CCCF, CFCF, and CFSF. In this case, a 16 × 16 regular mesh is adopted for square plates with different boundary conditions. The square of the nondimensional frequency parameterλ corresponding to the first four lowest frequencies is shown in Table 6. As a result, the GS-DSG method is almost superior to the other three methods and is consistent with the exact solution of all frequencies in this problem.Table 6
A nondimensional frequency parameterλ2 of the square plate (t/a) = 0.005) with various boundary conditions.
Plate typeMethodMode sequence number1234SSSFFEM-DSG11.706627.937642.134760.7093NS-FEM11.508423.168823.216826.8704NS-DSG11.396026.612840.293556.4514GS-DSG11.549827.091940.714557.6083Reference [74]11.68527.75641.19759.066SFSFFEM-DSG9.679716.111237.038439.8117NS-FEM9.549115.728823.123023.1655NS-DSG9.575715.298934.826438.2995GS-DSG9.599615.746835.494338.5924Reference [74]9.63116.13536.72638.945CCCFFEM-DSG24.306541.262165.672280.5606NS-FEM23.426325.586525.658430.8248NS-DSG23.485839.034160.859973.7164GS-DSG23.698139.485261.986074.9583Reference [74]24.02040.03963.49376.761CFCFFEM-DSG22.444426.989045.168663.3315NS-FEM21.754824.781124.821725.7548NS-DSG21.776825.856942.446558.8471GS-DSG21.963126.131742.909759.9521Reference [74]22.27226.52943.66464.466CFSFFEM-DSG15.325120.797740.543350.8558NS-FEM14.995820.082623.908923.9517NS-DSG15.026719.896538.104848.1404GS-DSG15.103920.250338.707148.7522Reference [74]15.28520.67339.88249.500


## 4.2.2. Elliptical Plate

In this section, a simply supported elliptical plate is considered. The geometric parameters of the plate are shown in Figure11 with thickness t=0.5. The material properties are E=2.0×1011, v=0.3, and ρ=8000. Since the analytical solution is unavailable for this problem, the result calculated by the ABAQUS software with a very fine mesh (33345 nodes) is employed as a reference.Figure 11
Geometric illustration.To illustrate the benefits of triangular grids, we use an unstructured mesh layout with 446 nodes, as shown in Figure12. The first 12 natural modes solved by the proposed method are plotted in Table 7. The spurious nonzero energy modes are marked by the blue border, which indicates that the NS-FEM method has the time instability problem. The natural frequencies corresponding to the modes together with other three method solutions are listed in Table 8. Natural frequencies marked by a black border in Table 8 denote the spurious nonzero modes. The relative errors of natural frequencies solved by GS-DSG together with solutions obtained using other three methods are plotted in Figure 13. The following conclusions can be drawn from the results: (1) the results obtained by GS-DSG calculation have no spurious nonzero energy modes and therefore no temporal instability problem; (2) the natural frequencies obtained using the FEM-DSG method are all higher than the reference solutions, with accuracy severely decreasing with increasing frequencies; (3) NS-DSG has been improved in the mode calculation of this problem, but its frequency calculation accuracy is still slightly lower than that of GS-DSG. Especially with increasing frequencies, the accuracy of NS-DSG decreases, whereas the presented method maintains accuracy not only for lower frequencies, but also for higher ones; (4) the presented method provides more accurate natural frequencies than FEM-DSG, NS-FEM, and NS-DSG and can effectively eliminate singular modes.Figure 12
Domain discretization.Table 7
First modes of the elliptical plate solved by different methods.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGReference123456789101112Table 8
First 12 parameterized natural frequenciesω of the elliptical plate with 446 nodes.
ModesFEM-DSGNS-FEMNS-DSGGS-DSGReference16.60936.39516.4636.48576.3595211.9388.58711.51311.57211.369319.6119.15118.49618.66218.425423.9849.57922.04222.38822.203529.91610.67127.36727.78227.628633.24410.90329.87330.40230.158743.22411.27338.24139.08939.012844.83311.66139.06439.93439.874954.67213.26845.67447.26947.5001059.02613.96349.71951.11951.4211159.87115.24650.81852.40552.5821269.46515.65556.69058.61858.565Figure 13
Relative errors for the first 12 parameterized natural frequencies of the elliptical plate.


## 4.3. Buckling Analysis

In the following examples, we use the proposed method to study the critical buckling load of rectangular plates with different length-width ratios and different edge loads, as shown in Figure14. For all cases considered here, the nondimensional buckling load factor is defined as K=λcrb2/π2D [35], where the edge width of the plate is expressed as b, the critical buckling load is expressed as λcr, and the bending stiffness is expressed as D. For the material parameters, E=2.0×1011 and v=0.3.Figure 14
Rectangular plates: (a) model; (b) axial compression; (c) biaxial compression; (d) shear in-plane.
(a)(b)(c)(d)


### 4.3.1. Rectangular Plates Subjected to Different Edge Loading

Firstly, a square plate with thicknesst under a simply supported boundary condition is considered. Uniaxial compression (UC), biaxial compression (BC), and shear in-plane (SP) are studied. The problem domain is discretized with a 16 × 16 uniformly distributed triangular mesh.The critical buckling load factorK solved using the GS-DSG together with the solutions calculated by FEM-DSG, NS-FEM, and NS-DSG is listed in Table 9. The reference solutions are given by Timoshenko [72], Tham [75], and Wang [40]. From the table, it can be clearly seen that (1) the FEM-DSG and NS-FEM are slightly less accurate than the NS-DSG and GS-DSG due to the resulting overly stiff and soft system, respectively; (2) numerical tests demonstrate that GS-DSG can provide a relatively good accuracy of the critical buckling load factor compared with other three methods. Figure 15 shows the buckling modes of simply supported square plates under different edge loading. It is clear that the results of the GS-DSG can provide very stable buckling modes.Table 9
The buckling load factorK of a simply supported rectangular plate with the length-to-width ratio a/b=1, thickness-to-width ratio t/b=0.01, and various loading types.
LoadFEM-DSGNS-FEMNS-DSGGS-DSGWang (α=0.2) [40]Tham [75]Timoshenko [72]UC3.983.913.883.983.934.004.00BC1.991.961.941.991.962.002.00SP9.298.959.029.379.429.409.33Figure 15
The buckling modes of the simply supported square with different types of edge loading: (a) axial compression; (b) biaxial compression; (c) shear in-plane.
(a)(b)(c)


### 4.3.2. Rectangular Plates with Different Length-to-Width and Thickness-to-Width Ratios

In this section, uniaxial compression rectangular plates with different length-to-width ratios and thickness-to-width ratios are considered. Simply supported boundary conditions are assumed. Four types of length-to-width ratios,a/b=1.0, 1.5, 2.0, and 2.5, and three types of thickness-to-width ratios, t/b = 0.05, 0.1, and 0.2, are investigated. The problem domain is discretized using a uniformly distributed triangular mesh with 16 elements along the y-axis.The critical buckling load factors solved by different schemes are given in Table10. The reference solutions are given by Liew [73], Kitipornchai [74], and Nguyen-Xuan [37]. As in the previous section, the GS-DSG method can obtain very good results compared with other methods. The axial buckling modes of simply supported rectangular plates with thickness-to-width ratios a/b=1.0; 1.5; 2.0; 2.5 are shown in Figure 16. Again, very stable buckling modes can be observed.Table 10
The factorK of axial buckling loads along the x axis of rectangular plates with various length-to-width ratios and various thickness-to-width ratios.
a/bt/bFEM-DSGNS-FEMNS-DSGGS-DSGMeshfree [68]Pb-2 Ritz [75]ES-DSG [37]1.00.053.83343.71673.76533.84143.92933.94443.94120.13.50123.35843.46373.51843.72703.78653.74020.22.63792.56622.65262.68393.14713.26373.12631.50.054.19984.04224.15794.23374.21164.25704.28520.13.80133.63083.78713.84193.89824.02503.98440.22.78792.70252.81762.84903.10323.30483.14612.00.053.88663.75973.85343.90883.86573.94443.98110.13.58673.43773.57173.61313.67973.78653.77110.22.76722.68432.79682.81433.07833.26373.14152.50.054.02733.88014.00494.06023.96004.06454.16910.13.68853.52553.68273.72523.73113.86833.89240.22.79182.70192.81672.84453.03063.24213.1234Figure 16
Axial buckling modes of simply supported rectangular plates with various length-to-width ratios: (a)a/b=1.0; (b) a/b=1.5; (c) a/b=2.0; (d) a/b=2.5.
(a)(b)(c)(d)


## 4.3.1. Rectangular Plates Subjected to Different Edge Loading

Firstly, a square plate with thicknesst under a simply supported boundary condition is considered. Uniaxial compression (UC), biaxial compression (BC), and shear in-plane (SP) are studied. The problem domain is discretized with a 16 × 16 uniformly distributed triangular mesh.The critical buckling load factorK solved using the GS-DSG together with the solutions calculated by FEM-DSG, NS-FEM, and NS-DSG is listed in Table 9. The reference solutions are given by Timoshenko [72], Tham [75], and Wang [40]. From the table, it can be clearly seen that (1) the FEM-DSG and NS-FEM are slightly less accurate than the NS-DSG and GS-DSG due to the resulting overly stiff and soft system, respectively; (2) numerical tests demonstrate that GS-DSG can provide a relatively good accuracy of the critical buckling load factor compared with other three methods. Figure 15 shows the buckling modes of simply supported square plates under different edge loading. It is clear that the results of the GS-DSG can provide very stable buckling modes.Table 9
The buckling load factorK of a simply supported rectangular plate with the length-to-width ratio a/b=1, thickness-to-width ratio t/b=0.01, and various loading types.
LoadFEM-DSGNS-FEMNS-DSGGS-DSGWang (α=0.2) [40]Tham [75]Timoshenko [72]UC3.983.913.883.983.934.004.00BC1.991.961.941.991.962.002.00SP9.298.959.029.379.429.409.33Figure 15
The buckling modes of the simply supported square with different types of edge loading: (a) axial compression; (b) biaxial compression; (c) shear in-plane.
(a)(b)(c)


## 4.3.2. Rectangular Plates with Different Length-to-Width and Thickness-to-Width Ratios

In this section, uniaxial compression rectangular plates with different length-to-width ratios and thickness-to-width ratios are considered. Simply supported boundary conditions are assumed. Four types of length-to-width ratios,a/b=1.0, 1.5, 2.0, and 2.5, and three types of thickness-to-width ratios, t/b = 0.05, 0.1, and 0.2, are investigated. The problem domain is discretized using a uniformly distributed triangular mesh with 16 elements along the y-axis.The critical buckling load factors solved by different schemes are given in Table10. The reference solutions are given by Liew [73], Kitipornchai [74], and Nguyen-Xuan [37]. As in the previous section, the GS-DSG method can obtain very good results compared with other methods. The axial buckling modes of simply supported rectangular plates with thickness-to-width ratios a/b=1.0; 1.5; 2.0; 2.5 are shown in Figure 16. Again, very stable buckling modes can be observed.Table 10
The factorK of axial buckling loads along the x axis of rectangular plates with various length-to-width ratios and various thickness-to-width ratios.
a/bt/bFEM-DSGNS-FEMNS-DSGGS-DSGMeshfree [68]Pb-2 Ritz [75]ES-DSG [37]1.00.053.83343.71673.76533.84143.92933.94443.94120.13.50123.35843.46373.51843.72703.78653.74020.22.63792.56622.65262.68393.14713.26373.12631.50.054.19984.04224.15794.23374.21164.25704.28520.13.80133.63083.78713.84193.89824.02503.98440.22.78792.70252.81762.84903.10323.30483.14612.00.053.88663.75973.85343.90883.86573.94443.98110.13.58673.43773.57173.61313.67973.78653.77110.22.76722.68432.79682.81433.07833.26373.14152.50.054.02733.88014.00494.06023.96004.06454.16910.13.68853.52553.68273.72523.73113.86833.89240.22.79182.70192.81672.84453.03063.24213.1234Figure 16
Axial buckling modes of simply supported rectangular plates with various length-to-width ratios: (a)a/b=1.0; (b) a/b=1.5; (c) a/b=2.0; (d) a/b=2.5.
(a)(b)(c)(d)


## 5. Conclusions

In this work, a GS-DSG method is formulated for Reissner–Mindlin plate problems in elastic-static, free vibration, and buckling analyses using 3-node triangular elements. The stabilization term is added to the discretized system equations by applying the smoothed Galerkin weak form. Through the formulations and numerical examples, the accuracy of the proposed method is demonstrated. Some concluding remarks can be drawn as follows:(1)
Several numerical examples indicate that GS-DSG is temporal stable for both free vibration and buckling problems(2)
In elastic-static analysis, the GS-DSG is free of shear locking and has higher accuracy in the displacement field compared with the FEM-DSG, NS-FEM, and NS-DSG methods(3)
In free vibration and buckling analyses, the GS-DSG effectively eliminates the spurious nonzero energy modes produced by the original NS-FEM and NS-DSG and provides a relatively accurate prediction of natural frequencies compared with other methods.



---

*Source: 8881983-2020-08-27.xml*