Abstract
The meshless Shepard and least-squares (MSLS) interpolation is a newly developed partition of unity- (PU-) based method which removes the difficulties with many other meshless methods such as the lack of the Kronecker delta property. The MSLS interpolation is efficient to compute and retain compatibility for any basis function used. In this paper, we extend the MSLS interpolation to the local Petrov-Galerkin weak form and adopt the duo nodal support domain. In the new formulation, there is no need for employing singular weight functions as is required in the original MSLS and also no need for background mesh for integration. Numerical examples demonstrate the effectiveness and robustness of the present method.
1. Introduction
Meshless methods have prospered both in theory and application in engineering problems in the past two decades as they offer the possibility of a discretised approach without the occurrence of mesh entanglement requiring remeshing. A wide range of meshless methods have been proposed as outlined in recent surveys [1–3]. Remarkable successes have been reported in applying these methods for analyzing challenging engineering problems, namely, fracture modeling [4–7], plate and shell analysis [8–15], three-dimensional problems [16–18], fluid structure interaction analysis [19], strain localization problems [20], large deformation problems [21], and other applications [22–29].
Some currently popular meshless approximations are the moving least-squares (MLS) approximation, Shepard shape functions, partition of unity (PU), radial basis functions (RBF), reproducing kernel particle method (RKPM), point interpolation (PI), and Kriging interpolation (KI). Among them, the MLS approximation [30] is one of the most widely used approximations at present due to its global continuity, completeness, and robustness. However, the MLS approximation suffers from a number of problems that practically limit its application, namely, the high computational cost in obtaining the shape functions and their derivatives, difficulty in retaining accuracy with respect to nodal arrangement, and also the difficulty with which essential boundary conditions can be imposed due to the lack of the Kronecker delta property. Efforts have been made to address these problems by various means. Breitkopf et al. [31] developed the analytical forms for computing shape functions and diffuse derivatives of shape functions by assuming that some terms are constant and complete derivatives of shape functions. However, these formulations are dependent on the number of nodes and the formulation grows unwieldy when there are a large number of nodes in the support domain. Orthogonal basis functions pertaining to MLS for efficient and accurate computation of shape functions are investigated in [32, 33]. Singular weight functions are introduced by Kaljević and Saigal [34] to produce an interpolatory MLS approximation for the direct imposition of essential boundary conditions. Chen and Wang [35] developed a matrix transformation method for the imposition of boundary conditions; however, the formulation is complicated for implementation. The Shepard and least square (MSLS) interpolation developed by the authors [36] satisfactorily keeps the consistency of the Shepard shape functions up to the order of basis function and satisfies the delta property. However, singular weight functions have to be used to enforce the interpolating property of shape functions, which results in the loss of the smoothness of the interpolation as well as local oscillation. To eliminate this problem, an improved MSLS interpolation possessing the delta property without using singular weight functions as PU was proposed in [37].
Apart from the interpolation, the integration scheme of the weak form is also an important factor affecting the solution accuracy, which has been the problem for many meshless methods. Background cells have usually been used to integrate the weak form as has been the case with the EFGM, RKPM, and PIM [38]. Due to the complexity of the shape functions, a large number of integration points are needed to avoid the underestimation of the weak form, which is computationally expensive but still not adequate to give accurate solutions. Furthermore when an irregular nodal arrangement is used, the background mesh has to be refined where nodes are densely distributed [39]. To remove this difficulty, the nodal integration scheme which is free from the background cell was proposed in [40] and later used in [41, 42]. However, the performance of this scheme is unstable and also reduces the accuracy of the results. The MLPG-type meshless methods well solved this problem in a natural way by partitioning the local domain into a number of subdomains that may or may not overlap. In this way, integration of high order accuracy can be obtained for the global stiffness matrix [43] without background mesh. A similar but conceptually different approach was developed in [44] using the partition of unity quadrature scheme. However, the algorithms are complicated and more computationally demanding compared to the classical quadrature on subcells. Although many different approaches have been carried out for the weak form integration, finding a simple, efficient, and accurate integration scheme for meshless methods remains an open question. In this paper, the local weak form along with numerical integration over local subdomains is used to derive the discrete equations.
The content of the paper is outlined as follows. In Section 2, formulation of the improved MSLS interpolation is described in detail including the local approximation and support domain with dual definitions. The Kronecker delta property of the interpolation is also proved. In the Section 3, the discretised formulation of the present interpolation is derived using the local Petrov-Galerkin weak form. It is followed by numerical tests demonstrating the convergence characteristic and accuracy of the present interpolation in Section 4. Discussions are given at the end highlighting the features of the present method and suggestions on further studies.
2. Formulation of the MSLS Based on Duo Nodal Supports
In this section, the MSLS based on duo nodal supports is described in detail. We start the description of the formulation using a 2D problem domain of arbitrary shape as shown in Figure 1. The formulation is described in the context of elastostatics, with the fundamental field variable being a displacement. For an arbitrary node , its displacement vector in 2D is , where and are the nodal displacements in the and directions, respectively (the following formulation is derived only for in the direction but an identical process can be used for in the direction). The interpolation at an arbitrary point inside the domain in the direction is expressed as where is a set of shape functions that forms a partition of unity (PU); that is, , is the node index, and is the number of the nodes whose supports include point as shown in Figure 2. here is not the nodal displacement in the FEM or the “fictitious” nodal values in the EFGM but the local approximation of node at where the superscript indicates local. Shepard shape functions used as PU are given by which is the same as in the original MSLS. The construction of the MSLS interpolation takes the following steps: firstly, construct the local approximation at each node and secondly apply the PU approximation to the local approximation to get the interpolation. The definition of the nodal support domain will be given in detail and the local approximations at a node will be described.
2.1. Local Approximation at a Node
The local approximation at an arbitrary node is given by where is the nodal displacement for the th node in support of node , is the total number of nodes in the local cover of node as shown in Figure 1, and is the modified least square shape function given by where are the modified least square shape functions of node and are determined by the following equations: is a polynomial basis, and is the number of monomials in the basis. In the development of the MSLS interpolation, we use a bilinear basis throughout in 2D such that , and It can be seen from (4) that , for and . Thus It has been proved in [36, 37] that a singular weight function used as will enforce the interpolation, the equivalent equation (1) here, satisfying the delta property. A similar approach has been used by [45] to produce interpolatory MLS approximation. However, the use of singular weight function will bring some other problems such as the loss of smoothness of interpolation.
2.2. Duo Support Domain of a Node
The support domain of a node is the area where a node exerts influence on the field variable. In this paper it is defined as a circle centered on that node although it may take other shapes such as a rectangle. Here, two support domains are defined at each node, one is used in the construction of local approximation and the other in the PU approximation. In Figure 2, for example, node has two support domains associated with it, namely, with radius and with radius . If a node, for example, node in Figure 1 falls inside , then node will be used in constructing the local approximation at node . Similarly, if a point, point in Figure 2, for example, is contained in , then the local approximation will contribute to the PU approximation at . For an arbitrary node, for example, node in Figure 1, the size of is defined by where is a scale factor that ranges between 1.0 and 2.0, is a coefficient such that for a node lying on the boundary and for all other nodes, and is the distance between and the fifth nearest neighbor node to . If there is a predefined triangular background mesh, can be defined as the maximum distance between and the nodes of triangles which are connected to the node .
For a node having its local support domain completely inside the domain, for example, the subdomain of node in Figure 2, the size of is the same as : For a node having its local support domain close to or intersecting the boundary, for example, node shown in Figure 2, the definition of subdomain follows these steps. Firstly, find the nearest boundary node to among the neighbor nodes which have been used as nodes in defining , and secondly calculate the distance between the nearest boundary node and , denoted as , and then the size of is set by If there is an interior node where prescribed values needed to be applied, the procedure described above for setting the support domain of near boundary nodes can be repeated to that node using (10) by assigning this interior node as a boundary node. If we want all nodes to take nodal values at the nodes, the size of the can be taken as the distance between the and its nearest node for every node . The following quadratic spline function is used as the weight function over support domain in (2): where is the distance between the point and node and is the coordinate of node . The aim of separately defining local domain and support domain is to produce MSLS interpolations having the delta property without using a singular weight, so that the difficulties associated with the use of singular weight function can be removed. Indeed this aim is achieved here if the domain for local approximation and domain for PU are defined by the method described above as will be proved later in Section 2.3.
2.3. Delta Property at a Node
Consider a boundary node to be applied with boundary conditions. If the support domain of the nodes is set according to (9) and (10), then the will be the only node contained in . Thus the MSLS interpolation equation (1) at becomes As there is only one node in the PU, then (2) becomes It is known by (7) that the local approximation at node satisfies Substituting (13) and (14) into (12) gives Hence, the present MSLS interpolation takes nodal value at boundary nodes and the essential boundary conditions or point load conditions can be directly imposed as in the FEM. Also, the present interpolation preserves the consistency up to the order of the basis function, which is a necessary requirement of accuracy. The proof is the same as has been presented in [36, 37].
3. The Local Petrov-Galerkin Weak Form
Let be the total number of nodes associated with the given point ; then (3) can be rewritten as where is the vector of Shepard shape function, is a matrix comprising the modified least square point interpolation (LSPI) shape functions, and is the MSLS shape function.
For domain bounded by (Figure 3), the equilibrium equations and boundary conditions of linear elasticity are given by where is the stress tensor, are the body forces, are the unit normal to the domain, and and are the global boundaries with prescribed displacements and tractions, respectively. Similar to [36], the local polygonal subdomains are constructed for the purpose of simplifying the integration and the discrete equations. For example, based on Delaunay algorithm, a local polygon is constructed by using the nodes in local cover , as shown in Figure 3. A generalized local weak form of the equilibrium equation in (17) is written as where is the integration domain or subdomain for node and is the test function. Using the divergence theorem in (18), we obtain the following local weak form: where is the outward unit normal to the boundary . The boundary for the subdomain is usually composed of three parts: the internal boundary , the boundary , and , over which the essential and natural boundary conditions are specified. Substituting in (19), the following is obtained: In order to simplify (20), we can deliberately select the three-node triangular FEM shape functions , which correspond to the node of the triangles constructing the polygonal subdomain , as test functions , such that they vanish over . Substituting shape functions for in (20), we obtain the following local weak form: For a local polygonal subdomain located entirely in the global domain , there is no intersection between and the global boundary , and the integrals over and in (21) vanish. For a local polygonal subdomain near the boundary, the first item of (21) can also be omitted because of the properties of the test functions . Substituting the MSLS approximation in (21) into the above equation leads to the following discretised system of linear equations: denoted as where is the elasticity matrix: Equation (23) can be individually integrated over each triangle constructing the local subdomain , as shown in Figure 3. In the present work, seven Gaussian points are used in each triangle.
4. Numerical Examples
The proposed MSLSM interpolation has been coded in C++. In this section, we show the performance of the present interpolation on a range of test problems. The results are compared with the exact solutions, the MLPG solutions, and the linear FEM solutions. The weight functions used in the MLPG for testing purpose are the Gaussian type weight functions given by where is defined by (8) and is used for all test examples. The scale factor in (8) is set to be 1.5 and the linear bases are used in MLPG and MSLSM. To study the convergence behavior we define the following error norms in displacement and energy, respectively: where is a vector collecting nodal displacement results and where is the infinitesimal strain tensor and is the Cauchy stress tensor. The relative displacement error and energy error are calculated by where the superscripts and refer to numerical solutions and exact (or reference) solutions, respectively.
4.1. A Constant Strain Patch Test
A constant strain patch test [46] using three distributions of 7, 28, and 126 irregular nodes is shown in Figure 4. Young’s modulus and Poission’s ratio of the material are 1000 and 0.25, respectively. The thickness of the plate is taken as a unit following plane stress assumption. Since the exact solution is linear, a linear basis for the MSLS interpolation is able to represent this solution. The computational results in Table 1 show that the present MSLSM passes the patch tests exactly up to the double precision of the computer.
4.2. A Cantilever Beam
A cantilever beam problem with dimensions of m and m, as shown in Figure 5, is tested first. The beam is subjected to parabolically distributed downward traction equivalent to a unit load at the right-hand end and is constrained at the left-hand end as shown in Figure 5. The elastic material properties used are Pa and , and the problem is solved under a plane strain assumption. We refer to the analytical solution of the problem given in [47]. The convergence of the present method is studied using four nodal arrangements with 50, 138, 402, and 965 nodes.
(a) Cantilever beam
(b) Nodal arrangement
The convergence rate is compared among FEM, MLPG, and the present MSLSM in Figures 6 and 7. It can be found that MSLSM shows a good accuracy and convergence rate. Figures 8 and 9 show the vertical displacement and the normal stress along the indicating accuracy of results using irregular 138 nodes and the proposed formulation.
As has been highlighted in [36, 37], the computational cost in obtaining the shape functions and its derivatives is much lower by the present LS interpolation than by the MLS approximation. This is shown in Table 2. It can also be observed from the table that the difference in computational efficiency between the two interpolations gets bigger when the number of nodes increases.
4.3. An Infinite Plate with a Circular Hole
The second example is an infinite plate with a circular hole of radius m. The plate is subjected to a far field traction Pa in the direction. A finite portion of the plate is considered for analysis and, due to the symmetry of the problem, only a quarter of the portion requires modeling, as shown in Figure 10. The elastic material properties used are Pa and and plane stress conditions are assumed. The stresses and displacements for this are given in an analytical solution in [47] as where is the shear modulus and is the Kolosov constant where for the plane strain assumption.
(a) A small portion taken for analysis
(b) Boundary conditions applied
Traction-prescribed boundary conditions consistent with the exact solution in (32) are applied at the top and right edges. Four distributions of 53, 188, 564, and 1012 nodes, which are shown in Figure 11, are employed for the convergence studies. Figures 12 and 13 show that the MSLS has a good convergence rate in the displacement and energy norm. The convergence slope of the present method is similar to that of MLPG though the latter is seen to be more accurate. The displacement obtained using MSLS and MLPG is shown in Figure 14. It is seen that the displacement values given by MSLS as well as MLPG are very close to the exact solution.
5. Conclusions
In this paper, we proposed a local weak form meshless Shepard and least-squares interpolation. The interpolation features the use of duo nodal supports for local approximation and global approximation, respectively. The present formulation satisfies the delta property at desired nodes without using singular weight functions. The shape functions constructed conform to the PU property. The local Petrov-Galerkin weak form is used so that there is no need for background mesh. It has been tested that the present formulation is more efficient when compared to the currently widely used MLS approximation. For the 2D example of 965 nodes, the computing time of the present MSLS is close to half of the computing time of MLS. The accuracy of the present method is slightly lower than the MLPG but higher than the FEM when the same number of nodes and the same order of basis function, for example, linear or quadratic bases, are used (linear basis corresponding to triangular element in the FEM). The formulation here is derived for 2D analysis but is readily extendable to 3D and the essential ideas are the same. Further development of the present interpolation such as application to problems of changing geometry/moving boundaries, finite deformation, elastoplasticity, and three-dimensional cracking is ongoing.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the support of the NSFC Programme (41130751, 51109162), 973 Program (2011CB013800), Shanghai Pujiang Talent Program (12PJ1409100), SRF for ROCS of State Education Ministry, and China Central University Funding.