A priori and a posteriori analysis of the meshless Galerkin boundary node method for three-dimensional elasticity

Engineering Analysis with Boundary Elements 36 (2012) 993–1004

Contents lists available at SciVerse ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A priori and a posteriori analysis of the meshless Galerkin boundary node method for three-dimensional elasticity Xiaolin Li n College of Mathematics Science, Chongqing Normal University, Chongqing 400047, PR China

a r t i c l e i n f o

abstract

Article history: Received 11 July 2011 Accepted 24 December 2011 Available online 1 February 2012

The meshless Galerkin boundary node method is presented in this paper for boundary-only analysis of three-dimensional elasticity problems. In this method, boundary conditions can be implemented directly and easily despite the employed moving least-squares shape functions lack the delta function property, and the resulting system matrices are symmetric and positive definite. A priori error estimates and the consequent rate of convergence are presented. A posteriori error estimates are also provided. Reliable and efficient error estimators and an efficient and convergent adaptive meshless algorithm are then derived. Numerical examples showing the efficiency of the method, confirming the theoretical properties of the error estimates, and illustrating the capability of the adaptive algorithm, are reported. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Meshless method Galerkin boundary node method Elasticity A priori error estimates A posteriori error estimates Adaptive algorithm

1. Introduction Meshless (or meshfree) methods to obtain numerical solutions for partial differential equations without resorting to an element frame have attracted much attention and gained great success in the field of computational science and engineering in the past few decades [1,2]. Compared to the finite element method (FEM) and the boundary element method (BEM), the core of this type of method is to alleviate the difficulty of meshing and remeshing the entire structure by simply adding or deleting nodes. Although many kinds of meshless methods have been proposed, these methods can be simply divided into the domain-type and the boundary-type. Several domain-type meshless methods, such as the element free Galerkin (EFG) method [1–3], the reproducing kernel particle method (RKPM) [1], the point interpolation method (PIM) [2], the generalized FEM [4,5], the h–p meshless method [1], the smoothed FEM [6] and the finite point method (FPM) [7] have been proposed and gained noticeable progress in solving a wide range of boundary value problems. Boundary integral equations (BIEs) are attractive computational techniques for linear and exterior problems as they can reduce the dimensionality of the original problem by one. Boundary-type meshless methods are developed by the combination of the meshless idea with BIEs, such as the boundary node method (BNM) [8,9], the hybrid BNM [10–13] and the boundary

n

Tel.: þ86 13527466263. E-mail address: [email protected]

0955-7997/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.12.015

face method [14]. In the three methods, the moving least-squares (MLS) approximations are used to generate the shape functions on the boundary of a domain. However, because MLS shape functions lack the delta function property, boundary conditions in these meshless methods cannot be implemented as easily as in the FEM or the BEM. The strategy employed in them to impose boundary conditions involves a new definition of the discrete norm used for the construction of the MLS approximations, which adds to the number of system equations. In order to construct meshless shape functions with delta function properties, Liu [2] has introduced the PIM into BIEs and produced boundary PIMs. Besides, Li et al. [15] have introduced the radial basis point interpolation into the hybrid displacement variational formulation and produced the hybrid radial BNM. Moreover, Peng and Cheng [16] have developed an improved MLS approximation that uses weighted orthogonal polynomials as basis functions to restore the delta function property of the MLS. The improved MLS scheme has been inserted into BIEs to propose a boundary-type meshless method called the boundary element-free method [16]. Li and Zhu [17] have recently proposed a boundary-type meshless method called the Galerkin boundary node method (GBNM). The GBNM combines a variational version of BIEs for governing equations with the MLS approximations for generation of the trial and test functions. Unlike other MLS-based methods mentioned above, boundary conditions in the GBNM do not present any difficulty and can be implemented with ease via multiplying the MLS shape function and integrating on the boundary. Another outstanding feature of the GBNM is the conservation of the symmetry and positive definiteness of

994

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

the variational problems in the process of numerical implementation. This method has been successfully tried for 2D problems in potential theory [17,18], elasticity [19] and fluid mechanics [20,21]. Very recently, the method has been extended to solve 3D problems in potential theory [22] and fluid mechanics [23]. The present paper extends the frontiers of the GBNM into solving problems in 3D elasticity. A priori error estimates, which ensure convergence of numerical methods, are crucial in meshless research. The associated mathematical proofs guarantee that meshless methods will converge to the true solution. During the past two decades, a large amount of research has been devoted to deriving a priori error estimation for domain-type meshless methods such as the h–p meshless method [1], the RKPM [1], the EFG [3], the generalized FEM [4], the smoothed FEM [2,6] and the FPM [7]. Nevertheless, although boundary-type meshless methods perform very well in practice, not much is rigorously known on the a priori error analysis of these schemes. Until now, a first a priori error analysis of the boundary-type meshless methods was given by Li and Zhu for the GBNM for 2D problems [17–21] and for 3D problems in potential theory [22] and fluid mechanics [23]. Hence, one goal of this paper is to establish a rigorous a priori error analysis of the GBNM for 3D elasticity problems. The optimal asymptotic convergence rates are given in Sobolev spaces. In meshless methods, since no predefined nodal connectivity or mesh is employed for field variable approximation, nodes can be inserted or removed conveniently. This prominent feature makes meshless methods especially suited for self-adaptive techniques. Indeed, the subject of a posteriori error estimates and corresponding adaptive procedures is central to the effective application of meshless methods for practical engineering computation. In recent years, a large amount of work has been performed concerning adaptive analysis based on a posteriori error estimation for domain-type meshless methods such as the h–p meshless method [24], the EFG method [2,25], the RKPM [26], the FPM [27] and the PIM [28]. Significant progress has been achieved in the theory and implementation of the adaptive procedures for these meshless methods. For boundary-type meshless methods, Chati et al. [9,29] have pioneered error indicators and an adaptive algorithm for the BNM using hypersingular residual techniques similar to those used in the BEM [9]. Besides, Guo and Chen [30] have developed an adaptive algorithm for the meshless local BIE method based on the dual error indicators. Very recently, the GBNM has been extended for a posteriori error estimate and adaptivity for 2D potential problems based on the difference between numerical solutions obtained using two successive nodal arrangements [31]. Another aim of this paper is to extend the a posteriori error results to 3D elasticity. The formulation of an accurate and reliable a posteriori error estimate is presented. Then, a posteriori error estimator for the error control of numerical solutions is derived. This error estimator has an upper and a lower bound by the constant multiples of the exact error in the energy norm. That is, this estimator is reliable and efficient. Finally, based on the a posteriori error estimation and a localization technique, computable local error indicators and an efficient and convergent adaptive meshless algorithm for h-adaptivity are established. The rest of this paper is organized as follows. In Section 2, detailed formulations of the GBNM for 3D elasticity problems are provided. Section 3 contains a priori error analysis and numerical examples showing the performance of the GBNM. In Section 4, a posteriori error analysis and adaptive algorithm are given. Numerical examples illustrating the capability of the adaptive algorithm are also presented in this section. Finally, conclusions are summarized in Section 5.

2. The GBNM for 3D elasticity Consider the following 3D elasticity problem: (

rr ¼ 0

in O on G

u9G ¼ u

ð1Þ

where r is the gradient operator, r is the stress tensor, O is a bounded or unbounded domain in R3 with boundary surface G, u ¼ ðu1 ,u2 ,u3 Þ is the displacement field, and u ¼ ðu 1 ,u 2 ,u 3 Þ is a given function on G. A general point of O is denoted by x ¼ ðx1 ,x2 ,x3 Þ. If G is a smooth open surface piece with a piecewise smooth boundary curve and O ¼ R3 \G , then problem (1) is the socalled screen crack problem. Let u be given satisfying u A H1=2 ðGÞ :¼ ðH1=2 ðGÞÞ3 , then problem (1) has a unique solution u. The solution can be expressed as [32,33] uj ðyÞ ¼

3 Z X i¼1

G

qi ðxÞU ij ðx,yÞ dSx ,

j ¼ 1; 2,3,

yAO

ð2Þ

in which q ¼ ðq1 ,q2 ,q3 Þ is the jump of the boundary traction n  r across G, n ¼ ðn1 ,n2 ,n3 Þ is the normal exterior to G, and U ij is the singular Kelvin fundamental solution U ij ðx,yÞ :¼

1 ½ð34nÞdij þ r ,i r ,j , 16pGð1nÞr

i,j ¼ 1; 2,3

ð3Þ

where G is the shear modulus, n is the Poisson ratio, dij is the Kronecker symbol, r ¼ 9xy9 and r ,i ¼ ðxi yi Þ=r. Eq. (2) gives the indirect integral equations of 3D elasticity. By direct method, we can also get the direct integral equations and dual integral [34,35]. Using Eq. (2) and the boundary condition of problem (1) leads to the following BIEs: ðAqÞj ðyÞ ¼

3 Z X

G

i¼1

qi ðxÞU ij ðx,yÞ dSx ¼ u j ðyÞ,

j ¼ 1; 2,3,

yAG

ð4Þ

which are suitable for the solution of the exterior as well as the interior problem. Here the boundary integral operator A : H1=2 ðGÞ-H1=2 ðGÞ is continuous and bijective. These BIEs have a unique solution in H1=2 ðGÞ and admit the variational problem: Find q A H1=2 ðGÞ such that /Aq,q0 SL2 ðGÞ ¼ /u,q0 SL2 ðGÞ ,

8q0 ¼ ðq01 ,q02 ,q03 Þ A H1=2 ðGÞ

ð5Þ

where /Aq,q0 SL2 ðGÞ :¼

3 Z Z X i,j ¼ 1

/u,q0 SL2 ðGÞ :¼

3 Z X j¼1

G

G G

qi ðxÞU ij ðx,yÞq0j ðyÞ dSx dSy

u j ðyÞq0j ðyÞ dSy

ð6Þ

To carry out integrations in the variational problem (5), the boundary surface G is discretized into cells. Assume that G is piecewise smooth and can be partitioned into finitely many smooth pieces Uk . Each Uk can be considered to be the image of Gk by a smooth bijection jk , i.e., Uk ¼ jk ðGk Þ, where Gk is a bounded polygonal domain in R2 . Let T kh be a triangulation of Gk by triangles. On each triangulation T kh , we construct the b-degree interpolant function of jk denoted by jkh . Then the image of Gk by the mapping jkh constitutes one piece Ukh of the surface Gh which we take as the background cell. Consequently, Gh is a connected parametric surface of degree b. After triangulation, we assume that Gh contains N cells Gi with an associated triangulation T h :¼ fG1 , G2 , . . . , GN g, where the parameter h denotes the nodal spacing.

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

In general, Gh is an approximation of G. However, when G is a polygonal surface, or, more generally, when G can be expressed in parametric form exactly as the boundary representation data structure in standard solid modeling packages [14], Gh and G can coincide with each other. It should be pointed out that cells in the GBNM are only used so that integrals in the variational problem (5) can be calculated numerically, and do not need to satisfy any aspect ratio or compatibility requirements such as in the BEM or the FEM. This salient feature makes the GBNM especially suited for adaptive procedures. In meshless approaches, the number of nodes corresponding to each cell is arbitrary. Nevertheless, to carry out accurate integration via Gaussian quadrature, it is recommended to introduce a small number of nodes per cell [8,9,23,24]. It has been shown that one node per cell and the location of the node at the centroid of the cell can yield excellent results [9,23,29]. This nodal placement is employed in this research. More specifically, one boundary node xi is selected for each cell Gi and is located at its centroid. Let Q h ¼ fxi gN i ¼ 1 denote the set of these boundary nodes, we then define the approximate spaces as W h ðGh Þ :¼ spanfFi ,1 r ir Ng,

Vh ðGh Þ :¼ ðW h ðGh ÞÞ3

ð7Þ

where Fi is the MLS shape function generated by Q h only [22]. Assume that Fi is g-times continuously differentiable, then Fi has minðg, bÞ compact support, i.e., Fi ðxÞ A C 0 ðRi Þ, where Ri  Gh is the influence domain of xi [22]. To compare a function defined on G with a function defined on Gh , we define a mapping C of G onto Gh . For x A G, CðxÞ is defined as the orthogonal projection of x onto Gh . Obviously, C is the identical mapping when Gh coincides with G. Then for any x,y A G, we have [32,33] C 1 9CðxÞCðyÞ9 r 9xy9 r C 2 9CðxÞCðyÞ9 2

2

99xy9 9CðxÞCðyÞ9 9 rCh

bþ1

ð8Þ 2

9CðxÞCðyÞ9

bþ1

9xCðxÞ9 r Ch

91Jð@CÞðxÞ9r Ch

bþ1

ð9Þ

3 Z X j¼1

Gh

Gh Gh

ð13Þ

b 9/Aq~ h , q~ 0h SL2 ðGÞ /Ah qh ,q0h SL2 ðGh Þ 9r Ch Jq~ h JH1=2 ðGÞ Jq~ 0h JH1=2 ðGÞ

Proof. According to Eqs. (8)–(10), it can be proved that 9U ij ðx,yÞU ij ðCðxÞ, CðyÞÞ9 rCh

mk

JuJHm ðGÞ ,

0 r kr m r b þ1

where u~ h ðxÞ ¼ ðu h JC1 ÞðxÞ  Jð@CÞðxÞ.

bþ1

1

9xy9

,

i,j ¼ 1; 2,3

Besides, applying Eq. (11) yields 9Jð@CÞðxÞ9 r C

91Jð@CÞðxÞ  Jð@CÞðyÞ9 rCh

and

bþ1

Therefore 9U ij ðx,yÞU ij ðCðxÞ, CðyÞÞ  Jð@CÞðxÞJð@CÞðyÞ9 r9U ij ðx,yÞ ½1Jð@CÞðxÞ  Jð@CÞðyÞ9 þ 9½U ij ðx,yÞU ij ðCðxÞ, CðyÞÞ Jð@CÞðxÞJð@CÞðyÞ9r Ch

bþ1

1

9xy9

,

i,j ¼ 1; 2,3

Consequently, via using Eq. (6) and changing the integrals in Eq. (13) into integrals on G by the mapping C, we obtain 9/Aq~ h , q~ 0h SL2 ðGÞ /Ah qh ,q0h SL2 ðGh Þ 9  X  3 Z Z q~ ih ðxÞ½U ij ðx,yÞU ij ðCðxÞ, CðyÞÞ ¼  i,j ¼ 1 G G    0 ~ Jð@CÞðxÞJð@CÞðyÞq jh ðyÞ dSx dSy   3 Z Z X

bþ1

r Ch

i,j ¼ 1 G G

bþ1

r Ch

9q~ ih ðxÞ99xy9

1

9q~ 0jh ðyÞ9 dSx dSy

b þ 1=2

Jq~ h JH0 ðGÞ Jq~ 0h JH0 ðGÞ r Ch

Jq~ h JH0 ðGÞ Jq~ 0h JH1=2 ðGÞ

ð15Þ

where we have used the inverse estimate Jq~ 0h JH0 ðGÞ r 1=2 Ch Jq~ 0h JH1=2 ðGÞ [22]. Finally, using an analogous inverse esti1=2 mate Jq~ h JH0 ðGÞ r Ch Jq~ h JH1=2 ðGÞ yields the desired result. &

b

/Ah qh ,qh SL2 ðGh Þ Z /Aq~ h , q~ h SL2 ðGÞ Ch Jq~ h J2H1=2 ðGÞ b

¼ Cð1h ÞJq~ h J2H1=2 ðGÞ

ð16Þ

Thus by the Lax–Milgram theorem, the proof is complete.

&

Now, the discrete equations of the GBNM for 3D elasticity can be deduced using the variational problem (12). On Vh ðGh Þ, the Galerkin approximation qh ¼ ðq1h ,q2h ,q3h Þ of the real solution q may be read as N X

qih ðxÞ ¼

Fm ðxÞqðmÞ , i ¼ 1; 2,3 i

ð17Þ

m¼1

Substituting Eq. (17) into the variational problem (12), one gets a 3Nn3N linear system 918 ðmÞ 9 8 1 9 08 11 13 ½akm  ½a12 > fq1 g > > > ff k g > > > > > km  ½akm  > = > =C> = < < < B ðmÞ 2 22 23 B ½a21 C  ½a  ½a  g fq ¼ ff k g ð18Þ km km 2 @> km A > > > > > > ; > ; > ; : fqðmÞ g > : ½a31  ½a32  ½a33  > : ff 3 g > km

in which u h ¼ ðu 1h ,u 2h ,u 3h Þ is an approximation of u restricted on Gh . In practical computations, u h can be determined as an interpolation of u with basis functions of piecewise polynomial degree b. If u A Hm ðGÞ, then we have [33] Juu~ h JHk ðGÞ r Ch

and

Proof. According to Lemma 2.1 and the coerciveness of the bilinear form /A  ,  SL2 ðGÞ , we have

8q0h ¼ ðq01h ,q02h ,q03h Þ A Vh ðGh Þ

u jh ðyÞq0jh ðyÞ dSy

q~ h ¼ qh JC1

ð11Þ

where the operator Ah is an approximation of A restricted on Gh , and 3 Z Z X /Ah qh ,q0h SL2 ðGh Þ :¼ qih ðxÞU ij ðx,yÞq0jh ðyÞ dSx dSy

/u h ,q0h SL2 ðGh Þ :¼

let

Theorem 2.1. The variational problem (12) has a unique solution.

ð12Þ

i,j ¼ 1

qh ,q0h A Vh ðGh Þ,

ð10Þ

where C 1 , C 2 and C are constants independent of h, and J is the Jacobian by the mapping C. In what follows, by C we will denote a general constant which may have different values at different occurrences. On Vh ðGh Þ, we may consider, instead of variational problem (5), an approximate variational problem to find qh ¼ ðq1h , q2h ,q3h Þ A Vh ðGh Þ such that /Ah qh ,q0h SL2 ðGh Þ ¼ /u h ,q0h SL2 ðGh Þ ,

Lemma 2.1. For any q~ 0h ¼ q0h JC1 , then

995

km

km

3

in which k,m ¼ 1; 2, . . . ,N, and Z Z aijkm ¼ Fm ðxÞU ij ðx,yÞFk ðyÞ dSx dSy ,

k

i,j ¼ 1; 2,3

ð19Þ

Rk Rm

ð14Þ j

fk ¼

Z Rk

u jh ðyÞFk ðyÞ dSy ,

j ¼ 1; 2,3

ð20Þ

996

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

In Eq. (19), as x-y the kernel U ij ðx,yÞ becomes weakly singular. The regularization procedure suggested in Ref. [23] is very attractive and is employed in this research. From Eq. (19) we can observe that ½aijkm  ¼ ½ajikm  and each submatrix ½aijkm  is symmetric, hence there are only 3NðN þ1Þ entries aijkm in the coefficient matrix that need to be computed. Therefore, the coefficient matrix of the system Eq. (18) is not only symmetric but also positive definite. Once the unknowns qiðmÞ are found from Eq. (18), the approximate solution uh of the problem (1) can be computed from an approximation form of the formula (2), i.e. 3 Z X ujh ðyÞ ¼ qih ðxÞU ij ðx,yÞ dSx i¼1

Gh

3 X N X

¼

qðmÞ i

i¼1m¼1

Z

On the other hand, with the operator A defined by Eq. (4) one gets Jqq~ h JHg2 ðGÞ ¼

/qq~ h ,gSL2 ðGÞ

sup rC

sup

The Galerkin scheme makes use of an additional integration over the boundary, which results in a heavy computational overload. However, the symmetry of the system matrix leads to substantial reductions in both memory space and solution time required by the method. Besides, the system matrix is not only symmetric but also positive definite, a property that makes the method an ideal choice for coupling the FEM. The positive definiteness also enables the use of more efficient equation solvers.

Set x ¼ A1 g, then according to the coerciveness of /A  ,  SL2 ðGÞ and Lemma 2.1 there holds /Aðqq~ h Þ, xSL2 ðGÞ ¼ /Aðqq~ h Þ, xSh xSL2 ðGÞ þ/Aq,Sh xSL2 ðGÞ /Aq~ h ,Sh xSL2 ðGÞ ¼ /Aðqq~ h Þ, xSh xSL2 ðGÞ þ/uu~ h ,Sh xSL2 ðGÞ

bþ1

þCJuu~ h JHg1 ðGÞ JSh xJHg þ 1 ðGÞ þ Ch

3. A priori error estimates

JxSh xJH1=2 ðGÞ rCh

g þ 3=2

b þ 1=2

JqJHm þ 1 ðGÞ þh

bþ1

JqJHm þ 1 ðGÞ þ h

ðJqJH0 ðGÞ þJuJHb þ 1 ðGÞ Þg

Proof. Let VðGÞ be the image of Vh ðGh Þ by the inverse mapping C1 and let Sh be the L2-projection onto VðGÞ, then for any q A Hm þ 1 ðGÞ we have [22] ð22Þ

JqJHm þ 1 ðGÞ

Besides, since Sh q A VðGÞ and Sh qJC A Vh ðGh Þ, we have

Therefore, inserting Eqs. (26) and (27) into Eq. (25) provides g þ 3=2

Jqq~ h JH1=2 ðGÞ

ðJuJHb þ 1 ðGÞ þ Jq~ h JH0 ðGÞ ÞgJxJHg þ 1 ðGÞ

Theorem 3.2. For any xA O with dx :¼ miny A G 9xy9 Z d 40, the following estimate holds for 9k9 Z0 and 0 r m r g r b 9@k uðxÞ@k uh ðxÞ9 rC

gX þ2

mþgþ3

ðdx Þj9k91 fh

bþ1

þh

ðJqJH0 ðGÞ þ JuJHb þ 1 ðGÞ Þg

ð29Þ

Proof. From Eqs. (2) and (21) it follows that   3 Z X    9uj ðxÞujh ðxÞ9 ¼  ðqi ðyÞU ij ðx,yÞq~ ih ðyÞU ij ðx, CðyÞÞ  Jð@CÞðyÞÞ dSy    G i¼1

rC

3 X

fJqi ðyÞq~ ih ðyÞJHg2 ðGÞ JU ij ðx,yÞJHg þ 2 ðGÞ

i¼1

þJU ij ðx,yÞU ij ðx, CðyÞÞJL2 ðGÞ Jq~ ih ðyÞJL2 ðGÞ g þJU ij ðx, CðyÞÞð1Jð@CÞðyÞÞJL2 ðGÞ Jq~ ih ðyÞJL2 ðGÞ g

gX þ2

ðdx Þj1

ð31Þ

Besides, using Eqs. (8)–(10), it can be proved that

Then according to Eqs. (15) and (16), and the continuity of the bilinear form /A  ,  SL2 ðGÞ , we obtain

bþ1

9U ij ðx,yÞU ij ðx, CðyÞÞ9 r Ch

9xy9

2

rCh

9U ij ðx, CðyÞÞ  ð1Jð@CÞðxÞÞ9 rCh Jqq~ h JH1=2 ðGÞ r JqSh qJH1=2 ðGÞ þ JSh qq~ h JH1=2 ðGÞ JqJHm þ 1 ðGÞ þ h

JqJH0 ðGÞ g

ðdx Þ2 ,

i,j ¼ 1; 2,3

Moreover, using Eq. (11) leads to

Therefore, using triangle inequality, Eqs. (14) and (22) yields

b þ 1=2

bþ1

ð32Þ

b þ 1=2 Jq~ h Sh qJH1=2 ðGÞ r CfJuu~ h JH1=2 ðGÞ þ JqSh qJH1=2 ðGÞ þ h JSh qJH0 ðGÞ g

JuJHb þ 1 ðGÞ þ h

ð30Þ

j¼1

/Ah ðSh qJCÞ,qh Sh qJCSL2 ðGh Þ

m þ 3=2

JqJHm þ 1 ðGÞ

j¼0

JU ij ðx,yÞJHg þ 2 ðGÞ rC

þ /AðqSh qÞ, q~ h Sh qSL2 ðGÞ þ /AðSh qÞ, q~ h Sh qSL2 ðGÞ

r Cfh

ð28Þ

Because of dx Z d 40 we get

/Ah ðqh Sh qJCÞ,qh Sh qJCSL2 ðGh Þ ¼ /u~ h u, q~ h Sh qSL2 ðGÞ

b þ 1=2

ð27Þ

ðJqJH0 ðGÞ þJuJHb þ 1 ðGÞ Þg

where 0 rm r g r b and q~ h ¼ qh JC1 is the counterpart of qh defined on G.

m þ 3=2

ð26Þ

Finally, by substituting Eq. (28) into Eq. (24) and using Eq. (23) we can prove the desired estimates. &

Theorem 3.1. Let q and qh be the solutions of variational problems (5) and (12), respectively. Assume that q A Hm þ 1 ðGÞ and u A Hb þ 1 ðGÞ, then

JqSh qJH1=2 ðGÞ rCh

JxJHg þ 1 ðGÞ

bþ1 JuJHb þ 1 ðGÞ Juu~ h JHg1 ðGÞ r Juu~ h JH0 ðGÞ r Ch

þh

In this subsection, we will prove that the numerical result obtained using the GBNM converges to the exact solution of the elasticity problem (1) gradually.

mþgþ3

ð25Þ

Besides, recalling again Eq. (14) yields

bþ1

3.1. A priori error analysis

Jq~ h JH0 ðGÞ JxJH0 ðGÞ

From Eq. (22) we have

/Aðqq~ h Þ, xSL2 ðGÞ rCfh

Jqq~ h JHg2 ðGÞ rCfh

ð24Þ

JA1 gJHg þ 1 ðGÞ

g A Hg þ 2 ðGÞ

r CJqq~ h JH1=2 ðGÞ JxSh xJH1=2 ðGÞ

Rm

m þ 3=2

9/Aðqq~ h Þ,A1 gSL2 ðGÞ 9

þ/Ah qh ,Sh xJCSL2 ðGh Þ /Aq~ h ,Sh xSL2 ðGÞ

Fm ðxÞU ij ðx,yÞ dSx , j ¼ 1; 2,3, y A O ð21Þ

Jqq~ h JH1=2 ðGÞ r Cfh

JgJHg þ 2 ðGÞ

g A Hg þ 2 ðGÞ

ð23Þ

bþ1

ðdx Þ1

ð33Þ

Hence, via using Theorem 3.1 and substituting Eqs. (31)–(33) into Eq. (30) we can prove Eq. (29) for 9k9 ¼ 0. Other cases are similar, and so are omitted. &

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

997

Theorem 3.2 gives the error of the displacement field u and its derivatives. The error result shows extremely high accuracy can be achieved not only for u but also for its derivatives. As a direct consequence of Theorem 3.2, we can easily show that the errors of stress r and displacement u in the GBNM are all of the same order. It is clear from Theorems 3.1 and 3.2 that the errors of the GBNM come from two terms: one from the approximation of the boundary function by the MLS approximations, and one from the approximation of the boundary surface G by the integration cell Gh . When Gh coincides with G, the boundary integration can be performed exactly on the boundary surface. In this case, the approximate spaces given by Eq. (7) can be rewritten as W h ðGÞ :¼ spanfFi ,1 ri r Ng,

Vh ðGÞ :¼ ðW h ðGÞÞ3

ð34Þ

Fig. 1. Results of u1 , s22 and s12 along the x1-axis of the cube for the patch test.

where Fi is the MLS shape function defined on G. And the approximate solution uh will be determined by uj ðyÞ ¼

3 Z X i¼1

G

qih ðxÞU ij ðx,yÞdSx ,

j ¼ 1; 2,3,

yAO

ð35Þ

where qh ¼ ðq1h ,q2h ,q3h Þ A Vh ðGÞ. Now, the results of Theorems 3.1 and 3.2 can be simplified as Jqqh JH1=2 ðGÞ r Ch

m þ 3=2

9@k uðxÞ@k uh ðxÞ9 r C

JqJHm þ 1 ðGÞ

gX þ2

ðdx Þj9k91 h

mþgþ3

JqJHm þ 1 ðGÞ

ð36Þ

j¼0

Therefore, in the case of the integration cell Gh coincides with the boundary surface G, the geometric error resulting from the approximation of the boundary by cells vanishes. 3.2. Numerical examples In this subsection, some numerical examples are performed to show the performance of the GBNM model established and to verify the validity of the a priori error estimates.

3.2.1. Patch test For working a numerical method for solid mechanics problems, the sufficient requirement for convergence is to pass the standard patch test. Therefore, the first numerical example is the standard patch test using the presented GBNM. The problem is studied in a cubic domain the dimensions of 2  2  2. The origin of the coordinate system is located on the centroid of the cube. The following analytical solution is employed for the patch test: u1 ¼

2x1 þx2 þx3 , 2

u2 ¼

x1 þ 2x2 þ x3 , 2

u3 ¼

x1 þx2 þ2x3 2

The calculated x1-component displacement u1 , the normal stress s22 and the shearing stress s12 along the x1-axis of the cube are plotted in Fig. 1 together with the analytical solutions. In this analysis, the cubic surface is discretized using 48 regularly distributed nodes. The material constants are chosen to be E¼1.0 and n ¼ 0:25. It can be found from Fig. 1 that the numerical results agree very well with the analytical ones. Therefore, we can conclude that the GBNM passes the patch test successfully.

3.2.2. Hollow sphere problem This example involves an object with curved surfaces—a hollow sphere. The sketched geometric configuration of this problem is shown in Fig. 2. For this benchmark problem, the analytical solutions for the radial displacement, and radial and

Fig. 2. Schematic diagram for the hollow sphere problem.

tangential stresses are available in the polar co-ordinate system as " # 3 Pa3 r b ur ¼ ð12nÞ þ ð1þ nÞ 3 3 2r Eðb a3 Þ 3

srr ¼

Pa3 ðb r 3 Þ 3 r 3 ða3 b Þ 3

syy ¼

Pa3 ðb þ 2r 3 Þ 3

2r 3 ðb a3 Þ

where r is the radial distance form 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 taken as P¼1, a¼1.0, b¼4.0, E¼1.0 and n ¼ 0:25 in the computation. The numerical results for ur , srr and syy are shown in Fig. 3 together with the analytical solutions. In this analysis, we employed 48 regularly distributed nodes on each sphere surface. Besides, quadratic triangles are used to describe the geometry of the sphere surfaces. From this figure, it can be observed that the numerical solutions are in excellent agreement with the analytical solutions. For investigating the influence of the integration cells to the solution accuracy, the boundary surface is approximated by two types of cells, i.e. linear triangles, b ¼ 1, and quadratic triangles, b ¼ 2. The values of the numerical approximations of ur , srr and

998

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

syy in both cases are listed in Table 1. The corresponding experimental convergence rate is also given. As we expected, the results from the proposed meshless method gradually converge to the analytical values along with the increase of the boundary nodes, and stress r and displacement u have almost the same experimental convergence rates. Comparing the numerical results, we can observe that the precision and convergence rates of b ¼ 2 are higher than those of b ¼ 1. The reason for this phenomenon is that the representation of the integration cells is quite a severe approximation on the original boundary surface and thus some corresponding loss of accuracy should be expected in the numerical solutions. These observations are consistent with the theoretical error analysis. Besides, according to the a priori error analysis presented in Section 3.1, the theoretical convergence rates are 2 for b ¼ 1 and 3 for b ¼ 2. From Table 1 we can observe that the experimental convergence rates approximate the theoretical rates. Consequently, the numerical results for this example all match the theoretical error estimates. 3.2.3. Linear field problem on a horn domain For the validation of the proposed method for irregular computational geometry, an elastic problem on a more complicated geometry is also solved. In this example, the domain is a closed horn with boundary surface G ¼ G1 [ G2 [ G3 . Here

and

G3 ¼ fðx1 ,x2 ,x3 Þ A R3 : x1 ¼ R, x2 ¼ r sin a, x3 ¼ 0:5r cos a, a A ½p, pg are bottom surfaces of the horn. Fig. 4 depicts the sketched geometric configuration and a discretization consisting of 288 integration cells with one node per cell. Because no analytical solution can be found for a practical load case on this complex structure, to assess the accuracy of the GBNM of this problem, we use the following linear displacement fields [36]: u1 ¼ 2x21 þ 3x22 þ 3x23 ,

u2 ¼ 3x21 2x22 þ 3x23 ,

u3 ¼ 3x21 þ 3x22 2x23

Displacements are imposed on all faces of the horn according to the analytical solutions. The parameters are taken as R¼9, r ¼1, E¼1.0 and n ¼ 0:25 in the computation. Fig. 5 displays a comparison between the present numerical results with the analytical solutions for displacement u and stress r along the arc given by the formulas x1 ¼ Rð1 þ sin bÞ, x2 ¼ 0, x3 ¼ Rð1 þcos bÞ, b A ½p,0:5p. It is clearly shown that the numerical solutions are seen to capture the behavior of the analytical solutions very well. Thus, we can conclude that the presented method free from mesh generation can be applied

G1 ¼ fðx1 ,x2 ,x3 Þ A R3 : x1 ¼ Rþ ðR þ2r cos aÞsin b,x2 ¼ r sin a, x3 ¼ R þ ðR þ 0:5r cos aÞcos b, a A ½p, p, b A ½p,0:5pg is the side surface of the horn, while

G2 ¼ fðx1 ,x2 ,x3 Þ A R3 : x1 ¼ 2r cos a, x2 ¼ r sin a, x3 ¼ R, a A ½p, pg

Fig. 3. ur , srr and syy for the hollow sphere problem.

Fig. 4. Schematic diagram and integration cells for a horn domain (the dot represents boundary nodes).

Table 1 Approximations and convergence rates for the hollow sphere problem (upper low: b ¼ 1; lower row: b ¼ 2). r

1.5

Numerical solutions

ur

srr syy 2.5

ur

srr syy

Exact solutions

N ¼96

N¼ 384

N ¼864

N¼ 1536

0.272150 0.292216  0.263515  0.282835 0.150931 0.164148

0.288535 0.293790  0.279118  0.284807 0.161813 0.166027

0.291226 0.294007  0.281923  0.285017 0.164343 0.166275

0.292724 0.294055  0.283583  0.285074 0.165244 0.166325

0.112070 0.120909 0.046572  0.049007 0.045863 0.048238

0.119021 0.121341 0.048254  0.049120 0.047610 0.048355

0.120411 0.121408 0.048733  0.049136 0.048064 0.048371

0.120812 0.121418 0.048953  0.049140 0.048190 0.048378

0.294092  0.285126 0.166373 0.121429  0.049143 0.048381

Experimental rates

Theoretical rates

1.96 2.84 1.86 2.73 1.88 2.80

2.00 3.00 2.00 3.00 2.00 3.00

1.98 2.86 1.84 2.73 1.88 2.62

2.00 3.00 2.00 3.00 2.00 3.00

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

999

Fig. 5. Results of displacement u and stress r for the linear field problem on a horn domain.

to linear elastic problems with complex geometry under easy implementation.

can define another approximate space Vh=2 ðGÞ and obtain another approximate solution qh=2 A Vh=2 ðGÞ. Then we have Jqqh=2 JH1=2 ðGÞ rJqqh JH1=2 ðGÞ

4. A posteriori error estimates and adaptive procedures

Let

In this section, a posteriori error estimates and adaptive refinement processes are developed for the meshless GBNM for 3D elasticity. For brevity, we consider here only screen crack problems in R3 , i.e. boundary value problems exterior to an open surface G.

In the GBNM, background cells are introduced for numerical integration over the boundary surface. In order to minimize the integration error, the integration order should match the MLS approximation order [25]. Consequently, one should subdivide cells or increase the integration order for cells with newly added nodes. Since cells can be divided into smaller ones without affecting their neighbors in any way, the cells are refined to keep the ratio of nodes per cell small and roughly the same in many adaptive meshless methods [24–26,29,31]. In the developed adaptive procedure, this scheme is more convenient to adopt to retain the integration accuracy as the nodal density in the local area increases. As before, let T h denote a triangulation of the boundary surface G into N cells and let Q h ¼ fxi gN i ¼ 1 denote the set of the associated boundary nodes. For simplicity, assume that the integration cell coincides with the boundary surface. Then based on Q h only, the approximate space Vh ðGÞ can be defined by Eq. (34). The weak solution of the considered screen elasticity problems is determined by the variational problem to find q A H1=2 ðGÞ such that 8q0 A H1=2 ðGÞ

ð37Þ

1=2

ðGÞ, the corresponding approximate variaSince Vh ðGÞ  H tional problem reads as: find qh A Vh ðGÞ such that /Aqh ,q SL2 ðGÞ ¼ /u,q SL2 ðGÞ , 0

0

0

8q A Vh ðGÞ

ð38Þ

Subtracting Eq. (38) from Eq. (37) we obtain the following Galerkin orthogonality /Aðqqh Þ,q0 SL2 ðGÞ ¼ 0,

e :¼ Jqqh JH1=2 ðGÞ

ð41Þ

be the exact error of the approximate solution qh , then using the Galerkin orthogonality (39) yields e2 ¼ Jqqh=2 J2H1=2 ðGÞ þJqh=2 qh J2H1=2 ðGÞ

ð42Þ

It is natural that the difference of the two numerical solutions qh and qh=2 can be employed to estimate the error e. Therefore, we define the error estimator as

4.1. A posteriori error analysis

/Aq,q0 SL2 ðGÞ ¼ /u,q0 SL2 ðGÞ ,

ð40Þ

8q0 A Vh ðGÞ

ð39Þ

In order to deduce a posteriori error estimates, let T h=2 be a second triangulation of G gained from a refinement of T h and let Q h=2 be the associated boundary nodes. Based on Q h=2 only, we

Z :¼ Jqh=2 qh JH1=2 ðGÞ

ð43Þ

From Eq. (42) it follows in a straightforward manner that Z r e, which shows the efficiency of Z with a fixed efficiency constant one. The reliability of Z is commonly proven with the aid of the saturation assumption Jqqh=2 JH1=2 ðGÞ rC s Jqqh JH1=2 ðGÞ ¼ C s e

ð44Þ

with a constant C s A ð0; 1Þ. Then according to Eq. (42) we get e2 rC 2s e2 þ Z2 and therefore reliability er ð1C 2s Þ1=2 Z

ð45Þ

On the other hand, if there exists a constant C r independent of h such that e rC r Z, then using again Eq. (42) we obtain 2 Jqqh=2 J2H1=2 ðGÞ ¼ e2 Z2 r ð1C 2 and thus the saturation r Þe qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi assumption (44) with C s ¼ 1C 2 r .

Given the above facts, we arrive at the following a posteriori error estimates. Theorem 4.1. The error estimator Z is always efficient and is reliable under the saturation assumption (44), i.e.

Z r e rð1C 2s Þ1=2 Z

ð46Þ

Besides, the saturation assumption (44) is equivalent to the reliability of Z. A posteriori error estimates of the form Eq. (46) are very important in practice since they are used to justify the effectiveness of the resulting adaptive procedure. Note that the only unproven assumption that we make for our error analysis is the saturation assumption. This assumption is very moderate and natural in the following sense: First, it is a strengthened version of

1000

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

Eq. (40). Second, it essentially shows that the approximate solution in the larger meshless space is a better approximation to the exact solution. However, it is difficult to provide a theoretical proof for the saturation assumption, but the validity of this assumption is observed numerically for 2D potential problems [31] and will be confirmed in subsequent numerical examples for 3D elasticity problems. In order to determine which cells must be refined and then to generate the optimal nodal configuration, the approximation error needs to be estimated locally. However, since J  J2H1=2 ðGÞ a PN 2 i ¼ 1 J  JH1=2 ðG Þ , the error estimator Z does not provide informai

tion for a local refinement directly. Fortunately, we have P 2 J  J2L2 ðGÞ ¼ N i ¼ 1 J  JL2 ðG Þ , hence error estimators based on the i

L2-norm can be adopted for the indicator-based adaptive refinement. We therefore define the following error estimator: pffiffiffi m :¼ J hðqh=2 qh ÞJL2 ðGÞ ð47Þ where h A L1 ðGÞ, h9Gi :¼ hi :¼ diamðGi Þ for Gi A T h , denotes the local cell-size function. As in the 2D potential problems [31], it can be verified that a posteriori error estimators Z and m are equivalent. Then using Theorem 4.1 we can immediately obtain the efficiency and reliability of m. 4.2. Adaptive procedures Based on the above a posteriori error estimates, we introduce an adaptive meshless procedure for the numerical solution of 3D elasticity. Due to the local property of J  JL2 ðGÞ , the error estimator qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN 2 m given in Eq. (47) can be rewritten as m ¼ i ¼ 1 mi , where the local error indicators pffiffiffiffi mi :¼ hi Jqh=2 qh JL2 ðGi Þ ,

i ¼ 1; 2, . . . ,N

ð48Þ

can be computed in cells and can be employed for cell marking in the following adaptive cell-refining algorithm. Algorithm 4.1. Choose an initial triangulation T hð0Þ , tolerance

algorithm should return qðkÞ instead of qðkÞ . On the other hand, h=2 h the computation of qðkÞ is much cheaper than that of qðkÞ . h h=2 Therefore, compared to the cost of computing the approximate , the cost of computing the estimated errors ZðkÞ and solution qðkÞ h=2

mðkÞ is fractional. The local error indicator mi given in Eq. (48) is computed for each individual quadrature cell. It is therefore necessary to determine the size of the influence domain locally. In the MLS approximations, to ensure regularity of the moment matrix and to improve the approximation precision, the size of the influence domain should be tuned such that the number of nodes contained in an influence domain is reasonable. Until now, there is no proper way to decide optimal sizes theoretically. In practical computations, this size is usually defined as a multiple of the nodal spacing for uniformly distributed nodes. For non-uniformly distributed nodes, there are also some recommendations, but they depend on experience [2,9,25,26,29,31]. Similar to the strategies used in Refs. [2,9,26,29,31], the strategy used in this study is that each evaluation point contains roughly the same number of boundary nodes in its influence domain. 4.3. Numerical examples To demonstrate the application and validity of the developed a posteriori error estimates and adaptive refinement procedure for solving 3D elasticity screen crack problems, some numerical examples are considered in this subsection. According to the Galerkin orthogonality (39) and variational problems (37) and (38), the error e can be computed as e ¼ ð/f,qSL2 ðGÞ /f,qh SL2 ðGÞ Þ1=2 where /f,qh SL2 ðGÞ ¼

P3

j¼1

ð49Þ

PN

j ðmÞ m ¼ 1 f m qjh

with

j fm

¼

R

Rm u j ðyÞFm

ðyÞ dSy , while /f,qSL2 ðGÞ is either evaluated analytically or obtained by Aitkin’s D2 extrapolation of the sequence of values for numerical solutions on uniformly refined nodes. Besides, from Eqs. (39) and (49), the error estimator Z can be computed as

e 4 0, adaptive parameter y A ð0; 1Þ, and set k¼0.

Z ¼ ð/f,qh=2 SL2 ðGÞ /f,qh SL2 ðGÞ Þ1=2

ðkÞ ðkÞ from a refinement of T ðkÞ ¼ fGðkÞ 1. Obtain T ðkÞ 1 , G2 , . . . , GN g. h=2 h

4.3.1. Example on the square screen In this example, we consider the screen given by the square plate G ¼ ½1; 12 . The boundary datum is u ¼ ðx2 ,x1 ,0Þ. For this problem, the analytical solution q is unknown, so the energy norm of the analytical solution is gained from Aitkin’s D2 extrapolation, /f,qSL2 ðGÞ ¼ 13,233:1872. In order to illustrate the singular behavior at edges and corners, Fig. 6 plots the computed numerical solution q1h on the square screen with uniformly distributed nodes N ¼32 and N ¼128, while Fig. 7 plots the computed numerical solution q2h . Since the computed numerical solution q3h approximates 0, the corresponding figure is omitted here. Fig. 6 shows that q1h yields strong singularities at edges x2 ¼ 71 and all corners of the screen, while Fig. 7 shows that q2h yields strong singularities at edges x1 ¼ 71 and all corners. Therefore, although there is no analytical solution available in the literature to verify the singularity for such a problem, Figs. 6 and 7 predict that generic singular behavior might occur along all edges of the square screen and the singularity at corners of the square is very severe. Hence adaptive procedure should yield a better convergence. Fig. 8 plots the error e and error estimator Z in the energy norm for the uniform and the adaptive refinement with respect to the number of boundary nodes. For uniform refinement, we get experimentally a convergence of the form OðN a Þ with the number 2 of degrees of freedom N ¼ Oð1=h Þ and a  1=4, where h is the

2. Compute approximate solutions qðkÞ A VðkÞ ðGÞ for T ðkÞ and h h h qðkÞ A VðkÞ ðGÞ for T ðkÞ . h=2 h=2 h=2 3. Compute error estimator ZðkÞ ¼ JqðkÞ qðkÞ J and stop the h=2 h H1=2 ðGÞ ðkÞ algorithm if Z r e. ðkÞ 4. For each cell GðkÞ compute local error indicators i AT h qffiffiffiffiffiffiffi ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ mi ¼ hi Jqh=2 qh JL2 ðGðkÞ Þ with hðkÞ i :¼ diamðGi Þ. i qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN ðkÞ 2 5. Compute mðkÞ ¼ and choose a minimal set i ¼ 1 ðmi Þ P ðkÞ ðkÞ Þ2 Z yðmðkÞ Þ2 . M D T h such that GðkÞ A MðkÞ ðmðkÞ i i

6. Refine all marked cells GðkÞ A MðkÞ and so generate a new i þ 1Þ triangulation T ðk , update k, and go to 1. h Ref. [31] presents a first convergence result for adaptive meshless algorithms. Similarly, it can be proven that the sequence of approximate solutions generated by Algorithm 4.1 is convergent, i.e. lim JqqðkÞ J ¼ lim JqqðkÞ J ¼ lim ZðkÞ ¼ lim mðkÞ ¼ 0 h H1=2 ðGÞ h=2 H1=2 ðGÞ

k-1

k-1

Since Eq. (40) implies that

k-1

qðkÞ h=2

k-1

is a better approximation of the

exact solution q than qðkÞ when Algorithm 4.1 stops, this h

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

1001

Fig. 6. Numerical solutions q1h for the square screen problem with (a) N ¼32 and (b) N ¼ 128.

Fig. 7. Numerical solutions q2h for the square screen problem with (a) N ¼32 and (b) N ¼ 128.

Fig. 8. The comparison of the error e and the error estimator Z for uniform and adaptive refinement for the square screen problem.

cell size. This is due to singularities of the unknown q at the edges of the square screen. In the case of a regular solution q A H1 ðGÞ, from Eq. (36) we can expect a convergence rate OðN3=4 Þ, which 3=2 corresponds to Oðh Þ. This motivates and leads to the implementation of robust adaptive procedures to improve the convergence rate. Numerical results in Fig. 8 indicate that the adaptive refinement accelerates the convergence up to a convergence rate of about OðN1 Þ. Therefore, the proposed adaptive algorithm recovers the optimal convergence rate. Besides, as expected from Theorem 4.1, this figure shows clearly that the curve of the error

estimator Z is always below the curve of the error e and both curves are almost parallel for both uniform and adaptive refinement. These observations strongly underline the empirical efficiency and reliability of the error estimator Z. Since Theorem 4.1 indicates that the reliability is equivalent to the saturation assumption, the saturation assumption holds for this example. The corresponding adaptively generated cell triangulations are depicted in Fig. 9. It can be found from this figure that the adaptive process has the same direction on all edges of the screen G. This fact matches the symmetric character of the numerical solutions shown in Figs. 6 and 7. Besides, we observe a high refinement towards all edges and corners. This observation is in good consistency with the expected singularities shown in Figs. 6 and 7. The refinement is also in agreement with our expectation as 3D edge singularities dominate over corner singularities. Moreover, Fig. 9 underlines that cells in the presented adaptive GBNM do not need to satisfy any aspect ratio or compatibility requirements such as in the BEM or the FEM. Therefore, adaptive procedures in the meshless GBNM can be implemented more easily than these in the BEM and the FEM.

4.3.2. Example on the disk screen In order to verify that the developed adaptive procedure has good capability to tackle problems with curved surfaces, we consider a penny-shaped crack with unit radius in an unbounded elastic solid and choose the boundary datum as u ¼ ð1; 1,1Þ. In Fig. 10, we report the computed numerical solution q3h on the screen with uniformly distributed nodes N ¼36 and N¼ 144. The figures for numerical solutions q1h and q2h are similar, and so are

1002

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

Fig. 9. Adaptively generated cells for the square screen problem: (a) initial triangulation with N ¼ 32; (b) step 2 with N ¼104; (c) step 5 with N ¼ 594; (d) step 7 with N ¼1208.

Fig. 10. Numerical solutions q3h for the disk screen problem with (a) N ¼ 36 and (b) N ¼144.

omitted. The numerical solutions suggest strong singularity at the circle of the disk screen. Some resulting refined cell configurations are depicted in Fig. 11. This figure indicates expected refinement towards the circle. This is in good agreement with the expected singularity shown in Fig. 10.

5. Conclusions In this study, the GBNM has been extended to solve 3D problems in elasticity. It is a boundary only meshless method that combines an equivalent variational formulation of BIEs and

MLS approximations. In this approach, boundary conditions can be implemented directly and easily despite the fact that the MLS shape functions lack the delta function property. In addition, the resulting formulation is not only symmetric but also positive definite. A priori error estimate of the GBNM for 3D elasticity problems has been derived. The error analysis shows that the error bound of the numerical solution is directly related to the nodal spacing. Besides, the error results mainly from the approximation of the boundary function by the MLS approximations and the approximation of the boundary by the integration cell structure. Moreover, displacements and stresses hold the same rate of

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

1003

Fig. 11. Adaptively generated cells for the disk screen problem: (a) initial triangulation with N ¼ 16; (b) step 3 with N ¼94; (c) step 6 with N¼ 562; (d) step 8 with N ¼1654.

convergence. Numerical experiments have been given and the numerical results are accurate and are in good consistency with the theoretical analysis. A posteriori error estimates and error estimators of the GBNM for 3D elasticity have also been developed. The estimators estimate the unknown error both reliably and efficiently. Combining the GBNM with a posteriori error estimators, an effective and convergent adaptive meshless algorithm has been established. Some examples with high singularities have been presented to conform the theoretical results and to show the effectiveness of the GBNM for adaptive analysis. As shown in the numerical examples, the a posteriori error estimates are reliable and the developed adaptive algorithm is attractive for obtaining adaptive meshless solutions of 3D elasticity problems. The MLS approximation, used in the GBNM, is computationally expensive, since for each integration sampling point, a linear system has to be solved. This makes the GBNM and other MLSbased meshless methods too expensive for problems requiring a dense distribution of points. Some techniques have been developed for the alleviation of this problem, such as fast algorithm [13,37], coupling algorithm [1,2] and parallel computing [9], but how to improve the efficiency of meshless methods is still an open problem. Since the system matrix in the GBNM is not only symmetric but also positive definite, this method is an ideal choice for coupling the FEM or the Galerkin BEM to improve the computational efficiency. Besides, fast algorithms based on multipole expansions, fast Fourier transform or singular value decomposition can also be used for acceleration of this method. Furthermore, the presented adaptive algorithm can be extended to solve crack problems. These are subjects for future research.

Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant nos. 11101454 and 11026198) and the Natural Science Foundation Project of CQ CSTC (Grant no. cstcjjA30002). The author is very grateful to the anonymous reviewers for their invaluable comments and suggestions which have greatly improved the quality and the presentation of this paper. References [1] Li S, Liu WK. Meshfree particle methods. Berlin: Springer; 2004. [2] Liu GR. Meshfree methods: moving beyond the finite element method.second ed. Boca Raton: CRC Press; 2009. [3] Krysl P, Belytschko T. Element-free Galerkin method: convergence of the continuous and discontinuous shape functions. Comput Methods Appl Mech Eng 1997;148:257–77. [4] Babuˇska I, Nistor V, Tarfulea N. Generalized finite element method for second-order elliptic operators with Dirichlet boundary conditions. J Comput Appl Math 2008;218:175–83. [5] Fries TP, Belytschko T. The extended/generalized finite element method: an overview of the method and its applications. Int J Numer Methods Eng 2010;84:253–304. [6] Liu GR, Nguyen-Xuan H, Nguyen-Thoi T. A theoretical study on the smoothed FEM (S-FEM) models: properties, accuracy and convergence rates. Int J Numer Methods Eng 2010;84:1222–56. [7] Cheng R, Cheng Y. Error estimates for the finite point method. Appl Numer Math 2008;58:884–98. [8] Mukherjee YX, Mukherjee S. The boundary node method for potential problems. Int J Numer Methods Eng 1997;40:797–815. [9] Mukherjee S, Mukherjee YX. Boundary methods: elements, contours, and nodes. Boca Raton: CRC; 2005. [10] Zhang J, Yao Z, Li H. A hybrid boundary node method. Int J Numer Methods Eng 2002;53:751–63.

1004

X. Li / Engineering Analysis with Boundary Elements 36 (2012) 993–1004

[11] Wang Q, Zheng JJ, Miao Y, Lv JH. The multi-domain hybrid boundary node method for 3D elasticity. Eng Anal Boundary Elem 2011;35:803–10. [12] Miao Y, He T, Yang Q, Zheng J. Multi-domain hybrid boundary node method for evaluating top-down crack in Asphalt pavements. Eng Anal Boundary Elem 2010;34:755–60. [13] Wang Q, Miao Y, Zheng JJ. The hybrid boundary node method accelerated by fast multipole expansion technique for 3D elasticity. CMES: Comput Model Eng Sci 2010;70:123–51. [14] Zhang J, Qin X, Han X, Li G. A boundary face method for potential problems in three dimensions. Int J Numer Methods Eng 2009;80:320–37. [15] Li X, Zhu J, Zhang S. A hybrid radial boundary node method based on radial basis point interpolation. Eng Anal Boundary Elem 2009;33:1273–83. [16] Peng M, Cheng Y. A boundary element-free method (BEFM) for two-dimensional potential problems. Eng Anal Boundary Elem 2009;32:77–82. [17] Li X, Zhu J. A Galerkin boundary node method and its convergence analysis. J Comput Appl Math 2009;230:314–28. [18] Li X, Zhu J. A Galerkin boundary node method for biharmonic problems. Eng Anal Boundary Elem 2009;33:858–65. [19] Li X, Zhu J. A Galerkin boundary node method for two-dimensional linear elasticity. CMES: Comput Model Eng Sci 2009;45:1–29. [20] Li X, Zhu J. A meshless Galerkin method for Stokes problems using boundary integral equations. Comput Methods Appl Mech Eng 2009;198:2874–85. [21] Li X, Zhu J. Meshless analysis of two-dimensional Stokes flows with the Galerkin boundary node method. Eng Anal Boundary Elem 2010;34:79–91. [22] Li X. Meshless Galerkin algorithms for boundary integral equations with moving least square approximations. Appl Numer Math 2011;61:1237–56. [23] Li X. The meshless Galerkin boundary node method for Stokes problems in three dimensions. Int J Numer Methods Eng 2011;88:442–72. [24] Duarte AC, Oden JT. An hp adaptive method using clouds. Comput Methods Appl Mech Eng 1996;139:237–62. [25] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. Int J Numer Methods Eng 2004;63:1559–82.

[26] Gan N, Li G, Long S. 3D adaptive RKPM method for contact problems with elastic–plastic dynamic large deformation. Eng Anal Boundary Elem 2009; 33:1211–22. [27] Angulo A, Pe´rez Pozo L, Perazzo F. A posteriori error estimator and an adaptive technique in meshless finite points method. Eng Anal Boundary Elem 2009;33:1322–38. [28] Tang Q, Zhang GY, Liu GR, Zhong ZH, He ZC. A three-dimensional adaptive analysis using the meshfree node-based smoothed point interpolation method (NS-PIM). Eng Anal Boundary Elem 2011;35:1123–35. [29] Chati MK, Paulino GH, Mukherjee S. The meshless standard and hypersingular boundary node methods—applications to error estimation and adaptivity in three-dimensional problems. Int J Numer Methods Eng 2001;50: 2233–69. [30] Guo XF, Chen HB. Dual error indicators for the local boundary integral equation method in 2D potential problems. Eng Anal Boundary Elem 2006; 30:702–8. [31] Li X. Adaptive methodology for the meshless Galerkin boundary node method. Eng Anal Boundary Elem 2011;35:750–60. [32] Dautray R, Lious JL. Mathematical analysis and numerical methods for science and technology. In: Integral equations and numerical methods, vol. 4. Berlin: Springer; 2000. [33] Zhu J, Yuan Z. Boundary element analysis. Beijing: Science Press; 2009. [34] Hong HK, Chen JT. Derivations of integral equations of elasticity. ASCE J Eng Mech 1988;114:1028–44. [35] Chen JT, Hong HK. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. Appl Mech Rev 1999;52: 17–33. [36] Gu J, Zhang J, Sheng X, Li G. B-spline approximation in boundary face method for three-dimensional linear elasticity. Eng Anal Boundary Elem 2011;35: 1159–67. [37] Zhang J, Zhuang C, Qin X, Li G, Sheng X. FMM-accelerated hybrid boundary node method for multi-domain problems. Eng Anal Boundary Elem 2010; 34:433–9.