Modeling Structural Dynamics Using FE-Meshfree QUAD4 Element with Radial-Polynomial Basis Functions

The PU (partition-of-unity) based FE-RPIM QUAD4 (4-node quadrilateral) element was proposed for statics problems. In this element, hybrid shape functions are constructed through multiplying QUAD4 shape function with radial point interpolation method (RPIM). In the present work, the FE-RPIM QUAD4 element is further applied for structural dynamics. Numerical examples regarding to free and forced vibration analyses are presented. The numerical results show that: (1) If CMM (consistent mass matrix) is employed, the FE-RPIM QUAD4 element has better performance than QUAD4 element under both regular and distorted meshes; (2) The DLMM (diagonally lumped mass matrix) can supersede the CMM in the context of the FE-RPIM QUAD4 element even for the scheme of implicit time integration.


Introduction
The FEM (finite element method) [1] has been widely adopted to model structural dynamics. The problem domain to be considered is discretized by a series of elements of simple shapes. In a two-dimensional problem domain, finite elements are either triangles or quadrilaterals [2]. TRIG3 (3-node triangular) element is much easier than QUAD4 element in mesh generation but has poorer accuracy. The QUAD4 element can generally obtain satisfactory results for regular mesh but performs badly for distorted mesh. Often in practice, the generation of a high-quality quadrilateral mesh is a time-consuming task for problems with complex geometric boundary.
A lot of work has been done in the past several decades to develop more powerful numerical methods than the FEM for the analysis of time-dependent problems. Due to the use of higher-order interpolation with specific quadrature formulae [3,4], the spectral finite element method (SFEM) is capable of simulating dynamic problems, such as wave propagation, and has better convergent rate than FEM [5]. The meshless methods, was also expected to become an effective procedure since the mesh is not needed [6][7][8][9][10][11][12][13][14][15][16]. However, some meshless methods' shape functions have no Kronecker-delta character, which means special treatment is required for essential boundary condition implementation. Moreover, the computational cost to form the trial functions of meshless methods cannot be neglected [17]. Hence, several schemes were proposed to optimize meshless methods [18,19].
In other front, a class of PU [20] based methods were proposed, such as the PUFEM (PU FEM) [21], the GFEM (generalized FEM) [22] and the NMM (numerical manifold method) . In these methods, global approximations with high order are usually built by using high order local approximations. However, the LD (linear dependence) problem [20] will arise if the local approximations and the PU functions were simultaneously constructed using explicit polynomials [21]. For the purpose of eliminating the LD problem, a lot of efforts have been made [46][47][48].
Recently, a series of FE-Meshfree elements including FE-LSPIM QUAD4 element [46], were proposed. In the FE-LSPIM QUAD4 element, shape functions are usually built through multiplying the FE shape function (QUAD4) and the meshfree shape functions (LSPIM). As a result, good properties of meshless method and FEM are inherited [46]. Compared to FEM, higher accuracy can be obtained by FE-LSPIM QUAD4 element, because global approximations with high order are constructed. In contrast to meshless method, FE-LSPIM QUAD4 shape functions have Kronecker-delta character, which means special treatment is not needed for essential boundary condition implementation. Moreover, the FE-LSPIM QUAD4 element is immune to LD problem. However, since a pure polynomial basis is adopted in the LSPIM, the moment matrix singularity may arise if the polynomial basis functions were inappropriately employed [14,49].
To avoid the drawback of LSPIM, the RPIM (radial point interpolation method) [14] has been employed to replace it, resulting in a new FE-RPIM QUAD4 element [49]. According to the report from Xu and Rajendran [49], FE-RPIM QUAD4 element has higher accuracy than FE-LSPIM QUAD4 element for linear and geometry nonlinear static problems if the same number of polynomial terms are employed. Apart from FE-LSPIM QUAD4 and FE-RPIM QUAD4 elements, there are also other types of FE-Meshfree elements [50][51][52][53].
The FE-RPIM QUAD4 element is further applied for structural dynamics in the present work. After the statement of the formulations related to FE-RPIM QUAD4 element in Section 2, equations related to elastodynamic problems are presented in Section 3. Besides, the expressions for the stiffness matrix, the CMM and the DLMM are given. In Section 4, five numerical examples are investigated using the FE-RPIM QUAD4 element. At the last section, some conclusions will be presented.

Formulation of Shape Functions
Before expounding formulation of the shape functions, two important concepts, namely, the nodal support domain (Ω i ) and the element support domain (Ω e ), are introduced. The nodal support domain is adopted to determine all the support nodes of a given node (also named as central node) to construct its local approximation. According to the scope of the support domain, we can define different order of nodal support domain (Ω i ). The nodal support domain of first order is determined through the nodal connectivity of first order. Figure 1 shows an example for the nodal support domain of first order for node 1, where 9 support nodes are obtained in Figure 1a and Figure 4 support nodes are determined in Figure 1b. Similarly, the nodal support domain of second order or third order can be defined. As shown in Figure 2,Ω is the union of 4 nodal support domains (Ω i ), which is defined asΩ  Let domain Ω e be defined with nodes {P 1 P 2 P 3 P 4 }. The global approximation u h (x) of FE-RPIM QUAD4 element can be simply written as [49] where w i (x) represents PU function, while u i (x) represents local approximation function related to node i. Since quadrilateral mesh is employed, w i (x) is equivalent to QUAD4 shape function. It is noticed that the global approximation of QUAD4 element can be considered as a special case of FE-RPIM QUAD4 element. If u i (x) is set as a constant, then u h (x) in Equation (1) becomes the global approximation for QUAD4 element. The local approximation functions are obtained by the RPIM [14], which is defined with: where M and n [i] represent the specified polynomial term number and the number of nodes within Ω i . a and b are two unknown vectors. r represents the radial basis functions, while p represents the polynomial functions formulated as: p(x) = {1 x y xy x 2 y 2 } for M = 6.
According to the discussion in [49], the FE-RPIM QUAD4 element with three-or four-term basis can give as accurate results as with six-term basis for static problems. Furthermore, the computational cost of six-term basis or four-term basis is more than that of three-term basis. Hence, in the present paper, we set M = 3, namely, a three-term basis is used; r(x, y) can be formulated with [14,16]: in which The numerical values of q and c have an influence on the accuracy of RPIM, which have been discussed in great details in [49]. In the present work, q and c are separately set to 2.01 and 0.0001 according to the recommendation in [49]. Enforcing Equation (2) to pass through all the nodes in within Ω i , the following equations are obtained: where u i is a vector of corresponding nodal displacements of all the nodes in Ω i , and Obviously, there are totally (M + n [i] ) parameters in Equation (8). However, only n [i] equations are available. Nevertheless, according to the work finished in [14,54], vectors of a and b can be eliminated as proposed in their work and local approximation function in Equation (2) is eventually expressed as [14]: Note that Wendland [55] has presented the evident of the existence of R −1 for any scattered nodes.

Properties of Shape Functions
Through substituting Equation (11) into Equation (1), u h (x) for the FE-RPIM QUAD4 element can then be formulated with a more compact form: where φ i (x) and d i represent a shape function and nodal displacement associated with node i. p represents node number withinΩ e .
The attractive properties of φ i (x) in Equation (19) can be listed as follows [49]: (ii) Compatibility property at the interface of elements.
(iii) High order completeness, in other words, reproducibility of all the assumed Cartesian terms (Equation (3)).
To investigate the changes of shape functions within element, we consider a twodimensional domain discretized with four quadrilateral elements and 9 nodes ( Figure 3). The 3D graphics of shape functions of FE-RPIM QUAD4 element for different nodes (node 1, node 2, node 3 and node 5) are plotted in Figure 4. For the purpose of com-parison, the 3D graphics of shape functions of QUAD4 element are also plotted, as shown in Figure 5. As can be seen in Figures 4 and 5, the shape function of FE-RPIM QUAD4 element is smoother than that of Quad4 element.

FE-RPIM QUAD4 for Dynamic Analysis
Let the problem domain Ω discretized using a series of QUAD4 elements: In the context of the FE-RPIM QUAD4 element, the discretized form of the system equations about dynamic analysis can be expressed as [56]: (21) in which M is the global mass matrix, K is the global stiffness matrix, and f is the global load vector, which can be computed as: where and in which ρ is the material density, D is the matrix of the elastic constants of the material, t is the specified traction vector applied on stress boundary Γ e σ , b is the body force per unit volume, N represents the matrix of shape function for Ω e , expressed as: and B is the gradient matrix expressed as: where φ e i represent a shape function associated with node i. p represents the node number withinΩ e .
The Rayleigh damping is adopted in the present work. Hence, the damping matrix C can be formulated into, in which β 1 , β 2 represent the coefficients of Rayleigh damping.

Time Integration Scheme
The Newmark method [56] will be employed in this study to solve Equation (21). Assuming that d n , v n and a n separately represent approximation values of d(t n ), then the approximations of d(t n+1 ) and . d(t n+1 ) are expressed as: Here, a n+1 represents the solution of Equation (32). with a 0 is computed with Equation (35).
It is noticed that M obtained from Equations (22) and (24) is named as CMM (consistent mass matrix). If M is a diagonally lumped mass matrix (DLMM), solving a 0 should be a very easy task. Moreover, if the damping effect is ignored, as is done frequently and β = 0, then M reduces to M. Consequently, if M is a positive DLMM whose inverse M −1 is easy to calculate, solving Equation (32) for a n+1 would be very easy.

Generalized Eigenvalue Problem
If the terms associated to damping matrix and force vector are neglected, Equation (21) will reduce to the following form [14]: The general solution for Equation (36) is formulated as: (37) in which t represents time. ω represents the natural frequency, while d represents the eigenvector. Based on Equations (36) and (37), ω is obtained from Equation (38).
Here, d determines the mode shape related to ω. Note that it is time-consuming to compute all the modes. Hence, the subspace iteration procedure is usually adopted to obtain only those modes with small values of ω.
However, if M is a DLMM, Equation (38) becomes a standard eigenvalue problem, where K is defined as: Since M is a DLMM, the inversion of M, namely, M −1 , can be easily obtained. There are more powerful algorithms available to solve the standard eigenvalue problem (Equation (39)) [57].

Diagonally Lumped Mass Matrix
Here, we propose to use the "special lumping technique" introduced by Hinton et al. [58] to obtain the diagonally lumped mass matrix (DLMM). This procedure always leads to positive lumped masses at the nodes. Moreover, the requirement for mass conservation can be ensured. With little change, this procedure can be used for FE-RPIM QUAD4 element. For the FE-RPIM QUAD4 element, we obtain: where M ce ii is the i-th diagonal entry of the element consistent mass matrix (CMM) obtained through Equation (24), M le ii is the i-th diagonal entry of the DLMM, α is a constant defined as: where ρ represents the density of material.

Numerical Examples
In this section, a static example and five examples regarding to free and forced vibration analyses will be presented to investigate the performance of FE-RPIM QUAD4 element. For comparison, other well-known element types, such as the TRIG3, QUAD4 and QUAD8 (8-node quadrilateral element) will also be used to conduct these tests. In the computation for dynamic problems, both the CMM and the DLMM will be employed in the context of FE-RPIM QUAD4 element, while only the CMM will be employed in the context of QUAD4, TRIG3, and QUAD8.
Physical units in this section will be according to the unit system of international standard without specification. n is the node number within discretized model. Natural frequency's relative error expressed in Equation (44) is adopted to evaluate the accuracy of the numerical models.
in which "num" and "ref" separately represent the numerical solution and reference solution.

Cook's Skew Beam
Before proceeding with dynamic analysis, the well-known example of Cook's skew beam (Figure 6a) is adopted to test the performance of the FE-RPIM QUAD4 element for static problem. A mesh with 10 × 10 quadrilateral elements and 121 nodes is plotted in Figure 6b. The reference solution for the vertical displacement of point A for Cook's skew beam is 23.96 [49]. With the mesh presented in Figure 6b, the solutions for the vertical displacement of point A based on the FE-RPIM QUAD4 element is 23.8170, which is better than that based on the QUAD4 element (22.6965). Note that this problem was also solved in [59] by the spectral finite element (SFEM) [60]. If the problem domain is discretized into 10 × 10 elements, where each element has 2 × 2 nodes, the solution based on the SFEM is 22.6940, which is inferior to that based on the FE-RPIM QUAD4 element (23.8170). However, If the problem domain is discretized into 1 element with 11 × 11 nodes, the solution based on the SFEM is 23.9414, which is even slightly better than that based on the FE-RPIM QUAD4 element (23.8170).

A Slender Rod
Free vibration analysis for a slender rod [61] is conducted, as shown in Figure 7. The ith natural frequency from the analytical solution is [62]: in which L represents the length of the rod. Geometric and mechanical parameters used in this example are L = 100, thickness t = 1, D = 1, Poisson's ratio v = 0, Young's modulus E = 72 × 10 9 , mass density ρ = 2700. Since L/D = 100 > 10, the geometric parameters used assure that the one-dimensional slender rod can be well represented by the two-dimensional model without causing unacceptable error. In the computation, both upside and down sides of the model are fixed in the normal direction.
For investigating the convergence of the solution, three regular discretized models (Figure 8) are constructed. Table 1 lists the numerical results from different numerical models. According to Table 1, the results from all the listed numerical models approach the analytical solution gradually, as the mesh density increases.  By using the consistent mass scheme, the result from FE-RPIM QUAD4 element is better than that from QUAD4 element and TRIG3 element, but slightly inferior to that from QUAD8 element. Note that QUAD8 element requires more nodes than FE-RPIM QUAD4 element to discretize the problem domain. In addition, results from the CMM and the DLMM in the context of FE-RPIM QUAD4 element are close to each other.

An Annulus
As the second example, an annulus without constraint shown in Figure 9 is employed to validate the FE-RPIM QUAD4 element. The geometric and mechanical parameters used in this example are R a = 0.4, R b = 0.5, v = 0.33, E = 72 × 10 9 , t = 1, ρ = 2700.
For the purpose of investigating the convergence of the solution, several regular discretized models shown in Figures 10 and 11 are constructed. Table 2 lists the frequencies assessed by different numerical models. The reference solution for this problem, which is listed in the last column of Table 2, can be found in [61]. According to Table 2, the result from FE-RPIM QUAD4 element is much better than that from QUAD4 and TRIG3 elements, but slightly inferior to that from QUAD8 element, if the CMM was employed.   In addition, FE-RPIM QUAD4 element can obtain better result if the DLMM is employed instead of the CMM.

Mesh Distortion Test
In the previous two examples, regular meshes are adopted. In this section, mesh distortion test based on the cantilever beam ( Figure 12) is conducted. The parameters used are D = 10 mm, L = 100 mm, v = 0.3, E = 2.1 × 10 4 kg/mm 2 , ρ = 8.0 × 10 −10 kg fs 2 /mm 4 , t = 1 mm. As shown in Figure 13a,b, the beam is meshed with two elements for QUAD4 element, FE-RPIM QUAD4 element and QUAD8 element. For TRIG3 element, the problem domain is meshed with four elements, as shown in Figure 13c. A distortion parameter, 2d/D, is adopted to control mesh distortion.  Table 3 lists the fundamental natural frequency assessed by different numerical models. The corresponding relative errors are plotted in Figure 14. To better see the reached accuracy of FE-RPIM QUAD4(CMM) and FE-RPIM QUAD4(DLMM), Figure 15 is plotted. According to Table 3, Figures 14 and 15, we can draw the following conclusions: (1) First, as distortion parameter's value increases, the errors based on FE-RPIM QUAD4 element do not change appreciably, while those based on QUAD4 element, TRIG3 element and QUAD8 elements change rapidly. The FE-RPIM QUAD4 element is immune to mesh distortion. (2) Second, accuracy of FE-RPIM QUAD4 element is always much higher than QUAD4 and TRIG3 elements. (3) Third, when 2d/D < 0.2, QUAD8 element's accuracy is higher than QUAD4, FE-RPIM QUAD4 and TRIG3 elements. However, as the value of 2d/D increases, accuracy through QUAD8 element deteriorates quickly. If meshes used are distorted severely, QUAD8 element's accuracy is much lower than FE-RPIM QUAD4 element. (4) Fourth, compared to CMM, FE-RPIM QUAD4 element can achieve better results if DLMM is employed.

A Plate with Four Holes
In this section, a plate with four holes shown in Figure 16 is considered. The mechanical parameters used are v = 0.3, E = 72 × 10 9 and ρ = 2700. Left side of the plate is fixed. For the purpose of investigating the convergence of the solution, four discretized models shown in Figures 17 and 18 are constructed. In order to get a reference solution, the QUAD4 element is adopted with a very dense discretized model (43,704 nodes, 43,058 elements). Table 4 lists the results from different numerical models. The results from different numerical models all approach the reference solution gradually, as the mesh density increases. In addition, by using the CMM, FE-RPIM QUAD4 element's accuracy is much higher than QUAD4 and TRIG3 elements.    The first six mode shapes obtained through the FE-RPIM QUAD4 element under CMM are presented in Figure 19. Mode shapes through CMM agree well with those through DLMM in the context of FE-RPIM QUAD4 element.

A Cantilever Beam under Harmonic Load
In the context of FE-RPIM QUAD4 element, accuracy for natural frequencies obtained through the DLMM is as good as or even better than those through the CMM. This conclusion was also hold FEM [61,63]. However, the CMM is considered to be a better choice to calculate deformation and mode shapes [61,63]. In this section, following the same purpose as in [61], an example is employed to show that the FE-RPIM QUAD4 element is able to obtain satisfactory results by employing DLMM even for the scheme of implicit time integration.
A cantilever beam ( Figure 20) under harmonic load ( f (t)) on the right end is considered. The parameters used are ρ = 1 kg/m 3   In the computation, scheme of implicit time integration with γ = 1 and β = 0.5 is used. Time step size ∆t = 1.57s is adopted. Shown in Figure 21 is the discretized model used by the FE-RPIM QUAD4 element. Apart from the FE-RPIM QUAD4 element, this dynamic problem is also investigated using QUAD4 element (27 nodes, 16 elements), TRIG3 element (27 nodes, 32 elements) and QUAD8 element (69 nodes, 16 elements). In order to compute a reference solution, the QUAD4 element is adopted with a very dense discretized model (6601 nodes, 6400 elements). As can be seen from Figure 22, if CMM scheme is employed, FE-RPIM QUAD4 element's accuracy is higher than QUAD4 and TRIG3 elements. More importantly, the result of FE-RPIM QUAD4 (DLMM) almost coincides with FE-RPIM QUAD4 (CMM), which means DLMM can supersede the CMM in the context of the FE-RPIM QUAD4 element.

Conclusions
FE-RPIM QUAD4 element is further extended for problems of structural dynamics. Some important conclusions, which can be drawn from this work, are as follows: (1) Based on 4-node quadrilateral mesh, FE-RPIM QUAD4 element's accuracy is much higher than QUAD4 and TRIG3 elements ( Table 2). (2) Although FE-RPIM QUAD4 element's accuracy is slightly inferior to QUAD8 element, QUAD8 element requires more nodes than FE-RPIM QUAD4 element to discretize the problem domain. In addition, FE-RPIM QUAD4 element can achieve results closing to the reference solution, even for coarse mesh ( Figure 22). (3) For distorted meshes, FE-RPIM QUAD4 element's accuracy is always much higher than QUAD4 and TRIG3 elements. Moreover, FE-RPIM QUAD4 element is immune to mesh distortion, but TRIG3, QUAD4 and QUAD8 elements give very bad results as the mesh quality deteriorates ( Figure 14). (4) In the tests associated to the analysis of free vibration, the result based on the DLMM are very close to those based on the CMM in the context of FE-RPIM QUAD4 element.
In the test on forced vibration analysis, the result from the DLMM also agrees well with that from the CMM, which means DLMM can supersede the CMM in the context of the FE-RPIM QUAD4 element even for the scheme of implicit time integration. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.