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

Authors: Xiaolin Li
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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, theh-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 Ls 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 Section2. 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 ℝd(d=2,3) with boundary Γ=Γu+Γt. A general point in ℝd is denoted by x=(x1,x2,…,xd) or y=(y1,y2,…,yd).Consider the following two- and three-dimensional problem of solid mechanics(2.1)μΔu+(λ+μ)∇(∇⋅u)=0,inΩ,u=u¯,on⁡Γu,t:=Tu=t¯,on⁡Γt,
where u=(u1,u2,…,ud) is the unknown displacement field; μ>0 and λ>-μ are given Lamé constants; Δ, ∇ and ∇· stand for the Laplacian, gradient, and divergence operators, respectively; u¯=(u¯1,u¯2,…,u¯d) and t¯=(t¯1,t¯2,…,t¯d) are prescribed boundary displacement and traction, respectively; and T is the following conformal derivative operator:(2.2)Tu=λ(∇⋅u)⋅n+2μ∂u∂n+μn×curl⁡u,
where n=(n1,n2,…,nd) 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]:(2.3)Π=12∫ΩuI,JCIJKLuK,LdΩ-∫Γt̃I(uI-ũI)dΓ-∫Γtt¯IũIdΓ,
where t̃I is the boundary traction, ũI is the boundary displacement satisfying ũI=u¯I on Γu, CIJKL=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 ɛ(2.4)σIJ(u)=λδIJ∑K=1dɛKK(u)+2μɛIJ(u),ɛIJ(u)=12(∂uI∂xJ+∂uJ∂xI),I,J=1,2,…,d,
where δIJ is the Kronecker delta symbol. Then, via carrying out the first variation of (2.3) one gets(2.5)δΠ=-∫ΩσIJ,JδuIdΩ+∫Γ(tI-t̃I)δuIdΓ-∫Γ(uI-ũI)δt̃IdΓ-∫Γt(t̃I-t¯I)δũIdΓ.
Thus, by setting δΠ to zero we obtain the following equivalent integral equations(2.6)∫Γ(tI-t̃I)δuIdΓ-∫ΩσIJ,JδuIdΩ=0,(2.7)∫Γ(uI-ũI)δt̃IdΓ=0,(2.8)∫Γt(t̃I-t¯I)δũIdΓ=0.
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 Ωsi, which is bounded by Γsi and Lsi and contains the boundary node yi. 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(2.9)∫Γsi+Lsi(tI-t̃I)widΓ-∫ΩsiσIJ,JwidΩ=0,(2.10)∫Γsi+Lsi(uI-ũI)widΓ=0,
where wi is a weight function or a test function. We emphasize that in the above equations the shape and dimension of the subdomains Ωsi 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 ℝ2 (or ℝ3). Hence, the subdomain Ωsi is chosen as the intersection of the domain Ω and a circle (or a sphere) centered at the boundary node yi (see Figure 1).SubdomainΩsi for boundary node yi for (a) two-dimensional case and (b) three-dimensional case.
(a)(b)


## 3. 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 uI and boundary traction tI 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(3.1)v(s)≈vh(s)=∑j=1NsajBj(s)+∑i=1mbipi(s)=BT(s)a+pT(s)b,
where s is a parametric coordinate on Γi; Bj(s) is a radial basis function (RBF); Ns is the node number of the interpolation domain; pi(s) is a monomials; m is the number of monomials; B(s)=[B1(s),B2(s),…,BNs(s)]T; a=[a1,a2,…,aNs]T; b=[b1,b2,…,bm]T; p(s)=[p1(s),p2(s),…,pm(s)]T; aj and bi 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)(3.2)Bj(s)=[|s-sj|2+c2]q,
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 Ns boundary nodes, we get(3.3)v̂=B0a+p0b,
in which v̂=[v(s1),v(s2),…,v(sNs)]T is a vector, while B0=[BT(s1),BT(s2),…,BT(sNs)]T and p0=[pT(s1),pT(s2),…,pT(sNs)]T are matrices. Besides, in order to guarantee unique solution, the following constraints should also be satisfied(3.4)p0Ta=0.
Then, according to (3.3) and (3.4), we obtain a=Sav̂ and b=Sbv̂, where Sb=[p0TB0-1p0]-1p0TB0-1 and Sa=B0-1[I-p0Sb]. Consequently, substituting a and b into (3.1) yields(3.5)vh(s)=[BT(s)Sa+pT(s)Sb]v̂=Φ(s)v̂,
where the shape function vector is(3.6)Φ(s)=BT(s)Sa+pT(s)Sb:=[ϕ1(s),ϕ2(s),…,ϕi(s),…,ϕNs(s)].Now, in (2.9) and (2.10), ũI and t̃I on Γsi can be represented as(3.7)ũI(s)=∑j=1Nsϕj(s)uIj,t̃I(s)=∑j=1Nsϕj(s)tIj,
where uIj and tIj are the nodal values uI(yj) and tI(yj), respectively; ϕj(s) is the contributions from the node yj 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 Lsi are not defined yet. However, this problem can be tackled by selecting the weight function wi such that the size of the support of wi is less than the radius of the subdomain Ωsi, then all integrals over Lsi 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(3.8)wi(x)={exp⁡[-(di/ci)2]-exp⁡[-(ri/ci)2]1-exp⁡[-(ri/ci)2],0≤di≤ri,0,di≥ri,
where ri is the radius of Ωsi; di is the distance between a sampling point x, in domain Ω, and the nodal point yi; and ci is a constant controlling the shape of the weight function wi. As a result, wi(x) vanishes on Lsi.Finally, the domain variablesuI and tI in (2.9) and (2.10) are interpolated by the fundamental solution as(3.9)uI(x)=∑j=1N∑K=1duIK*(x,yj)γKj,tI(x)=∑j=1N∑K=1dtIK*(x,yj)γKj,
where γKj is the unknown parameter; N is the total number of boundary nodes; tIK*=TuIK*; uIK* is the fundamental solution of the Lamé system(3.10)uIK*(x,y)=λ+3μ4π(d-1)μ(λ+2μ)[γd(x,y)δIK+λ+μλ+3μ(xI-yI)(xK-yK)T|x-y|d],
with γd(x,y)=-ln⁡|x-y| for d=2 and γd(x,y)=1/|x-y| for d=3.


## 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(4.1)∑j=1N∫Γsi(∑K=1dtIK*(x,yj)γKj-ϕj(x(s))tIj)wi(x)dΓ=0,∑j=1N∫Γsi(∑K=1duIK*(x,yj)γKj-ϕj(x(s))uIj)wi(x)dΓ=0,
where I=1,2,…,d. Then applying the above equations to all boundary nodes provides the final system of equations(4.2)Tγ=Ht̃,(4.3)Uγ=Hũ.
In the case d=2, γ=[γ11,γ21,γ12,γ22,…,γ1N,γ2N]T, t̃=[t11,t21,t12,t22,…,t1N,t2N]T, ũ=[u11,u21,u12,u22,…,u1N,u2N]T, and
(4.4)TIJ=∫Γsi[t11*(x,yj)t12*(x,yj)t21*(x,yj)t22*(x,yj)]wi(x)dΓ,UIJ=∫Γsi[u11*(x,yj)u12*(x,yj)u21*(x,yj)u22*(x,yj)]wi(x)dΓ,HIJ=∫Γsi[ϕj(x(s))00ϕj(x(s))]wi(x)dΓ.
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-1Hũ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 off-diagonal terms. In the same way we tackle for d=3.Once the unknownst̃ 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 displacementu and the stress σ at an internal point x are evaluated via the following traditional BIEs,(4.5)u(x)=∫Γu*(x,y)t̃(y)dΓy-∫Γt*(x,y)ũ(y)dΓy=∑i=1NΓ∫Γiu*(x,y)t̃(y)dΓy-∑i=1NΓ∫Γit*(x,y)ũ(y)dΓy,(4.6)σ(x)=∫Γ-σ*(x,y)t̃(y)dΓy-∫Γt*(x,y)ũ(y)dΓy=∑i=1NΓ∫Γi-σ*(x,y)t̃(y)dΓy-∑i=1NΓ∫Γit*(x,y)ũ(y)dΓy,
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.Subdividing of the segmentΓi at an evaluation point x for (a) two-dimensional case and (b) three-dimensional case.
(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 ri, 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 ci 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.


### 5.1. Internal Pressurized Hollow Cylinder

A hollow cylinder under unit internal pressure is shown in Figure3. The radii of the inner and outer cylinders are R1 and R2, 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(5.1)σrr(r,θ)=R12R22-R12(1-R22r2),σθθ(r,θ)=R12R22-R12(1+R22r2),σrθ(r,θ)=0.
The corresponding displacements are(5.2)ur(r,θ)=R12E(R22-R12)[(1-ν)r+(1+ν)R22r],uθ(r,θ)=0,
where ur and uθ are the radial and tangential component of the displacement.Figure 3
Schematic diagram for the internal pressurized hollow cylinder.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 asR1=1 and R2=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.Figure 4
Radial displacementur along the radius for the internal pressurized hollow cylinder.Figure 5
Stressσ along the radius for the internal pressurized hollow cylinder.


### 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.Figure 6
Schematic diagram for a rigid flat punch on a semi-infinite foundation.Consider a rigid punch of lengthL=12 subjected to a uniform pressure of p=100 on the top. The parameters of the soil foundation are taken as E=3.0×104 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 asur=0.5(uc+ue), where ur is the vertical displacement of the rigid area in contact with the rigid punch, and uc and ue are vertical displacements at the center and edge, respectively. The analytical solution of uc and ue can be obtained in [21]. The analytical solution of the contact stress along the contact surface is(5.3)σx2=2pL2πL2-4x12.
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.Figure 7
Contact stress on the contact surface between the punch and foundation.


### 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 Figure8. For this benchmark problem, the analytical solutions for the radial displacement, and radial and tangential stresses are available in the polar coordinate system as(5.4)ur=Pa3rE(b3-a3)[(1-2ν)+(1+ν)b32r3],σrr=Pa3(b3-r3)r3(a3-b3),σθθ=Pa3(b3+2r3)2r3(b3-a3),
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 8
Schematic diagram for the three-dimensional Lamé problem.Figure9 exhibits the comparison between the present numerical results with the analytical solutions for ur, σ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.Figure 9
Results ofur, σrr, and σθθ for the three-dimensional Lamé problem.


### 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σ0 in the x3-axis direction as shown in Figure 10. The analytical solution for the normal stress (σ33), for points in the plane x3=0, is(5.5)σ33=σ0[1+4-5ν2(7-5ν)(ar)3+92(7-5ν)(ar)5],
where r=x12+x22.Figure 10
Schematic diagram for the three-dimensional Kirsch problem.The problem is solved her fora=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 x1-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.Figure 11
Stressσ33 for the three-dimensional Kirsch problem along the x1-axis.


## 5.1. Internal Pressurized Hollow Cylinder

A hollow cylinder under unit internal pressure is shown in Figure3. The radii of the inner and outer cylinders are R1 and R2, 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(5.1)σrr(r,θ)=R12R22-R12(1-R22r2),σθθ(r,θ)=R12R22-R12(1+R22r2),σrθ(r,θ)=0.
The corresponding displacements are(5.2)ur(r,θ)=R12E(R22-R12)[(1-ν)r+(1+ν)R22r],uθ(r,θ)=0,
where ur and uθ are the radial and tangential component of the displacement.Figure 3
Schematic diagram for the internal pressurized hollow cylinder.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 asR1=1 and R2=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.Figure 4
Radial displacementur along the radius for the internal pressurized hollow cylinder.Figure 5
Stressσ along the radius for the internal pressurized hollow cylinder.


## 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.Figure 6
Schematic diagram for a rigid flat punch on a semi-infinite foundation.Consider a rigid punch of lengthL=12 subjected to a uniform pressure of p=100 on the top. The parameters of the soil foundation are taken as E=3.0×104 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 asur=0.5(uc+ue), where ur is the vertical displacement of the rigid area in contact with the rigid punch, and uc and ue are vertical displacements at the center and edge, respectively. The analytical solution of uc and ue can be obtained in [21]. The analytical solution of the contact stress along the contact surface is(5.3)σx2=2pL2πL2-4x12.
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.Figure 7
Contact stress on the contact surface between the punch and foundation.


## 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 Figure8. For this benchmark problem, the analytical solutions for the radial displacement, and radial and tangential stresses are available in the polar coordinate system as(5.4)ur=Pa3rE(b3-a3)[(1-2ν)+(1+ν)b32r3],σrr=Pa3(b3-r3)r3(a3-b3),σθθ=Pa3(b3+2r3)2r3(b3-a3),
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 8
Schematic diagram for the three-dimensional Lamé problem.Figure9 exhibits the comparison between the present numerical results with the analytical solutions for ur, σ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.Figure 9
Results ofur, σrr, and σθθ for the three-dimensional Lamé problem.


## 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σ0 in the x3-axis direction as shown in Figure 10. The analytical solution for the normal stress (σ33), for points in the plane x3=0, is(5.5)σ33=σ0[1+4-5ν2(7-5ν)(ar)3+92(7-5ν)(ar)5],
where r=x12+x22.Figure 10
Schematic diagram for the three-dimensional Kirsch problem.The problem is solved her fora=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 x1-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.Figure 11
Stressσ33 for the three-dimensional Kirsch problem along the x1-axis.


## 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.



---

*Source: 298903-2012-05-30.xml*