FE-Meshfree QUAD 4 Element with Modified Radial Point Interpolation Function for Structural Dynamic Analysis

+e 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. +is method yields high accuracy and convergence rate without necessitating extra nodes or DOFs. In this study, the FE-Meshfreemethod 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. +e 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.


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][2][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.ey are suitable for solving complex practical problems such as large deformation, fracture propagation simulation, and impact-induced failure [6][7][8][9][10][11][12][13][14][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][19][20][21][22][23][24][25][26][27].
Another common and extensively studied hybrid method is the partition-of-unity method (PUM) [28,29].
e 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][41][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.e 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 twodimensional solids.e rest of this paper is organized as follows.Section 2 briefly reviews the formulation of the composite shape function and the modified radial point interpolation function.is 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.

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 e j can be written as follows: where N i (x)(i � 1, 2, . . ., 4) are isoparametric shape functions used in the classic QUAD4 element and u i (x) 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 p z (x) and radial basis function r j (x), which can be expressed as follows: In the above, m is the number of polynomial terms and n denotes the number of nodes in the node support domain Ω i .
e vectors a and b in equation (2) are vectors yet to be determined.p(x) and r(x) denote the corresponding polynomial and radial basis functions, respectively.If we consider three respective terms for 2D problems, p(x) can be written as follows: e radial basis function vector r(x) is expressed as where r k (d k ) is the function of distance d k , which denotes the space distance between the selected interpolation point x and a node x k .Also for the 2D problem, d k can be expressed as Enforcing equation (2) to pass through all the nodes in the nodal support domain yields the following equation: where u i 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: 2 Shock and Vibration e related matrices R and P are R e following constraint condition can be imposed to guarantee a unique displacement approximation as follows [54]: From equation ( 6), we can obtain: Substituting equation (10) into equation ( 9), the vector b can be solved: Substituting equation (11) into equation (10), the vector a can be obtained: Eventually the nodal displacement function can be expressed by substituting equations (11) and (13) into equation (2): where e 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 de ned by the element including nodes in element j.As shown in Figure 2, the element support domain 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.
Substituting equations (10) and (11) into equation (1) allows us to rewrite the element displacement eld as follows: In the above, u L×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 Similarly, the element displacement eld in y-direction is written as follows: where v L×1 denotes the y-direction displacements of all nodes in Ω e .
e Cartesian derivatives of shape function are obtained by e composite shape function ψ is characterized by the Kronecker delta property, interelement compatibility property, and higher-order completeness properties, i.e., the

Shock and Vibration
reproducibility of all the Cartesian terms appearing under the assumed basis of equation (1).

Modi ed Radial Based Function.
During the numerical analysis process, many speci c radial basis functions are used in radial point interpolation method (RPIM) to ensure favorable performances.
e typical radial basis functions, usually employed in engineering problems, are tabulated in Table 1.However, the tuning parameters show signi cant in uences 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.
In this study, we used the modi ed basis function proposed by Do [55] to construct the composite shape functions without any supporting tuning parameters.
e radial basis function can be taken as a compactly supported form: In the above, s(d k ) denotes a compactly supported function, M(d k ) represents a radial basis function, and d s is the size of the support domain.Generally, s(d k ) can be chosen as a Heaviside step function, if However, as the nodal support domain and the corresponding nodes have been de ned, s(d k ) does not need to be used, and equation ( 21) can be rewritten as e following radial basis function in a compactly supported form based on the original kernel function of the quartic spline (QS) [6] is applied: e 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.
e 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 de ned as follows: where κ is the scale factor and d c 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.d c is determined as follows for the irregular nodal distribution [56]: where A s denotes the area of the nodal support domain and n A s is the corresponding number of nodes in the nodal support domain.Substituting equation (24) into equation (23), equation ( 23) can be rewritten as From equation ( 26), θ/κ can be regarded as one variable; therefore, κ 1 can be adopted in the presented formulation.

Name
Expression Parameters Gaussian (exp)  Shock and Vibration

Structural Dynamic Analysis
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 where T is the nodal displacement vector and B represents the standard displacement gradient matrix.For 2D problems, it can be determined as follows: e stress can be obtained by the following constitutive relation: where D is the matrix of the elastic constants of the material.Substituting equations ( 28) and ( 29)into equation ( 27) yields According to the principle of virtual work, the discrete governing equation can be expressed as or in which

Free Vibration Analysis.
e 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]: A general solution to this system is 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 ω: e natural frequencies and corresponding mode shapes of a given structure are often referred to as the dynamic characteristics of the structure.

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 where α and β are the Rayleigh damping coefficients.ere 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][59][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: Substituting equation (40) into equation ( 33) yields the following algebraic equation:

Shock and Vibration
5 In the second substep, the following equations are expressed by the three-point Euler backward method: Similarly, the control equation ( 34) at time t + h can also be expressed as 9 Once d t+h/2 and d t+h is solved by equations ( 41) and ( 43

Numerical Examples
To evaluate the accuracy of the method, the following displacement error norm and the strain energy norm are defined: where superscript "exact" denotes the exact value of the problem, "num" denotes result solved by numerical solutions, and N d is the number of total field nodes.e strain energy of the numerical solution E num and the total strain energy of the exact solution E exact are defined as follows: where ε i is the strain of the exact solution.In the actual computation on equation ( 45), we used a very fine mesh (N e ⟶ ∞) to calculate the "exact" strain energy E exact .
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: 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.e 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 where p(x) 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 For the distributed nodes in the patch domain, a pointwise error vector can be evaluated as where U ex n and U n denote the numerical solution and the corresponding exact one, respectively.
Summing the resulted errors for each monomial, we can then obtain the error e as e Euclidian error can be evaluated as where Q represents a positive semidefinite matrix, and the maximum eigenvalue is derived as 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 taking m � 1 in equation ( 47) leads to 6 Shock and Vibration We choose a [−1,1] × [−1,1] square patch domain using 10 × 10 regular nodes and 100 irregular nodes as shown in Figure 3.For irregular interior nodes, the corresponding coordinates can be obtained as where Δx and Δy are the initial regular meshes' sizes along x and y axes, respectively.r c 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.

Example 2: Cantilever Beam Subjected to Tip-Shear Force.
A cantilever beam with length L 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. e beam is assumed to have a unit thickness satisfying the plane stress condition.e analytical displacement solution is described in detail by Timoshenko and Goodier [63]: where the moment of inertia I of the beam is given by e corresponding stresses are e parameters are E � 3.0 × 10 7 kPa, v � 0.3, D � 12 m, L � 48 m, and P � −1000 N.
During the computational process, the nodes on the left boundary are constrained using the exact displacements obtained from equations ( 55) and ( 56).
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.Figure 5 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.e convergence of the strain energy is shown in Figure 7. e 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.
e FE-Meshfree (QS) shows a very close-to-exact stiffness and hence very accurate results.e strain energy, displacement, and energy error norms of the FE-Meshfree (QS) are all better than those of the FEM-QUAD4.
Figure 10 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.

Example 3: Cantilever Beam Free Vibration Analysis.
As shown in Figure 13, 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.e plane stress condition is     Shock and Vibration considered.e relative geometrical and material parameters can be summarized as follows: thickness t 1.0 m, Young's modulus E 2.1 × 10 4 N/m 2 , Poisson's ratio v 0.3, and mass density ρ 8.0 × 10 −10 kg/m 3 .Based on the Euler-Bernoulli beam theory, the fundamental frequency f 1 0.08276 × 10 4 Hz can be obtained as a reference.

E ect of the Correlation Parameter θ.
e e ects of the correlation parameter θ in the modi ed RPIM on the solution accuracy and stability property are investigated here.
e 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.e rst 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.
is 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.

Convergence Study.
For the convergence rates, several types of regular meshes using FEM-QUAD4, FE-Meshfree

Shock and Vibration
(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. e rst 12 natural frequencies computed for the four meshes are tabulated in Tables 4-7.e 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.
Figure 15 shows an error plot of the rst 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 rst 12 eigenmodes obtained for FE-Meshfree (QS) using the 20 × 2 meshes.ese modeshape plots are in accordance with previously published results for ES-FEM [21].

Sensitivity to Distorted Meshes.
e cantilever beam with irregular mesh shown in Figure 13 is also employed to verify the e ectiveness of the proposed method.As shown in Figure 17, four distorted meshes with 50 × 5 elements are used to assess the in uence of mesh quality.Table 8 gives the comparison results of the computed natural frequencies for the rst 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.

Example 4: Shear Wall Free Vibration Analysis.
In this part, a shear wall with four openings, as shown in Figure 19 is investigated.e parameters in the computation include Young's modulus E 10000 N/m 2 , Poisson's ratio v 0.2, thickness t 1.0 m, and mass density ρ 1.0 kg/m 3 .e bottom edge is fully clamped, and a plane stress case is considered here.e mesh for computation with 476 elements and 559 nodes is shown in Figure 20. is 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.
e commercial software applications, ANSYS and ABAQUS, are employed here for comparison.
e natural frequencies of the rst eight modes are tabulated in Table 9. e results show great agreement with those obtained by other methods.Figure 21 displays the rst 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].

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 Figure 22. e rod is discretized with 723 irregularly distributed nodes, as shown in Figure 23.e plane stress problem is considered with E � 10 GPa, v � 0.3, and ρ � 7.8 × 10 3 kg/m 3 .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).e corresponding dynamic responses can be shown in Figures 26 and 27 with ω f � 0.01 and 0.05, respectively.en, a constant step load f(t) � 1 is added at the apex from t � 0. e 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.
Per the dynamic responses in Figures 26-28, both methods tested provided nearly identical vibration frequency

Shock and Vibration
with the amplitude of FE-Meshfree (QS) being closer to the reference solution.is 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 sti ness.

Example 7: Spherical Shell Implicit Forced Vibration
Analysis.As shown in Figure 29, 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/m 2 , v 0.

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.erefore, 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) e 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 (Section 4.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) e proposed method enables to obtain more accurate free vibration analysis results at a higher convergence rate.In addition, the vibration period

Figure 14 :
Figure 14: Domain discretization of beam using four-node element with regular elements.

Figure 18 :Figure 19 :
Figure 18: Error in the computed rst two frequencies for distortion sensitivity test.

Figure 20 :Figure 21
Figure 20: Model of shear wall with four openings.

Figure 30 :Figure 31 :Figure 32 :
Figure 30: Domain discretization of half of the spherical shell.

Figure 33 :
Figure 33: Transient response of spherical shell subjected to constant step loading: (a) without damping; (b) with damping.

Table 2 :
Displacement error norm of numerical results for standard patch test on different meshes.