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.

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 -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 [25], strain (ANS) methods [6, 7], and also the methods in [816]. 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 [2044] 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) [4547] 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, 6267].

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:where , , and are the displacements of the plate midplane in the , , and directions and and denote the rotations with respect to and axes, respectively, as shown in Figure 1.

The relevant strain vector can be written in terms of the midplane deformations of equations (1)–(3), which giveswhere , , and are the membrane strain, the bending strain (curvature), and the shear strain, respectively:

By applying the principle of virtual work, the weak form can be stated as follows:where is the membrane stiffness constitutive coefficients, represents the bending stiffness constitutive coefficients, and denotes the transverse shear stiffness constitutive coefficients defined asin which is the shear modulus, is the shear correction factor, and the matrix contains the constitutive coefficients:where and 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 is the variation of the displacement field and is the inertia matrix containing the mass density and thickness :

For the buckling analysis, when the plate is subjected to in-plane prebuckling stresses , the corresponding weak form can be reformulated aswhere

Equation (13) can be rewritten in a compact form asin whichwhere “” represents the derivative of with respect to x.

2.2. Discrete Formulation for Reissner–Mindlin Plates

The bounded domain is discretized into triangular elements such that and , . 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 in the element is interpolated. The same shape functions are used for both displacements and rotations aswhere is the total number of nodes, is the generalized nodal displacement at node , and is a diagonal matrix of shape functions given by

Substituting equation (17) into (5), the membrane strain , the bending strain , and the shear strain can be written asin whichwhere the subscript .

From equations (14)–(16), the geometrical strain can be written asin which

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 denotes the vector of global nodal displacement at all of the nodes and is the force vector (including forces and torques) defined aswhere and denote the distributed load and prescribed boundary load, respectively.

In equation (27), the global stiffness matrix can be expressed asin which

For free vibration, we havewhere denotes the natural frequency and is the global inertia matrix:

For the buckling analysis, we havewherewhich is the geometrical stiffness matrix, and 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 in Section 2.2, the structural stiffness matrix is composed of three parts, namely, , , and . On the one hand, for and , we take the method of gradient stable node-based smoothed integration to avoid temporal instability and spatial instability; for , 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 that 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

Clearly, if the interpolated strain can be somehow constructed from the displacement and meet the following orthogonality condition:

Then, 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 and expressed as a finite element form:where is the shape function corresponding to strain in element I and 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 is the space base which is one order lower than the space for the displacement .

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 . Then, the following equations can be obtained:

In the nodal integration scheme, node is the only point in , and there exists only one unknown . Hence, equations (42) and (43) cannot be satisfied at the same time. To this end, we introduce the derivatives of the function by means of Taylor’s expansion such that and are introduced and can serve as the other two unknowns. Taylor’s expansion for can be formulated aswhere H.O.T means higher-order terms.

Substitution of equations (56)–(58) into (42) and (43) leads towhere is the area of and

By solving equation (47), the corrected nodal derivative and its derivatives, i.e., and , are obtained. Following the same derivation, Taylor’s expansion for is

The equation for y-derivatives can be written as

To simplify the calculation, the equivalent circle domain can be assumed for the subdomain . The area , first moments and , and second moments of inertia , , and 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.

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 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 and for , in which is the total number of nodes of the problem domain.

From equations (44)–(46) and equations (49)–(51), the smoothed strain and in surrounding node can be expressed aswherein which

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 as i, j, and k, respectively. The shear strain can be given aswhere and are the discrete shear gaps at the node given bywhere

From equations (64)–(70), the shear strain in each element can be rewritten asin whichwhere 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 should be corrective as :where is a positive constant [69, 71] and the characteristic length can be estimated as the diameter of the equation circle domain:

3.4. Discrete Formulation for GS-DSG Method

We now seek for a weak form solution of the generalized displacement field 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 is the global smoothed stiffness matrix assembled in the form of

The summation in equation (79) means an assembly process same as the practice in the FEM, and is the number of the nodes of the whole problem domain . is given from equation (32) by using the stabilized discrete shear gap method to eliminate the shear locking. and are the stiffness matrices associated with node given as

For free vibration, we have

For the buckling analysis, we havewhere denotes the geometrical stiffness matrix assembled in the form of

The nodal geometrical stiffness matrix in equation (84) can be calculated bywith

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 aswhere 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 , 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 load as shown in Figure 3. The geometric and material parameters are length and thickness ; and . Due to symmetry, only a quarter of the plate is modeled to reduce the computation cost, and uniform meshes are employed. The center deflection is normalized as , where 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.

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 , where 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.

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 Figure 6. Two thickness are considered, and . Material parameters are and . The uniform load applied to the entire plate is given by , and the concentrated load applied to point A is taken as . 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. 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.

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 [35], where is the circle frequency value, is the geometry size given in each problem, is the mass density, is the thickness, and is the bending stiffness.

4.2.1. Square Plates

Square plates of length , width , and thickness are considered. The material parameters are , , and . 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 ratios is considered. The thickness-to-length for the thin plate is and for the thick plate is , 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 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.

The 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. The nondimensional frequency parameter corresponding to the first six frequencies of the CCCC plate is shown in Table 4. The corresponding modes under 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. 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.

4.2.2. Elliptical Plate

In this section, a simply supported elliptical plate is considered. The geometric parameters of the plate are shown in Figure 11 with thickness . The material properties are , , and . 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. 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.

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 Figure 14. For all cases considered here, the nondimensional buckling load factor is defined as [35], where the edge width of the plate is expressed as , the critical buckling load is expressed as , and the bending stiffness is expressed as . For the material parameters, and .

4.3.1. 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 in-plane (SP) are studied. The problem domain is discretized with a 16 × 16 uniformly distributed triangular mesh.

The critical buckling load factor 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.

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, , 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 Table 10. 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 ; 1.5; 2.0; 2.5 are shown in Figure 16. Again, very stable buckling modes can be observed.

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.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This work was supported by the National Natural Science Foundation of China under Grant no. 11472137 and the Fundamental Research Funds for the Central Universities under Grant nos. 309181A8801 and 30919011204.