# FE-Meshfree QUAD4 Element with Modified Radial Point Interpolation Function for Structural Dynamic Analysis

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

---

## Abstract

The partition-of-unity method based on FE-Meshfree QUAD4 element synthesizes the respective advantages of meshfree and finite element methods by exploiting composite shape functions to obtain high-order global approximations. This method yields high accuracy and convergence rate without necessitating extra nodes or DOFs. In this study, the FE-Meshfree method is extended to the free and forced vibration analysis of two-dimensional solids. A modified radial point interpolation function without any supporting tuning parameters is applied to construct the composite shape functions. The governing equations of elastodynamic problem are transformed into a standard weak formulation and then discretized into time-dependent equations which are solved via Bathe time integration scheme to conduct the forced vibration analysis. Several numerical test problems are solved and compared against previously published numerical solutions. Results show that the proposed FE-Meshfree QUAD4 element owns greater tolerance for mesh distortion and provides more accurate solutions.


---

## Body

## 1. Introduction

Analysis of structural dynamics problems is a crucial step in many engineering applications. However, existing analytical methods are limited, as it is very challenging to find an exact solution to a class of dynamic problems. In addition, such solutions are reachable only with simple geometrical configurations and boundary conditions. Accordingly, numerical methods have emerged as an alternative approach to gain approximate solutions, which now represent the most widely used computational engineering tools.Finite element method (FEM), serving as one of the most popular numerical approaches, is employed to analyze the structural vibration characteristics [1–3]. It is accurate, convenient, and flexible and has been successfully applied in solid mechanics, fluid mechanics, heat transfer, and other fields. However, the gradient of lower-order elements often exhibits discontinuity at element boundaries and nodes, accordingly rendering it unreliable in certain situations. Some isoparametric elements are also highly sensitive to mesh distortion, which decreases the accuracy of some FEM-based solutions [4, 5]. Meshfree methods only employ nodes to discretize the problem domain and therefore are immune to mesh distortion problems. They are suitable for solving complex practical problems such as large deformation, fracture propagation simulation, and impact-induced failure [6–15]. However, there exist important drawbacks for meshfree methods, such as the essential boundary condition implementation handicap, high computational cost, and complex trial function construction process, thereby giving rise to several hybrid schemes to enhance the applicability of meshfree methods. Liu et al. [16], for example, developed the smoothed finite element method (SFEM) by combining the existing FEM technology and the strain-smoothing technique of meshfree methodology [17]. In SFEM, the strain is smoothed based on nodes or the edges of elements. SFEM provides an upper bound solution for problems of solid mechanics making it complementary to FEM. Other hybrid methods have also been proposed over the past decades [18–27].Another common and extensively studied hybrid method is the partition-of-unity method (PUM) [28, 29]. The PUM can also be utilized to develop HP-clouds [30], generalized finite element method (GFEM) [31], particle-PUM [32], numerical manifold method [33], and extended FEM [34, 35]. PUM also serves as a convenient approach to construct the higher-order global approximations without adding external nodes, thereby achieving higher accuracy and convergence rate than other methods. However, PUM has a serious disadvantage known as the “linear dependence” (LD) problem [36], wherein the global stiffness matrix remains singular after applying the basic boundary condition to eliminate the rigid body displacement. LD occurs when both the PUM functions and the local approximation functions are taken as explicit polynomials. Tian et al. [37] conducts a comprehensive review about some effective approaches to eliminate the LD problem. Recently, a new family of elements called “FE-Meshfree” elements based on the PUM has been developed as a combination of the respective strengths of meshfree and FEM techniques such as the PU-based FE-Meshfree QUAD4 element by Rajendran et al. [38, 39], FE-Meshfree RRIA3 element by Xu [36], FE-Meshfree elements with continuous nodal stress for 2D and 3D problems by Tang and Yang [40–42], FE-Meshfree QUAD4 element with radial-polynomial basis functions by Xu [43], PU-based FE-Meshfree QUAD4 and triangular with continuous nodal stress by Yang [44, 45], and others [46, 47]. FE-Meshfree features the composite shape functions, which are constructed by combining the meshfree and finite element shape functions. The shape functions of FE-Meshfree elements possess the much desired Kronecker delta property, which is crucial in the implementation process of essential boundary conditions under FEM, interelement compatibility properties, and all the completeness properties that are necessary to ensure the reproducibility of not only the linear polynomial terms but also the higher-order terms in the local approximation. FE-Meshfree elements, as opposed to FEM elements, inherit some advantages, like better accuracy, higher convergence rate, and greater robustness to mesh distortion. FE-Meshfree elements are also free from the LD problem which exists in many of the PUM-based finite elements. PU-based FE-Meshfree method has already been applied successfully within statics, free vibration [48, 49], forced vibration [50], geometric nonlinearity [51, 52], acoustic engineering [53], and other fields.In the present study, the FE-Meshfree method with modified radial point interpolation function is to be applied in the free and forced vibration analysis of two-dimensional solids. The rest of this paper is organized as follows. Section2 briefly reviews the formulation of the composite shape function and the modified radial point interpolation function. This is followed by the equations for free and forced vibration analysis in Section 3. Section 4 is devoted to the numerical tests, in which the convergence and accuracy of the proposed method is verified. Finally, conclusions are drawn in Section 5.


## 2. Composite Shape Functions Based on Modified Radial Based Function

### 2.1. Composite Shape Function

Regarding the composite shape functions for the FE-Meshfree QUAD4 element, constituted by the combination of radial-polynomial basis functions proposed by Xu [43], the displacement interpolation in j-th element ej can be written as follows:(1)ux=∑i=14Nixuix,where Nixi=1,2,…,4 are isoparametric shape functions used in the classic QUAD4 element and uix denotes the nodal displacement function vector which can be constructed by the nodes in the node support domain Ωi which is defined by the element including node i. As shown in Figure 1, the node support domain Ω1 for node i=1 is marked by a dashed line, which includes nodes 1,2,3,4,5,6,7,8,9. Regarding the unknown nodal displacement functions, they can be interpolated by the polynomial basis function pzx and radial basis function rjx, which can be expressed as follows:(2)uix=∑j=1nrjxaj+∑z=1mpzxbz=rTxa+pTxb.Figure 1
Nodal support domain.In the above,m is the number of polynomial terms and n denotes the number of nodes in the node support domain Ωi.The vectorsa and b in equation (2) are vectors yet to be determined. px and rx denote the corresponding polynomial and radial basis functions, respectively. If we consider three respective terms for 2D problems, px can be written as follows:(3)px=1xyT.The radial basis function vectorrx is expressed as(4)rx=r1d1r2d2…rkdk…rnx,yT,where rkdk is the function of distance dk, which denotes the space distance between the selected interpolation point x and a node xk. Also for the 2D problem, dk can be expressed as(5)dk=dkx,y=x−xk2+y−yk2.Enforcing equation (2) to pass through all the nodes in the nodal support domain yields the following equation:(6)ui=Ra+Pb,where ui is a nodal displacement vector of all the nodes in the nodal support domain and a and b are, respectively, the vector of unknown coefficients for radial basis functions and polynomial basis functions:(7)aT=a1a2…an,bT=b1b2…bm.The related matricesR and P are(8)R=r1x1,y1r2x1,y1…rnx1,y1r1x2,y2r2x2,y2…rnx2,y2⋮⋮⋱⋮r1xn,ynr2xn,yn…rnxn,yn,P=p1x1,y1p2x1,y1…pmx1,y1p1x2,y2p2x2,y2…pmx2,y2⋮⋮⋱⋮p1xn,ynp2xn,yn…pmxn,yn.The following constraint condition can be imposed to guarantee a unique displacement approximation as follows [54]:(9)PTa=0.From equation (6), we can obtain:(10)a=R−1ui−R−1Pb.Substituting equation (10) into equation (9), the vector b can be solved:(11)b=Sbui,(12)Sb=PTR−1P−1PTR−1.Substituting equation (11) into equation (10), the vector a can be obtained:(13)a=Saui,(14)Sa=R−1I−PSb=R−1−R−1PSb.Eventually the nodal displacement function can be expressed by substituting equations (11) and (13) into equation (2):(15)uix,y=Φiui,i=1,2,3,4,Φi=rTx,ySa+pTx,ySb,where(16)Φi=Φ1iΦ2i…ΦNi,ui=u1u2…uNT.The arraysΦi, i=1,2,3,4 are to be assembled into a matrix Φ4×L, in which L is the number of nodes in the j-th element support domain Ωj which is defined by the element including nodes in element j. As shown in Figure 2, the element support domain Ω1=Ω1∪Ω2∪Ω3∪Ω4=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 for element j=1 is marked by a dot dash line.Figure 2
Element support domain.Substituting equations (10) and (11) into equation (1) allows us to rewrite the element displacement field as follows:(17)ux=N1×4Φ4×LuL×1=N1×4Φ4×LuL×1=ψ1×LuL×1.In the above,uL×1 denotes the x-direction displacement of all the nodes in an element support domain Ωe, and the composite shape function array, ψ1×L, can be represented as(18)ψ1×L=ψ1,ψ2,⋯,ψL=N1×4Φ4×L.Similarly, the element displacement field iny-direction is written as follows:(19)vx=ψ1×LvL×1,where vL×1 denotes the y-direction displacements of all nodes in Ωe.The Cartesian derivatives of shape function are obtained by(20)ψ,x1×L=ψ1,xψ2,x⋯ψL,x=N,x1×4Φ4×L+N1×4Φ,x4×L,ψ,y1×L=ψ1,yψ2,y⋯ψL,y=N,y1×4Φ4×L+N1×4Φ,y4×L.The composite shape functionψ is characterized by the Kronecker delta property, interelement compatibility property, and higher-order completeness properties, i.e., the reproducibility of all the Cartesian terms appearing under the assumed basis of equation (1).


### 2.2. Modified Radial Based Function

During the numerical analysis process, many specific radial basis functions are used in radial point interpolation method (RPIM) to ensure favorable performances. The typical radial basis functions, usually employed in engineering problems, are tabulated in Table1. However, the tuning parameters show significant influences on the quality of the basis function during the shape functions construction process to achieve a smooth solution. So far, the optimal values of the shape parameters are still being explored and remain questionable.Table 1
Typical radial basis functions.
Name
Expression
Parameters
Gaussian (exp)
r
j
x
,
y
=
exp
−
b
d
j
/
d
c
2
b
Multiquadrics (MQ)
r
j
x
,
y
=
d
j
2
+
b
d
c
2
q
b, q
Inverse multiquadrics (IMQ)
r
j
x
,
y
=
1
/
d
j
2
+
b
d
c
2
q
b, qIn this study, we used the modified basis function proposed by Do [55] to construct the composite shape functions without any supporting tuning parameters. The radial basis function can be taken as a compactly supported form:(21)rkx=sd¯kMd¯k,d¯k=dkds.In the above,sd¯k denotes a compactly supported function, Md¯k represents a radial basis function, and ds is the size of the support domain. Generally, sd¯k can be chosen as a Heaviside step function, if d¯k≤1, sd¯k=1, otherwise, sd¯k=0. However, as the nodal support domain and the corresponding nodes have been defined, sd¯k does not need to be used, and equation (21) can be rewritten as(22)rkx=Md¯k,d¯k=dkds.The following radial basis function in a compactly supported form based on the original kernel function of the quartic spline (QS) [6] is applied:(23)Md¯k=1−6×θd¯k2+8×θd¯k3−3×θd¯k4.The correlation parameterθ can be set to any value due to the normalized distance, which is demonstrated by numerical examples below, and we revisit this in Section 4.3.1. The circular or rectangular support domain can be applied in the meshfree method, but in this work, the node support domain is inherently suitable as a support domain and its equivalent radius is defined as follows:(24)ds=κdc,where κ is the scale factor and dc is a characteristic length related to the nodal distance, which can be estimated as the average distance between a group of nodes in the support domain. dc is determined as follows for the irregular nodal distribution [56]:(25)dc=AsnAs−1,where As denotes the area of the nodal support domain and nAs is the corresponding number of nodes in the nodal support domain.Substituting equation (24) into equation (23), equation (23) can be rewritten as(26)Mdk=1−6×θκdkdc2+8×θκdkdc3−3×θκdkdc4.From equation (26), θ/κ can be regarded as one variable; therefore, κ=1 can be adopted in the presented formulation.


## 2.1. Composite Shape Function

Regarding the composite shape functions for the FE-Meshfree QUAD4 element, constituted by the combination of radial-polynomial basis functions proposed by Xu [43], the displacement interpolation in j-th element ej can be written as follows:(1)ux=∑i=14Nixuix,where Nixi=1,2,…,4 are isoparametric shape functions used in the classic QUAD4 element and uix denotes the nodal displacement function vector which can be constructed by the nodes in the node support domain Ωi which is defined by the element including node i. As shown in Figure 1, the node support domain Ω1 for node i=1 is marked by a dashed line, which includes nodes 1,2,3,4,5,6,7,8,9. Regarding the unknown nodal displacement functions, they can be interpolated by the polynomial basis function pzx and radial basis function rjx, which can be expressed as follows:(2)uix=∑j=1nrjxaj+∑z=1mpzxbz=rTxa+pTxb.Figure 1
Nodal support domain.In the above,m is the number of polynomial terms and n denotes the number of nodes in the node support domain Ωi.The vectorsa and b in equation (2) are vectors yet to be determined. px and rx denote the corresponding polynomial and radial basis functions, respectively. If we consider three respective terms for 2D problems, px can be written as follows:(3)px=1xyT.The radial basis function vectorrx is expressed as(4)rx=r1d1r2d2…rkdk…rnx,yT,where rkdk is the function of distance dk, which denotes the space distance between the selected interpolation point x and a node xk. Also for the 2D problem, dk can be expressed as(5)dk=dkx,y=x−xk2+y−yk2.Enforcing equation (2) to pass through all the nodes in the nodal support domain yields the following equation:(6)ui=Ra+Pb,where ui is a nodal displacement vector of all the nodes in the nodal support domain and a and b are, respectively, the vector of unknown coefficients for radial basis functions and polynomial basis functions:(7)aT=a1a2…an,bT=b1b2…bm.The related matricesR and P are(8)R=r1x1,y1r2x1,y1…rnx1,y1r1x2,y2r2x2,y2…rnx2,y2⋮⋮⋱⋮r1xn,ynr2xn,yn…rnxn,yn,P=p1x1,y1p2x1,y1…pmx1,y1p1x2,y2p2x2,y2…pmx2,y2⋮⋮⋱⋮p1xn,ynp2xn,yn…pmxn,yn.The following constraint condition can be imposed to guarantee a unique displacement approximation as follows [54]:(9)PTa=0.From equation (6), we can obtain:(10)a=R−1ui−R−1Pb.Substituting equation (10) into equation (9), the vector b can be solved:(11)b=Sbui,(12)Sb=PTR−1P−1PTR−1.Substituting equation (11) into equation (10), the vector a can be obtained:(13)a=Saui,(14)Sa=R−1I−PSb=R−1−R−1PSb.Eventually the nodal displacement function can be expressed by substituting equations (11) and (13) into equation (2):(15)uix,y=Φiui,i=1,2,3,4,Φi=rTx,ySa+pTx,ySb,where(16)Φi=Φ1iΦ2i…ΦNi,ui=u1u2…uNT.The arraysΦi, i=1,2,3,4 are to be assembled into a matrix Φ4×L, in which L is the number of nodes in the j-th element support domain Ωj which is defined by the element including nodes in element j. As shown in Figure 2, the element support domain Ω1=Ω1∪Ω2∪Ω3∪Ω4=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 for element j=1 is marked by a dot dash line.Figure 2
Element support domain.Substituting equations (10) and (11) into equation (1) allows us to rewrite the element displacement field as follows:(17)ux=N1×4Φ4×LuL×1=N1×4Φ4×LuL×1=ψ1×LuL×1.In the above,uL×1 denotes the x-direction displacement of all the nodes in an element support domain Ωe, and the composite shape function array, ψ1×L, can be represented as(18)ψ1×L=ψ1,ψ2,⋯,ψL=N1×4Φ4×L.Similarly, the element displacement field iny-direction is written as follows:(19)vx=ψ1×LvL×1,where vL×1 denotes the y-direction displacements of all nodes in Ωe.The Cartesian derivatives of shape function are obtained by(20)ψ,x1×L=ψ1,xψ2,x⋯ψL,x=N,x1×4Φ4×L+N1×4Φ,x4×L,ψ,y1×L=ψ1,yψ2,y⋯ψL,y=N,y1×4Φ4×L+N1×4Φ,y4×L.The composite shape functionψ is characterized by the Kronecker delta property, interelement compatibility property, and higher-order completeness properties, i.e., the reproducibility of all the Cartesian terms appearing under the assumed basis of equation (1).


## 2.2. Modified Radial Based Function

During the numerical analysis process, many specific radial basis functions are used in radial point interpolation method (RPIM) to ensure favorable performances. The typical radial basis functions, usually employed in engineering problems, are tabulated in Table1. However, the tuning parameters show significant influences on the quality of the basis function during the shape functions construction process to achieve a smooth solution. So far, the optimal values of the shape parameters are still being explored and remain questionable.Table 1
Typical radial basis functions.
Name
Expression
Parameters
Gaussian (exp)
r
j
x
,
y
=
exp
−
b
d
j
/
d
c
2
b
Multiquadrics (MQ)
r
j
x
,
y
=
d
j
2
+
b
d
c
2
q
b, q
Inverse multiquadrics (IMQ)
r
j
x
,
y
=
1
/
d
j
2
+
b
d
c
2
q
b, qIn this study, we used the modified basis function proposed by Do [55] to construct the composite shape functions without any supporting tuning parameters. The radial basis function can be taken as a compactly supported form:(21)rkx=sd¯kMd¯k,d¯k=dkds.In the above,sd¯k denotes a compactly supported function, Md¯k represents a radial basis function, and ds is the size of the support domain. Generally, sd¯k can be chosen as a Heaviside step function, if d¯k≤1, sd¯k=1, otherwise, sd¯k=0. However, as the nodal support domain and the corresponding nodes have been defined, sd¯k does not need to be used, and equation (21) can be rewritten as(22)rkx=Md¯k,d¯k=dkds.The following radial basis function in a compactly supported form based on the original kernel function of the quartic spline (QS) [6] is applied:(23)Md¯k=1−6×θd¯k2+8×θd¯k3−3×θd¯k4.The correlation parameterθ can be set to any value due to the normalized distance, which is demonstrated by numerical examples below, and we revisit this in Section 4.3.1. The circular or rectangular support domain can be applied in the meshfree method, but in this work, the node support domain is inherently suitable as a support domain and its equivalent radius is defined as follows:(24)ds=κdc,where κ is the scale factor and dc is a characteristic length related to the nodal distance, which can be estimated as the average distance between a group of nodes in the support domain. dc is determined as follows for the irregular nodal distribution [56]:(25)dc=AsnAs−1,where As denotes the area of the nodal support domain and nAs is the corresponding number of nodes in the nodal support domain.Substituting equation (24) into equation (23), equation (23) can be rewritten as(26)Mdk=1−6×θκdkdc2+8×θκdkdc3−3×θκdkdc4.From equation (26), θ/κ can be regarded as one variable; therefore, κ=1 can be adopted in the presented formulation.


## 3. Structural Dynamic Analysis

### 3.1. Discrete Governing Equations

Here, we consider a deformable bodyΩ subjected to body force f, traction t on boundary Γt, and displacement boundary conditions u=u¯ on Γu, where Γ=Γt∪Γu, Γt∩Γu=∅. The virtual weak form of the dynamic equation can be expressed as(27)∫ΩδεTσdΩ−∫ΩδuTf−ρu¨−cu˙dΩ−∫ΓtδuTtdΓ=0,where ρ and c are the density and the coefficient of damping. ε and σ denote the strain and stress in the element. As the spatial discretization procedure in the FEM, the virtual displacements and strains within element can be expressed as(28)u=∑INIdI,δu=∑INIδdI,(29)ε=∑IBIdI,δε=∑IBIδdI,where dI=uIvIT is the nodal displacement vector and B represents the standard displacement gradient matrix. For 2D problems, it can be determined as follows:(30)BI=∇SNIx=NI,x00NI,yNI,yNI,x.The stress can be obtained by the following constitutive relation:(31)σ=Dε,where D is the matrix of the elastic constants of the material.Substituting equations (28) and (29)into equation (27) yields(32)∫ΩδdTBTσdΩ−∫ΩδdTNTf−ρu¨−cu˙dΩ−∫ΓtδdTNTtdΓ=0.According to the principle of virtual work, the discrete governing equation can be expressed as(33)∫ΩBTσdΩ−∫ΩNTf−ρu¨−cu˙dΩ−∫ΓtNTtdΓ=0,or(34)Md¨+Cd˙+Kd=fext,in which(35)fext=∫ΩNTfdΩ+∫ΓtNTtdΓ,M=∫ΩNTρNdΩ,K=∫ΩBTDBdΩ,C=∫ΩNTcNdΩ.


### 3.2. Free Vibration Analysis

The damping and external forces are not taken into account during free vibration analysis. Accordingly, equation (34) can be reduced to a system of homogeneous equations [1]:(36)Md¨+Kd=0.A general solution to this system is(37)d=d¯expiωt,where i denotes the imaginary unit, t is time, d¯ indicates the eigenvector, and ω is the natural frequency. Substituting equation (37) into equation (36) leads to the following eigenvalue equation for natural frequency ω:(38)−ω2M+Kd¯=0.The natural frequencies and corresponding mode shapes of a given structure are often referred to as the dynamic characteristics of the structure.


### 3.3. Forced Vibration Analysis

In the forced vibration analysis, equation (34) is a function of both space and time. For simplicity, a constant damping matrix C which is assumed to be a linear combination of M and K can be utilized as(39)C=αM+βK,where α and β are the Rayleigh damping coefficients.There exist many numerical integration methods, such as Newmark method, Wilson method, generalized-α method [57], and so on. Here, the Bathe method, as briefly discussed as below [58–60], is used in this work. Two equal substeps are included in the complete time step h in Bathe method. In the first substep, the following equations are obtained by using the trapezoidal rule:(40)d˙t+h/2=4hdt+h/2−dt−d˙t,d¨t+h/2=16h2dt+h/2−dt−h2d˙t−d¨t.Substituting equation (40) into equation (33) yields the following algebraic equation:(41)16h2M+4hC+Kdt+h/2−M16h2dt+8hd˙t+d¨t−C4hdt+d˙t=fext.In the second substep, the following equations are expressed by the three-point Euler backward method:(42)d˙t+h=1hdt−4hdt+h/2+3hdt+h,d¨t+h=1hd˙t−4hd˙t+h/2+3hd˙t+h.Similarly, the control equation (34) at time t+h can also be expressed as(43)9h2M+3hC+Kdt+h+M1hd˙t−4hd˙t+h/2+3h2dt−12h2dt+h/2+C1hdt−4hdt+h/2=Q.Oncedt+h/2 and dt+h is solved by equations (41) and (43), d˙t+h/2, d¨t+h/2, d˙t+h, and d¨t+h can be obtained via equations (40) and (42).


## 3.1. Discrete Governing Equations

Here, we consider a deformable bodyΩ subjected to body force f, traction t on boundary Γt, and displacement boundary conditions u=u¯ on Γu, where Γ=Γt∪Γu, Γt∩Γu=∅. The virtual weak form of the dynamic equation can be expressed as(27)∫ΩδεTσdΩ−∫ΩδuTf−ρu¨−cu˙dΩ−∫ΓtδuTtdΓ=0,where ρ and c are the density and the coefficient of damping. ε and σ denote the strain and stress in the element. As the spatial discretization procedure in the FEM, the virtual displacements and strains within element can be expressed as(28)u=∑INIdI,δu=∑INIδdI,(29)ε=∑IBIdI,δε=∑IBIδdI,where dI=uIvIT is the nodal displacement vector and B represents the standard displacement gradient matrix. For 2D problems, it can be determined as follows:(30)BI=∇SNIx=NI,x00NI,yNI,yNI,x.The stress can be obtained by the following constitutive relation:(31)σ=Dε,where D is the matrix of the elastic constants of the material.Substituting equations (28) and (29)into equation (27) yields(32)∫ΩδdTBTσdΩ−∫ΩδdTNTf−ρu¨−cu˙dΩ−∫ΓtδdTNTtdΓ=0.According to the principle of virtual work, the discrete governing equation can be expressed as(33)∫ΩBTσdΩ−∫ΩNTf−ρu¨−cu˙dΩ−∫ΓtNTtdΓ=0,or(34)Md¨+Cd˙+Kd=fext,in which(35)fext=∫ΩNTfdΩ+∫ΓtNTtdΓ,M=∫ΩNTρNdΩ,K=∫ΩBTDBdΩ,C=∫ΩNTcNdΩ.


## 3.2. Free Vibration Analysis

The damping and external forces are not taken into account during free vibration analysis. Accordingly, equation (34) can be reduced to a system of homogeneous equations [1]:(36)Md¨+Kd=0.A general solution to this system is(37)d=d¯expiωt,where i denotes the imaginary unit, t is time, d¯ indicates the eigenvector, and ω is the natural frequency. Substituting equation (37) into equation (36) leads to the following eigenvalue equation for natural frequency ω:(38)−ω2M+Kd¯=0.The natural frequencies and corresponding mode shapes of a given structure are often referred to as the dynamic characteristics of the structure.


## 3.3. Forced Vibration Analysis

In the forced vibration analysis, equation (34) is a function of both space and time. For simplicity, a constant damping matrix C which is assumed to be a linear combination of M and K can be utilized as(39)C=αM+βK,where α and β are the Rayleigh damping coefficients.There exist many numerical integration methods, such as Newmark method, Wilson method, generalized-α method [57], and so on. Here, the Bathe method, as briefly discussed as below [58–60], is used in this work. Two equal substeps are included in the complete time step h in Bathe method. In the first substep, the following equations are obtained by using the trapezoidal rule:(40)d˙t+h/2=4hdt+h/2−dt−d˙t,d¨t+h/2=16h2dt+h/2−dt−h2d˙t−d¨t.Substituting equation (40) into equation (33) yields the following algebraic equation:(41)16h2M+4hC+Kdt+h/2−M16h2dt+8hd˙t+d¨t−C4hdt+d˙t=fext.In the second substep, the following equations are expressed by the three-point Euler backward method:(42)d˙t+h=1hdt−4hdt+h/2+3hdt+h,d¨t+h=1hd˙t−4hd˙t+h/2+3hd˙t+h.Similarly, the control equation (34) at time t+h can also be expressed as(43)9h2M+3hC+Kdt+h+M1hd˙t−4hd˙t+h/2+3h2dt−12h2dt+h/2+C1hdt−4hdt+h/2=Q.Oncedt+h/2 and dt+h is solved by equations (41) and (43), d˙t+h/2, d¨t+h/2, d˙t+h, and d¨t+h can be obtained via equations (40) and (42).


## 4. Numerical Examples

To evaluate the accuracy of the method, the following displacement error norm and the strain energy norm are defined:(44)ed=∑i=1Nduiexact−uinum2∑i=1Nduiexact2,ee=Enum−EexactEexact,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 strain energy of the numerical solutionEnum and the total strain energy of the exact solution Eexact are defined as follows:(45)Enum=12uTKu,Eexact=12∑i=1Ne∫ΩiεiTDεidΩ,where εi is the strain of the exact solution. In the actual computation on equation (45), we used a very fine mesh (Ne⟶∞) to calculate the “exact” strain energy Eexact.In order to assess accuracy and convergence of each method for free vibration analysis, the relative error in the natural frequency is defined as follows:(46)Re=fnum−freffref,where the superscript “ref” represents the reference solution.


### 4.1. Example 1: Patch Test

Patch tests often serve as useful tools to evaluate the convergence of a numerical method based on the Galerkin weak form. Here, we use the generalized patch test introduced in [61], in which the coefficients of polynomial are not specifically predefined. The patch test is regarded as pass in the case that the maximum eigenvalue of the resulted matrix tends to zero.As listed in [62], for the patch test, we can obtain(47)uex=∑aijpijx=pxa,pijx=xiyj,i+j≤m,where px and a represent the row vector of monomials and the corresponding vector of coefficients, respectively, and m denotes the arbitrary maximum order of the monomials. Regarding each of these monomials in equation (47) as separate exact solution, a set of M numerical solutions can be obtained as(48)unex=pnx,n=1,2,…,M,M=m+1m+22.For the distributed nodes in the patch domain, a pointwise error vector can be evaluated as(49)en=Unex−Un,n=1,2,…,M,where Unex and Un denote the numerical solution and the corresponding exact one, respectively.Summing the resulted errors for each monomial, we can then obtain the errore as(50)e=∑n=1Manen=Ea,E=e1e2⋯eM.The Euclidian error can be evaluated as(51)e2=1NdeTe=1NdaTQa,Q=ETE,where Q represents a positive semidefinite matrix, and the maximum eigenvalue is derived as(52)Q−λmaxI=0,where λmax is the maximum eigenvalue and I is the unit matrix. On condition that the patch test is passed, the error norm and λmax both will be zero.In this study, we conduct a patch test based on polynomials of first order, hence takingm=1 in equation (47) leads to(53)uex=a0+a1x+a2y.We choose a [−1,1] × [−1,1] square patch domain using 10 × 10 regular nodes and 100 irregular nodes as shown in Figure3. For irregular interior nodes, the corresponding coordinates can be obtained as(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. Table 2 shows that the present method can produce the liner displacement field.Figure 3
Regular and irregular meshes for patch test: (a) regular; (b) irregularity coefficientαr=0.1; (c) irregularity coefficient αr=0.2; (d) irregularity coefficient αr=0.3; (e) irregularity coefficient αr=0.4.
(a)
(b)
(c)
(d)
(e)Table 2
Displacement error norm of numerical results for standard patch test on different meshes.
Solution
Regular
Irregularαr=0.1
Irregular,αr=0.2
Irregular,αr=0.3
Irregular,αr=0.4
λ
max
5.45 × 10−5
8.72 × 10−5
1.72 × 10−4
2.59 × 10−4
3.45 × 10−4


### 4.2. Example 2: Cantilever Beam Subjected to Tip-Shear Force

A cantilever beam with lengthL and height D is investigated as a benchmark problem. In this test, the system is subjected to a parabolic traction at the free end as shown in Figure 4. 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 [63]:(55)ux=Py6EI6L−3xx+2+vy2−D24,(56)uy=−P6EI3vy2L−x+144+5vD2x+x23L−x,where the moment of inertia I of the beam is given by I=D3/12.Figure 4
Cantilever beam problem.The corresponding stresses are(57)σx=PL−xyI,(58)σy=0,(59)τxy=−P2ID24−y2.The parameters areE=3.0×107kPa, 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 (55) and (56). The distributed parabolic shear stresses, expressed as equations (57)–(59), acted on the right boundary.To compare with other methods for convenience in the following examples, we use FE-Meshfree (QS) to denote the presented method, i.e., FE-Meshfree with QS kernel function. Figure5 shows the steady condition number for different θ, especially θ>10. As shown in Figure 6, if the irregularity factor αr<0.2, the errors of the deflection of point A are almost same. If the irregularity factor αr>0.2, the error increases as the factor increases. 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-QUAD4 modes are overly stiff and hence produced lower bounds. The FE-Meshfree (QS) shows a very close-to-exact stiffness and hence very accurate results. The strain energy, displacement, and energy error norms of the FE-Meshfree (QS) are all better than those of the FEM-QUAD4.Figure 5
Condition number.Figure 6
Deflection of point A for irregular mesh.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. FE-Meshfree (QS) takes longer to compute than FEM-QUAD4. 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 in FE-Meshfree (QS) compared to FEM-QUAD4.Figure 10
Computation time.Figure 11
Energy error norm with respect to computation time.Figure 12
Displacement error norm with respect to computation time.


### 4.3. Example 3: Cantilever Beam Free Vibration Analysis

As shown in Figure13, a cantilever beam fixed on the left end with length L = 100 m and height H = 10 m is conducted here as a benchmark problem. The plane stress condition is considered. The relative geometrical and material parameters can be summarized as follows: thickness t = 1.0 m, Young’s modulus E=2.1×104N/m2, Poisson’s ratio v=0.3, and mass density ρ=8.0×10−10kg/m3. Based on the Euler–Bernoulli beam theory, the fundamental frequency f1=0.08276×104Hz can be obtained as a reference.Figure 13
Cantilever beam problem.


#### 4.3.1. Effect of the Correlation Parameterθ

The effects of the correlation parameterθ in the modified RPIM on the solution accuracy and stability property are investigated here. The parameter sensitivity analysis of the free vibration of a cantilever beam is conducted by employing 50 × 5 meshes (Figure 14), in which the correlation parameter θ is set to be ranged from 1 to 1,000. The first 12 frequencies of the problems determined using various θ, for comparison against the reference fundamental frequency solutions are tabulated in Table 3. Results show great agreement with the exact solutions and high stability regardless of the correlation parameter. Furthermore, nearly the same results are obtained for θ≥10. This unique feature of the presented radial basis function enables constant results for any correlation parameter. Accordingly, correlation parameter θ=10 is used for all subsequent computations.Figure 14
Domain discretization of beam using four-node element with regular elements.Table 3
First 12 natural frequencies (×104 Hz) for different correlation parameters θ.
θ
=
1
θ
=
5
θ
=
10
θ
=
50
θ
=
100
θ
=
500
θ
=
1000
f
1
0.0821
0.0822
0.0823
0.0823
0.0823
0.0823
0.0823
f
2
0.4937
0.4936
0.4937
0.4938
0.4938
0.4938
0.4938
f
3
1.2826
1.2826
1.2826
1.2826
1.2826
1.2826
1.2826
f
4
1.3038
1.3008
1.3010
1.3011
1.3012
1.3012
1.3012
f
5
2.3770
2.3653
2.3655
2.3656
2.3657
2.3657
2.3657
f
6
3.6374
3.6093
3.6094
3.6094
3.6094
3.6095
3.6095
f
7
3.8458
3.8450
3.8450
3.8450
3.8450
3.8450
3.8450
f
8
5.0252
4.9719
4.9716
4.9715
4.9715
4.9714
4.9714
f
9
6.4024
6.3985
6.3984
6.3983
6.3983
6.3983
6.3983
f
10
6.5005
6.4124
6.4116
6.4112
6.4111
6.4111
6.4111
f
11
8.0361
7.9033
7.9019
7.9009
7.9008
7.9008
7.9008
f
12
8.9455
8.9350
8.9347
8.9345
8.9345
8.9345
8.9345


#### 4.3.2. Convergence Study

For the convergence rates, several types of regular meshes using FEM-QUAD4, FE-Meshfree (MQ), FE-Meshfree (Exp), are illustrated against the proposed method in this part. Four discrete mesh conditions with regular grids are considered as follows: 10 × 1, 20 × 2, 50 × 5, and 100 × 10. The first 12 natural frequencies computed for the four meshes are tabulated in Tables4–7. The reference results shown in the last column of Table 4 are gathered in the ABAQUS QUAD8 with 200 × 20 elements. Results show that the proposed method outperforms the other methods. In addition, the performance of FE-Meshfree (MQ) and FE-Meshfree (Exp) is also sensitive to the corresponding relation parameters.Table 4
Natural frequencies (×104 Hz) of a cantilever beam for 10 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(10 × 1) 10
22
0.1000
0.0824
0.0825
0.0909
0.1130
0.0848
0.0822
0.6077
0.5058
0.5050
0.5455
0.6732
0.5120
0.4931
1.2863
1.2846
1.2844
1.2849
1.2858
1.2839
1.2823
1.6422
1.3672
1.3621
1.4281
1.7140
1.3539
1.2990
3.0958
2.5583
2.5437
2.5729
2.9437
2.4719
2.3607
3.8918
3.8421
3.8430
3.8384
3.8596
3.8071
3.6002
4.9311
4.0208
3.9947
3.9079
4.2015
3.8405
3.8442
6.5947
5.6927
5.6693
5.4218
5.4614
5.3572
4.9566
7.1110
6.3692
6.3762
6.3422
6.4073
6.3656
6.3897
9.4517
7.5079
7.5313
7.1592
6.8653
7.1909
6.3972
9.5747
8.8600
8.8821
8.7702
8.6496
8.8518
7.8719
12.188
9.4823
9.6067
9.1641
9.8540
9.4100
8.9331Table 5
Natural frequencies (×104 Hz) of a cantilever beam for 40 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(20 × 2) 40
63
0.0872
0.0824
0.0825
0.0875
0.0981
0.0849
0.0822
0.5263
0.4981
0.4986
0.5249
0.5821
0.5112
0.4931
1.2837
1.2834
1.2834
1.2836
1.2839
1.2831
1.2823
1.4010
1.3233
1.3255
1.3826
1.5161
1.3523
1.2990
2.5816
2.4265
2.4325
2.5107
2.7232
2.4680
2.3607
3.8573
3.7297
3.7423
3.8232
3.8472
3.7767
3.6002
4.0002
3.8469
3.8472
3.8453
4.1064
3.8463
3.8442
5.6040
5.1688
5.1903
5.2518
5.5903
5.2123
4.9566
6.4492
6.3995
6.4018
6.3901
6.3954
6.3995
6.3897
7.3601
6.7004
6.7325
6.7524
7.1240
6.7308
6.3972
9.0702
8.2958
8.3397
8.2973
8.6741
8.3043
7.8719
9.2470
8.9314
8.9379
8.9054
8.9147
8.9338
8.9331Table 6
Natural frequencies (×104 Hz) of a cantilever beam for 250 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(50 × 5) 250
306
0.0830
0.0823
0.0822
0.0831
0.0852
0.0827
0.0822
0.4989
0.4937
0.4936
0.4986
0.5104
0.4962
0.4931
1.2827
1.2826
1.2826
1.2827
1.2828
1.2825
1.2823
1.3168
1.3010
1.3005
1.3127
1.3421
1.3063
1.2990
2.3992
2.3655
2.3640
2.3839
2.4343
2.3726
2.3607
3.6701
3.6094
3.6064
3.6331
3.7057
3.6165
3.6002
3.8467
3.8450
3.8449
3.8448
3.8449
3.8445
3.8442
5.0697
4.9716
4.9667
4.9984
5.0934
4.9765
4.9566
6.4062
6.3984
6.3983
6.3965
6.3955
6.3971
6.3897
6.5590
6.4116
6.4043
6.4391
6.5558
6.4121
6.3972
8.1119
7.9019
7.8918
7.9274
8.0652
7.8957
7.8719
8.9562
8.9347
8.9348
8.9290
8.9251
8.9320
8.9331Table 7
Natural frequencies (×104 Hz) of a cantilever beam for 1000 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(100 × 10) 1000
1111
0.0824
0.0822
0.0822
0.0824
0.0830
0.0824
0.0822
0.4947
0.4933
0.4933
0.4943
0.4978
0.4940
0.4931
1.2825
1.2824
1.2824
1.2824
1.2825
1.2824
1.2823
1.3037
1.2995
1.2994
1.3019
1.3110
1.3011
1.2990
2.3707
2.3618
2.3614
2.3654
2.3817
2.3641
2.3607
3.6183
3.6023
3.6015
3.6068
3.6312
3.6050
3.6002
3.8450
3.8445
3.8444
3.8445
3.8446
3.8443
3.8442
4.9858
4.9600
4.9586
4.9649
4.9981
4.9627
4.9566
6.3996
6.3944
6.3923
6.3972
6.3970
6.3968
6.3897
6.4332
6.3976
6.3975
6.3992
6.4415
6.3972
6.3972
7.9335
7.8783
7.8753
7.8822
7.9341
7.8797
7.8719
8.9391
8.9337
8.9337
8.9326
8.9313
8.9328
8.9331Figure15 shows an error plot of the first two natural frequencies obtained using the proposed method as well as FEM-QUAD4, FE-Meshfree (MQ), and FE-Meshfree (Exp). Results show that the corresponding error of the proposed method is generally lower than that of other methods. Even for the coarse mesh, the results of the proposed method still show better consistence with the reference solution. Figure 16 displays the first 12 eigenmodes obtained for FE-Meshfree (QS) using the 20 × 2 meshes. These modeshape plots are in accordance with previously published results for ES-FEM [21].Figure 15
Computed first two natural frequencies. (a) Mode 1. (b) Mode 2.
(a)
(b)Figure 16
First 12 modeshapes of cantilever beam.


#### 4.3.3. Sensitivity to Distorted Meshes

The cantilever beam with irregular mesh shown in Figure13 is also employed to verify the effectiveness of the proposed method. As shown in Figure 17, four distorted meshes with 50 × 5 elements are used to assess the influence of mesh quality. Table 8 gives the comparison results of the computed natural frequencies for the first 12 natural frequencies. A comparison of relative error in the computed natural frequencies is shown in Figure 18. As αr increases, the accuracy of FEM-QUAD4 decreases, whereas the FE-Meshfree (QS) does not show markedly deviation from their original (high) accuracy, indicating great robustness of the proposed method to mesh distortion phenomena.Figure 17
Beam domain discretization using four-node element with irregular elements. (a)αr=0.1. (b) αr=0.2. (c) αr=0.3. (d) αr=0.4.
(a)
(b)
(c)
(d)Table 8
First 12 natural frequencies (×104 Hz) of a cantilever beam using irregular nodal distribution (50 × 5 elements).
Method
α
r
=
0
α
r
=
0.1
α
r
=
0.2
α
r
=
0.3
α
r
=
0.4
Reference (ABAQUS QUAD8)
FEM-QUAD4
0.0830
0.0831
0.0834
0.0835
0.0838
0.0822
0.4989
0.4992
0.5016
0.5019
0.5036
0.4931
1.2827
1.2827
1.2827
1.2827
1.2827
1.2823
1.3168
1.3176
1.3243
1.3245
1.3297
1.2990
2.3992
2.4007
2.4114
2.4139
2.4242
2.3607
3.6701
3.6725
3.6914
3.6963
3.7096
3.6002
3.8467
3.8468
3.8470
3.8470
3.8471
3.8442
5.0697
5.0728
5.0974
5.1059
5.1249
4.9566
6.4062
6.4063
6.4072
6.4072
6.4074
6.3897
6.5590
6.5643
6.5999
6.6089
6.6360
6.3972
8.1119
8.1186
8.1651
8.1733
8.2060
7.8719
8.9562
8.9565
8.9528
8.9585
8.9606
8.9331
FE-Meshfree (QS,θ=10)
0.0823
0.0822
0.0822
0.0822
0.0822
0.0822
0.4937
0.4937
0.4936
0.4936
0.4936
0.4931
1.2826
1.2826
1.2825
1.2824
1.2822
1.2823
1.3010
1.3010
1.3008
1.3009
1.3010
1.2990
2.3655
2.3655
2.3655
2.3655
2.3661
2.3607
3.6094
3.6092
3.6098
3.6095
3.6113
3.6002
3.8450
3.8450
3.8445
3.8444
3.8441
3.8442
4.9716
4.9715
4.9722
4.9728
4.9755
4.9566
6.3984
6.3984
6.3978
6.3974
6.3966
6.3897
6.4116
6.4114
6.4126
6.4142
6.4179
6.3972
7.9019
7.9017
7.9043
7.9039
7.9110
7.8719
8.9347
8.9346
8.9337
8.9333
8.9331
8.9331Figure 18
Error in the computed first two frequencies for distortion sensitivity test.
(a)
(b)


### 4.4. Example 4: Shear Wall Free Vibration Analysis

In this part, a shear wall with four openings, as shown in Figure19 is investigated. The parameters in the computation include Young’s modulus E=10000 N/m2, Poisson’s ratio v=0.2, thickness t=1.0 m, and mass density ρ=1.0 kg/m3. The bottom edge is fully clamped, and a plane stress case is considered here. The mesh for computation with 476 elements and 559 nodes is shown in Figure 20. This problem was analyzed previously using a boundary element method by Brabbia et al. [64] and also by Liu and Gu [13] using a local radial point interpolation method with the multiquadrics radial function. The commercial software applications, ANSYS and ABAQUS, are employed here for comparison.Figure 19
Shear wall with four openings.Figure 20
Model of shear wall with four openings.The natural frequencies of the first eight modes are tabulated in Table9. The results show great agreement with those obtained by other methods. Figure 21 displays the first eight eigenmodes obtained by the proposed method, from which it can be observed that the modeshapes plot match well with those of ES-FEM [21].Table 9
Natural frequencies (rad/s) of shear wall.
Model
FEM-QUAD4
FE-Meshfree (QS)
LRPIM (MQ) [13]
Brebbia (BEM) [64]
ANSYS [49]
ABAQUS
QUAD4-CNS [50]
1
2.074
2.085
2.086
2.079
2.057
2.073
2.094
2
7.114
7.117
7.152
7.181
7.067
7.096
7.124
3
7.627
7.632
7.647
7.644
7.620
7.625
7.590
4
11.978
12.063
12.019
11.833
11.840
11.938
12.192
5
15.444
15.410
15.628
15.947
15.313
15.341
15.451
6
18.429
18.351
18.548
18.644
18.342
18.345
18.330
7
19.991
19.882
20.085
20.268
19.887
19.876
19.858
8
22.348
22.258
22.564
22.765
22.236
22.210
22.225Figure 21
First 8 modeshapes of shear wall with four openings. (a) Mode 1. (b) Mode 2. (c) Mode 3. (d) Mode 4. (e) Mode 5. (f) Mode 6. (g) Mode 7. (h) Mode 8.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


### 4.5. Example 5: Connecting Rod Free Vibration Analysis

Here, we perform free vibration analysis of a connecting rod with a relatively complex configuration in which the boundary conditions can be clearly observed in Figure22. The rod is discretized with 723 irregularly distributed nodes, as shown in Figure 23. The plane stress problem is considered with E = 10 GPa, v = 0.3, and ρ = 7.8 × 103 kg/m3. A reference solution is also obtained in ABAQUS QUAD8 with 24,149 elements. Table 10 gives the first 10 natural frequencies, where FE-Meshfree (QS) yields comparable results to ABAQUS (similar to the previous example).Figure 22
Geometric model and boundary conditions of automotive connecting rod.Figure 23
Domain discretization of connecting bar using four-node elements.Table 10
First ten natural frequencies (×106 rad/s) of connecting rod.
Mode
FEM-QUAD4
FE-Meshfree (QS)
Reference (ABAQUS QUAD8)
1
1.3673
1.3592
1.3418
2
2.8093
2.7995
2.7042
3
3.9773
3.9289
3.7932
4
4.1908
4.1890
4.1078
5
5.7172
5.6162
5.5291
6
7.3185
7.3137
7.0938
7
7.8352
7.6765
7.5696
8
8.2621
8.2325
8.1509
9
8.4762
8.4387
8.3444
10
8.6747
8.6310
8.5351


### 4.6. Example 6: Cantilever Beam Implicit Forced Vibration Analysis

The benchmark cantilever beam and mesh shown in Figures24 and 25 are investigated via the Bathe method. The plane strain problem is considered with numerical parameters L=4.0 m, H=1.0 m, t=1.0 m, E=1.0 N/m2, v=0.3, ρ=1.0 kg/m3, α=0.005, and β=0.272. The domain is represented by 10 × 2 elements. The first loading is in the harmonic form with ft=cosωft. The corresponding dynamic responses can be shown in Figures 26 and 27 with ωf=0.01 and 0.05, respectively. Then, a constant step load ft=1 is added at the apex from t=0. The time step Δt=2.0 is used for time integration. Both damping effect and no damping effect are considered here. A reference solution is obtained in ABAQUS QUAD8 with 2,600 elements.Figure 24
Cantilever beam.Figure 25
Domain discretization of cantilever beam.Figure 26
Transient response of cantilever beam subjected to harmonic loading withωf=0.01: (a) without damping; (b) with damping.
(a)
(b)Figure 27
Transient response of cantilever beam subjected to harmonic loading withωf=0.05: (a) without damping; (b) with damping.
(a)
(b)Per the dynamic responses in Figures26–28, both methods tested provided nearly identical vibration frequency with the amplitude of FE-Meshfree (QS) being closer to the reference solution. This indicates that the FE-Meshfree (QS) can be readily applied to actual vibration analyses with excellent accuracy, which may be partially due to the fact that the FE-Meshfree (QS) model owns very close-to-exact stiffness.Figure 28
Transient response of cantilever beam subjected to constant step loading: (a) without damping; (b) with damping.
(a)
(b)


### 4.7. Example 7: Spherical Shell Implicit Forced Vibration Analysis

As shown in Figure29, a spherical shell subjected to a concentrated loading at its apex is examined. In total, 30 × 3 asymmetric elements are used (Figure 30) with the numerical parameters R=12 m, t=0.1 m, ϕ=10.9°, E=1.0 N/m2, v=0.3, and ρ=1.0 kg/m3. A reference solution is obtained by ABAQUS QUAD8 with 2,000 elements. The first loading is applied in the harmonic form ft=cosωft. Its dynamic responses can be observed from Figures 31 and 32 with ωf=0.01 and 0.05, respectively. A constant step load ft=1 is then added at the apex from t=0. The time step Δt=2.0 is used for time integration. The dynamic responses shown in Figures 31–33 indicate that FE-Meshfree (QS) is more accurate than FEM-QUAD4. It yields results very close to the reference one for both damping and no damping effect situations.Figure 29
Spherical shell.Figure 30
Domain discretization of half of the spherical shell.Figure 31
Transient response of spherical shell subjected to harmonic loading withωf=0.01: (a) without damping; (b) with damping.
(a)
(b)Figure 32
Transient response of spherical shell subjected to harmonic loading withωf=0.05: (a) without damping; (b) with damping.
(a)
(b)Figure 33
Transient response of spherical shell subjected to constant step loading: (a) without damping; (b) with damping.
(a)
(b)


## 4.1. Example 1: Patch Test

Patch tests often serve as useful tools to evaluate the convergence of a numerical method based on the Galerkin weak form. Here, we use the generalized patch test introduced in [61], in which the coefficients of polynomial are not specifically predefined. The patch test is regarded as pass in the case that the maximum eigenvalue of the resulted matrix tends to zero.As listed in [62], for the patch test, we can obtain(47)uex=∑aijpijx=pxa,pijx=xiyj,i+j≤m,where px and a represent the row vector of monomials and the corresponding vector of coefficients, respectively, and m denotes the arbitrary maximum order of the monomials. Regarding each of these monomials in equation (47) as separate exact solution, a set of M numerical solutions can be obtained as(48)unex=pnx,n=1,2,…,M,M=m+1m+22.For the distributed nodes in the patch domain, a pointwise error vector can be evaluated as(49)en=Unex−Un,n=1,2,…,M,where Unex and Un denote the numerical solution and the corresponding exact one, respectively.Summing the resulted errors for each monomial, we can then obtain the errore as(50)e=∑n=1Manen=Ea,E=e1e2⋯eM.The Euclidian error can be evaluated as(51)e2=1NdeTe=1NdaTQa,Q=ETE,where Q represents a positive semidefinite matrix, and the maximum eigenvalue is derived as(52)Q−λmaxI=0,where λmax is the maximum eigenvalue and I is the unit matrix. On condition that the patch test is passed, the error norm and λmax both will be zero.In this study, we conduct a patch test based on polynomials of first order, hence takingm=1 in equation (47) leads to(53)uex=a0+a1x+a2y.We choose a [−1,1] × [−1,1] square patch domain using 10 × 10 regular nodes and 100 irregular nodes as shown in Figure3. For irregular interior nodes, the corresponding coordinates can be obtained as(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. Table 2 shows that the present method can produce the liner displacement field.Figure 3
Regular and irregular meshes for patch test: (a) regular; (b) irregularity coefficientαr=0.1; (c) irregularity coefficient αr=0.2; (d) irregularity coefficient αr=0.3; (e) irregularity coefficient αr=0.4.
(a)
(b)
(c)
(d)
(e)Table 2
Displacement error norm of numerical results for standard patch test on different meshes.
Solution
Regular
Irregularαr=0.1
Irregular,αr=0.2
Irregular,αr=0.3
Irregular,αr=0.4
λ
max
5.45 × 10−5
8.72 × 10−5
1.72 × 10−4
2.59 × 10−4
3.45 × 10−4


## 4.2. Example 2: Cantilever Beam Subjected to Tip-Shear Force

A cantilever beam with lengthL and height D is investigated as a benchmark problem. In this test, the system is subjected to a parabolic traction at the free end as shown in Figure 4. 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 [63]:(55)ux=Py6EI6L−3xx+2+vy2−D24,(56)uy=−P6EI3vy2L−x+144+5vD2x+x23L−x,where the moment of inertia I of the beam is given by I=D3/12.Figure 4
Cantilever beam problem.The corresponding stresses are(57)σx=PL−xyI,(58)σy=0,(59)τxy=−P2ID24−y2.The parameters areE=3.0×107kPa, 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 (55) and (56). The distributed parabolic shear stresses, expressed as equations (57)–(59), acted on the right boundary.To compare with other methods for convenience in the following examples, we use FE-Meshfree (QS) to denote the presented method, i.e., FE-Meshfree with QS kernel function. Figure5 shows the steady condition number for different θ, especially θ>10. As shown in Figure 6, if the irregularity factor αr<0.2, the errors of the deflection of point A are almost same. If the irregularity factor αr>0.2, the error increases as the factor increases. 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-QUAD4 modes are overly stiff and hence produced lower bounds. The FE-Meshfree (QS) shows a very close-to-exact stiffness and hence very accurate results. The strain energy, displacement, and energy error norms of the FE-Meshfree (QS) are all better than those of the FEM-QUAD4.Figure 5
Condition number.Figure 6
Deflection of point A for irregular mesh.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. FE-Meshfree (QS) takes longer to compute than FEM-QUAD4. 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 in FE-Meshfree (QS) compared to FEM-QUAD4.Figure 10
Computation time.Figure 11
Energy error norm with respect to computation time.Figure 12
Displacement error norm with respect to computation time.


## 4.3. Example 3: Cantilever Beam Free Vibration Analysis

As shown in Figure13, a cantilever beam fixed on the left end with length L = 100 m and height H = 10 m is conducted here as a benchmark problem. The plane stress condition is considered. The relative geometrical and material parameters can be summarized as follows: thickness t = 1.0 m, Young’s modulus E=2.1×104N/m2, Poisson’s ratio v=0.3, and mass density ρ=8.0×10−10kg/m3. Based on the Euler–Bernoulli beam theory, the fundamental frequency f1=0.08276×104Hz can be obtained as a reference.Figure 13
Cantilever beam problem.


### 4.3.1. Effect of the Correlation Parameterθ

The effects of the correlation parameterθ in the modified RPIM on the solution accuracy and stability property are investigated here. The parameter sensitivity analysis of the free vibration of a cantilever beam is conducted by employing 50 × 5 meshes (Figure 14), in which the correlation parameter θ is set to be ranged from 1 to 1,000. The first 12 frequencies of the problems determined using various θ, for comparison against the reference fundamental frequency solutions are tabulated in Table 3. Results show great agreement with the exact solutions and high stability regardless of the correlation parameter. Furthermore, nearly the same results are obtained for θ≥10. This unique feature of the presented radial basis function enables constant results for any correlation parameter. Accordingly, correlation parameter θ=10 is used for all subsequent computations.Figure 14
Domain discretization of beam using four-node element with regular elements.Table 3
First 12 natural frequencies (×104 Hz) for different correlation parameters θ.
θ
=
1
θ
=
5
θ
=
10
θ
=
50
θ
=
100
θ
=
500
θ
=
1000
f
1
0.0821
0.0822
0.0823
0.0823
0.0823
0.0823
0.0823
f
2
0.4937
0.4936
0.4937
0.4938
0.4938
0.4938
0.4938
f
3
1.2826
1.2826
1.2826
1.2826
1.2826
1.2826
1.2826
f
4
1.3038
1.3008
1.3010
1.3011
1.3012
1.3012
1.3012
f
5
2.3770
2.3653
2.3655
2.3656
2.3657
2.3657
2.3657
f
6
3.6374
3.6093
3.6094
3.6094
3.6094
3.6095
3.6095
f
7
3.8458
3.8450
3.8450
3.8450
3.8450
3.8450
3.8450
f
8
5.0252
4.9719
4.9716
4.9715
4.9715
4.9714
4.9714
f
9
6.4024
6.3985
6.3984
6.3983
6.3983
6.3983
6.3983
f
10
6.5005
6.4124
6.4116
6.4112
6.4111
6.4111
6.4111
f
11
8.0361
7.9033
7.9019
7.9009
7.9008
7.9008
7.9008
f
12
8.9455
8.9350
8.9347
8.9345
8.9345
8.9345
8.9345


### 4.3.2. Convergence Study

For the convergence rates, several types of regular meshes using FEM-QUAD4, FE-Meshfree (MQ), FE-Meshfree (Exp), are illustrated against the proposed method in this part. Four discrete mesh conditions with regular grids are considered as follows: 10 × 1, 20 × 2, 50 × 5, and 100 × 10. The first 12 natural frequencies computed for the four meshes are tabulated in Tables4–7. The reference results shown in the last column of Table 4 are gathered in the ABAQUS QUAD8 with 200 × 20 elements. Results show that the proposed method outperforms the other methods. In addition, the performance of FE-Meshfree (MQ) and FE-Meshfree (Exp) is also sensitive to the corresponding relation parameters.Table 4
Natural frequencies (×104 Hz) of a cantilever beam for 10 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(10 × 1) 10
22
0.1000
0.0824
0.0825
0.0909
0.1130
0.0848
0.0822
0.6077
0.5058
0.5050
0.5455
0.6732
0.5120
0.4931
1.2863
1.2846
1.2844
1.2849
1.2858
1.2839
1.2823
1.6422
1.3672
1.3621
1.4281
1.7140
1.3539
1.2990
3.0958
2.5583
2.5437
2.5729
2.9437
2.4719
2.3607
3.8918
3.8421
3.8430
3.8384
3.8596
3.8071
3.6002
4.9311
4.0208
3.9947
3.9079
4.2015
3.8405
3.8442
6.5947
5.6927
5.6693
5.4218
5.4614
5.3572
4.9566
7.1110
6.3692
6.3762
6.3422
6.4073
6.3656
6.3897
9.4517
7.5079
7.5313
7.1592
6.8653
7.1909
6.3972
9.5747
8.8600
8.8821
8.7702
8.6496
8.8518
7.8719
12.188
9.4823
9.6067
9.1641
9.8540
9.4100
8.9331Table 5
Natural frequencies (×104 Hz) of a cantilever beam for 40 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(20 × 2) 40
63
0.0872
0.0824
0.0825
0.0875
0.0981
0.0849
0.0822
0.5263
0.4981
0.4986
0.5249
0.5821
0.5112
0.4931
1.2837
1.2834
1.2834
1.2836
1.2839
1.2831
1.2823
1.4010
1.3233
1.3255
1.3826
1.5161
1.3523
1.2990
2.5816
2.4265
2.4325
2.5107
2.7232
2.4680
2.3607
3.8573
3.7297
3.7423
3.8232
3.8472
3.7767
3.6002
4.0002
3.8469
3.8472
3.8453
4.1064
3.8463
3.8442
5.6040
5.1688
5.1903
5.2518
5.5903
5.2123
4.9566
6.4492
6.3995
6.4018
6.3901
6.3954
6.3995
6.3897
7.3601
6.7004
6.7325
6.7524
7.1240
6.7308
6.3972
9.0702
8.2958
8.3397
8.2973
8.6741
8.3043
7.8719
9.2470
8.9314
8.9379
8.9054
8.9147
8.9338
8.9331Table 6
Natural frequencies (×104 Hz) of a cantilever beam for 250 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(50 × 5) 250
306
0.0830
0.0823
0.0822
0.0831
0.0852
0.0827
0.0822
0.4989
0.4937
0.4936
0.4986
0.5104
0.4962
0.4931
1.2827
1.2826
1.2826
1.2827
1.2828
1.2825
1.2823
1.3168
1.3010
1.3005
1.3127
1.3421
1.3063
1.2990
2.3992
2.3655
2.3640
2.3839
2.4343
2.3726
2.3607
3.6701
3.6094
3.6064
3.6331
3.7057
3.6165
3.6002
3.8467
3.8450
3.8449
3.8448
3.8449
3.8445
3.8442
5.0697
4.9716
4.9667
4.9984
5.0934
4.9765
4.9566
6.4062
6.3984
6.3983
6.3965
6.3955
6.3971
6.3897
6.5590
6.4116
6.4043
6.4391
6.5558
6.4121
6.3972
8.1119
7.9019
7.8918
7.9274
8.0652
7.8957
7.8719
8.9562
8.9347
8.9348
8.9290
8.9251
8.9320
8.9331Table 7
Natural frequencies (×104 Hz) of a cantilever beam for 1000 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(100 × 10) 1000
1111
0.0824
0.0822
0.0822
0.0824
0.0830
0.0824
0.0822
0.4947
0.4933
0.4933
0.4943
0.4978
0.4940
0.4931
1.2825
1.2824
1.2824
1.2824
1.2825
1.2824
1.2823
1.3037
1.2995
1.2994
1.3019
1.3110
1.3011
1.2990
2.3707
2.3618
2.3614
2.3654
2.3817
2.3641
2.3607
3.6183
3.6023
3.6015
3.6068
3.6312
3.6050
3.6002
3.8450
3.8445
3.8444
3.8445
3.8446
3.8443
3.8442
4.9858
4.9600
4.9586
4.9649
4.9981
4.9627
4.9566
6.3996
6.3944
6.3923
6.3972
6.3970
6.3968
6.3897
6.4332
6.3976
6.3975
6.3992
6.4415
6.3972
6.3972
7.9335
7.8783
7.8753
7.8822
7.9341
7.8797
7.8719
8.9391
8.9337
8.9337
8.9326
8.9313
8.9328
8.9331Figure15 shows an error plot of the first two natural frequencies obtained using the proposed method as well as FEM-QUAD4, FE-Meshfree (MQ), and FE-Meshfree (Exp). Results show that the corresponding error of the proposed method is generally lower than that of other methods. Even for the coarse mesh, the results of the proposed method still show better consistence with the reference solution. Figure 16 displays the first 12 eigenmodes obtained for FE-Meshfree (QS) using the 20 × 2 meshes. These modeshape plots are in accordance with previously published results for ES-FEM [21].Figure 15
Computed first two natural frequencies. (a) Mode 1. (b) Mode 2.
(a)
(b)Figure 16
First 12 modeshapes of cantilever beam.


### 4.3.3. Sensitivity to Distorted Meshes

The cantilever beam with irregular mesh shown in Figure13 is also employed to verify the effectiveness of the proposed method. As shown in Figure 17, four distorted meshes with 50 × 5 elements are used to assess the influence of mesh quality. Table 8 gives the comparison results of the computed natural frequencies for the first 12 natural frequencies. A comparison of relative error in the computed natural frequencies is shown in Figure 18. As αr increases, the accuracy of FEM-QUAD4 decreases, whereas the FE-Meshfree (QS) does not show markedly deviation from their original (high) accuracy, indicating great robustness of the proposed method to mesh distortion phenomena.Figure 17
Beam domain discretization using four-node element with irregular elements. (a)αr=0.1. (b) αr=0.2. (c) αr=0.3. (d) αr=0.4.
(a)
(b)
(c)
(d)Table 8
First 12 natural frequencies (×104 Hz) of a cantilever beam using irregular nodal distribution (50 × 5 elements).
Method
α
r
=
0
α
r
=
0.1
α
r
=
0.2
α
r
=
0.3
α
r
=
0.4
Reference (ABAQUS QUAD8)
FEM-QUAD4
0.0830
0.0831
0.0834
0.0835
0.0838
0.0822
0.4989
0.4992
0.5016
0.5019
0.5036
0.4931
1.2827
1.2827
1.2827
1.2827
1.2827
1.2823
1.3168
1.3176
1.3243
1.3245
1.3297
1.2990
2.3992
2.4007
2.4114
2.4139
2.4242
2.3607
3.6701
3.6725
3.6914
3.6963
3.7096
3.6002
3.8467
3.8468
3.8470
3.8470
3.8471
3.8442
5.0697
5.0728
5.0974
5.1059
5.1249
4.9566
6.4062
6.4063
6.4072
6.4072
6.4074
6.3897
6.5590
6.5643
6.5999
6.6089
6.6360
6.3972
8.1119
8.1186
8.1651
8.1733
8.2060
7.8719
8.9562
8.9565
8.9528
8.9585
8.9606
8.9331
FE-Meshfree (QS,θ=10)
0.0823
0.0822
0.0822
0.0822
0.0822
0.0822
0.4937
0.4937
0.4936
0.4936
0.4936
0.4931
1.2826
1.2826
1.2825
1.2824
1.2822
1.2823
1.3010
1.3010
1.3008
1.3009
1.3010
1.2990
2.3655
2.3655
2.3655
2.3655
2.3661
2.3607
3.6094
3.6092
3.6098
3.6095
3.6113
3.6002
3.8450
3.8450
3.8445
3.8444
3.8441
3.8442
4.9716
4.9715
4.9722
4.9728
4.9755
4.9566
6.3984
6.3984
6.3978
6.3974
6.3966
6.3897
6.4116
6.4114
6.4126
6.4142
6.4179
6.3972
7.9019
7.9017
7.9043
7.9039
7.9110
7.8719
8.9347
8.9346
8.9337
8.9333
8.9331
8.9331Figure 18
Error in the computed first two frequencies for distortion sensitivity test.
(a)
(b)


## 4.3.1. Effect of the Correlation Parameterθ

The effects of the correlation parameterθ in the modified RPIM on the solution accuracy and stability property are investigated here. The parameter sensitivity analysis of the free vibration of a cantilever beam is conducted by employing 50 × 5 meshes (Figure 14), in which the correlation parameter θ is set to be ranged from 1 to 1,000. The first 12 frequencies of the problems determined using various θ, for comparison against the reference fundamental frequency solutions are tabulated in Table 3. Results show great agreement with the exact solutions and high stability regardless of the correlation parameter. Furthermore, nearly the same results are obtained for θ≥10. This unique feature of the presented radial basis function enables constant results for any correlation parameter. Accordingly, correlation parameter θ=10 is used for all subsequent computations.Figure 14
Domain discretization of beam using four-node element with regular elements.Table 3
First 12 natural frequencies (×104 Hz) for different correlation parameters θ.
θ
=
1
θ
=
5
θ
=
10
θ
=
50
θ
=
100
θ
=
500
θ
=
1000
f
1
0.0821
0.0822
0.0823
0.0823
0.0823
0.0823
0.0823
f
2
0.4937
0.4936
0.4937
0.4938
0.4938
0.4938
0.4938
f
3
1.2826
1.2826
1.2826
1.2826
1.2826
1.2826
1.2826
f
4
1.3038
1.3008
1.3010
1.3011
1.3012
1.3012
1.3012
f
5
2.3770
2.3653
2.3655
2.3656
2.3657
2.3657
2.3657
f
6
3.6374
3.6093
3.6094
3.6094
3.6094
3.6095
3.6095
f
7
3.8458
3.8450
3.8450
3.8450
3.8450
3.8450
3.8450
f
8
5.0252
4.9719
4.9716
4.9715
4.9715
4.9714
4.9714
f
9
6.4024
6.3985
6.3984
6.3983
6.3983
6.3983
6.3983
f
10
6.5005
6.4124
6.4116
6.4112
6.4111
6.4111
6.4111
f
11
8.0361
7.9033
7.9019
7.9009
7.9008
7.9008
7.9008
f
12
8.9455
8.9350
8.9347
8.9345
8.9345
8.9345
8.9345


## 4.3.2. Convergence Study

For the convergence rates, several types of regular meshes using FEM-QUAD4, FE-Meshfree (MQ), FE-Meshfree (Exp), are illustrated against the proposed method in this part. Four discrete mesh conditions with regular grids are considered as follows: 10 × 1, 20 × 2, 50 × 5, and 100 × 10. The first 12 natural frequencies computed for the four meshes are tabulated in Tables4–7. The reference results shown in the last column of Table 4 are gathered in the ABAQUS QUAD8 with 200 × 20 elements. Results show that the proposed method outperforms the other methods. In addition, the performance of FE-Meshfree (MQ) and FE-Meshfree (Exp) is also sensitive to the corresponding relation parameters.Table 4
Natural frequencies (×104 Hz) of a cantilever beam for 10 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(10 × 1) 10
22
0.1000
0.0824
0.0825
0.0909
0.1130
0.0848
0.0822
0.6077
0.5058
0.5050
0.5455
0.6732
0.5120
0.4931
1.2863
1.2846
1.2844
1.2849
1.2858
1.2839
1.2823
1.6422
1.3672
1.3621
1.4281
1.7140
1.3539
1.2990
3.0958
2.5583
2.5437
2.5729
2.9437
2.4719
2.3607
3.8918
3.8421
3.8430
3.8384
3.8596
3.8071
3.6002
4.9311
4.0208
3.9947
3.9079
4.2015
3.8405
3.8442
6.5947
5.6927
5.6693
5.4218
5.4614
5.3572
4.9566
7.1110
6.3692
6.3762
6.3422
6.4073
6.3656
6.3897
9.4517
7.5079
7.5313
7.1592
6.8653
7.1909
6.3972
9.5747
8.8600
8.8821
8.7702
8.6496
8.8518
7.8719
12.188
9.4823
9.6067
9.1641
9.8540
9.4100
8.9331Table 5
Natural frequencies (×104 Hz) of a cantilever beam for 40 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(20 × 2) 40
63
0.0872
0.0824
0.0825
0.0875
0.0981
0.0849
0.0822
0.5263
0.4981
0.4986
0.5249
0.5821
0.5112
0.4931
1.2837
1.2834
1.2834
1.2836
1.2839
1.2831
1.2823
1.4010
1.3233
1.3255
1.3826
1.5161
1.3523
1.2990
2.5816
2.4265
2.4325
2.5107
2.7232
2.4680
2.3607
3.8573
3.7297
3.7423
3.8232
3.8472
3.7767
3.6002
4.0002
3.8469
3.8472
3.8453
4.1064
3.8463
3.8442
5.6040
5.1688
5.1903
5.2518
5.5903
5.2123
4.9566
6.4492
6.3995
6.4018
6.3901
6.3954
6.3995
6.3897
7.3601
6.7004
6.7325
6.7524
7.1240
6.7308
6.3972
9.0702
8.2958
8.3397
8.2973
8.6741
8.3043
7.8719
9.2470
8.9314
8.9379
8.9054
8.9147
8.9338
8.9331Table 6
Natural frequencies (×104 Hz) of a cantilever beam for 250 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(50 × 5) 250
306
0.0830
0.0823
0.0822
0.0831
0.0852
0.0827
0.0822
0.4989
0.4937
0.4936
0.4986
0.5104
0.4962
0.4931
1.2827
1.2826
1.2826
1.2827
1.2828
1.2825
1.2823
1.3168
1.3010
1.3005
1.3127
1.3421
1.3063
1.2990
2.3992
2.3655
2.3640
2.3839
2.4343
2.3726
2.3607
3.6701
3.6094
3.6064
3.6331
3.7057
3.6165
3.6002
3.8467
3.8450
3.8449
3.8448
3.8449
3.8445
3.8442
5.0697
4.9716
4.9667
4.9984
5.0934
4.9765
4.9566
6.4062
6.3984
6.3983
6.3965
6.3955
6.3971
6.3897
6.5590
6.4116
6.4043
6.4391
6.5558
6.4121
6.3972
8.1119
7.9019
7.8918
7.9274
8.0652
7.8957
7.8719
8.9562
8.9347
8.9348
8.9290
8.9251
8.9320
8.9331Table 7
Natural frequencies (×104 Hz) of a cantilever beam for 1000 elements.
Number of elements
Number of nodes
FEM-QUAD4
FE-Meshfree (QS,θ=10)
FE-Meshfree (MQ,b=1×10−4,q=2.01)
FE-Meshfree (MQ,b=1.0,q=2.01)
FE-Meshfree (Exp,b=1.0)
FE-Meshfree (Exp,b=0.1)
Reference (ABAQUS Quad8)
(100 × 10) 1000
1111
0.0824
0.0822
0.0822
0.0824
0.0830
0.0824
0.0822
0.4947
0.4933
0.4933
0.4943
0.4978
0.4940
0.4931
1.2825
1.2824
1.2824
1.2824
1.2825
1.2824
1.2823
1.3037
1.2995
1.2994
1.3019
1.3110
1.3011
1.2990
2.3707
2.3618
2.3614
2.3654
2.3817
2.3641
2.3607
3.6183
3.6023
3.6015
3.6068
3.6312
3.6050
3.6002
3.8450
3.8445
3.8444
3.8445
3.8446
3.8443
3.8442
4.9858
4.9600
4.9586
4.9649
4.9981
4.9627
4.9566
6.3996
6.3944
6.3923
6.3972
6.3970
6.3968
6.3897
6.4332
6.3976
6.3975
6.3992
6.4415
6.3972
6.3972
7.9335
7.8783
7.8753
7.8822
7.9341
7.8797
7.8719
8.9391
8.9337
8.9337
8.9326
8.9313
8.9328
8.9331Figure15 shows an error plot of the first two natural frequencies obtained using the proposed method as well as FEM-QUAD4, FE-Meshfree (MQ), and FE-Meshfree (Exp). Results show that the corresponding error of the proposed method is generally lower than that of other methods. Even for the coarse mesh, the results of the proposed method still show better consistence with the reference solution. Figure 16 displays the first 12 eigenmodes obtained for FE-Meshfree (QS) using the 20 × 2 meshes. These modeshape plots are in accordance with previously published results for ES-FEM [21].Figure 15
Computed first two natural frequencies. (a) Mode 1. (b) Mode 2.
(a)
(b)Figure 16
First 12 modeshapes of cantilever beam.


## 4.3.3. Sensitivity to Distorted Meshes

The cantilever beam with irregular mesh shown in Figure13 is also employed to verify the effectiveness of the proposed method. As shown in Figure 17, four distorted meshes with 50 × 5 elements are used to assess the influence of mesh quality. Table 8 gives the comparison results of the computed natural frequencies for the first 12 natural frequencies. A comparison of relative error in the computed natural frequencies is shown in Figure 18. As αr increases, the accuracy of FEM-QUAD4 decreases, whereas the FE-Meshfree (QS) does not show markedly deviation from their original (high) accuracy, indicating great robustness of the proposed method to mesh distortion phenomena.Figure 17
Beam domain discretization using four-node element with irregular elements. (a)αr=0.1. (b) αr=0.2. (c) αr=0.3. (d) αr=0.4.
(a)
(b)
(c)
(d)Table 8
First 12 natural frequencies (×104 Hz) of a cantilever beam using irregular nodal distribution (50 × 5 elements).
Method
α
r
=
0
α
r
=
0.1
α
r
=
0.2
α
r
=
0.3
α
r
=
0.4
Reference (ABAQUS QUAD8)
FEM-QUAD4
0.0830
0.0831
0.0834
0.0835
0.0838
0.0822
0.4989
0.4992
0.5016
0.5019
0.5036
0.4931
1.2827
1.2827
1.2827
1.2827
1.2827
1.2823
1.3168
1.3176
1.3243
1.3245
1.3297
1.2990
2.3992
2.4007
2.4114
2.4139
2.4242
2.3607
3.6701
3.6725
3.6914
3.6963
3.7096
3.6002
3.8467
3.8468
3.8470
3.8470
3.8471
3.8442
5.0697
5.0728
5.0974
5.1059
5.1249
4.9566
6.4062
6.4063
6.4072
6.4072
6.4074
6.3897
6.5590
6.5643
6.5999
6.6089
6.6360
6.3972
8.1119
8.1186
8.1651
8.1733
8.2060
7.8719
8.9562
8.9565
8.9528
8.9585
8.9606
8.9331
FE-Meshfree (QS,θ=10)
0.0823
0.0822
0.0822
0.0822
0.0822
0.0822
0.4937
0.4937
0.4936
0.4936
0.4936
0.4931
1.2826
1.2826
1.2825
1.2824
1.2822
1.2823
1.3010
1.3010
1.3008
1.3009
1.3010
1.2990
2.3655
2.3655
2.3655
2.3655
2.3661
2.3607
3.6094
3.6092
3.6098
3.6095
3.6113
3.6002
3.8450
3.8450
3.8445
3.8444
3.8441
3.8442
4.9716
4.9715
4.9722
4.9728
4.9755
4.9566
6.3984
6.3984
6.3978
6.3974
6.3966
6.3897
6.4116
6.4114
6.4126
6.4142
6.4179
6.3972
7.9019
7.9017
7.9043
7.9039
7.9110
7.8719
8.9347
8.9346
8.9337
8.9333
8.9331
8.9331Figure 18
Error in the computed first two frequencies for distortion sensitivity test.
(a)
(b)


## 4.4. Example 4: Shear Wall Free Vibration Analysis

In this part, a shear wall with four openings, as shown in Figure19 is investigated. The parameters in the computation include Young’s modulus E=10000 N/m2, Poisson’s ratio v=0.2, thickness t=1.0 m, and mass density ρ=1.0 kg/m3. The bottom edge is fully clamped, and a plane stress case is considered here. The mesh for computation with 476 elements and 559 nodes is shown in Figure 20. This problem was analyzed previously using a boundary element method by Brabbia et al. [64] and also by Liu and Gu [13] using a local radial point interpolation method with the multiquadrics radial function. The commercial software applications, ANSYS and ABAQUS, are employed here for comparison.Figure 19
Shear wall with four openings.Figure 20
Model of shear wall with four openings.The natural frequencies of the first eight modes are tabulated in Table9. The results show great agreement with those obtained by other methods. Figure 21 displays the first eight eigenmodes obtained by the proposed method, from which it can be observed that the modeshapes plot match well with those of ES-FEM [21].Table 9
Natural frequencies (rad/s) of shear wall.
Model
FEM-QUAD4
FE-Meshfree (QS)
LRPIM (MQ) [13]
Brebbia (BEM) [64]
ANSYS [49]
ABAQUS
QUAD4-CNS [50]
1
2.074
2.085
2.086
2.079
2.057
2.073
2.094
2
7.114
7.117
7.152
7.181
7.067
7.096
7.124
3
7.627
7.632
7.647
7.644
7.620
7.625
7.590
4
11.978
12.063
12.019
11.833
11.840
11.938
12.192
5
15.444
15.410
15.628
15.947
15.313
15.341
15.451
6
18.429
18.351
18.548
18.644
18.342
18.345
18.330
7
19.991
19.882
20.085
20.268
19.887
19.876
19.858
8
22.348
22.258
22.564
22.765
22.236
22.210
22.225Figure 21
First 8 modeshapes of shear wall with four openings. (a) Mode 1. (b) Mode 2. (c) Mode 3. (d) Mode 4. (e) Mode 5. (f) Mode 6. (g) Mode 7. (h) Mode 8.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 4.5. Example 5: Connecting Rod Free Vibration Analysis

Here, we perform free vibration analysis of a connecting rod with a relatively complex configuration in which the boundary conditions can be clearly observed in Figure22. The rod is discretized with 723 irregularly distributed nodes, as shown in Figure 23. The plane stress problem is considered with E = 10 GPa, v = 0.3, and ρ = 7.8 × 103 kg/m3. A reference solution is also obtained in ABAQUS QUAD8 with 24,149 elements. Table 10 gives the first 10 natural frequencies, where FE-Meshfree (QS) yields comparable results to ABAQUS (similar to the previous example).Figure 22
Geometric model and boundary conditions of automotive connecting rod.Figure 23
Domain discretization of connecting bar using four-node elements.Table 10
First ten natural frequencies (×106 rad/s) of connecting rod.
Mode
FEM-QUAD4
FE-Meshfree (QS)
Reference (ABAQUS QUAD8)
1
1.3673
1.3592
1.3418
2
2.8093
2.7995
2.7042
3
3.9773
3.9289
3.7932
4
4.1908
4.1890
4.1078
5
5.7172
5.6162
5.5291
6
7.3185
7.3137
7.0938
7
7.8352
7.6765
7.5696
8
8.2621
8.2325
8.1509
9
8.4762
8.4387
8.3444
10
8.6747
8.6310
8.5351


## 4.6. Example 6: Cantilever Beam Implicit Forced Vibration Analysis

The benchmark cantilever beam and mesh shown in Figures24 and 25 are investigated via the Bathe method. The plane strain problem is considered with numerical parameters L=4.0 m, H=1.0 m, t=1.0 m, E=1.0 N/m2, v=0.3, ρ=1.0 kg/m3, α=0.005, and β=0.272. The domain is represented by 10 × 2 elements. The first loading is in the harmonic form with ft=cosωft. The corresponding dynamic responses can be shown in Figures 26 and 27 with ωf=0.01 and 0.05, respectively. Then, a constant step load ft=1 is added at the apex from t=0. The time step Δt=2.0 is used for time integration. Both damping effect and no damping effect are considered here. A reference solution is obtained in ABAQUS QUAD8 with 2,600 elements.Figure 24
Cantilever beam.Figure 25
Domain discretization of cantilever beam.Figure 26
Transient response of cantilever beam subjected to harmonic loading withωf=0.01: (a) without damping; (b) with damping.
(a)
(b)Figure 27
Transient response of cantilever beam subjected to harmonic loading withωf=0.05: (a) without damping; (b) with damping.
(a)
(b)Per the dynamic responses in Figures26–28, both methods tested provided nearly identical vibration frequency with the amplitude of FE-Meshfree (QS) being closer to the reference solution. This indicates that the FE-Meshfree (QS) can be readily applied to actual vibration analyses with excellent accuracy, which may be partially due to the fact that the FE-Meshfree (QS) model owns very close-to-exact stiffness.Figure 28
Transient response of cantilever beam subjected to constant step loading: (a) without damping; (b) with damping.
(a)
(b)


## 4.7. Example 7: Spherical Shell Implicit Forced Vibration Analysis

As shown in Figure29, a spherical shell subjected to a concentrated loading at its apex is examined. In total, 30 × 3 asymmetric elements are used (Figure 30) with the numerical parameters R=12 m, t=0.1 m, ϕ=10.9°, E=1.0 N/m2, v=0.3, and ρ=1.0 kg/m3. A reference solution is obtained by ABAQUS QUAD8 with 2,000 elements. The first loading is applied in the harmonic form ft=cosωft. Its dynamic responses can be observed from Figures 31 and 32 with ωf=0.01 and 0.05, respectively. A constant step load ft=1 is then added at the apex from t=0. The time step Δt=2.0 is used for time integration. The dynamic responses shown in Figures 31–33 indicate that FE-Meshfree (QS) is more accurate than FEM-QUAD4. It yields results very close to the reference one for both damping and no damping effect situations.Figure 29
Spherical shell.Figure 30
Domain discretization of half of the spherical shell.Figure 31
Transient response of spherical shell subjected to harmonic loading withωf=0.01: (a) without damping; (b) with damping.
(a)
(b)Figure 32
Transient response of spherical shell subjected to harmonic loading withωf=0.05: (a) without damping; (b) with damping.
(a)
(b)Figure 33
Transient response of spherical shell subjected to constant step loading: (a) without damping; (b) with damping.
(a)
(b)


## 5. Conclusions

In this study, FE-Meshfree QUAD4 element with modified radial point interpolation function has been extended to the free vibration and forced vibration analysis of 2D solids. Based on the above simulations, we can make the following conclusions:(1)
Compared to meshfree methods, the composite shape function of FE-Meshfree QUAD4 element with modified radial point interpolation function possess the much desired Kronecker delta property. Therefore, the essential boundary conditions can be easily applied as in the FEM.(2)
As in the classical FEM, the shape function of FE-Meshfree possesses the much desired Kronecker delta property, and accordingly the essential boundary conditions can be easily applied. Moreover, higher-order global approximations can be readily achieved without the addition of extra nodes or DOFs.(3)
The normalized radial basis function yields composite shape functions without any tuning parameters and thus can be used to approximate displacement fields.(4)
Compared to the FEM-QUAD4, the mesh distortion test (Section4.3.3) shows that the error of the proposed method is insensitive to the increase of the distortion parameters, indicating better robustness to mesh distortion of the proposed method.(5)
The proposed method enables to obtain more accurate free vibration analysis results at a higher convergence rate. In addition, the vibration period obtained for forced vibration analysis is very close to the reference ones.Although the presented method works well in 2D linear elastic, free vibration, and forced vibration problems, most problems are 3D (solid, plate, shell, etc.) in reality. Formulations and numerical tests of these elements will be presented and applied for more complex problems in our following work.



---

*Source: 3269276-2019-01-08.xml*