Numerical Solution of Solid Mechanics Problems Using a Boundary-Only and Truly Meshless Method

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.


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 Mathematical Problems in Engineering 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 h − p 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 L s 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 boundaryonly 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.

Boundary Variational Principle for Solid Mechanics Problems
Let Ω be a bounded domain in where n n 1 , n 2 , . . ., n d is the unit outward normal on Γ.Based on the modified variational principle, the solution of problem 2.1 is the function u which minimizes the following functional 19 : where t I is the boundary traction, u I is the boundary displacement satisfying u I u I on Γ u , C IJKL 2Gνδ IJ δ KL / 1 − 2ν Gδ IL δ JK , G is the shear modulus and ν is the Poisson's ratio.We introduce the stress tensor σ and strain tensor ε where δ IJ 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 t I t I 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 Ω i s , which is bounded by Γ i s and L i s and contains the boundary node y i .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 w i is a weight function or a test function.We emphasize that in the above equations the shape and dimension of the subdomains Ω i s 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 R 2 or R 3 .Hence, the subdomain Ω i s is chosen as the intersection of the domain Ω and a circle or a sphere centered at the boundary node y i see Figure 1 .
Figure 1: Subdomain Ω i s for boundary node y i for a two-dimensional case and b three-dimensional case.

Variables Interpolation
Assume that the boundary Γ is the union of piecewise smooth segments Γ i , i 1, 2, . . ., N Γ .To avoid the discontinuity at corners and edges, the point interpolation for displacement u I and boundary traction t I on each Γ i is constructed independently.So in the following interpolation scheme, let the variable v denote u I or t I for simplicity.Then the radial basis point interpolation for v can be defined as 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 c and q are shape parameters.In this situation, the required polynomial basis is linear as p s 1, s T .
Letting the approximation represented by 3.1 pass through all the N s boundary nodes, we get v B 0 a p 0 b, 3.3 in which v v s 1 , v s 2 , . . ., v s N s T is a vector, while B 0 B T s 1 , B T s 2 , . . ., B T s N s T and p 0 p T s 1 , p T s 2 , . . ., p T s N s T are matrices.Besides, in order to guarantee unique solution, the following constraints should also be satisfied Then, according to 3.
where the shape function vector is are the nodal values u I y j and t I y j , respectively; φ j s is the contributions from the node y j to the evaluation point x s , which is not equal to zero in the interpolation domain of the jth node only.
In 2.9 and 2.10 , u I and t I on L i s are not defined yet.However, this problem can be tackled by selecting the weight function w i such that the size of the support of w i is less than the radius of the subdomain Ω i s , then all integrals over L i s 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 r i is the radius of Ω i s ; d i is the distance between a sampling point x, in domain Ω, and the nodal point y i ; and c i is a constant controlling the shape of the weight function w i .As a result, w i x vanishes on L i s .
Finally, the domain variables u I and t I in 2.9 and 2.10 are interpolated by the fundamental solution as

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 I 1, 2, . . ., d.Then applying the above equations to all boundary nodes provides the final system of equations

Mathematical Problems in Engineering
The evaluation of the main diagonal terms of matrix U involves only weak singularities, while the main diagonal terms of matrix T are strongly singular ones.To avoid direct numerical integration of these terms, a uniform displacement field is assumed as u 1  1, 0 T and u 2 0, 1 T without any traction on the boundary.Substituting them into 3.9 yields γ k U −1 H u k with k 1, 2. Inserting γ k into 4.2 leads to Tγ k H0, where 0 is a column vector.Consequently, the main diagonal terms of matrix T can be computed by the offdiagonal terms.In the same way we tackle for d 3.
Once the unknowns t and u are found, the values of the displacement u and the traction t at any boundary point are computed using 3.7 .The displacement u 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 u and the stress σ at an internal point x are evaluated via the following traditional BIEs, where u * x, y , t * x, y , and σ * x, y are the fundamental solution with y and x being the field point and source point, respectively; N Γ 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, l, and the distance between the evaluation point and the center of the quarter, d.If l < d, 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.
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 boundaryonly.

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 r i , of local integration subdomain is taken to be 0.9h, while the size of interpolation domain is chosen as 3.0h, where h is the minimum distance between the adjacent nodes.Besides, the parameters in 3.2 are taken as c 2.0h and q 0.5, and the parameter c i is taken to be 1.5h.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.

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 R 1 and R 2 , respectively.The plane stress condition is assumed with the Young's modulus E 10 and the Poisson's ratio ν 0.3.In the polar coordinate system r, θ , the analytical solution for stresses is The corresponding displacements are BC .In the computation, the geometry is chosen as R 1 1 and R 2 2. 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.

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 L 12 subjected to a uniform pressure of p 100 on the top.The parameters of the soil foundation are taken as E 3.0 × 10 4 and ν 0.3.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 u r 0.5 u c u e , where u r is the vertical displacement of the rigid area in contact with the rigid punch, and u c and u e are vertical displacements at the center and edge, respectively.The analytical solution of u c and u e 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.

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 r is the radial distance from the centroid of the sphere to the point of interest in the sphere, P is the internal pressure, a is the inner radius of sphere, and b is the outer radius.The parameters are chosen as P 1, a 2.0, b 4.0, E 1.0, and ν 0.3 in the computation.
Figure 9 exhibits the comparison between the present numerical results with the analytical solutions for u r , σ rr , 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.

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 σ 0 in the x 3 -axis direction as shown in Figure 10.The analytical solution for the normal stress σ 33 , for points in the plane x 3 0, is x 2 2 .The problem is solved her for a 1.0, b 10.0, and σ 0 1.0.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 x 1 -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.

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.

3 and 3 . 4 ,
we obtain a S a v and b S b v, where S b p T 0 B −1 0 p 0 −1 p T 0 B −1 0 and S a B −1 0 I − p 0 S b .Consequently, substituting a and b into 3.1 yields

y u y dΓ y , 4
x, y t y dΓ y − Γ t * x, y u y dΓ y N Γ i 1 Γ i −σ * x, y t y dΓ y − N Γ i 1 Γ i t * x, y u y dΓ y , 4.6

Figure 2 :
Figure 2: Subdividing of the segment Γ i at an evaluation point x for a two-dimensional case and b three-dimensional case.

Figure 5 : 2 Figure 6 :
Figure 5: Stress σ along the radius for the internal pressurized hollow cylinder.

1 Figure 7 :
Figure 7: Contact stress on the contact surface between the punch and foundation.

Figure 9 :
Figure 9: Results of u r , σ rr , and σ θθ for the three-dimensional Lamé problem.

Figure 10 :
Figure 10: Schematic diagram for the three-dimensional Kirsch problem.