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.

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 [13]. 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 [615]. 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 [1827].

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 [4042], 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. Section 2 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 can be written as follows:where are isoparametric shape functions used in the classic QUAD4 element and denotes the nodal displacement function vector which can be constructed by the nodes in the node support domain which is defined by the element including node . As shown in Figure 1, the node support domain for node is marked by a dashed line, which includes nodes . Regarding the unknown nodal displacement functions, they can be interpolated by the polynomial basis function and radial basis function , which can be expressed as follows:

In the above, is the number of polynomial terms and denotes the number of nodes in the node support domain .

The vectors and in equation (2) are vectors yet to be determined. and denote the corresponding polynomial and radial basis functions, respectively. If we consider three respective terms for 2D problems, can be written as follows:

The radial basis function vector is expressed aswhere is the function of distance , which denotes the space distance between the selected interpolation point and a node . Also for the 2D problem, can be expressed as

Enforcing equation (2) to pass through all the nodes in the nodal support domain yields the following equation:where is a nodal displacement vector of all the nodes in the nodal support domain and and are, respectively, the vector of unknown coefficients for radial basis functions and polynomial basis functions:

The related matrices and are

The 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 can be solved:

Substituting equation (11) into equation (10), the vector can be obtained:

Eventually the nodal displacement function can be expressed by substituting equations (11) and (13) into equation (2):where

The arrays , are to be assembled into a matrix , in which is the number of nodes in the j-th element support domain which is defined by the element including nodes in element . As shown in Figure 2, the element support domain for element is marked by a dot dash line.

Substituting equations (10) and (11) into equation (1) allows us to rewrite the element displacement field as follows:

In the above, denotes the x-direction displacement of all the nodes in an element support domain , and the composite shape function array, , can be represented as

Similarly, the element displacement field in y-direction is written as follows:where denotes the y-direction displacements of all nodes in .

The Cartesian derivatives of shape function are obtained by

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 Table 1. 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.

In 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:

In the above, denotes a compactly supported function, represents a radial basis function, and is the size of the support domain. Generally, can be chosen as a Heaviside step function, if , , otherwise, . However, as the nodal support domain and the corresponding nodes have been defined, does not need to be used, and equation (21) can be rewritten as

The following radial basis function in a compactly supported form based on the original kernel function of the quartic spline (QS) [6] is applied:

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:where is the scale factor and 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. is determined as follows for the irregular nodal distribution [56]:where denotes the area of the nodal support domain and 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, 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 , traction on boundary , and displacement boundary conditions on , where , . The virtual weak form of the dynamic equation can be expressed aswhere and 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 aswhere is the nodal displacement vector and represents the standard displacement gradient matrix. For 2D problems, it can be determined as follows:

The stress can be obtained by the following constitutive relation:where 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 asorin which

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]:

A general solution to this system iswhere denotes the imaginary unit, is time, indicates the eigenvector, and is the natural frequency. Substituting equation (37) into equation (36) leads to the following eigenvalue equation for natural frequency :

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 which is assumed to be a linear combination of and can be utilized aswhere 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 [5860], is used in this work. Two equal substeps are included in the complete time step 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:

In the second substep, the following equations are expressed by the three-point Euler backward method:

Similarly, the control equation (34) at time can also be expressed as

Once and is solved by equations (41) and (43), , , , and 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:where superscript “exact” denotes the exact value of the problem, “num” denotes result solved by numerical solutions, and is the number of total field nodes.

The strain energy of the numerical solution and the total strain energy of the exact solution are defined as follows:where is the strain of the exact solution. In the actual computation on equation (45), we used a very fine mesh () to calculate the “exact” strain energy .

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. 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 obtainwhere and represent the row vector of monomials and the corresponding vector of coefficients, respectively, and denotes the arbitrary maximum order of the monomials. Regarding each of these monomials in equation (47) as separate exact solution, a set of numerical solutions can be obtained as

For the distributed nodes in the patch domain, a pointwise error vector can be evaluated aswhere and denote the numerical solution and the corresponding exact one, respectively.

Summing the resulted errors for each monomial, we can then obtain the error as

The Euclidian error can be evaluated aswhere represents a positive semidefinite matrix, and the maximum eigenvalue is derived aswhere is the maximum eigenvalue and is the unit matrix. On condition that the patch test is passed, the error norm and both will be zero.

In this study, we conduct a patch test based on polynomials of first order, hence taking in equation (47) leads to

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 aswhere and are the initial regular meshes’ sizes along x and y axes, respectively. denotes a random number between −1.0 and 1.0 and is a prescribed irregularity factor between 0.0 and 0.5. A larger 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.

4.2. 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. 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]:where the moment of inertia of the beam is given by .

The corresponding stresses are

The parameters are , , , , and .

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. Figure 5 shows the steady condition number for different , especially . As shown in Figure 6, if the irregularity factor , the errors of the deflection of point A are almost same. If the irregularity factor , 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 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.

4.3. Example 3: Cantilever Beam Free Vibration Analysis

As shown in Figure 13, a cantilever beam fixed on the left end with length = 100 m and height = 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 = 1.0 m, Young’s modulus , Poisson’s ratio , and mass density . Based on the Euler–Bernoulli beam theory, the fundamental frequency can be obtained as a reference.

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 . This unique feature of the presented radial basis function enables constant results for any correlation parameter. Accordingly, correlation parameter is used for all subsequent computations.

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 Tables 47. 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.

Figure 15 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].

4.3.3. Sensitivity to Distorted Meshes

The cantilever beam with irregular mesh shown in Figure 13 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 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.

4.4. Example 4: Shear Wall Free Vibration Analysis

In this part, a shear wall with four openings, as shown in Figure 19 is investigated. The parameters in the computation include Young’s modulus  N/m2, Poisson’s ratio , thickness  m, and mass density  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.

The natural frequencies of the first eight modes are tabulated in Table 9. 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].

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 Figure 22. The rod is discretized with 723 irregularly distributed nodes, as shown in Figure 23. The plane stress problem is considered with E = 10 GPa, = 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).

4.6. Example 6: Cantilever Beam Implicit Forced Vibration Analysis

The benchmark cantilever beam and mesh shown in Figures 24 and 25 are investigated via the Bathe method. The plane strain problem is considered with numerical parameters  m,  m,  m,  N/m2, ,  kg/m3, , and . The domain is represented by 10 × 2 elements. The first loading is in the harmonic form with . The corresponding dynamic responses can be shown in Figures 26 and 27 with and 0.05, respectively. Then, a constant step load is added at the apex from . The time step 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 2628, 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.

4.7. 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  m,  m, ,  N/m2, , and  kg/m3. A reference solution is obtained by ABAQUS QUAD8 with 2,000 elements. The first loading is applied in the harmonic form . Its dynamic responses can be observed from Figures 31 and 32 with and 0.05, respectively. A constant step load is then added at the apex from . The time step is used for time integration. The dynamic responses shown in Figures 3133 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.

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 (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)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.

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the National Natural Science Foundation of China under grant no. 11472137 and the Fundamental Research Funds for the Central Universities under grant no. 309181A8801.