A Modiﬁed Radial Point Interpolation Method (M-RPIM) for Free Vibration Analysis of Two-Dimensional Solids

: The classical radial point interpolation method (RPIM) is a powerful meshfree numerical technique for engineering computation. In the original RPIM, the moving support domain for the quadrature point is usually employed for the ﬁeld function approximation, but the local supports of the nodal shape functions are always not in alignment with the integration cells constructed for numerical integration. This misalignment can result in additional numerical integration error and lead to a loss in computation accuracy. In this work, a modiﬁed RPIM (M-RPIM) is proposed to address this issue. In the present M-RPIM, the misalignment between the constructed integration cells and the nodal shape function supports is successfully overcome by using a ﬁxed support domain that can be easily constructed by the geometrical center of the integration cell. Several numerical examples of free vibration analysis are conducted to evaluate the abilities of the present M-RPIM and it is found that the computation accuracy of the original RPIM can be markedly improved by the present M-RPIM.


Introduction
The classical finite element method (FEM), which is based on the weighted residual technique, is a versatile and well-developed numerical approach in the field of modern computation mechanics [1]. Many mature commercial software packages (such as ANSYS, ABAQUS and NASTRAN) based on the FE approach have been developed and used in various engineering applications. Though the standard FEM has achieved great success in practical engineering computation, the FEM still suffers from several inherent shortcomings compared to other advanced numerical techniques [2][3][4][5][6][7][8][9][10][11][12]. Among them, one important issue is that the FEM is essentially a mesh-based method and the involved problem domain should be firstly discretized into a series of elements that are connected by nodes for FE analysis. Therefore, the additional burdensome tasks for meshing operations cannot always be circumvented. Additionally, the solution accuracy of the FEM is usually sensitive to mesh qualities and the solutions from low-quality meshes are always not sufficiently accurate. To obtain sufficiently fine solutions, more attention should be given to obtain high-quality meshes. These issues will be greater when the standard FEM is employed to manage problems related to dynamic cracks and large deformation of complicated geometric shapes.
To alleviate the dependence of the conventional FE approach on predefined meshes, a series of smoothed FEMs [13][14][15][16][17] and meshfree techniques [18][19][20][21][22][23][24] have been proposed, function approximation u h (x) can be expressed in the following form by using the radial basis function (RBF) and polynomial basis function (PBF) [18]: (1) in which R i (x) stands for the RBF used and P j (x) represents the PBF used; n denotes the number of RBF used for interpolation, namely, there are n field nodes in the support domain of the sampling point x, m denotes the number of PBF used for interpolation, and the complete linear polynomial ([1 x y]) is used in this work, namely, m = 3; a i and b j are the unknown interpolation coefficients. There are many different types of RBF that can be used to formulate the RPIM, and different RBFs have different features [18]. In this work, the well-known multiquadrics (MQ) function is used to construct the required field function approximation owing to its several excellent characteristics. The expression of the MQ function is as follows [18,21]: in which r i denotes the distance from the field node to the sampling point, d c is the average nodal interval of the field nodes used, and α c and q denote two undetermined parameters that are closely related to the computation accuracy of the RPIM; q = 1.03 and α c = 1 are used in this work because very good numerical results can always be obtained for solid mechanics with these parameters.
With the aim to determine the coefficients a i and b j , Equation (1) should satisfy a series of reasonable constraint conditions. Firstly, it is usually assumed that the constructed field function approximation can exactly pass through the function values of all the nodes located in the support domain of the sampling point x; these constraints can be expressed by: R 1 (r 1 ) R 2 (r 1 ) · · · R n (r 1 ) R 1 (r 2 ) R 2 (r 2 ) · · · R n (r 2 ) . . . . . . . . . . . .
in which R 0 and P 0 are the so-called moment matrices corresponding to the RBF and PBF, respectively. To uniquely determine the unknown interpolation coefficients a i and b j , the following additional constraints should also be satisfied: The combination of all the constraining conditions shown in Equations (2) and (6) can result in the following matrix equation: Mathematics 2022, 10, 2889 4 of 20 Then, the undetermined interpolation coefficients can be calculated by Substituting the interpolation coefficients obtained into Equation (1) and following the standard formulation of the conventional RPIM, the required nodal interpolation shape function can be obtained by In the standard RPIM, the field nodes participating in building the field function approximation for the sampling point, which are usually quadrature points, are determined by a support domain. The shape of the support domain can be a square or a circle. The sampling point is usually the center of the defined support domain, while the background cells for numerical integration are always constructed independently of the support domain. As a result, the different sampling points (or quadrature points) in one integration cell may have different support domains, namely, the required field nodes to construct the field function approximation are different. In summary, since the moving support domain is used in the traditional RPIM, the support domain of the nodal shape functions always do not align with the background integration cells, which then leads to considerable numerical integration error and degrades the quality of the numerical solutions obtained.
To effectively overcome the abovementioned misalignment between the nodal shape function supports and the background integration cells, in this work a modified RPIM (M-RPIM) is employed to analyze the free vibration of two-dimensional solids. In this M-RPIM, a fixed support domain (as shown in Figure 1) rather than the moving support domain in the standard RPIM is used to select the required field nodes for the construction of the field function approximation. In other words, the identical field nodes are used for interpolation for any quadrature points in one background integration cell. The fixed support domain used can still be a square or a circle (the square support domain is used in this work); however, this fixed support domain is always centered by the geometrical center of the integration cell, not centered by the sampling points (which are usually the quadrature points) as in the conventional RPIM. The difference between the original RPIM and the present M-RPIM in constructing the field function approximation can be shown as follows: in which Ω Q stands for the moving support domains, which are centered by the quadrature points in one background cell, Ω * represents the fixed support domains that are directly centered by the centroids of the background integration cells.
Mathematics 2022, 10, x FOR PEER REVIEW in which Ω Q stands for the moving support domains, which are centered by quadrature points in one background cell, Ω * represents the fixed support dom that are directly centered by the centroids of the background integration cells.

Formulation of the Elastodynamics of Two-Dimensional Solids
Based on the small displacement assumption, the partial differential equation ( of the boundary-value problem for the elastodynamics of solids can be written by in which Ω denotes the problem domain considered, b stands for the body forc represents the stress tensor, ρ is the mass density, u is the displacement vector & & u signifies second derivatives of u .
As usual, the following two kinds of boundary conditions are always consid for the two-dimensional elastodynamics of solids: in which Γ E and Γ N denote the essential boundary condition and the na boundary condition, respectively; u and t are the imposed displacement vector traction vector on the corresponding boundary conditions.

Formulation of the Elastodynamics of Two-Dimensional Solids
Based on the small displacement assumption, the partial differential equation (PDE) of the boundary-value problem for the elastodynamics of solids can be written by in which Ω denotes the problem domain considered, b stands for the body force, σ represents the stress tensor, ρ is the mass density, u is the displacement vector and .. u signifies second derivatives of u.
As usual, the following two kinds of boundary conditions are always considered for the two-dimensional elastodynamics of solids: in which Γ E and Γ N denote the essential boundary condition and the natural boundary condition, respectively; u and t are the imposed displacement vector and traction vector on the corresponding boundary conditions. Using the boundary conditions shown in Equation (12) and following the virtual displacement principle, the weak form of Equation (11) for the elastodynamics of twodimensional solids can be obtained by in which δu and δε stand for the virtual displacement and strain, respectively. Using the Galerkin weighted residual techniques and the field function approximation shown in Equation (1), the matrix equation for the weak form shown in Equation (13) can be obtained [1,18] by the following relationship: in which M is the usual mass matrix, K is the usual stiffness matrix, F is the applied force vector and C is the matrix containing the damping effects. Without considering the damping effects and the external force, Equation (14) reduces to Equation (15) is the governing matrix equation obtained for the free vibration analysis of two-dimensional solids.
Assuming that the displacement solution to Equation (15) is time harmonic, namely, in which j = √ −1, ω denotes the angular frequency, and U is the amplitude of the displacement distribution.
Substituting Equation (16) into Equation (15), then Equation (15) can be rewritten as From Equation (17), we can observe that the typical eigenvalue problem should be solved to perform the analysis of free vibration problems.

Numerical Example
In this section, several typical numerical examples are considered to assess the capability of the proposed M-RPIM in free vibration analysis of the two-dimensional solids. For the convenience of discussion, the natural frequency values from the present M-RPIM are compared to those from the original RPIM and the standard finite element approach with bilinear quadrilateral elements (FEM-Q4). In all the numerical examples considered, identical node arrangements are employed for these three different numerical methods (M-RPIM, RPIM and FEM-Q4). For simplification, the quadrilateral meshes used are directly employed as the background cells to perform the numerical integration for the RPIM and M-RPIM, unless otherwise noted. To effectively examine and compare the accuracy and convergence of numerical solutions from the different numerical methods, the following relative error indicator is employed in this work: in which f num denotes the natural frequency results from the numerical methods (M-RPIM, RPIM and FEM-Q4) and f re f represents the reference natural frequency results, which are usually obtained from the commercial finite element software packages with a very refined mesh.

Free Vibration Analysis of the Cantilever Beam
Firstly, the free vibration of a cantilever beam is considered here. As shown in Figure 2, the geometric configuration of the cantilever beam has length L = 100 mm and height D = 10 mm. A unit thickness (t = 1 mm) is considered for this beam and, hence, this numerical example can be simplified as a plane stress problem. The material constants of this beam are taken as Young's modulus E = 2.1 × 10 11 Pa, Poisson's ratio v = 0.3 and mass density ρ = 8 × 10 3 kg/m 3 . The regular node arrangements are used to discretize the problem domain of this cantilever beam for the three different numerical methods. For a detailed analysis and discussion, a series of different node arrangement patterns with different nodal intervals are used here (see Figure 3).
Mathematics 2022, 10, x FOR PEER REVIEW 7 of 20 results, which are usually obtained from the commercial finite element software packages with a very refined mesh.

Free Vibration Analysis of the Cantilever Beam
Firstly, the free vibration of a cantilever beam is considered here. As shown in Figure 2, the geometric configuration of the cantilever beam has length L = 100 mm and height D = 10 mm. A unit thickness (t = 1 mm) is considered for this beam and, hence, this numerical example can be simplified as a plane stress problem. The material constants of this beam are taken as Young's modulus E = 2.1 × 10 11 Pa, Poisson's ratio v = 0.3 and mass density ρ = 8 × 10 3 kg/m 3 . The regular node arrangements are used to discretize the problem domain of this cantilever beam for the three different numerical methods. For a detailed analysis and discussion, a series of different node arrangement patterns with different nodal intervals are used here (see Figure 3).  . The different node arrangement patterns that are employed to discretize the cantilever beam for different numerical methods: (a) uniform mesh pattern, which is used to discretize the cantilever beam for the standard FEM-Q4; (b) node arrangement pattern used to discretize the cantilever beam for the RPIM and M-RPIM.

Computation Accuracy Study
Utilizing a series of different node arrangement patterns, the first twelve natural frequency solutions from the three numerical methods are listed in Tables 1-4. Among them, the corresponding RPIM and M-RPIM solutions are obtained when the size of the nodal support domain is taken as α s = 2.5h (h denotes average nodal interval of the meshes used). The reference solutions from eight-node quadrilateral element (FEM-Q8) with a very refined mesh pattern (average nodal interval h = 0.1 mm) are also provided in the tables for comparison.
From the results listed in the tables, we can observe that the original RPIM cannot always provide more accurate solutions than the standard FEM-Q4 in calculating the natural frequency values of this cantilever beam, although the higher order interpolation (not the bilinear interpolation in the FEM-Q4) is employed in the RPIM when the nodal support domain α s = 2.5h. This is mainly caused by the misalignment between the Mathematics 2022, 10, x FOR PEER REVIEW results, which are usually obtained from the commercial finite element so packages with a very refined mesh.

Free Vibration Analysis of the Cantilever Beam
Firstly, the free vibration of a cantilever beam is considered here. As sho Figure 2, the geometric configuration of the cantilever beam has length L = 100 m height D = 10 mm. A unit thickness (t = 1 mm) is considered for this beam and, this numerical example can be simplified as a plane stress problem. The m constants of this beam are taken as Young's modulus E = 2.1 × 10 11 Pa, Poisson's ra 0.3 and mass density ρ = 8 × 10 3 kg/m 3 . The regular node arrangements are u discretize the problem domain of this cantilever beam for the three different num methods. For a detailed analysis and discussion, a series of different node arrang patterns with different nodal intervals are used here (see Figure 3).  . The different node arrangement patterns that are employed to discretize the ca beam for different numerical methods: (a) uniform mesh pattern, which is used to discre cantilever beam for the standard FEM-Q4; (b) node arrangement pattern used to discre cantilever beam for the RPIM and M-RPIM.

Computation Accuracy Study
Utilizing a series of different node arrangement patterns, the first twelve n frequency solutions from the three numerical methods are listed in Tables 1-4. A them, the corresponding RPIM and M-RPIM solutions are obtained when the size nodal support domain is taken as α s = 2.5h (h denotes average nodal interval meshes used). The reference solutions from eight-node quadrilateral element (FE with a very refined mesh pattern (average nodal interval h = 0.1 mm) are also pro in the tables for comparison.
From the results listed in the tables, we can observe that the original RPIM always provide more accurate solutions than the standard FEM-Q4 in calculati natural frequency values of this cantilever beam, although the higher order interp (not the bilinear interpolation in the FEM-Q4) is employed in the RPIM when the support domain α s = 2.5h. This is mainly caused by the misalignment betwe . The different node arrangement patterns that are employed to discretize the cantilever beam for different numerical methods: (a) uniform mesh pattern, which is used to discretize the cantilever beam for the standard FEM-Q4; (b) node arrangement pattern used to discretize the cantilever beam for the RPIM and M-RPIM.

Computation Accuracy Study
Utilizing a series of different node arrangement patterns, the first twelve natural frequency solutions from the three numerical methods are listed in Tables 1-4. Among them, the corresponding RPIM and M-RPIM solutions are obtained when the size of the nodal support domain is taken as α s = 2.5h (h denotes average nodal interval of the meshes used). The reference solutions from eight-node quadrilateral element (FEM-Q8) with a very refined mesh pattern (average nodal interval h = 0.1 mm) are also provided in the tables for comparison. From the results listed in the tables, we can observe that the original RPIM cannot always provide more accurate solutions than the standard FEM-Q4 in calculating the natural frequency values of this cantilever beam, although the higher order interpolation (not the bilinear interpolation in the FEM-Q4) is employed in the RPIM when the nodal support domain α s = 2.5h. This is mainly caused by the misalignment between the constructed integration cells and the local support domains of the nodal interpolation functions. Owing to this misalignment, the integrands obtained in the original RPIM are not always continuously differentiable, then considerable numerical integration error is generated and leads to an additional loss in computation accuracy. However, from the tables we can observe that very good agreement between the M-RPIM solutions and the reference solutions can be achieved, and the M-RPIM solutions are much more accurate than the RPIM solutions. The main reason for this is that in the M-RPIM a fixed nodal support domain (not a moving support domain), which is built by the centroids of the integration cells, is directly used to perform the required numerical integration; then, the abovementioned misalignment between the integration cells and the local nodal support domains can be easily removed. As a result, the integrands obtained are completely continuously differentiable in the integration cells, so the numerical integration error can be markedly reduced and the computation accuracy can be significantly improved by the present M-RPIM for free vibration analysis. In addition, the vibration modes of the cantilever beam corresponding to the first twelve natural frequency values from the present M-RPIM are plotted in Figure 4; we can observe that the vibration modes obtained are quite stable and the physical mode shapes can be accurately achieved.

Convergence Study
In this subsection, the convergence performance of the numerical solutions from different numerical approaches is investigated in great detail. As shown in Figure 5, the comparison of the relative error (Re) results of the computed natural frequency values from different numerical methods versus the nodal interval (1/h) are given; the sign R in the legend of Figure 5 denotes the convergence rate of different numerical techniques. For simplicity, only the first two natural frequency values (Mode 1 and Mode 2) are considered here. From Figure 5, it can be observed that the convergence rate of the original RPIM is unexpectedly lower than the standard FEM-Q4 when the size of the nodal interpolation function support domain is taken as α s = 2.5h. This observation indicates that the misalignment between the integration cells and the local support domain of the nodal shape function in the original RPIM indeed can result in considerable numerical integration error; thus, the convergence rate can be markedly reduced.

Convergence Study
In this subsection, the convergence performance of the numerical solutions f different numerical approaches is investigated in great detail. As shown in Figure 5 comparison of the relative error (Re) results of the computed natural frequency va from different numerical methods versus the nodal interval (1/h) are given; the sign the legend of Figure 5 denotes the convergence rate of different numerical techniq For simplicity, only the first two natural frequency values (Mode 1 and Mode 2) considered here. From Figure 5, it can be observed that the convergence rate of original RPIM is unexpectedly lower than the standard FEM-Q4 when the size of nodal interpolation function support domain is taken as α s = 2.5h. This observa indicates that the misalignment between the integration cells and the local sup domain of the nodal shape function in the original RPIM indeed can resul considerable numerical integration error; thus, the convergence rate can be mark reduced.
However, from Figure 5 we also can see that the present M-RPIM is able to ach a higher convergence rate than the original RPIM and standard FEM-Q4. These find again demonstrate that the proposed program in this paper overcomes the misalignm between the constructed integration cells, and that the nodal shape function ind effectively supports suppression of possible numerical integration error. For vibration analysis of solids, therefore, the present M-RPIM has a higher convergence than the original RPIM and standard FEM-Q4. However, from Figure 5 we also can see that the present M-RPIM is able to achieve a higher convergence rate than the original RPIM and standard FEM-Q4. These findings again demonstrate that the proposed program in this paper overcomes the misalignment between the constructed integration cells, and that the nodal shape function indeed effectively supports suppression of possible numerical integration error. For free vibration analysis of solids, therefore, the present M-RPIM has a higher convergence rate than the original RPIM and standard FEM-Q4.

Computation Efficiency Study
From the analysis and discussion above, it can be observed that the present M-RPIM behaves better than the original RPIM and standard FEM-Q4 in terms of computation accuracy and convergence properties. However, the computation efficiency of the proposed M-RPIM is not yet studied. Note that computation efficiency is also a crucial index to assess the capabilities of numerical methods in engineering computation; comparison of the computation efficiency for the three different numerical methods is performed here. To analyze the computation efficiency, a series of different node arrangement schemes shown in Figure 3 are again employed. Figure 6 gives the comparison of the relative error results (Re) of the natural frequency values versus the computation cost for the three numerical methods. For simplicity, we still only consider the first two modes. From Figure 6, we can observe that the required computational cost for the standard FEM-Q4 is much less than for the original RPIM and the present M-RPIM when the identical node arrangement scheme is employed. This is because many more quadrature points are used to perform the numerical integration in RPIM and M-RPIM compared to standard FEM-Q4.

Computation Efficiency Study
From the analysis and discussion above, it can be observed that the present M-RPIM behaves better than the original RPIM and standard FEM-Q4 in terms of computation accuracy and convergence properties. However, the computation efficiency of the proposed M-RPIM is not yet studied. Note that computation efficiency is also a crucial index to assess the capabilities of numerical methods in engineering computation; comparison of the computation efficiency for the three different numerical methods is performed here. To analyze the computation efficiency, a series of different node arrangement schemes shown in Figure 3 are again employed. Figure 6 gives the comparison of the relative error results (Re) of the natural frequency values versus the computation cost for the three numerical methods. For simplicity, we still only consider the first two modes. From Figure 6, we can observe that the required computational cost for the standard FEM-Q4 is much less than for the original RPIM and the present M-RPIM when the identical node arrangement scheme is employed. This is because many more quadrature points are used to perform the numerical integration in RPIM and M-RPIM compared to standard FEM-Q4. Nevertheless, the computation accuracy of the standard FEM-Q4 cannot surpass the M-RPIM because a higher local approximation is used in this meshless numerical technique. From Figure 6, we also can observe that the original RPIM is actually numerically more expensive than the M-RPIM for the identical node arrangement scheme. This is because a moving support domain is used for different quadrature points in the original RPIM. In other words, for each quadrature point, the related operation in determining the support domain (namely the node selection for interpolation) should be performed once, while in the present M-RPIM, a fixed support domain is employed for any quadrature points within one integration cell; hence, the required operation in determining the support domain only should be performed once for each integration cell. Thus, in the M-RPIM, less computational cost is required in the node selection compared to the RPIM. Note that the present M-RPIM also has higher computation accuracy than the original RPIM; hence, the present M-RPIM also possesses higher computation efficiency than the original RPIM in engineering computation. This point can be clearly seen in Figure 6. From Figure 6, we also can observe that the original RPIM is actually numerically more expensive than the M-RPIM for the identical node arrangement scheme. This is because a moving support domain is used for different quadrature points in the original RPIM. In other words, for each quadrature point, the related operation in determining the support domain (namely the node selection for interpolation) should be performed once, while in the present M-RPIM, a fixed support domain is employed for any quadrature points within one integration cell; hence, the required operation in determining the support domain only should be performed once for each integration cell. Thus, in the M-RPIM, less computational cost is required in the node selection compared to the RPIM. Note that the present M-RPIM also has higher computation accuracy than the original RPIM; hence, the present M-RPIM also possesses higher computation efficiency than the original RPIM in engineering computation. This point can be clearly seen in Figure 6.

Free Vibration Analysis of the Cantilever Beam with Variable Cross-Section
The second numerical example considered here is a cantilever beam with variable crosssection. The geometric configuration of the variable cross-section beam is shown in Figure 7 and the related material constants are taken as Young's modulus E = 3 × 10 7 Pa, Poisson's ratio v = 0.3 and mass density ρ = 1 kg/m 3 . The regular node arrangement scheme is used to discretize this variable cross-section beam and the corresponding node distributions for the standard FEM-Q4 and the two meshless methods (RPIM and M-RPIM) are given in Figure 8. The first twelve natural frequency values computed using different numerical methods are listed in Table 5. Similar to the first numerical example, the corresponding natural frequency results from the eight-node quadrilateral element (FEM-Q8) with a very refined mesh pattern (5151 nodes and 5000 elements) are also provided as the reference solutions. It is clearly seen that the accuracy of FEM-Q4 results is worse than the original RPIM and the present M-RPIM results. However, the RPIM results are not more accurate than the M-RPIM ones, and the most accurate natural frequency solutions of this variable cross-section cantilever beam can be provided by the present M-RPIM. In addition, the first twelve mode shapes of this variable cross-section cantilever beam obtained from the proposed M-RPIM are depicted in Figure 9. It is easy to find that the eigenmode of this variable cross-section cantilever beam can be accurately predicted by the present M-RPIM. This numerical example demonstrates that the abilities of the original RPIM in engineering computation can be markedly improved by the present M-RPIM.

Free Vibration Analysis of the Cantilever Beam with Variable Cross-Section
The second numerical example considered here is a cantilever beam with variable cross-section. The geometric configuration of the variable cross-section beam is shown in Figure 7 and the related material constants are taken as Young's modulus E = 3 × 10 7 Pa, Poisson's ratio v = 0.3 and mass density ρ = 1 kg/m 3 . The regular node arrangement scheme is used to discretize this variable cross-section beam and the corresponding node distributions for the standard FEM-Q4 and the two meshless methods (RPIM and M-RPIM) are given in Figure 8. The first twelve natural frequency values computed using different numerical methods are listed in Table 5. Similar to the first numerical example, the corresponding natural frequency results from the eight-node quadrilateral element (FEM-Q8) with a very refined mesh pattern (5151 nodes and 5000 elements) are also provided as the reference solutions. It is clearly seen that the accuracy of FEM-Q4 results is worse than the original RPIM and the present M-RPIM results. However, the RPIM results are not more accurate than the M-RPIM ones, and the most accurate natural frequency solutions of this variable cross-section cantilever beam can be provided by the present M-RPIM. In addition, the first twelve mode shapes of this variable cross-section cantilever beam obtained from the proposed M-RPIM are depicted in Figure 9. It is easy to find that the eigenmode of this variable cross-section cantilever beam can be accurately predicted by the present M-RPIM. This numerical example demonstrates that the abilities of the original RPIM in engineering computation can be markedly improved by the present M-RPIM.

Free Vibration Analysis of the Cantilever Beam with Holes
The last numerical example is also a cantilever beam in plane stress condition. Unlike the previous numerical examples, the cantilever beam considered here has three identical holes (see Figure 10). The geometric parameters of this beam are given in Figure 10 and the material constants are taken as Young's modulus E = 2.1 × 10 11 Pa, Poisson's ratio v = 0.3 and mass density ρ = 8 × 10 3 kg/m 3 . The node arrangement scheme for the different numerical methods are plotted in Figure 11, and the average nodal interval h = 0.002 m. Similar to the previous two numerical examples, the first twelve natural frequency values from the different numerical methods are listed in Table 6, and the corresponding mode shapes from the present M-RPIM are given in Figure 12. In Table 6, the reference solutions are also computed from the eight-node quadrilateral element (FEM-Q8) with a very refined mesh pattern (average node interval h = 0.0001 m). Similarly, Table 6 and Figure 12 show that we obtain results similar to those in the previous two numerical examples, namely, the present M-RPIM can generate much more accurate numerical solutions than the original RPIM and FEM-Q4 for free vibration analysis; the present method has great potential for more complicated engineering computation.

Free Vibration Analysis of the Cantilever Beam with Holes
The last numerical example is also a cantilever beam in plane stress condition. Unlike the previous numerical examples, the cantilever beam considered here has three identical holes (see Figure 10). The geometric parameters of this beam are given in element (FEM-Q8) with a very refined mesh pattern (average node interval h = 0.0001 m). Similarly, Table 6 and Figure 12 show that we obtain results similar to those in the previous two numerical examples, namely, the present M-RPIM can generate much more accurate numerical solutions than the original RPIM and FEM-Q4 for free vibration analysis; the present method has great potential for more complicated engineering computation.  Figure 11. The node arrangement patterns employed to discretize the cantilever beam with three identical holes for the different numerical methods: (a) mesh pattern used to discretize the cantilever beam with three identical holes for the standard FEM-Q4; (b) node arrangement pattern used to discretize the cantilever beam with three identical holes for the RPIM and M-RPIM.  element (FEM-Q8) with a very refined mesh pattern (average node interval h = 0.0001 m). Similarly, Table 6 and Figure 12 show that we obtain results similar to those in the previous two numerical examples, namely, the present M-RPIM can generate much more accurate numerical solutions than the original RPIM and FEM-Q4 for free vibration analysis; the present method has great potential for more complicated engineering computation. (a) (b) Figure 11. The node arrangement patterns employed to discretize the cantilever beam with three identical holes for the different numerical methods: (a) mesh pattern used to discretize the cantilever beam with three identical holes for the standard FEM-Q4; (b) node arrangement pattern used to discretize the cantilever beam with three identical holes for the RPIM and M-RPIM.  Figure 11. The node arrangement patterns employed to discretize the cantilever beam with three identical holes for the different numerical methods: (a) mesh pattern used to discretize the cantilever beam with three identical holes for the standard FEM-Q4; (b) node arrangement pattern used to discretize the cantilever beam with three identical holes for the RPIM and M-RPIM.

Conclusions
In this work, a modified radial point interpolation method (M-RPIM) is proposed to enhance the capacities of the original RPIM for the free vibration analysis of two-dimensional solids. In the present M-RPIM, the numerical approximation established in integration cells is continuously differentiable while the corresponding numerical approximation in the original RPIM is always not continuously differentiable. Therefore, the possible numerical integration error in the original RPIM can be markedly reduced by the present M-RPIM.
Several supporting numerical examples are employed to investigate fully and in detail the performance of the proposed M-RPIM in solving free vibration problems. It is demonstrated that the proposed M-RPIM not only is able to surpass the original RPIM and the standard FEM-Q4 in terms of computation accuracy and convergence properties when the identical node arrangement scheme is employed, but the proposed method also has higher computation efficiency. This is because the fixed support domain is employed for any quadrature points in the integration cells; hence, the additional operations to determine the support domain for each quadrature point are not required. Owing to these excellent features, the present M-RPIM has great potential for solving more complex problems in practical engineering application.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.