# A Gradient Stable Node-Based Smoothed Finite Element Method for Solid Mechanics Problems

Authors: Guangsong Chen; Linfang Qian; Jia Ma
Journal: Shock and Vibration
 (2019)

---

## 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.


---

## Body

## 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 [11–13], 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 [47–51].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 Section2. 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 fieldu for node xL in a nonoverlapping smoothing domain ΩLs (Figure 1) can be expressed as follows:(1)u¯i,jxL=∫ΩLsui,jxWx,x−xLdΩ,where ui,j=∂ui/∂xj is the gradient of the field u and u¯i,j is the smoothed gradient. Wx,x−xL is a smoothing function satisfying the following partition of unity:(2)∫ΩLsWx,x−xLdΩ=1.Figure 1
Triangular elements and smoothing cells associated with nodes.The Heaviside-type piecewise constant function can be employed as follows:(3)Wx,x−xL=1AL,x∈ΩLs,0,x∉ΩLs,where AL is the area of smoothing domain ΩLs.Introducing the divergence theorem to equation (1) yields the following:(4)u¯i,jxL=1AL∫ΓLsuixnjdΓ,where nj is the component of the surface outward normal n of boundary ΓLs.Substituting the smoothed gradientu¯i,j in equation (4) into the following smoothing operation of strain tensor εij yields a smoothed strain:(5)ε¯ijxL=∫ΩLsεijWx,x−xLdΩ=12u¯i,j+u¯j,i.The gradient smoothing operation on the second-order derivatives∂2uk/∂xi∂xj of the displacement field u is performed as follows:(6)u¯k,ij=∫ΩLsuk,ijWx,x−xLdΩ=1AL∫ΓLsuk,injdΓ.The derivative of the strainεkl can be expressed as follows:(7)εkl,m=12∂∂xm∂uk∂xl+∂∂xm∂ul∂xk.Substituting equation (6) into equation (7), the corresponding smoothed derivative of strain tensor εkl can then be given as follows:(8)ε¯kl,m=12u¯k,lm+u¯l,km.


## 3. Node-Based Smoothed FEM

The strong form of the governing equation of linear elastostatics problems can be expressed as follows:(9)divσ+b=0inΩ,u=uΓonΓu,σ⋅n=τonΓτ,where b denotes the body force vector, σ is the Cauchy stress tensor, uΓ represents the prescribed displacement on the essential boundary Γu, τ is the prescribed traction on the natural boundary Γτ, and n 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:(10)∫Ωδε¯uσ¯u−uTbdΩ−∫ΓtδuTtdΓ=0.Assuming that the problem domainΩ is discretized by Ne elements with Nd nodes in total and Ns smoothing domains ΩLs, L=1,2,…,Nd are constructed. In the smoothing domain ΩLs, the displacement field uix and its gradient ui,jx can be interpolated using the following equations:(11)uix=∑I∈GLΦIxuIi,(12)ui,jx=∑I∈GL∂ΦIx∂xjuIi=∑I∈GLΦI,jxuIi,where ΦIx is the shape function of the node I that is continuous, uIi represents the displacement component of node I, and GL is a group of nodes associated with the smoothing domain ΩLs or so-called supporting nodes of xL.Substituting equation (11) into equation (4) leads to(13)u¯i,jx=1AL∫ΓLs∑I∈GLΦIxuIinjdΓ=∑I∈GLΦ¯I,jxLuIi,where the so-called smoothed first-order derivatives of shape function Φ¯I,jxL can be denoted as(14)Φ¯I,jxL=1AL∫ΓLsΦIxnjdΓ.The matrix form of the smoothed strain tensorε¯ij is(15)ε¯xL=∑I∈GLB¯IxLuI,where(16)uI=uI1uI2,ε¯=ε¯11ε¯22ε¯12,B¯IxL=Φ¯I,100Φ¯I,2Φ¯I,2Φ¯I,1.The Cauchy stress tensor can then be immediately calculated as follows:(17)σ¯xL=Dε¯xL=∑I∈GLDB¯IxLuI,where D is the elastic matrix of the isotropic linear elastic material and σ¯=σ¯11σ¯22σ¯12T.Substituting equations (15) and (17) into equation (10) yields the following:(18)δuT∫ΩB¯TDB¯udΩ−δuT∫ΩΦTbdΩ−δuT∫ΓtΦTtdΓ=0.Eliminating the arbitraryδu yields the following discrete equilibrium equation:(19)K¯u=fext,where(20)K¯IJ=∫ΩB¯ITDB¯JdΩ=∑LB¯xLITDB¯xLJAL,(21)fIext=∫ΩΦITbIdΩ+∫ΓtΦITtIΓdΓ=∑LΦITbIAL+∫ΓtΦITtIdΓ.


## 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ΦIx and corresponding derivatives ΦI,xx and ΦI,yx can be, respectively, approximately expanded as follows:(22)ΦIx≈ΦIxL+x−xLΦI,xxL+y−yLΦI,yxL,ΦI,xx≈ΦI,xxL+x−xLΦI,xxxL+y−yLΦI,xyxL,ΦI,yx≈ΦI,yxL+x−xLΦI,xyxL+y−yLΦI,yyxL,where xL=xLyL is the center node of the cell.The displacements field inΩL can also be evaluated via Taylor expansion:(23)uihx≈uihxL+ui,jhxLxj−xLj=∑I∈GLΦIxLuIi+∑I∈GLΦI,jxLxj−xLjuIi.Considering a gradient expansion of the form in equation (24) for the strain εx in ΩL with xL gives(24)εx=Lux≈LuxL+L∇uxLx−xL=LuxL︸a+L∇xuxLx−xL︸b+L∇yuxLy−yL︸c,where L is the differential operator matrix given by(25)L=∂∂x00∂∂y∂∂y∂∂x,∇xuxL=∑I∈GLΦI,1xLuI,∇yuxL=∑I∈GLΦI,2xLuI.Assume(26)εL=LuxL,εLx=L∇xuxL,εLy=L∇yuxL.Equation (24) can be rewritten in the following compact form:(27)εx≈εL+εLxx−xL+εLyy−yL.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., ΦI,ijx≡0; hence, terms (b) and (c) of equation (24) do not contribute to the stabilization if ΦI,ijx is calculated directly. Here, we replace equation (27) with the gradient smoothing technique presented in Section 2, that is,(28)ε¯x≈ε¯L+ε¯Lxx−xL+ε¯Lyy−yL.The smoothed divergence of strain tensors in equation (26) can be truncated in the following vector form:(29)ε¯L=∑I∈GLB¯IuI,ε¯Lx=∑I∈GLB¯IxuI,ε¯Ly=∑I∈GLB¯IyuI,in which(30)B¯I=Φ¯I,100Φ¯I,2Φ¯I,2Φ¯I,1,B¯Ix=Φ¯I,1100Φ¯I,12Φ¯I,12Φ¯I,11,B¯Iy=Φ¯I,1200Φ¯I,22Φ¯I,22Φ¯I,12.Eventually, the stiffness matrix in equation (20) can be calculated as follows:(31)K¯IJ=∑LK¯IJxL,where(32)K¯IJxL=B¯ITDB¯JAL+B¯IxTDB¯J+B¯ITDB¯JxSLx+B¯IyTDB¯J+B¯ITDB¯JySLy+B¯IxTDB¯Jy+B¯IyTDB¯JxJLxy+B¯IxTDB¯JxJLx+B¯IyTDB¯JyJLy,in which AL, SIx, SIy, JLx, JLy, and JLxy in equation (32) are the area, first, and second moments of inertia of each nodal domain, respectively.(33)AL=∫ΩLdΩ,SIx=∫ΩLx−xLdΩ,SIy=∫ΩLy−yLdΩ,JLx=∫ΩLx−xL2dΩ,JLy=∫ΩLy−yL2dΩ,JLxy=∫ΩLx−xLy−yLdΩ,under the following assumption:(34)∫ΩIx−xLdΩ=∫ΩIy−yLdΩ=∫ΩIx−xLy−yLdΩ=0.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 ΩL 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:(35)K¯IJxL=B¯ITDB¯JAL+B¯IxTDB¯JxJLx+B¯IyTDB¯JyJLy.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 AL, JLx, and JLy cannot be directly transformed into a boundary integral as there is no gradient integration. For boundary integral treatment of AL, JLx, and JLy, a technique taken from the literature [52] is used.Define a vectorF^x1,x2 for two-dimensional problems which satisfy the following condition:(36)f=divF^x1,x2,where F^ can be expressed as follows:(37)F^=F^1i1+F^2i2,in which i1 and i2 are the unit vectors of axis of coordination x1 and x2 and F^1 and F^2 satisfy the integration condition. Substituting equation (37) into equation (36) and then applying the divergence theorem can lead to(38)∫ΩLfxdΩ=∫ΓLF^⋅ndΓ,where n denotes the outward normal vector of the smoothed subdomain ΓL and can be obtained by the addition of n1i1 and n2i2, in which n1 and n2 are the cosine of the angles between x1, x2, and n, respectively.To satisfy the condition in equation (36), we assume F^2=0 while F^ is the indefinite integral of function f expressed as follows:(39)F^1x1,x2=∫fx1,x2dx1+cx2,∫ΓLcx2n1dΓ=0.Similarly, assumingF^1=0, F^ is the indefinite integral of function f expressed as follows:(40)F^2x1,x2=∫fx1,x2dx2+cx1,∫ΓLcx1n2dΓ=0.The formulation and mathematical proof of equations (39) and (40) can be found in [53]. The effects of terms cx1 and cx2 can be neglected.If functionf is set to be 1, x−xL2, and y−yL2, respectively, WL, MLx, and MLy can be expressed as(41)AL=∫ΓLF^A⋅ndΓ,JLx=∫ΓLF^Jx⋅ndΓ,JLy=∫ΓLF^Jy⋅ndΓ,where(42)F^A=F^A,0T,F^Jx=F^Jx,0T,F^Jy=0,F^JyT,in which(43)F^A=∫1dx=x,F^Jx=∫x−xL2dx=13x−xL3,F^Jy=∫y−yL2dy=13y−yL3.Finally,AL, JLx, and JLy can be calculated as follows:(44)AL=∫ΓLF^An1dΓ,JLx=∫ΓLF^Jxn1dΓ,JLy=∫ΓLF^Jyn2dΓ.


## 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ΦIx and corresponding derivatives ΦI,xx and ΦI,yx can be, respectively, approximately expanded as follows:(22)ΦIx≈ΦIxL+x−xLΦI,xxL+y−yLΦI,yxL,ΦI,xx≈ΦI,xxL+x−xLΦI,xxxL+y−yLΦI,xyxL,ΦI,yx≈ΦI,yxL+x−xLΦI,xyxL+y−yLΦI,yyxL,where xL=xLyL is the center node of the cell.The displacements field inΩL can also be evaluated via Taylor expansion:(23)uihx≈uihxL+ui,jhxLxj−xLj=∑I∈GLΦIxLuIi+∑I∈GLΦI,jxLxj−xLjuIi.Considering a gradient expansion of the form in equation (24) for the strain εx in ΩL with xL gives(24)εx=Lux≈LuxL+L∇uxLx−xL=LuxL︸a+L∇xuxLx−xL︸b+L∇yuxLy−yL︸c,where L is the differential operator matrix given by(25)L=∂∂x00∂∂y∂∂y∂∂x,∇xuxL=∑I∈GLΦI,1xLuI,∇yuxL=∑I∈GLΦI,2xLuI.Assume(26)εL=LuxL,εLx=L∇xuxL,εLy=L∇yuxL.Equation (24) can be rewritten in the following compact form:(27)εx≈εL+εLxx−xL+εLyy−yL.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., ΦI,ijx≡0; hence, terms (b) and (c) of equation (24) do not contribute to the stabilization if ΦI,ijx is calculated directly. Here, we replace equation (27) with the gradient smoothing technique presented in Section 2, that is,(28)ε¯x≈ε¯L+ε¯Lxx−xL+ε¯Lyy−yL.The smoothed divergence of strain tensors in equation (26) can be truncated in the following vector form:(29)ε¯L=∑I∈GLB¯IuI,ε¯Lx=∑I∈GLB¯IxuI,ε¯Ly=∑I∈GLB¯IyuI,in which(30)B¯I=Φ¯I,100Φ¯I,2Φ¯I,2Φ¯I,1,B¯Ix=Φ¯I,1100Φ¯I,12Φ¯I,12Φ¯I,11,B¯Iy=Φ¯I,1200Φ¯I,22Φ¯I,22Φ¯I,12.Eventually, the stiffness matrix in equation (20) can be calculated as follows:(31)K¯IJ=∑LK¯IJxL,where(32)K¯IJxL=B¯ITDB¯JAL+B¯IxTDB¯J+B¯ITDB¯JxSLx+B¯IyTDB¯J+B¯ITDB¯JySLy+B¯IxTDB¯Jy+B¯IyTDB¯JxJLxy+B¯IxTDB¯JxJLx+B¯IyTDB¯JyJLy,in which AL, SIx, SIy, JLx, JLy, and JLxy in equation (32) are the area, first, and second moments of inertia of each nodal domain, respectively.(33)AL=∫ΩLdΩ,SIx=∫ΩLx−xLdΩ,SIy=∫ΩLy−yLdΩ,JLx=∫ΩLx−xL2dΩ,JLy=∫ΩLy−yL2dΩ,JLxy=∫ΩLx−xLy−yLdΩ,under the following assumption:(34)∫ΩIx−xLdΩ=∫ΩIy−yLdΩ=∫ΩIx−xLy−yLdΩ=0.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 ΩL 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:(35)K¯IJxL=B¯ITDB¯JAL+B¯IxTDB¯JxJLx+B¯IyTDB¯JyJLy.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 AL, JLx, and JLy cannot be directly transformed into a boundary integral as there is no gradient integration. For boundary integral treatment of AL, JLx, and JLy, a technique taken from the literature [52] is used.Define a vectorF^x1,x2 for two-dimensional problems which satisfy the following condition:(36)f=divF^x1,x2,where F^ can be expressed as follows:(37)F^=F^1i1+F^2i2,in which i1 and i2 are the unit vectors of axis of coordination x1 and x2 and F^1 and F^2 satisfy the integration condition. Substituting equation (37) into equation (36) and then applying the divergence theorem can lead to(38)∫ΩLfxdΩ=∫ΓLF^⋅ndΓ,where n denotes the outward normal vector of the smoothed subdomain ΓL and can be obtained by the addition of n1i1 and n2i2, in which n1 and n2 are the cosine of the angles between x1, x2, and n, respectively.To satisfy the condition in equation (36), we assume F^2=0 while F^ is the indefinite integral of function f expressed as follows:(39)F^1x1,x2=∫fx1,x2dx1+cx2,∫ΓLcx2n1dΓ=0.Similarly, assumingF^1=0, F^ is the indefinite integral of function f expressed as follows:(40)F^2x1,x2=∫fx1,x2dx2+cx1,∫ΓLcx1n2dΓ=0.The formulation and mathematical proof of equations (39) and (40) can be found in [53]. The effects of terms cx1 and cx2 can be neglected.If functionf is set to be 1, x−xL2, and y−yL2, respectively, WL, MLx, and MLy can be expressed as(41)AL=∫ΓLF^A⋅ndΓ,JLx=∫ΓLF^Jx⋅ndΓ,JLy=∫ΓLF^Jy⋅ndΓ,where(42)F^A=F^A,0T,F^Jx=F^Jx,0T,F^Jy=0,F^JyT,in which(43)F^A=∫1dx=x,F^Jx=∫x−xL2dx=13x−xL3,F^Jy=∫y−yL2dy=13y−yL3.Finally,AL, JLx, and JLy can be calculated as follows:(44)AL=∫ΓLF^An1dΓ,JLx=∫ΓLF^Jxn1dΓ,JLy=∫ΓLF^Jyn2dΓ.


## 5. Free and Forced Vibration Analyses

### 5.1. Governing Equations

Here, we consider a deformable bodyΩ subjected to body force b, traction t on boundary Γt, and displacement boundary conditions u=u¯ on Γu, where Γ=Γt∪Γu and Γt∩Γu=∅. The discretized equation of dynamic analysis can be denoted as follows:(45)Mu¨+Cu˙+Ku=fext,where u and fext represent the displacement vector and load vector, respectively. M, C, and K are the global mass matrix, damp matrix, and stiffness matrix, respectively, which can be defined as(46)MIJ=∑LMIJL,KIJ=∑LK¯IJL,CIJ=∑LCIJL,where K¯IJL 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:(47)Me=diagm1m2m3,where mi=ρhAL/3 is the mass of the ith cell corresponding to local node L and ρ and h 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]:(48)Mu¨+Ku=0.A general solution to this system can be expressed as(49)u=u¯expiωt,where i denotes the imaginary unit, t is the time, u¯ indicates the eigenvector, and ω is the natural frequency. Substituting equation (49) into equation (48) leads to the following eigenvalue equation for natural frequency ω:(50)−ω2M+Ku¯=0.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 C which is assumed to be a linear combination of M and K can be utilized:(51)C=αM+βK,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 t+h can be expressed as(52)a0M+a1C+Kut+h=ft+hext+Ma0ut+a2u˙t+a3u¨t+Ca1ut+a4u˙t+a5u¨t,u¨t+h=a0ut+h−ut−a2u˙t−a3u¨t,u˙t+h=u˙t+a6u¨t+a7u¨t+h,where(53)a0=1δΔt2,a1=γδΔt,a2=1δΔt,a3=12δ−1,a3=12δ−1,a4=γδ−1,a5=Δt2γδ−2,a6=Δt1−γ,a7=γΔt,in which δ and γ are the parameters which can control the precision and stability of the algorithm. Generally, δ=0.25 and γ=0.5 are used.


## 5.1. Governing Equations

Here, we consider a deformable bodyΩ subjected to body force b, traction t on boundary Γt, and displacement boundary conditions u=u¯ on Γu, where Γ=Γt∪Γu and Γt∩Γu=∅. The discretized equation of dynamic analysis can be denoted as follows:(45)Mu¨+Cu˙+Ku=fext,where u and fext represent the displacement vector and load vector, respectively. M, C, and K are the global mass matrix, damp matrix, and stiffness matrix, respectively, which can be defined as(46)MIJ=∑LMIJL,KIJ=∑LK¯IJL,CIJ=∑LCIJL,where K¯IJL 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:(47)Me=diagm1m2m3,where mi=ρhAL/3 is the mass of the ith cell corresponding to local node L and ρ and h 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]:(48)Mu¨+Ku=0.A general solution to this system can be expressed as(49)u=u¯expiωt,where i denotes the imaginary unit, t is the time, u¯ indicates the eigenvector, and ω is the natural frequency. Substituting equation (49) into equation (48) leads to the following eigenvalue equation for natural frequency ω:(50)−ω2M+Ku¯=0.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 C which is assumed to be a linear combination of M and K can be utilized:(51)C=αM+βK,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 t+h can be expressed as(52)a0M+a1C+Kut+h=ft+hext+Ma0ut+a2u˙t+a3u¨t+Ca1ut+a4u˙t+a5u¨t,u¨t+h=a0ut+h−ut−a2u˙t−a3u¨t,u˙t+h=u˙t+a6u¨t+a7u¨t+h,where(53)a0=1δΔt2,a1=γδΔt,a2=1δΔt,a3=12δ−1,a3=12δ−1,a4=γδ−1,a5=Δt2γδ−2,a6=Δt1−γ,a7=γΔt,in which δ and γ are the parameters which can control the precision and stability of the algorithm. Generally, δ=0.25 and γ=0.5 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 Figure2. For irregular interior nodes, the corresponding coordinates are(54)x′=x+Δx⋅rc⋅αr,y′=y+Δy⋅rc⋅αr,where Δx and Δy are the initial regular meshes sizes along x and y axes, respectively. rc denotes a random number between −1.0 and 1.0, and αr is a prescribed irregularity factor between 0.0 and 0.5. A larger αr value means a more irregular shape of generated meshes in the patch.Figure 2
Regular and irregular meshes for the patch test: (a) regular; (b) irregularity coefficientαr=0.2; (c) irregularity coefficient αr=0.3; (d) irregularity coefficient αr=0.4.
(a)
(b)
(c)
(d)Displacement is prescribed at the boundaries as follows:(55)ux,y=0.1+0.1x+0.2y0.05+0.15x+0.1y.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:(56)ed=∑i=1Nduiexact−uinum2∑i=1Nduiexact2,where superscript “exact” denotes the exact value of the problem, “num” denotes result solved by numerical solutions, and Nd 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.Table 1
Displacement error norm of numerical results for the standard patch test on different meshes.
Solution
FEM
NS-FEM
GS-FEM
Regular
7.5426E − 15
4.7573E − 15
1.2906E − 14
Irregular,αr=0.2
9.0044E − 15
4.7056E − 15
8.9583E − 15
Irregular,αr=0.3
2.3392E − 14
4.6646E − 15
9.1383E − 15
Irregular,αr=0.4
1.5240E − 14
5.5588E − 15
1.7212E − 14


## 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 Figure2. For irregular interior nodes, the corresponding coordinates are(54)x′=x+Δx⋅rc⋅αr,y′=y+Δy⋅rc⋅αr,where Δx and Δy are the initial regular meshes sizes along x and y axes, respectively. rc denotes a random number between −1.0 and 1.0, and αr is a prescribed irregularity factor between 0.0 and 0.5. A larger αr value means a more irregular shape of generated meshes in the patch.Figure 2
Regular and irregular meshes for the patch test: (a) regular; (b) irregularity coefficientαr=0.2; (c) irregularity coefficient αr=0.3; (d) irregularity coefficient αr=0.4.
(a)
(b)
(c)
(d)Displacement is prescribed at the boundaries as follows:(55)ux,y=0.1+0.1x+0.2y0.05+0.15x+0.1y.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:(56)ed=∑i=1Nduiexact−uinum2∑i=1Nduiexact2,where superscript “exact” denotes the exact value of the problem, “num” denotes result solved by numerical solutions, and Nd 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.Table 1
Displacement error norm of numerical results for the standard patch test on different meshes.
Solution
FEM
NS-FEM
GS-FEM
Regular
7.5426E − 15
4.7573E − 15
1.2906E − 14
Irregular,αr=0.2
9.0044E − 15
4.7056E − 15
8.9583E − 15
Irregular,αr=0.3
2.3392E − 14
4.6646E − 15
9.1383E − 15
Irregular,αr=0.4
1.5240E − 14
5.5588E − 15
1.7212E − 14


## 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:(57)ee=ψnum−ψexactψexact.The strain energy of the numerical solutionψnum and the total strain energy of the exact solution ψexact is defined as follows:(58)ψnum=12uTK¯u,(59)ψexact=12∑i=1Ne∫ΩiεiTDεidΩ,where εi is the strain of the exact solution. In the actual computation on equation (59), we used a very fine mesh (Ne⟶∞) to calculate the “exact” strain energy ψexact.


### 7.1. Cantilever Beam Subjected to Tip-Shear Force

A cantilever beam with lengthl 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]:(60)ux=Py6EI6l−3xx+2+vy2−d24,(61)uy=−P6EI3vy2l−x+144+5vd2x+x23l−x,where the moment of inertia I of the beam is given by I=d3/12.Figure 3
Cantilever beam problem.The corresponding stresses are(62)σx=Pl−xyI,(63)σy=0,(64)τxy=−P2Id24−y2.The parameters areE=3.0×107N/m2, v=0.3, d=12m, l=48m, and P=−1000N.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.Figure 4
Beam domain discretization using triangular and quadrilateral elements.
(a)
(b)The distribution vertical displacement alongx (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.Figure 5
Displacementv along y = 0.Figure 6
Displacement relative errorv−v0 along y = 0.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 7
Strain energy convergence.Figure 8
Convergence rate of displacement error norms.Figure 9
Energy error norm convergence.Figure10 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.Figure 10
Computation time.Figure 11
Energy error norm with respect to computation time.Figure 12
Displacement error norm with respect to computation time.


### 7.2. Infinite Plate with a Circular Hole

Figure13 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 v 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]:(65)σ11=1−a2r232cos2θ+cos4θ+3a42r4cos4θ,(66)σ22=−a2r212cos2θ−cos4θ−3a42r4cos4θ,(67)σ12=−a2r212sin2θ+sin4θ+3a42r4sin4θ,where r,θ 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 are(68)u1=a8GraK+1cosθ+2arK+1cosθ+cos3θ−2a3r3cos3θ,u2=a8GraK−1sinθ+2ar1−Ksinθ+sin3θ−2a3r3sin3θ,where G=E/21+v and K are defined in terms of Poisson’s ratio by K=3−4v for plane strain cases. The domain is discretized using 288 triangular and 144 quadrilateral elements, as shown in Figure 14.Figure 13
Infinite plate with circular hole and its quarter model.
(a)
(b)Figure 14
Domain discretization of the infinite plate with a circular hole using triangular and quadrilateral elements.
(a)
(b)Figures15 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.Figure 15
Displacementu along the bottom boundary.Figure 16
Displacementv along the left boundary.Figure 17
Stress along the left boundary (x = 0).Figure 18
Stress along the bottom boundary (y = 0).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.Figure 19
Convergence of strain energy.Figure 20
Convergence rate of displacement error norm.Figure 21
Convergence rate of energy error norm.


### 7.3. Cantilever Beam Free Vibration Analysis

As shown in Figure3, a cantilever beam fixed on the left end with length L = 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 t = 1.0 mm, Young’s modulus E=2.1×1010N/m2, Poisson’s ratio v=0.3, and mass density ρ=8.0×10−4kg/m3. Based on the Euler–Bernoulli beam theory, the fundamental frequency f1=0.08276×104Hz is set as a reference.Table2 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.Table 2
First 12 natural frequencies (104 Hz) of a cantilever beam.
No. of elements
No. of nodes
FEM-T3
FEM-Q4
NS-FEM
GS-FEM
Reference FEM-Q4 (100 × 10)
(10 × 1) 10 4-node elements for FEM-Q4 and 20 triangular elements for other methods
22
0.1692
0.0992
0.0576
0.1302
0.0824
0.9163
0.5791
0.3243
0.7161
0.4944
1.2869
1.2834
0.7441
1.2846
1.2824
2.1843
1.4830
0.9875
1.7230
1.3022
3.5942
2.6183
1.0112
2.8548
2.3663
3.8338
3.8140
1.1346
3.7728
3.6085
5.0335
3.8824
1.2783
3.9961
3.8442
6.2421
5.1924
1.5712
5.0417
4.9674
6.4154
6.2345
2.3697
5.9011
6.3960
7.5940
6.4846
3.2685
6.0277
6.4023
8.4790
7.7039
3.7064
6.4965
7.8853
8.7033
8.4632
3.8642
6.8162
8.9290
(20 × 2) 40 4-node elements for FEM-Q4 and 80 triangular elements for other methods
63
0.1117
0.0870
0.0675
0.0847
0.0824
0.6539
0.5199
0.4032
0.5026
0.4944
1.2843
1.2830
1.0518
1.2822
1.2824
1.6748
1.3640
1.2810
1.3055
1.3022
2.9554
2.4685
1.6467
2.3274
2.3663
3.8424
3.7477
1.8786
3.4747
3.6085
4.3866
3.8378
2.7823
3.8241
3.8442
5.8836
5.1322
3.0926
4.6703
4.9674
6.3751
6.3585
3.6783
5.8601
6.3960
7.4046
6.5731
3.8089
6.2917
6.4023
8.8210
8.0342
4.0543
7.0120
7.8853
8.9411
8.8187
4.1605
8.0983
8.9290
(40 × 4) 160 4-node elements for FEM-Q4 and 320 triangular elements for other methods
205
0.0906
0.0835
0.0778
0.0814
0.0824
0.5409
0.5004
0.4654
0.4876
0.4944
1.2831
1.2827
1.2199
1.2801
1.2824
1.4161
1.3174
1.2818
1.2831
1.3022
2.5570
2.3926
1.6689
2.3175
2.3663
3.8433
3.6462
2.2012
3.5159
3.6085
3.8786
3.8431
3.2517
3.8391
3.8442
5.3087
5.0150
3.3270
4.8097
4.9674
6.3935
6.3883
3.8344
6.1531
6.3960
6.8093
6.4561
4.5248
6.3692
6.4023
8.3473
7.9398
4.6406
7.5141
7.8853
8.9183
8.9057
5.3275
8.8515
8.9290Figure 22
First 12 modes of the cantilever beam by NS-FEM.Figure 23
First 12 modes of the cantilever beam by GS-FEM.


### 7.4. Shear Wall Free Vibration Analysis

We next investigate a shear wall with four openings, as shown in Figure24. The parameters in the computation included Young’s modulus E=10000N/m2, Poisson’s ratio v=0.2, thickness t=1.0m, and mass density ρ=1.0kg/m3. 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].Figure 24
Shear wall with four openings.Figure 25
Domain discretization using triangular and quadrilateral elements.
(a)
(b)Table3 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.Table 3
First 12 natural frequencies (×104 Hz) of shear wall.
No. of elements
No. of nodes
FEM-T3
FEM-Q4
NS-FEM
GS-FEM
Reference FEM-Q8
Reference [31]
476 4-node elements for FEM-Q4 and 952 triangular elements for other methods
559
0.3417
0.3300
0.2907
0.3197
0.3200
0.3309
1.1625
1.1293
1.0365
1.1024
1.1064
1.1429
1.2191
1.2136
1.1973
1.2088
1.2096
1.2166
1.9957
1.9000
1.6163
1.8139
1.8256
1.8833
2.5406
2.4416
2.1838
2.3626
2.3828
2.5381
2.9819
2.9196
2.3398
2.8481
2.8753
2.9673
3.2451
3.1634
2.6970
3.0793
3.1164
3.2258
3.6095
3.5348
2.7149
3.4367
3.4810
3.6232
3.7654
3.6608
2.8882
3.5615
3.6026
—
3.8378
3.7483
2.9267
3.6656
3.7072
—
4.0628
4.0067
2.9663
3.9518
3.9817
—
4.2578
4.1494
3.0992
4.0544
4.1184
—
1904 4-node elements for FEM-Q4 and 80 triangular elements for other methods
2072
0.3285
0.3234
0.3073
0.3189
0.3200
0.3309
1.1270
1.1139
1.0779
1.1030
1.1064
1.1429
1.2129
1.2110
1.2042
1.2092
1.2096
1.2166
1.8887
1.8502
1.7298
1.8152
1.8256
1.8833
2.4385
2.4020
2.3024
2.3711
2.3828
2.5381
2.9146
2.8899
2.3583
2.8646
2.8753
2.9673
3.1606
3.1317
2.7046
3.1032
3.1164
3.2258
3.5337
3.4993
2.8148
3.4639
3.4810
3.6232
3.6593
3.6222
3.0478
3.5874
3.6026
—
3.7499
3.7210
3.2445
3.6918
3.7072
—
4.0109
3.9905
3.3095
3.9711
3.9817
—
4.1703
4.1311
3.3965
4.0947
4.1184
—Figure 26
First 12 modes of shear wall by NS-FEM.Figure 27
First 12 modes of shear wall by GS-FEM.


### 7.5. Cantilever Beam Forced Vibration Analysis

The benchmark cantilever beam and mesh shown in Figures28 and 29 via the Newmark method are investigated here for forced vibration analysis. The plane strain problem is considered with numerical parameters L=4.0m, H=1.0m, t=1.0m, E=1.0N/m2, v=0.3, ρ=1.0kg/m3, α=0.005, and β=0.272. The domain is represented by 10 × 2 elements, and the first loading is in the harmonic form with ft=cosωft. The corresponding dynamic responses are shown in Figures 30 and 31 with ωf=0.01 and 0.05, respectively. The time step Δt=2.0 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.Figure 28
Cantilever beam.Figure 29
Domain discretization using triangular and quadrilateral elements.
(a)
(b)Figure 30
Transient response of the cantilever beam subjected to harmonic loading withωf=0.01: (a) without damping; (b) with damping.
(a)
(b)Figure 31
Transient response of the cantilever beam subjected to harmonic loading withωf=0.05: (a) without damping; (b) with damping.
(a)
(b)As per the dynamic responses in Figures30–32, 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 ft=1 at the apex from t=0 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.Figure 32
Transient response of the cantilever beam subjected to constant step loading: (a) without damping; (b) with damping.
(a)
(b)


### 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 Figure33 which consists of multisegment variable cross-sectional structures. The plane stress problem is considered with numerical parameters t=1.0, E=45 GPa, v=0.3, and ρ=1800 kg/m3. A harmonic excitation is acted at end of the beam F=Asinωft with A=5000 and ωf=1000. The clearance δ=5×10−5m is considered, and the Hertz contact force model Fc=Kδn with n=1.2 and K=1×106 is applied. The time step Δt=5×10−5 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.Figure 33
Nonlinear clamped beam structure.Figure 34
Domain discretization using triangular and quadrilateral elements.
(a)
(b)As per the dynamic responses displayed in Figures35–37, 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.Figure 35
Transient response of point A.Figure 36
Transient response of point B.Figure 37
Transient response of point C.


## 7.1. Cantilever Beam Subjected to Tip-Shear Force

A cantilever beam with lengthl 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]:(60)ux=Py6EI6l−3xx+2+vy2−d24,(61)uy=−P6EI3vy2l−x+144+5vd2x+x23l−x,where the moment of inertia I of the beam is given by I=d3/12.Figure 3
Cantilever beam problem.The corresponding stresses are(62)σx=Pl−xyI,(63)σy=0,(64)τxy=−P2Id24−y2.The parameters areE=3.0×107N/m2, v=0.3, d=12m, l=48m, and P=−1000N.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.Figure 4
Beam domain discretization using triangular and quadrilateral elements.
(a)
(b)The distribution vertical displacement alongx (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.Figure 5
Displacementv along y = 0.Figure 6
Displacement relative errorv−v0 along y = 0.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 7
Strain energy convergence.Figure 8
Convergence rate of displacement error norms.Figure 9
Energy error norm convergence.Figure10 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.Figure 10
Computation time.Figure 11
Energy error norm with respect to computation time.Figure 12
Displacement error norm with respect to computation time.


## 7.2. Infinite Plate with a Circular Hole

Figure13 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 v 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]:(65)σ11=1−a2r232cos2θ+cos4θ+3a42r4cos4θ,(66)σ22=−a2r212cos2θ−cos4θ−3a42r4cos4θ,(67)σ12=−a2r212sin2θ+sin4θ+3a42r4sin4θ,where r,θ 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 are(68)u1=a8GraK+1cosθ+2arK+1cosθ+cos3θ−2a3r3cos3θ,u2=a8GraK−1sinθ+2ar1−Ksinθ+sin3θ−2a3r3sin3θ,where G=E/21+v and K are defined in terms of Poisson’s ratio by K=3−4v for plane strain cases. The domain is discretized using 288 triangular and 144 quadrilateral elements, as shown in Figure 14.Figure 13
Infinite plate with circular hole and its quarter model.
(a)
(b)Figure 14
Domain discretization of the infinite plate with a circular hole using triangular and quadrilateral elements.
(a)
(b)Figures15 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.Figure 15
Displacementu along the bottom boundary.Figure 16
Displacementv along the left boundary.Figure 17
Stress along the left boundary (x = 0).Figure 18
Stress along the bottom boundary (y = 0).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.Figure 19
Convergence of strain energy.Figure 20
Convergence rate of displacement error norm.Figure 21
Convergence rate of energy error norm.


## 7.3. Cantilever Beam Free Vibration Analysis

As shown in Figure3, a cantilever beam fixed on the left end with length L = 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 t = 1.0 mm, Young’s modulus E=2.1×1010N/m2, Poisson’s ratio v=0.3, and mass density ρ=8.0×10−4kg/m3. Based on the Euler–Bernoulli beam theory, the fundamental frequency f1=0.08276×104Hz is set as a reference.Table2 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.Table 2
First 12 natural frequencies (104 Hz) of a cantilever beam.
No. of elements
No. of nodes
FEM-T3
FEM-Q4
NS-FEM
GS-FEM
Reference FEM-Q4 (100 × 10)
(10 × 1) 10 4-node elements for FEM-Q4 and 20 triangular elements for other methods
22
0.1692
0.0992
0.0576
0.1302
0.0824
0.9163
0.5791
0.3243
0.7161
0.4944
1.2869
1.2834
0.7441
1.2846
1.2824
2.1843
1.4830
0.9875
1.7230
1.3022
3.5942
2.6183
1.0112
2.8548
2.3663
3.8338
3.8140
1.1346
3.7728
3.6085
5.0335
3.8824
1.2783
3.9961
3.8442
6.2421
5.1924
1.5712
5.0417
4.9674
6.4154
6.2345
2.3697
5.9011
6.3960
7.5940
6.4846
3.2685
6.0277
6.4023
8.4790
7.7039
3.7064
6.4965
7.8853
8.7033
8.4632
3.8642
6.8162
8.9290
(20 × 2) 40 4-node elements for FEM-Q4 and 80 triangular elements for other methods
63
0.1117
0.0870
0.0675
0.0847
0.0824
0.6539
0.5199
0.4032
0.5026
0.4944
1.2843
1.2830
1.0518
1.2822
1.2824
1.6748
1.3640
1.2810
1.3055
1.3022
2.9554
2.4685
1.6467
2.3274
2.3663
3.8424
3.7477
1.8786
3.4747
3.6085
4.3866
3.8378
2.7823
3.8241
3.8442
5.8836
5.1322
3.0926
4.6703
4.9674
6.3751
6.3585
3.6783
5.8601
6.3960
7.4046
6.5731
3.8089
6.2917
6.4023
8.8210
8.0342
4.0543
7.0120
7.8853
8.9411
8.8187
4.1605
8.0983
8.9290
(40 × 4) 160 4-node elements for FEM-Q4 and 320 triangular elements for other methods
205
0.0906
0.0835
0.0778
0.0814
0.0824
0.5409
0.5004
0.4654
0.4876
0.4944
1.2831
1.2827
1.2199
1.2801
1.2824
1.4161
1.3174
1.2818
1.2831
1.3022
2.5570
2.3926
1.6689
2.3175
2.3663
3.8433
3.6462
2.2012
3.5159
3.6085
3.8786
3.8431
3.2517
3.8391
3.8442
5.3087
5.0150
3.3270
4.8097
4.9674
6.3935
6.3883
3.8344
6.1531
6.3960
6.8093
6.4561
4.5248
6.3692
6.4023
8.3473
7.9398
4.6406
7.5141
7.8853
8.9183
8.9057
5.3275
8.8515
8.9290Figure 22
First 12 modes of the cantilever beam by NS-FEM.Figure 23
First 12 modes of the cantilever beam by GS-FEM.


## 7.4. Shear Wall Free Vibration Analysis

We next investigate a shear wall with four openings, as shown in Figure24. The parameters in the computation included Young’s modulus E=10000N/m2, Poisson’s ratio v=0.2, thickness t=1.0m, and mass density ρ=1.0kg/m3. 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].Figure 24
Shear wall with four openings.Figure 25
Domain discretization using triangular and quadrilateral elements.
(a)
(b)Table3 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.Table 3
First 12 natural frequencies (×104 Hz) of shear wall.
No. of elements
No. of nodes
FEM-T3
FEM-Q4
NS-FEM
GS-FEM
Reference FEM-Q8
Reference [31]
476 4-node elements for FEM-Q4 and 952 triangular elements for other methods
559
0.3417
0.3300
0.2907
0.3197
0.3200
0.3309
1.1625
1.1293
1.0365
1.1024
1.1064
1.1429
1.2191
1.2136
1.1973
1.2088
1.2096
1.2166
1.9957
1.9000
1.6163
1.8139
1.8256
1.8833
2.5406
2.4416
2.1838
2.3626
2.3828
2.5381
2.9819
2.9196
2.3398
2.8481
2.8753
2.9673
3.2451
3.1634
2.6970
3.0793
3.1164
3.2258
3.6095
3.5348
2.7149
3.4367
3.4810
3.6232
3.7654
3.6608
2.8882
3.5615
3.6026
—
3.8378
3.7483
2.9267
3.6656
3.7072
—
4.0628
4.0067
2.9663
3.9518
3.9817
—
4.2578
4.1494
3.0992
4.0544
4.1184
—
1904 4-node elements for FEM-Q4 and 80 triangular elements for other methods
2072
0.3285
0.3234
0.3073
0.3189
0.3200
0.3309
1.1270
1.1139
1.0779
1.1030
1.1064
1.1429
1.2129
1.2110
1.2042
1.2092
1.2096
1.2166
1.8887
1.8502
1.7298
1.8152
1.8256
1.8833
2.4385
2.4020
2.3024
2.3711
2.3828
2.5381
2.9146
2.8899
2.3583
2.8646
2.8753
2.9673
3.1606
3.1317
2.7046
3.1032
3.1164
3.2258
3.5337
3.4993
2.8148
3.4639
3.4810
3.6232
3.6593
3.6222
3.0478
3.5874
3.6026
—
3.7499
3.7210
3.2445
3.6918
3.7072
—
4.0109
3.9905
3.3095
3.9711
3.9817
—
4.1703
4.1311
3.3965
4.0947
4.1184
—Figure 26
First 12 modes of shear wall by NS-FEM.Figure 27
First 12 modes of shear wall by GS-FEM.


## 7.5. Cantilever Beam Forced Vibration Analysis

The benchmark cantilever beam and mesh shown in Figures28 and 29 via the Newmark method are investigated here for forced vibration analysis. The plane strain problem is considered with numerical parameters L=4.0m, H=1.0m, t=1.0m, E=1.0N/m2, v=0.3, ρ=1.0kg/m3, α=0.005, and β=0.272. The domain is represented by 10 × 2 elements, and the first loading is in the harmonic form with ft=cosωft. The corresponding dynamic responses are shown in Figures 30 and 31 with ωf=0.01 and 0.05, respectively. The time step Δt=2.0 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.Figure 28
Cantilever beam.Figure 29
Domain discretization using triangular and quadrilateral elements.
(a)
(b)Figure 30
Transient response of the cantilever beam subjected to harmonic loading withωf=0.01: (a) without damping; (b) with damping.
(a)
(b)Figure 31
Transient response of the cantilever beam subjected to harmonic loading withωf=0.05: (a) without damping; (b) with damping.
(a)
(b)As per the dynamic responses in Figures30–32, 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 ft=1 at the apex from t=0 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.Figure 32
Transient response of the cantilever beam subjected to constant step loading: (a) without damping; (b) with damping.
(a)
(b)


## 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 Figure33 which consists of multisegment variable cross-sectional structures. The plane stress problem is considered with numerical parameters t=1.0, E=45 GPa, v=0.3, and ρ=1800 kg/m3. A harmonic excitation is acted at end of the beam F=Asinωft with A=5000 and ωf=1000. The clearance δ=5×10−5m is considered, and the Hertz contact force model Fc=Kδn with n=1.2 and K=1×106 is applied. The time step Δt=5×10−5 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.Figure 33
Nonlinear clamped beam structure.Figure 34
Domain discretization using triangular and quadrilateral elements.
(a)
(b)As per the dynamic responses displayed in Figures35–37, 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.Figure 35
Transient response of point A.Figure 36
Transient response of point B.Figure 37
Transient response of point C.


## 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.



---

*Source: 8610790-2019-06-04.xml*