Abstract
Combining the hybrid displacement variational formulation and the radial basis point interpolation, a truly meshless and boundary-only method is developed in this paper for the numerical solution of solid mechanics problems in two and three dimensions. In this method, boundary conditions can be applied directly and easily. Besides, it is truly meshless, that is, it only requires nodes generated on the boundary of the domain, and does not require any element either for variable interpolation or for numerical integration. Some numerical examples are presented to demonstrate the efficiency of the method.
1. Introduction
Boundary integral equations (BIEs) have widely been used for the numerical solution of a variety of boundary value problems in solid mechanics as they can reduce the computational dimensions of the original problem by one and give a simple discretization of the exterior problems. The numerical discretization of BIEs is commonly known as the boundary element method (BEM) [1]. For many problems, the BEM is undoubtedly superior to the “domain discretization’’ types of methods, such as the finite element method (FEM) and the finite difference method. In the BEM, for example, only the two-dimensional bounding surface of a three-dimensional body needs to be discretized. However, as the FEM, the BEM depends on the generation of meshes, adapted or not. In some cases, this can be time-consuming and very difficult.
Meshless (or meshfree) methods for numerical solutions of boundary value problems have been developed for alleviating the meshing-related difficulties. Compared to the FEM and the BEM, the core of this type of method is to get rid of, or at least to alleviate, the difficulty of meshing and remeshing the entire structure by simply adding or deleting nodes. Some domain-type meshless methods, such as the element-free Galerkin method, the meshless method, the reproducing kernel particle method, the point interpolation method, and the meshless local Petrov-Galerkin (MLPG) methods, have been proposed and gained great success in solving a wide range of problems for solids. These meshless methods are domain based, as the FEM, in which the problem domain is discretized. For an extensive overview on the subject of meshless methods, containing most of the previously proposed methods, some monographs [2, 3] can be read.
The meshless idea has also been used in BIEs. The first BIEs-based meshless method was the boundary node method (BNM) [4, 5]. This approach takes the advantages of both BIEs in dimension reduction and moving least squares (MLSs) approximations in elements elimination. Nevertheless, background cells are required for numerical integration. In order to get rid of background cells, Atluri et al. proposed the meshless local boundary integral equation (LBIE) method [6]. Although absolutely no domain and boundary elements are required, the LBIE method is not strictly a boundary method since it requires evaluation of integrals over certain surfaces (called in [6]) that can be regarded as “closure surfaces’’ of boundary elements. To avoid this hindrance, Zhang et al. proposed the hybrid boundary node method (HBNM) that combines the MLS approximation with the hybrid displacement variational formulation [7]. In this method, the integration is limited to a fixed local region, thus no cells are needed either for interpolation or for integration.
However, because the MLS approximation lacks the delta function property of the usual FEM and BEM shape functions, it is difficult to impose boundary conditions in the BNM and the HBNM. This problem becomes even more serious in boundary-type meshless methods, since a large number of boundary conditions need to be satisfied [2]. The technique used in the BNM involves a new definition of the discrete norm used for the construction of the MLS approximation and thus doubles the number of system equations. This technique is also employed in the HBNM, together with the addition of a penalty formulation. In the BNM and the HBNM, the basic unknown quantities are approximations of the nodal variables. They are not the real nodal variables, and thus boundary conditions could not be directly applied. To restore the delta function property of the MLS, Liew et al. presented an improved MLS approximation that uses weighted orthogonal polynomials as basis functions [8]. The improved MLS has been inserted into BIEs to develop the boundary element-free method (BEFM). The BEFM is a direct boundary type meshless method, which has been used for problems in linear elasticity [8], elastodynamics [9], and potential theory [10]. Additional, via combining a variational form of BIEs and the MLS approximation, another technique is developed in the Galerkin BNM to impose boundary conditions [11, 12].
The radial basis point interpolation (RBPI) [2, 13] is another meshless interpolation technique that uses radial basis functions and polynomial basis to construct meshless shape functions. Compared to the MLS approximation, the shape functions so constructed have the delta function property. Consequently, the construction of an RBPI shape function is more efficient than the MLS procedures. There have been many meshless methods based on the RBPI for the numerical solution of mathematical problems in engineering. Typical of them are the radial point interpolation method [13], the local radial point interpolation method [14], and the local boundary integral equation method based on radial basis functions [15]. Besides, some boundary type meshless methods are developed by the combination of the RBPI with BIEs, such as the boundary radial point interpolation method (BRPIM) [16] and the hybrid BRPIM [17]. Since the RBPI shape functions possess Kronecker delta function properties, these BRPIMs have some advantages. It is easy to enforce boundary conditions, and it is computationally efficient. However, similar to the BNM and the BEFM, an underlying background cell structure is still used in these BRPIMs for numerical integration.
This paper extends the previous works to develop a truly meshless and boundary-only method for the numerical solution of two- and three-dimensional problems in solid mechanics by a combined use of the RBPI with the hybrid displacement variational formulation. This method is called the hybrid radial boundary node method. As has been illustrated in [18], the RBPI is used in this method to construct shape functions with delta function properties based on arbitrary distributed boundary nodes. So unlike the BNM and the HBNM, the present method is a direct numerical approach in which the basic unknown quantity is the real solution of nodal variables, and boundary conditions can be applied directly and easily, which leads to greater computational precision. Besides, with the help of the hybrid displacement variational formulation, the present method does not involve any mesh for either interpolation or integration. Thus, the inherent inefficiency of the BNM, the BEFM, and the BRPIMs due to the use of the background integration cell is alleviated in this novel method, which leads to tremendous improvement in computational efficiency.
The following discussions begin with a description of the boundary variational principle for solid mechanics problems in Section 2. Then, a detailed variables interpolation is provided in Section 3. Section 4 assembles the final system of equations. Numerical examples are given in Section 5. Section 6 contains some conclusions.
2. Boundary Variational Principle for Solid Mechanics Problems
Let be a bounded domain in with boundary . A general point in is denoted by or .
Consider the following two- and three-dimensional problem of solid mechanics where is the unknown displacement field; and are given Lamé constants; , and stand for the Laplacian, gradient, and divergence operators, respectively; and are prescribed boundary displacement and traction, respectively; and is the following conformal derivative operator: where is the unit outward normal on .
Based on the modified variational principle, the solution of problem (2.1) is the function which minimizes the following functional [19]: where is the boundary traction, is the boundary displacement satisfying on , , is the shear modulus and is the Poisson’s ratio.
We introduce the stress tensor and strain tensor where is the Kronecker delta symbol. Then, via carrying out the first variation of (2.3) one gets Thus, by setting to zero we obtain the following equivalent integral equations If the traction boundary condition is imposed, (2.8) will be satisfied and thus it can be ignored in what follows.
Since the variational principle is a universal theory, (2.6) and (2.7) should be satisfied in any subdomain , which is bounded by and and contains the boundary node . Figure 1 depicts the sketched geometric configuration for both two- and three-dimensional cases. Therefore, to obtain a truly boundary-only meshless method, (2.6) and (2.7) are replaced by where is a weight function or a test function. We emphasize that in the above equations the shape and dimension of the subdomains can be arbitrary. This observation forms the basis for the present truly meshless method. Obviously, a circle (or a sphere) is the simplest regularly shaped subdomains in (or . Hence, the subdomain is chosen as the intersection of the domain and a circle (or a sphere) centered at the boundary node (see Figure 1).
(a)
(b)
3. Variables Interpolation
Assume that the boundary is the union of piecewise smooth segments , . To avoid the discontinuity at corners and edges, the point interpolation for displacement and boundary traction on each is constructed independently. So in the following interpolation scheme, let the variable denote or for simplicity. Then the radial basis point interpolation for can be defined as where is a parametric coordinate on ; is a radial basis function (RBF); is the node number of the interpolation domain; is a monomials; is the number of monomials; ; ; ; ; and are unknown coefficients.
The effectiveness and accuracy of the interpolation depends on the choice of the RBFs. A number of different types of RBFs [20] such as linear distance functions, thin plates plines, multiquadrics, Gaussians and RBFs with compact supports have been proposed and may be used for this purpose. Characteristics of radial functions have been widely investigated. The variable in RBFs is only the distance. Hence, the forms of interpolation formulations are the same for both two-dimensional problems and three-dimensional problems. The following multiquadrics radial function is used in this work (other RBF can also be used similar) where and are shape parameters. In this situation, the required polynomial basis is linear as .
Letting the approximation represented by (3.1) pass through all the boundary nodes, we get in which is a vector, while and are matrices. Besides, in order to guarantee unique solution, the following constraints should also be satisfied Then, according to (3.3) and (3.4), we obtain and , where and . Consequently, substituting and into (3.1) yields where the shape function vector is
Now, in (2.9) and (2.10), and on can be represented as where and are the nodal values and , respectively; is the contributions from the node to the evaluation point , which is not equal to zero in the interpolation domain of the th node only.
In (2.9) and (2.10), and on are not defined yet. However, this problem can be tackled by selecting the weight function such that the size of the support of is less than the radius of the subdomain , then all integrals over vanish. A variety of weight functions have been investigated in the past for meshless methods. In this paper, Gaussian weight function is chosen and can be written as where is the radius of ; is the distance between a sampling point , in domain , and the nodal point ; and is a constant controlling the shape of the weight function . As a result, vanishes on .
Finally, the domain variables and in (2.9) and (2.10) are interpolated by the fundamental solution as where is the unknown parameter; is the total number of boundary nodes; ; is the fundamental solution of the Lamé system with for and for .
4. Meshless Formulations for Solids
From (3.9) it follows that the second integral term of (2.9) is only attributed to the principal diagonal of the matrix. This fact will be taken into account when calculating the boundary integrals. Thus substituting (3.7) and (3.9) into (2.9) and (2.10) yields where . Then applying the above equations to all boundary nodes provides the final system of equations In the case , , , , and The evaluation of the main diagonal terms of matrix involves only weak singularities, while the main diagonal terms of matrix are strongly singular ones. To avoid direct numerical integration of these terms, a uniform displacement field is assumed as and without any traction on the boundary. Substituting them into (3.9) yields with . Inserting into (4.2) leads to , where is a column vector. Consequently, the main diagonal terms of matrix can be computed by the off-diagonal terms. In the same way we tackle for .
Once the unknowns and are found, the values of the displacement and the traction at any boundary point are computed using (3.7). The displacement and the stress at an internal point may be computed simply using (3.9). Although this scheme avoids further integrations, it has the drawback of serious “boundary layer effect,’’ that is, the accuracy of the result near the boundary is very sensitive to the proximity of the interior points near the boundary. Similar to schemes proposed in [7, 18], an adaptive integration scheme is further developed to circumvent this problem.
The displacement and the stress at an internal point are evaluated via the following traditional BIEs, where , , and are the fundamental solution with and being the field point and source point, respectively; is the number of the segments which compose the whole boundary. Since every segment can be represented by a unit sector or square in the parametric space, the integrals on each segment in (4.5) and (4.6) can be computed easily. Here, an adaptive technique is developed to compute these integrals on the segment. In this technique, the unit sector is first divided into four equal quarters, as shown in Figure 2. Then, for each quarter, we compute the diagonal length, , and the distance between the evaluation point and the center of the quarter, . If , this quarter is considered as a regular integration segment. Otherwise, it will be further divided into subquarters, and this process goes on, until all segments become regular. Finally, Gaussian quadrature is applied for all segments. This adaptive scheme is very accurate even when the evaluation point is very close to the boundary. It should be pointed out that the segments are in the parametric space and are not like the boundary element in the BEM.
(a)
(b)
From the above discussion, it can be concluded that all integrations are computed along the boundary only and that no boundary elements are used both for interpolation and integration purposes. Thus, the present numerical method is truly meshless and boundary-only.
5. Numerical Experiments
In order to demonstrate the efficiency and accuracy of the present method, some numerical examples are considered and their results are compared with the analytical results. As in many meshless methods, there exist some parameters in the present method. For numerical computations, as in many works [2, 5, 7, 13, 16, 17], these parameters can be fixed. In all examples, the radius , of local integration subdomain is taken to be , while the size of interpolation domain is chosen as , where is the minimum distance between the adjacent nodes. Besides, the parameters in (3.2) are taken as and , and the parameter is taken to be . In order to deal with the traction discontinuities at corners and edges, the nodes are not arranged at these locations and the support domain for interpolation is truncated.
5.1. Internal Pressurized Hollow Cylinder
A hollow cylinder under unit internal pressure is shown in Figure 3. The radii of the inner and outer cylinders are and , respectively. The plane stress condition is assumed with the Young’s modulus and the Poisson’s ratio . In the polar coordinate system , the analytical solution for stresses is The corresponding displacements are where and are the radial and tangential component of the displacement.
Due to the symmetry of the problem, only one quarter of the cylinder is modeled. The boundary of the cylinder in the first quadrant is discretized by 60 boundary nodes (12 uniformly distributed nodes on AB, CD, and AD, and 24 uniformly distributed nodes on BC). In the computation, the geometry is chosen as and . Figure 4 plots the radial displacement versus the node position along the radius, while Figure 5 plots the stress results. It can be clearly seen that the numerical solution is in excellent agreement with the analytical ones.
5.2. Rigid Flat Punch on Semi-Infinite Foundation
As boundary type methods have a clear advantage over domain type methods for problems with infinite domain, the present truly meshless method is also employed to obtain a solution for an indentation produced by a rigid flat punch in a semi-infinite soil foundation [2], as shown in Figure 6. In this case, only the contact surface between the punch and the half-space needs to be discretized.
Consider a rigid punch of length subjected to a uniform pressure of on the top. The parameters of the soil foundation are taken as and . The punch is considered to be perfectly smooth, and there has not been any friction in the interface between the punch and the foundation. Due to symmetry, only the right half part is discretized by 31 uniformly boundary nodes.
The vertical displacement of the foundation is assumed to be the same of that of the punch. A prescribed vertical displacement of the punch is imposed on the contact surface as boundary constraints. This displacement can be expressed as , where is the vertical displacement of the rigid area in contact with the rigid punch, and and are vertical displacements at the center and edge, respectively. The analytical solution of and can be obtained in [21]. The analytical solution of the contact stress along the contact surface is Figure 7 plots the comparison between the contact stresses along the contact surface computed analytically and by the present meshless method. It can be found that excellent agreement between the analytical and numerical solutions is achieved.
5.3. Lamé Problem
The three-dimensional Lamé problem consists of a hollow sphere under internal pressure. The sketched geometric configuration of this problem is shown in Figure 8. For this benchmark problem, the analytical solutions for the radial displacement, and radial and tangential stresses are available in the polar coordinate system as where is the radial distance from the centroid of the sphere to the point of interest in the sphere, is the internal pressure, is the inner radius of sphere, and is the outer radius. The parameters are chosen as , , , , and in the computation.
Figure 9 exhibits the comparison between the present numerical results with the analytical solutions for , , and . In this analysis, only 48 regularly distributed nodes are used to discretize each surface of the hollow sphere. This figure indicates that the numerical solutions are seen to capture the behavior of the exact solutions very well.
5.4. Kirsch Problem
The three-dimensional Kirsch problem is a portion of an infinite cube with a small spherical cavity subjected to a unidirectional tensile load of in the -axis direction as shown in Figure 10. The analytical solution for the normal stress (), for points in the plane , is where .
The problem is solved her for , , and . Seventy-two uniformly spaced nodes are used on the inner spherical surface and ninety-six nodes on the outer cube boundary. Figure 11 plots a comparison between the numerical solution and the analytical solution for the normal stress along the -axis ahead of the cavity. From this figure it can be found that the numerical solution for this example is in good consistency with the analytical solution.
6. Conclusions
A novel numerical method has been developed in this paper for analysis of two- and three-dimensional solid mechanics problems. The present method inherently possesses some desirable numerical merits which include truly meshless and boundary-only. In this method, meshless shape functions constructed by the radial basis point interpolation possess delta function properties, thus the basic unknown quantities are the real solutions of the nodal variables and boundary conditions can be directly and easily implemented. Besides, it only uses a cluster of unorganized nodes on the bounding surface of the domain to be solved, and thus absolutely no discretization grids are required either for variable interpolation or for numerical integration. Some examples have been given and the numerical results are accurate.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant no. 11101454) and the Natural Science Foundation Project of CQ CSTC (Grant nos. cstcjjA30002 and cstc2011jjA30003).