Abstract

This paper presents a gradient stable node-based smoothed finite element method (GS-FEM) which resolves the temporal instability of the node-based smoothed finite element method (NS-FEM) while significantly improving its accuracy. In the GS-FEM, the strain is expanded at the first order by Taylor expansion in a node-supported domain, and the strain gradient is then smoothed within each smoothing domain. Therefore, the stiffness matrix includes stable terms derived by the gradient of the strain. The GS-FEM model is softer than the FEM but stiffer than the NS-FEM and yields far more accurate results than the FEM-T3 or NS-FEM. It even has comparative accuracy compared with those of the FEM-Q4. The GS-FEM owns no spurious nonzero-energy modes and is thus temporally stable and well-suited for dynamic analyses. Additionally, the GS-FEM is demonstrated on static, free, and forced vibration example analyses of solids.

1. Introduction

The finite element method (FEM) [1, 2] is proven to be reliable and robust and hence has been applied successfully to countless practical engineering and science problems related to solid mechanics, fluid mechanics, and heat transfer. However, inherent disadvantages such as the discontinuous stress field at the interelement boundaries, especially for lower-order elements, render it ineffective in certain situations. The meshfree method has emerged in recent years as a workable alternative to the traditional FEM. It has been applied to large deformation, fracture propagation simulation, and impact-induced failure problems. Major advancements in meshfree methods include smoothed particle hydrodynamics (SPH) [3], element-free Galerkin method (EFG) [4], reproducing kernel particle method (RKPM) [5], meshfree local Petrov–Galerkin method (MLPG) [6], the radial point interpolation method (RPIM) [7], Hp-clouds method [8], moving particle finite element method (MPFEM) [9], and finite point method (FPM) [10]. Although meshfree methods are free of several of the FEM drawbacks, they are not totally ideal—for example, they present difficulties in essential boundary condition implementation, as well as high computational cost and complex trial function construction processes.

Hybrid schemes composed of meshfree and FEM methodologies may encompass the advantages of both while mitigating their respective shortcomings [1113], such as the partition of the unity finite element method (PUFEM) [14], generalized finite element method (GFEM) [15], FE-meshfree [16, 17], meshfree-enriched FEM (ME-FEM) [18], and RKPM [5, 19]. Zeng and Liu [20] combined the strain-smoothing technique of meshfree methods [21] and the existing FEM technology to establish smoothed finite element methods (S-FEMs) including the CS-SFEM for both 2D and 3D problems [22, 23], node-based SFEM (NS-FEM) for both 2D and 3D [24, 25], edge-based SFEM (ES-FEN) for 2D and 3D [26, 27], face-based SFEM (FS-FEM) for 3D [28], and other hybrid schemes such as αFEM [29, 30], βFEM [31], and smoothed FE-meshfree [32]. The NS-FEM can be considered a variant of the FEM. It has properties that are complementary to the FEM and can be applied easily to triangular or tetrahedral elements without modification of the formulation or procedures [33]. NS-FEM has garnered a great deal of research interest in recent years for properties such as insensitivity to element distortion and robustness against volumetric locking [33]. The computational time and computational efficiency of NS-FEM have also been investigated in previous studies for linear elastostatics [34, 35]. Node-based solutions, however, tend to suffer spatial and temporal stability problems [35]. Spatial stability issues can be resolved through smoothing operations such as stabilized conforming nodal integration (SCNI) [21], natural stable nodal integration (NSNI) [36], and quadratically consistent one point (QC1) [37], while temporal stability remains problematic. The NS-FEM has been proven spatially stable but temporally instable. Even when an unconditionally stable time-integration scheme is applied to solve a transient dynamic problem, unphysical numerical responses still appear [26]. The hallmark of temporal instability is spurious nonzero-energy modes which often appears in free vibration analyses. Such modes are spatially stable and will not accompany zero-energy modes; however, when they are excited at higher energy levels, they behave unphysically [35].

Beissel and Belytschko [38] first proposed a stabilized nodal integration procedure to eliminate spurious nonzero-energy modes by adding the square of the residual of the equilibrium equation to the potential energy function. This solution has been further extended to 2D and 3D problems to form a stabilization procedure for NS-FEM [35, 39] with a recommended range for the stabilization parameter. Taylor series expansions of the displacement fields [40, 41] can also be used to reduce the instability in direct node integration (DNI), but high-order derivatives appear in underling formulations resulting in an increase in its computational cost. Other forms of stabilization consisting of the Taylor expansion and displacement smoothing have been proposed [42, 43], wherein the nodal integration technique is directly applied to obtain stable solutions. Liu et al. [29] proposed an -FEM combining NS-FEM and standard FEM which can be used to stabilize the NS-FEM by introducing stiffening effects from the standard FEM stiffness matrix with a small . The calculation complexity and uncertain values hinder the practical application of -FEM. Nguyen-Xuan et al. [44] combined the discrete shear gap (DSG) method with a stabilization technique into the NS-FEM to eliminate transverse shear locking and maintain the stability of the formulation but could not resolve the uncertain value problem. Puso et al. [45] developed an effective nodal integration technique by adding integration points, and their method is proved effective on both small and large deformation problems. Feng et al. [46] proposed a stable nodal integration method with strain gradient for dynamics analysis of solid structures based on NS-FEM. It can achieve appropriate system stiffness in strain energy between FEM and NS-FEM solutions and indeed provide temporally stable results. There exists a variety of additions of terms on gradient term constructions, which have been adopted by researchers to stabilize the performance and proved to be highly effective [4751].

This paper proposes a novel gradient stable node-based smoothed finite element method (GS-FEM) for static and dynamic analyses for solid mechanics problems. The strain is expanded at the first order by Taylor expansion in a node-supported domain, and then gradient smoothing is implemented for the strain within each smoothing domain. The resulting stiffness matrix includes stable terms derived by the strain gradient. As a result, a temporally stable result is obtained—there are no spurious nonzero-energy modes in the free vibration analysis. The accuracy of the GS-FEM for numerical solutions is investigated by comparing against FEM-T3, FEM-Q4, and NS-FEM with the same mesh. Results show that the GS-FEM provides accurate solutions using triangular elements. It has comparative accuracy with the standard FEM using quadrilateral elements for static, free, and forced vibration analyses of solid mechanics problems.

The rest of this paper is organized as follows: the gradient smoothing method is first reviewed in Section 2. Section 3 discusses the node-based smoothed FEM and the corresponding Galerkin weak form of elastostatic problems. The gradient stabilization of NS-FEM is presented in Section 4, and free and forced vibration analyses are described in Section 5. The standard patch test we conducted are reported in Section 6. Section 7 demonstrates the effectiveness of the proposed stabilizations by comparing against FEM-T3, FEM-Q4, and NS-FEM. Section 8 provides a brief summary and conclusion.

2. Gradient Smoothing Method

The smoothing operation on the gradient of displacement field for node in a nonoverlapping smoothing domain (Figure 1) can be expressed as follows:where is the gradient of the field and is the smoothed gradient. is a smoothing function satisfying the following partition of unity:

The Heaviside-type piecewise constant function can be employed as follows:where is the area of smoothing domain .

Introducing the divergence theorem to equation (1) yields the following:where is the component of the surface outward normal of boundary .

Substituting the smoothed gradient in equation (4) into the following smoothing operation of strain tensor yields a smoothed strain:

The gradient smoothing operation on the second-order derivatives of the displacement field is performed as follows:

The derivative of the strain can be expressed as follows:

Substituting equation (6) into equation (7), the corresponding smoothed derivative of strain tensor can then be given as follows:

3. Node-Based Smoothed FEM

The strong form of the governing equation of linear elastostatics problems can be expressed as follows:where denotes the body force vector, is the Cauchy stress tensor, represents the prescribed displacement on the essential boundary , is the prescribed traction on the natural boundary , and indicates the outward surface normal of .

Under the smoothed or weakened form with a displacement field satisfying the essential boundary conditions [34], we can obtain the following equation:

Assuming that the problem domain is discretized by elements with nodes in total and smoothing domains , are constructed. In the smoothing domain , the displacement field and its gradient can be interpolated using the following equations:where is the shape function of the node that is continuous, represents the displacement component of node , and is a group of nodes associated with the smoothing domain or so-called supporting nodes of .

Substituting equation (11) into equation (4) leads towhere the so-called smoothed first-order derivatives of shape function can be denoted as

The matrix form of the smoothed strain tensor iswhere

The Cauchy stress tensor can then be immediately calculated as follows:where is the elastic matrix of the isotropic linear elastic material and .

Substituting equations (15) and (17) into equation (10) yields the following:

Eliminating the arbitrary yields the following discrete equilibrium equation:where

4. Gradient Stabilization of NS-FEM

4.1. Gradient Stabilized Nodal Integration

The direct nodal integration leads to unstable solutions due to the rank deficiency of approximated matrices. Taylor expansion can be applied to resolve this problem. The shape function and corresponding derivatives and can be, respectively, approximately expanded as follows:where is the center node of the cell.

The displacements field in can also be evaluated via Taylor expansion:

Considering a gradient expansion of the form in equation (24) for the strain in with giveswhere is the differential operator matrix given by

Assume

Equation (24) can be rewritten in the following compact form:

Equations (26) and (27) show that the second-order derivatives of the shape functions are unknown. Note that the assumed displacement field used in this work does not have second-order derivatives over the whole problem domain as the FEM shape function is applied, i.e., ; hence, terms (b) and (c) of equation (24) do not contribute to the stabilization if is calculated directly. Here, we replace equation (27) with the gradient smoothing technique presented in Section 2, that is,

The smoothed divergence of strain tensors in equation (26) can be truncated in the following vector form:in which

Eventually, the stiffness matrix in equation (20) can be calculated as follows:wherein which , , , , , and in equation (32) are the area, first, and second moments of inertia of each nodal domain, respectively.under the following assumption:

The effect of the assumption in equation (34) is expected to be negligible since nodes are located at or near the centroid of the domain generally. Moreover, for nodes located on the edges of the domain, equation (34) will be fulfilled [36] in general.

The stiffness matrix in equation (32) can be calculated as follows:

Equation (35) contains stable terms introduced by the smoothed gradient expansion of the strain at the node. The additional terms in the stabilization are necessarily positive for nonzero strain and thus provide additional correctness insurance. Moreover, in contrast to other stabilized methods, the constants associated with the additional terms do not include the tuning parameter. Stabilization with equation (35) is referred to here as “gradient-stabilized nodal integration.”

4.2. Nongradient Term Boundary Integration

Equation (33) shows that , , and cannot be directly transformed into a boundary integral as there is no gradient integration. For boundary integral treatment of , , and , a technique taken from the literature [52] is used.

Define a vector for two-dimensional problems which satisfy the following condition:where can be expressed as follows:in which and are the unit vectors of axis of coordination and and and satisfy the integration condition. Substituting equation (37) into equation (36) and then applying the divergence theorem can lead towhere denotes the outward normal vector of the smoothed subdomain and can be obtained by the addition of and , in which and are the cosine of the angles between , , and , respectively.

To satisfy the condition in equation (36), we assume while is the indefinite integral of function expressed as follows:

Similarly, assuming , is the indefinite integral of function expressed as follows:

The formulation and mathematical proof of equations (39) and (40) can be found in [53]. The effects of terms and can be neglected.

If function f is set to be 1, , and , respectively, , , and can be expressed aswherein which

Finally, , , and can be calculated as follows:

5. Free and Forced Vibration Analyses

5.1. Governing Equations

Here, we consider a deformable body subjected to body force , traction on boundary , and displacement boundary conditions on , where and . The discretized equation of dynamic analysis can be denoted as follows:where and represent the displacement vector and load vector, respectively. , , and are the global mass matrix, damp matrix, and stiffness matrix, respectively, which can be defined aswhere can be obtained from equation (35).

5.2. Mass Matrix

Here, the lump mass matrix for the linear triangular element is used in the process of dynamic analysis:where is the mass of the ith cell corresponding to local node L and and denote the mass density and the thickness of the cell, respectively.

5.3. Free Vibration Analysis

Damping and external forces are not considered during free vibration analysis. Accordingly, equation (45) can be reduced to a system of homogeneous equations [1]:

A general solution to this system can be expressed aswhere denotes the imaginary unit, is the time, indicates the eigenvector, and is the natural frequency. Substituting equation (49) into equation (48) leads to the following eigenvalue equation for natural frequency :

The natural frequencies and corresponding mode shapes of a given structure are often referred to as the “dynamic characteristics” of the structure.

5.4. Forced Vibration Analysis

Equation (45) denotes a function of both space and time for a forced vibration analysis. For simplicity, a constant damping matrix which is assumed to be a linear combination of and can be utilized:where and are the Rayleigh damping coefficients.

There are many numerical integration methods available (e.g., Newmark method, Wilson method, and generalized- method [54]). Here, the Newmark method is used. The control equation of equation (45) at time can be expressed aswherein which and are the parameters which can control the precision and stability of the algorithm. Generally, and are used.

6. Numerical Implementation

6.1. Standard Patch Test

A patch test serves as a means of assessing the convergence of a numerical method based on the Galerkin weak form. We performed linear patch tests of the proposed solution using 5 × 5 regular nodes and 25 irregular nodes, as shown in Figure 2. For irregular interior nodes, the corresponding coordinates arewhere and are the initial regular meshes sizes along x and y axes, respectively. denotes a random number between −1.0 and 1.0, and is a prescribed irregularity factor between 0.0 and 0.5. A larger value means a more irregular shape of generated meshes in the patch.

Displacement is prescribed at the boundaries as follows:

The exact displacement solution at any point in the problem domain can be also calculated by equation (55). The following displacement error norm is defined as follows:where superscript “exact” denotes the exact value of the problem, “num” denotes result solved by numerical solutions, and is the number of total field nodes. The numerical errors in displacement calculated by FEM, NS-FEM, and GS-FEM solutions with regular and increasingly irregularly distributed nodes are listed in Table 1. It is shown that the GS-FEM accurately produces the liner displacement field.

7. Numerical Examples

This section provides examples of the GS-FEM in comparison with FEM-T3, FEM-Q4, NS-FEM, and analytical solutions for a series of numerical problems. All the numerical processes are conducted in MATLAB. The convergence rate of the proposed method is measured by two standards: the displacement error norm and the strain energy norm. The displacement error norm is defined in equation (56), and the strain energy norm is as follows:

The strain energy of the numerical solution and the total strain energy of the exact solution is defined as follows:where is the strain of the exact solution. In the actual computation on equation (59), we used a very fine mesh () to calculate the “exact” strain energy .

7.1. Cantilever Beam Subjected to Tip-Shear Force

A cantilever beam with length l and height d is conducted as a benchmark problem. In this test, the system is subjected to a parabolic traction at the free end, as shown in Figure 3. The beam is assumed to have a unit thickness satisfying the plane stress condition. The analytical displacement solution is described in detail by Timoshenko and Goodier [55]:where the moment of inertia of the beam is given by .

The corresponding stresses are

The parameters are , , , , and .

During the computational process, the nodes on the left boundary are constrained using the exact displacements obtained from equations (60) and (61). The distributed parabolic shear stresses, expressed as equations (62)–(64), act on the right boundary. Figure 4 gives an example of the discretization with 288 triangular and 144 quadrilateral elements.

The distribution vertical displacement along x (y = 0) and the corresponding errors are plotted in Figures 5 and 6, where the results are also compared with the FEM-T3, FEM-Q4, and NS-FEM. The GS-FEM is more accurate than FEM-T3 or NS-FEM, which is overstiff and oversoft compared to the exact solution. It yields results as accurate as those of the FEM-Q4.

The results obtained by FEM-Q8 with 10000 elements are used as the reference solution. According to equations (56)–(59), the displacement error norm and strain energy norm can be obtained for comparison. The convergence of the strain energy is shown in Figure 7. The convergence rates of displacement and energy error norms are shown in Figures 8 and 9. As expected, the FEM modes are overly stiff and hence produce lower bounds. The NS-FEM is overly soft and provides an upper bound. The GS-FEM shows a very close-to-exact stiffness and hence very accurate results. The strain energy, displacement, and energy error norms of the GS-FEM are all better than those of the FEM-T3 or NS-FEM, and even the FEM-Q4 apart from the large h. In effect, GS-FEM provides more accurate results in the case of small DOFs (large elements) in the system. Results reported in Section 7.2 further confirm this point.

Figure 10 shows the computation time of different methods using the same direct solver with the same sets of nodes. GS-FEM takes longer to converge than that of FEM-T3, FEM-Q4, or NS-FEM. However, Figures 11 and 12 show that the efficiency of computation (i.e., computation time necessary for the same accuracy) in terms of both energy and displacement error norms is superior for GS-FEM compared to FEM-T3 or NS-FEM, and even the FEM-Q4 apart from the large h.

7.2. Infinite Plate with a Circular Hole

Figure 13 shows a plate with a central circular hole of 1 m radius subjected to a unidirectional tensile load of 1.0 N/m2 to infinity in the x-direction. Due to its symmetry, only the upper right quadrant of the plate is modeled here. The plane strain condition is considered with E of 1.0 × 103 N/m2 and of 0.3. Symmetric conditions are imposed on the left and bottom edges, and the inner boundary of the hole is traction free as reference [26]. The exact solution of the stress can be given as follows [55]:where represent the polar coordinates and is measured counterclockwise from the positive x-axis. Traction boundary conditions are imposed on the right (x = 5.0) and top (y = 5.0) edges based on the exact solution (equations (65)(67)). The displacement components corresponding to the stresses listed above arewhere and are defined in terms of Poisson’s ratio by for plane strain cases. The domain is discretized using 288 triangular and 144 quadrilateral elements, as shown in Figure 14.

Figures 15 and 16 show the distribution displacement of the GS-FEM compared with the FEM-T3, FEM-Q4, and NS-FEM. The FEM-T3 is overstiff, while the NS-FEM is oversoft compared to the exact solution. The GS-FEM provides the most accurate results—even better than those of the FEM-Q4. Figures 17 and 18 show that all the computed stresses using the GS-FEM are in good agreement with the analytical solutions. The present stresses are also very smooth though they have not been postprocessed.

The results obtained by FEM-Q8 with 6400 elements are used as the reference solution. According to equation (56)–(59), the displacement error norm and strain energy norm can be obtained for comparison. The convergence of the strain energy is shown in Figure 19. The GS-FEM provides better results than FEM-T3, FEM-Q4, or NS-FEM. The convergence rates of the error norms in displacement and energy are plotted in Figures 20 and 21. As expected, the displacement and energy norms of the GS-FEM are superior to those of FEM-T3 and NS-FEM. They are even similar to those of the FEM-Q4.

7.3. Cantilever Beam Free Vibration Analysis

As shown in Figure 3, a cantilever beam fixed on the left end with length  = 100 mm and height D = 10 mm is used here as a benchmark problem. The plane stress condition is again imposed in this test. The geometrical and material parameters include thickness  = 1.0 mm, Young’s modulus , Poisson’s ratio , and mass density . Based on the Euler–Bernoulli beam theory, the fundamental frequency is set as a reference.

Table 2 lists the first 12 natural frequencies of the beam in bold for spurious nonzero-energy modes. The first 12 modes using the NS-FEM and GS-FEM are given in Figures 22 and 23. The natural frequencies corresponding to spurious nonzero-energy modes for the NS-FEM method are marked in bold. We find that spurious nonzero-energy modes which appear in the NS-FEM are resolved by the GS-FEM. And the natural frequencies obtained using the GS-FEM are much accurate than those of the FEM-T3 and NS-FEM and close to the reference solution. The GS-FEM provides similar accuracy and converge rate to those of the FEM-Q4.

7.4. Shear Wall Free Vibration Analysis

We next investigate a shear wall with four openings, as shown in Figure 24. The parameters in the computation included Young’s modulus , Poisson’s ratio , thickness , and mass density . The bottom edge is fully clamped, and a plane stress case is assumed. Two types of meshes using triangular and quadrilateral elements are used, as shown in Figure 25. Numerical results using the FEM-Q8 with 6104 nodes and 1922 elements for the same problem are computed and used as a reference solution. This problem is also analyzed previously using a boundary element method by Brebbia et al. [56].

Table 3 lists the first 12 natural frequencies with bold form for spurious nonzero-energy modes. The first 12 modes using the NS-FEM and GS-FEM are shown in Figures 26 and 27. Again, it shows that the GS-FEM does not have any of the spurious nonzero-energy modes which appear in the NS-FEM; the natural frequencies obtained using the GS-FEM are much closer to the reference solution obtained by FEM-Q8 with 6104 nodes and 1922 elements.

7.5. Cantilever Beam Forced Vibration Analysis

The benchmark cantilever beam and mesh shown in Figures 28 and 29 via the Newmark method are investigated here for forced vibration analysis. The plane strain problem is considered with numerical parameters , , , , , , , and . The domain is represented by 10 × 2 elements, and the first loading is in the harmonic form with . The corresponding dynamic responses are shown in Figures 30 and 31 with and 0.05, respectively. The time step is used for time integration. Both damping and no damping effects are investigated. The reference solution is obtained by using ABAQUS with 2,600 8-node biquadratic quadrilateral elements.

As per the dynamic responses in Figures 3032, again, the amplitude of the GS-FEM is much more accurate than that of FEM-T3 and comparable to that of the FEM-Q4. The amplitude and corresponding period obtained by NS-FEM are both largest among the methods we test due to the model’s overly soft properties (i.e., its so-called temporal instabilization). A constant step load at the apex from is also added. Without damping, Figure 32 shows where the deflection at the apex tends toward the constant value over time. With damping, the response is damped out more quickly.

7.6. Nonlinear Structure Analysis

This example is given for demonstrating the effectiveness of the presented method for complex nonlinear problems with clearance supported caused by system uncertainties (e.g., manufacturing tolerance and wear). The nonlinearity comes from impact caused by the clearance and the deformation of the structure. The nonlinear clamped beam structure is shown in Figure 33 which consists of multisegment variable cross-sectional structures. The plane stress problem is considered with numerical parameters ,  GPa, , and  kg/m3. A harmonic excitation is acted at end of the beam with and . The clearance is considered, and the Hertz contact force model with and is applied. The time step is used for time integration. The discretized domain is shown as Figure 34. The reference solution is obtained by using ABAQUS with 2,800 8-node biquadratic quadrilateral elements.

As per the dynamic responses displayed in Figures 3537, the amplitude of the GS-FEM is much more accurate than that of FEM-T3 and NS-FEM. It is comparable to that of the FEM-Q4 and close to the reference solution demonstrating the effectiveness of the presented method for complexity nonlinear problems.

8. Conclusions

In this paper, a gradient stable node-based smoothed finite element method (GS-FEM) is proposed for stable and accurate solutions to static, free, and forced vibration analyses of 2D linear and nonlinear solid mechanics problems. Some conclusions can be summarized as follows:(1)The GS-FEM using triangular elements is temporally stable and accurate without requiring any tuning parameters for stabilization. It consistently passes the standard patch test. The formulation is straightforward, and the implementation is as easy as the FEM without any increases in DOFs. The numerical results of the GS-FEM using triangular elements are more accurate than those of FEM-Q4 apart from a small number of elements.(2)Although the computation time of the GS-FEM is longer than that of the FEM-T4, FEM-Q4, or NS-FEM, it is more computationally efficient in terms of both energy and displacement error norms.(3)The GS-FEM with triangular elements gives highly accurate results on free vibration analysis problems. No spurious nonzero-energy modes appear in such vibration analyses, and hence, the GS-FEM is temporally stable.(4)The vibration amplitude and period obtained by the GS-FEM using triangular elements is as accurate as that of FEM-Q4 on forced vibration analysis for both linear and complex nonlinear problems. And its vibration amplitude is closer to that of the higher-order FEM-Q8 than FEM-T3 or NS-FEM.

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 they have no conflicts of interest.

Acknowledgments

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