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

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.*emodified smoothed strain with gradient information is derived from the Hu–Washizu three-field variational principle, resulting in the stabilization terms in the system equations. *e 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. *e results demonstrate that the presented method is free of shear locking and can overcome the temporal instability issues, simultaneously obtaining excellent solutions.


Introduction
in-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. e analyses of plate structures in elastic-static, free vibration, and buckling fields stand for the key three aspects in their engineering application. ere 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 C 1 -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 C 0 -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][3][4][5], strain (ANS) methods [6,7], and also the methods in [8][9][10][11][12][13][14][15][16]. All the methods show excellent performance in reducing the shear-locking deficiency and increasing the solution accuracy. e 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  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][46][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. e 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. e proposed method can achieve appropriate system stiffness in train energy between FEM and NS-FEM solutions and indeed provide temporally stable results.
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. e 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. e 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. e 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.

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: in which G is the shear modulus, κ � 5/6 is the shear correction factor, and the matrix D 0 contains the constitutive coefficients: 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., where δd is the variation of the displacement field and m is the inertia matrix containing the mass density ρ and thickness t: For the buckling analysis, when the plate is subjected to in-plane prebuckling stresses σ 0 , the corresponding weak form can be reformulated as Shock and Vibration 3 where Equation (13) can be rewritten in a compact form as where "u 0,x " represents the derivative of u 0 with respect to x.

Discrete Formulation for Reissner-Mindlin
Plates. e bounded domain Ω is discretized into N e triangular elements such that Ω � ∪ N e e�1 Ω 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. e same shape functions are used for both displacements and rotations as where N n is the total number of nodes, is the generalized nodal displacement at node i, and N i (x) is a diagonal matrix of shape functions given by Substituting equation (17) into (5), the membrane strain ε m , the bending strain ε b , and the shear strain ε s can be written as in which where the subscript i � 1, 2, 3. From equations (14)-(16), the geometrical strain ε g can be written as in which 4 Shock and Vibration 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: 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 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 in which For free vibration, we have where ω denotes the natural frequency and M is the global inertia matrix: For the buckling analysis, we have where 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.

The Formulation of Gradient Stable Node-Based Smoothed Integration and Discrete Shear Gap Technique
According to the introduction in Section 2.2, the structural stiffness matrix is composed of three parts, namely, K m , K b , and K s . On the one hand, for K m and K b , we take the method of gradient stable node-based smoothed integration to avoid temporal instability and spatial instability; for K s , on the other hand, the discrete shear gap technique is employed to avoid the shear-locking problem. e formula for the smoothed integral is derived from the Hu-Washizu threefield variational principle as shown in Section 3.1, and its discrete form is given in Section 3.2.
e discrete shear technology is introduced in Section 3.3.

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 that u and σ are the displacement and assumed Cauchy stress, respectively, ε is the interpolated strain or smoothed strain, and ε is the actual strain. e Hu-Washizu three-field weak form for the elastic-static problem can be written as Clearly, if the interpolated strain ε can be somehow constructed from the displacement u and meet the following orthogonality condition: en, a form containing only independent variables can be obtained as simple as the classical one: 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: Shock and Vibration 5 where N I 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: where q(x) 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 q(x) � 1. en, the following equations can be obtained: In the nodal integration scheme, node L is the only point in Ω L , and there exists only one unknown N I,x (x L ). 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,xx (x L ) and N I,xy (x L ) are introduced and can serve as the other two unknowns. Taylor's expansion for N I,x can be formulated as where H.O.T means higher-order terms. Substitution of equations (56)-(58) into (42) and (43) leads to where By solving equation (47), the corrected nodal derivative N I,x (x L ) and its derivatives, i.e., N I,xx (x L ) and N I,xy (x L ), are obtained. Following the same derivation, Taylor's expansion for N I,y is e equation for y-derivatives can be written as To simplify the calculation, the equivalent circle domain can be assumed for the subdomain Ω L . e area A L , first moments I Lx and I Ly , and second moments of inertia I Lxx , I Lxy , and I Lyy of each nodal domain can be expressed as 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.

Discrete Form of Gradient Stabilized Nodal Integration.
In this part, the nodal integration formulation will be introduced. We can discretize the problem domain Ω with N e 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 Ω � ∪ N n k�1 Ω k and Ω i ∩ Ω j � ϕ for i ≠ j, in which N n is the total number of nodes of the problem domain.
From equations (44)- (46) and equations (49)-(51), the smoothed strain ε m (x) and ε b (x) in Ω L surrounding node L can be expressed as 6 Shock and Vibration where in which

Formulation of the Stabilized Discrete Shear Gap
Technique. e discrete shear gap method is adopted here to eliminate the shear locking. In each triangular element, the nodes are denoted anticlockwise as i, j, and k, respectively. e shear strain can be given as where Δw xm and Δw ym are the discrete shear gaps at the node m given by where Shock and Vibration 7 From equations (64)- (70), the shear strain ε s in each element can be rewritten as in which where A e 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. erefore, the transverse shear stiffness constitutive coefficients D s should be corrective as D s : where α is a positive constant [69,71] and the characteristic length h e can be estimated as the diameter of the equation circle domain:

Discrete Formulation for GS-DSG Method.
We now seek for a weak form solution of the generalized displacement field d that satisfies the following smoothed Galerkin weak form: Substituting equations (74)-(76) into (77), a set of discretized algebraic system equations can be obtained in the following matrix form: where K is the global smoothed stiffness matrix assembled in the form of e summation in equation (79) means an assembly process same as the practice in the FEM, and N n is the number of the nodes of the whole problem domain Ω. K (e) s is given from equation (32) For free vibration, we have For the buckling analysis, we have where K g denotes the geometrical stiffness matrix assembled in the form of with

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 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 as E, Poisson's ratio is expressed as v, and mass density is expressed as ρ.   9 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. e center deflection w c is normalized as w � w c D/qL 4 , where D � Et 3 /(12(1 − v 2 )) 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. e result calculated by the ABAQUS software is used for reference, using S4R elements and a large number (37,249) of nodes. Table 1 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 e u , 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. Shock and Vibration 9

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. e geometric parameters are shown in Figure 6. Two thickness are considered, t � 0.01 and t � 0.5. Material parameters are E � 2.1 × 10 9 and v � 0.3. e 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 7 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. e 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. at 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. e 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.

Free Vibration Analysis.
In this section, numerical examples of free vibration for various plates are given. e nondimensional frequency parameter is normalized by λ � (ωa 2 ���� ρt/D ) 1/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 � Et 3 /(12(1 − v 2 )) is the bending stiffness.     First, the SSSS plate corresponding to length-to-width ratios a/b � 1 is considered. e 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. e 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. e 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. e second problem is that the mesh of CCCC square plate as shown in Figure 10(b) is the same as that of the SSSS plate.
e nondimensional frequency parameter λ corresponding to the first six frequencies of the CCCC plate is shown in Table 4. e 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.
In 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. e 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.

Elliptical Plate.
In this section, a simply supported elliptical plate is considered. e geometric parameters of the plate are shown in Figure 11 with thickness t � 0.5. e material properties are E � 2.0 × 10 11 , 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.
To illustrate the benefits of triangular grids, we use an unstructured mesh layout with 446 nodes, as shown in Figure 12. e first 12 natural modes solved by the proposed method are plotted in Table 7. e spurious nonzero energy modes are marked by the blue border, which indicates that the NS-FEM method has the time instability problem. e 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. e relative errors of natural frequencies solved by GS-DSG together with solutions obtained using other three methods are plotted in Figure 13. e 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    Shock and Vibration    Shock and Vibration 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.

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             Figure 14. For all cases considered here, the nondimensional buckling load factor is defined as K � λ cr b 2 /(π 2 D) [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 × 10 11 and v � 0.3.

Rectangular Plates Subjected to Different Edge Loading.
Firstly, a square plate with thickness t under a simply supported boundary condition is considered. Uniaxial compression (UC), biaxial compression (BC), and shear inplane (SP) are studied. e problem domain is discretized with a 16 × 16 uniformly distributed triangular mesh. e critical buckling load factor K solved using the GS-DSG together with the solutions calculated by FEM-DSG, NS-FEM, and NS-DSG is listed in Table 9. e reference solutions are given by Timoshenko [72], am [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.

Rectangular Plates with Different Length-to-Width and
ickness-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-towidth 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. e problem domain is discretized using a uniformly distributed triangular mesh with 16 elements along the y-axis. e critical buckling load factors solved by different schemes are given in Table 10. e 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. e 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.

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. e stabilization term is added to the discretized system equations by applying the smoothed Galerkin weak form.
rough 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.

Shock and Vibration
Data Availability e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.