- Email: [email protected]
A hybrid radial boundary node method based on radial basis point interpolation
ARTICLE IN PRESS Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A hybrid radial boundary node method based on radial basis point interpolation Xiaolin Li a,, Jialin Zhu a, Shougui Zhang b a b
College of Mathematics and Physics, Chongqing University, Chongqing 400044, PR China College of Mathematics and Computer Science, Chongqing Normal University, Chongqing 400047, PR China
a r t i c l e in fo
abstract
Article history: Received 14 January 2009 Accepted 4 June 2009 Available online 17 July 2009
The hybrid boundary node method (HBNM) retains the meshless attribute of the moving least squares (MLS) approximation and the reduced dimensionality advantages of the boundary element method. However, the HBNM inherits the deficiency of the MLS approximation, in which shape functions lack the delta function property. Thus in the HBNM, boundary conditions are implemented after they are transformed into their approximations on the boundary nodes with the MLS scheme. This paper combines the hybrid displacement variational formulation and the radial basis point interpolation to develop a direct boundary-type meshless method, the hybrid radial boundary node method (HRBNM) for two-dimensional potential problems. The HRBNM is truly meshless, i.e. absolutely no elements are required either for interpolation or for integration. The radial basis point interpolation is used to construct shape functions with delta function property. So unlike the HBNM, the HRBNM is a direct numerical method 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. Some selected numerical tests illustrate the efficiency of the method proposed. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Meshless Hybrid radial boundary node method Radial basis function Boundary integral equations
1. Introduction In recent years, the meshless or meshfree methods have generated considerable interest among researchers for solving boundary value problems [1,2]. Compared to the finite element method (FEM) and the boundary element method (BEM), the main feature of this type of method is the absence of an explicit mesh, and the approximate solutions are constructed entirely based on a cluster of scattered nodes. Meshless methods have been found to have special advantages on problems to which the conventional mesh-based methods are difficult to be applied. These include problems with complicated boundary, moving boundary, large deformations, mesh adaptively, and so on. The moving least squares (MLS) as an approximation method has been introduced by Lancaster and Salkauskas [3]. Since the numerical approximations of the MLS start from a cluster of scattered nodes instead of meshes, there have many meshless methods based on the MLS scheme for the numerical solution of differential equations. Typical of them are the diffuse element method [4], the element-free Galerkin (EFG) method [5] and the meshless local Petrov–Galerkin (MLPG) method [6]. The EFG method uses a global weak form. Although no mesh is required in
Corresponding author. Fax: +86 23 65120520.
E-mail address: [email protected] (X. Li). 0955-7997/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2009.06.003
this method for variable interpolation, background cells are necessary for ‘energy’ integration in general. The MLPG approach is formulated on a local weak form. The main advantage of the MLPG approach over the FEM and the EFG method is that it does not need a ‘finite element mesh’, either for purposes of interpolation of the solution variables or for the integration of the ‘energy’. The MLS approximation has also been used in boundary integral equations (BIEs). Among them are the meshless local boundary integral equation (LBIE) method [7], the boundary node method (BNM) [8], the hybrid boundary node method (HBNM) [9], the boundary element-free method (BEFM) [10,11] and the Galerkin boundary node method (GBNM) [12,13]. The first boundary-type meshless method was the BNM. This method takes the advantages of both the BIE in dimension reduction and the MLS approximation in elements elimination. Nevertheless, background cells are required for integration. To avoid this hindrance, Zhu et al. proposed the LBIE method. It is equivalent to a sort of MLPG approach, which needs absolutely no domain and boundary elements. The LBIE method, however, is not strictly a boundary method since it requires evaluation of integrals over certain surfaces (called Ls in [7]) that can be regarded as ‘closure surfaces’ of boundary elements. To obtain a truly boundary-type meshless method, Zhang et al. [9] proposed the HBNM that combines the MLS approximation with the hybrid displacement variational formulation. In this approach, the integration is limited to a fixed local region, thus no cells are needed either
ARTICLE IN PRESS 1274
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
for the interpolation purposes or for integration process. The HBNM is a very promising method, which has been used for problems in potential theory [9,14,15], linear elasticity [16–18] and elastodynamics [19]. These MLS-based meshless methods have gained remarkable progress in solving a broad class of boundary value problems. However, since the MLS approximations lack the delta function property of the usual FEM and BEM shape functions, boundary conditions in these meshless techniques are difficult to be enforced. This problem becomes even more serious in the boundary-type meshless methods, since a large number of boundary conditions need to be satisfied. The strategy proposed in the BNM [8] involves a new definition of the discrete norm used for the construction of the MLS approximations, making the BNM computationally much more expensive than the traditional BEM. 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. However, they are not the real nodal variables, and thus boundary conditions could not be directly applied. The BNM and the HBNM are the indirect boundary-type meshless methods. In addition to the abovementioned strategy, Kitipornchai et al. [10] presented an improved MLS approximation that uses weighted orthogonal polynomials as basis functions to restore the delta function property of the MLS. The BEFM is based on the improved MLS, thus it is a direct boundary-type meshless method. Recently, via using an equivalent variational form of BIEs, another effective technique is developed in the GBNM to impose boundary conditions [12]. The radial point interpolation method (RPIM) [2,20] is another meshless interpolation technique proposed by Liu et al. The RPIM approximates its field variable by letting the problem function pass through each scattered node in its local support domain and radial basis functions (RBFs) are used for point interpolation. Different from the MLS approximation, the shape functions so constructed have the delta function property. As a result, the construction of a RPIM shape function is more efficient than the MLS procedures. Some meshless methods using the RPIM shape functions have been developed. Typical of them are the RPIM [20], the local radial point interpolation method (LRPIM) [21,22] and the boundary radial point interpolation method (BRPIM) [23,24]. As the MLPG approach, the LRPIM is a truly domain-type meshless method. The BRPIM is a boundary-type method, but it requires background cells for numerical integration. In these RPIM-based methods, boundary conditions can be enforced as easily as in the conventional FEM and BEM, thus they are much more efficient than the methods using MLS shape functions. Up to date, these approaches have been successfully implemented for problems of solid [20–24], fluid [25] and elastodynamics [26]. An interesting question arises here—Is there possibly a boundary-type meshless method that only involves the real solution of nodal variables and does not need any element either for interpolation or for integration? Such a method is a direct numerical approach, in which boundary conditions can be applied directly and easily, compared with the HBNM and the BNM; and it is a truly meshless method, which require no mesh for either interpolation or integration purpose, compared with the BNM and the BRPIM. The answer to this question is positive. The new method is called the hybrid radial boundary node method (HRBNM) and is the subject of this paper. The HRBNM combines the RPIM with the hybrid displacement variational formulation. Hybrid boundary element formulations are based on generalized variational principles. In the late eighties, two hybrid models have been developed: the hybrid stress boundary element model [27] and the hybrid displacement boundary element model [28].
The first is based on the Hellinger–Reissner principle with the stresses in the domain and displacements on the boundary as independent functions, and the second on a modified functional using three independent variables, two of them defined on the boundary and the other inside the domain. Both models use the classical fundamental solution to approximate the domain variable and thus allowing for the transfer of the domain integrals to the boundary. Besides, the resulting coefficient matrices are symmetric. Recently, with combination of the RPIM and the hybrid displacement boundary element model, a hybrid BRPIM has been developed [29]. The hybrid BRPIM results in symmetric matrices, but an underlying cell structure is still used for integration. As in the HBNM, the aim of using the hybrid displacement variational formulation in the present HRBNM is to gain a truly boundary-type meshless method by localizing the integration domain to a regular sub-domain, other than to gain the symmetrical system of equations. The rest of this paper is organized as follows. Section 2 describes the radial basis point interpolation scheme. In Section 3, the HRBNM is developed. Then, a detailed numerical implementation is provided in the next section. Section 5 gives some numerical results. Finally, the conclusion is presented in Section 6.
2. Radial basis point interpolation Let O be a bounded domain of R2 with smooth boundary G. A generic point in R2 is denoted by x ¼ ðx1 ; x2 Þ or y ¼ ðy1 ; y2 Þ. Let n ¼ ðn1 ; n2 Þ denote the outward normal direction on G. It is assumed that G is the union of piecewise smooth segments Gi , i ¼ 1; 2; . . . ; NG . In contrast to the BRPIM, the point interpolation in the present method is independently constructed on each segment, respectively, so that the discontinuity at corners can be avoided. As in the traditional BEM formulation, the point interpolation for u and its normal derivative q ¼ @u=@n on each Gi can be constructed independently. In the following interpolation scheme, let the variable v denote u or q for simplicity. In an interpolation domain with Ns scattered nodes, the function v can be approximated as vðsÞ vh ðsÞ ¼
Ns X
aj Bj ðsÞ ¼ BT ðsÞa
ð1Þ
j¼1
Gi ; where s is a parametric coordinate on BT ðsÞ ¼ ½B1 ðsÞ; B2 ðsÞ; . . . ; BNs ðsÞ; Bj ðsÞ is a RBF; aT ¼ ½a1 ; a2 ; . . . ; aNs ; aj are normally obtained by collocation; i.e. via solving vðsk Þ ¼ BT ðsk Þa;
1rkrNs
ð2Þ
Notice that the theory of RBFs provides a solid theoretical foundation to approximate v by vh. If the condition vðsk Þ ¼ vh ðsk Þ, 1rkrNs , is imposed, then the system of Eq. (2) is uniquely solvable providing that the coefficient matrix BT0 ¼ ½BT ðs1 Þ; BT ðs2 Þ; . . . ; BT ðsNs Þ
ð3Þ
is non-singular. It is well known that positive definiteness of the matrix is sufficient to guarantee its invertibility. However, most of the globally defined RBFs are only conditionally positive definite [30]. To guarantee the unique solvability of the interpolation problem and achieve polynomial precision, a polynomial term has to be added to Eq. (1), that is vðsÞ vh ðsÞ ¼
Ns X j¼1
aj Bj ðsÞ þ
m X i¼1
bi pi ðsÞ ¼ BT ðsÞa þ pT ðsÞb
ð4Þ
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
where pi ðsÞ is a monomial; m is the number of monomials; T bi is an unknown coefficient; b ¼ ½b1 ; b2 ; . . . ; bm and pT ðsÞ ¼ ½p1 ðsÞ; p2 ðsÞ; . . . ; pm ðsÞ. The effectiveness and accuracy of the interpolation depends on the choice of the RBFs. A number of different types of RBFs have been proposed and may be used for this purpose. Three usually used RBFs are: (1) Multiquadric (MQ): Bj ðsÞ ¼ ½ðs sj Þ2 þ C 2 Z
ð5Þ
(2) Gaussian (EXP): Bj ðsÞ ¼ exp½bðs sj Þ2
ð6Þ
(3) Thin plate spline (TPS): Bj ðsÞ ¼ ðs sj Þ2 lnjs sj j
ð7Þ
where C, Z and b are shape parameters. For the MQ and TPS, the following linear polynomial basis is required: pT ðsÞ ¼ ½1; s;
m¼2
ð8Þ
But the EXP can be implemented without any additional polynomial terms as it is positive definite [31]. Let the approximation represented by Eq. (4) pass through all the Ns boundary nodes, namely v^ ¼ B0 a þ p0 b ð9Þ where v^ T ¼ ½vðs1 Þ; vðs2 Þ; . . . ; vðsNs Þ, B0 is given in Eq. (3) and pT0 ¼ ½pT ðs1 Þ; pT ðs2 Þ; . . . ; pT ðsNs Þ. In order to guarantee unique solutions, the following constraints are imposed:
the functional P given by Z Z 1 P¼ u;I u;J dO qu dG 2 Gq
1275
ð17Þ
O
with the boundary condition (15). The variational formulation (17) should also satisfy the subsidiary compatibility conditions ð18Þ u~ ¼ u on G where u~ denotes the boundary potential field and u is the potential in the domain but very close to the boundary. Introduce the compatibility condition (18) in the variational expression (17) in terms of a Lagrange multiplier l. It can be concluded that the Lagrange multiplier l represents the boundary normal flux q~ [28]. As a result, the modified variational principle can be written as Z Z Z 1 ~ uÞ ~ dG P¼ u;I u;J dO qðu q u~ dG ð19Þ 2 O G Gq By carrying out the first variation of Eq. (19), we obtain Z Z Z ~ du dG ðu uÞ ~ dq~ dG dP ¼ u;II du dO þ ðq qÞ O G G Z þ ðq~ qÞdu~ dG
ð20Þ
Gq
Then, one can get the following integral equations via setting dP to zero: Z Z ~ du dG u;II du dO ¼ 0 ðq qÞ ð21Þ G
Z
O
~ dq~ dG ¼ 0 ðu uÞ
ð22Þ
G
pT0 a ¼ 0
ð10Þ
Solving Eqs. (9) and (10), we obtain ^ a ¼ Sa v; b ¼ Sb v^
ð11Þ
ð12Þ
where the shape function vector is
UðsÞ ¼ BT ðsÞSa þ pT ðsÞSb ¼ ½j1 ðsÞ; j2 ðsÞ; . . . ; ji ðsÞ; . . . ; jNs ðsÞ
ðq~ qÞdu~ dG ¼ 0
ð23Þ
Gq
1 T 1 where Sb ¼ ½pT0 B1 and Sa ¼ B1 0 p 0 p 0 B0 0 ½I p0 Sb . Then substituting Eq. (11) into Eq. (4) yields
vh ðsÞ ¼ ½BT ðsÞSa þ pT ðsÞSb v^ ¼ UðsÞv^
Z
ð13Þ
Obviously, the shape function UðsÞ obtained through the above process possesses the property of a delta function.
Since Eq. (23) will be satisfied if the flux boundary condition, q~ ¼ q, is imposed, we will ignore it in what follows. We can show that Eqs. (21) and (22) should be satisfied in any sub-domain Os , which is bounded by Gs and Ls (see Fig. 1). As in the LBIE method and the LRPIM, Eqs. (21) and (22) can be replaced by the following weak formulations: Z Z ~ i dG ðq qÞw u;II wi dO ¼ 0 ð24Þ Gis þLis
Z Gis þLis
Ois
~ i dG ¼ 0 ðu uÞw
3. The hybrid radial boundary node method 3.1. The boundary variational principle We consider a mixed boundary value problem in potential theory
DuðxÞ ¼ 0 for x 2 O
ð14Þ
uðxÞ ¼ uðxÞ
ð15Þ
for x 2 Gu
qðxÞ ¼ u;I ðxÞnI ¼ qðxÞ
for x 2 Gq
ð16Þ
where ð Þ;I denotes @ð Þ=@xI ; u and q are the given potential and normal flux on the essential boundary Gu and on the flux boundary Gq ¼ G Gu, respectively. Based on the classical one-field variational principle, the solution of a potential problem is the function u which minimizes
Fig. 1. The sub-domain Ois for the boundary node yi .
ð25Þ
ARTICLE IN PRESS 1276
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
where wi is the weight function or the test function; Ois denotes a sub-domain, which is located entirely inside O and contains the boundary node yi ; Gis and Lis are the boundary of Ois . Notice that in the above equations, the shape and dimension of the sub-domains Ois could be arbitrary. Thus we can choose a simple regular shape for Ois. Since a circle is the simplest regularly shaped sub-domain in R2 , we choose the sub-domain Ois as the intersection of the domain O and a circle centered at the boundary node yi (see Fig. 1).
In Eqs. (24) and (25), u~ and q~ on Gis can be represented by the point interpolation, namely 8 Ns X > > ~ > uðsÞ jj ðsÞuj ; ¼ > > > < j¼1 xðsÞ 2 G ð26Þ Ns X > > > qðsÞ ~ j ðsÞq ; ¼ > j j > > : j¼1 where jj ðsÞ is the contribution from the boundary node yj to the evaluation point xðsÞ, which is given by Eq. (13) and is not equal to zero in the interpolation domain of the j-th node only. It is worth noting that in the HRBNM, since the shape functions given by Eq. (13) have the delta function property, uj and qj of Eq. (26) are equal to the nodal values uðyj Þ and qðyj Þ, respectively. In the HBNM where the MLS is used, the shape functions lack the delta function property, which may lead to difficulties in imposition of boundary conditions, so that uj and qj are approximate values of the nodal values uðyj Þ and qðyj Þ, respectively. This is the main difference between the HRBNM and the HBNM. Bear in mind that u~ and q~ on the local boundary Lis are not known a priori. To get rid of these unknowns in the integral over Lis , we should select the weight function wi such that no integration is needed along Lis . As the hybrid boundary element model [28], the domain variables u and q can be approximated as 8 N X > > > uðxÞ ¼ U ðx; yj Þsj ; > > > < j¼1 x2O ð27Þ N X > > > qðxÞ ¼ Q ðx; y Þ s ; > j j > > : j¼1 where sj is the unknown parameter; N is the total number of boundary nodes; U ðx; yÞ ¼ ð1=2pÞlnjx yj is the fundamental solution of the Laplace operator; and Q ¼ @U =@n. 3.3. System of equations As u is represented by Eq. (27), the last integral in Eq. (24) vanishes if we exclude the node yi from the sub-domain Ois at which the singularity occurs. Indeed, this integral term is only attributed to the principal diagonal of the matrix. And it will be taken into account when calculating the boundary integrals. Thus by substituting Eqs. (26) and (27) into Eqs. (24) and (25), we obtain N Z N Z X X Q ðx; yj Þsj wi ðxÞ dG ¼ jj ðxðsÞÞqj wi ðxÞ dG ð28Þ Gis
N Z X j¼1
Gis
Gis
j¼1
U ðx; yj Þsj wi ðxÞ dG ¼
Q r ¼ Hq
ð30Þ
Ur ¼ Hu
ð31Þ
where Z Q ðx; yj Þwi ðxÞ dG Qij ¼
ð32Þ
Gis
Uij ¼
3.2. Variables interpolation
j¼1
Applying the above equations to all boundary nodes yields
N Z X j¼1
Gis
jj ðxðsÞÞuj wi ðxÞ dG
ð29Þ
Hij ¼
Z Gis
Z Gis
U ðx; yj Þwi ðxÞ dG
ð33Þ
jj ðsÞwi ðxÞ dG
ð34Þ
rT ¼ ½s1 ; s2 ; . . . ; sN
ð35Þ
q T ¼ ½q1 ; q2 ; . . . ; qN
ð36Þ
u T ¼ ½u1 ; u2 ; . . . ; uN
ð37Þ
The matrix U is a non-singular square matrix, thus it is invertible and the unknown vector r can be obtained from Eq. (31) as
r ¼ U1 Hu
ð38Þ
then substituting Eq. (38) into Eq. (30) yields QU1 Hu ¼ Hq
ð39Þ
Since the point interpolation satisfies the delta function property, boundary conditions in the present HRBNM can be applied as easily as in the traditional BEM. For a well-defined problem, either u or q is known at each boundary node. Thus the system equation (39) can be solved in a standard way to get the unknown boundary values. Then the unknown vector r is obtained via Eq. (38).
4. Numerical aspects 4.1. Integration = Gis , and In Eqs. (32) and (33), the integrals are regular when yj2 these integrals can be evaluated by usual Gaussian quadrature. If yj 2 Gis , a special logarithmic weighted numerical integration formula can be employed to compute the integrals in Eq. (33). For Eq. (32), the integral is difficult to compute when yj 2 Gis . In the situation, if iaj, namely the region Gis includes more than one boundary node, it is very hard to obtain this 1=r type singular integrals. To avoid this drawback, we assume that the region Gis covers only one boundary node in this work. Consequently, only the main diagonal terms of the matrix Q involve a strongly singular integral. This integral can be obtained by applying a uniform potential. More specifically, from Eq. (39) we have QU1 HI ¼ 0
ð40Þ
where I is a unit vector. As a result, the value of the main diagonal terms of matrix Q can be easily evaluated by the off-diagonal terms using Eq. (40). It should be pointed out that the uniform potential can be utilized in this research, since the shape function obtained by the point interpolation possesses the delta function property.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
4.2. Computation of potentials and its derivatives Once the unknowns u and q are found, the values of potential u~ and its normal derivative q~ at any boundary point could be evaluated using Eq. (26). The potential and its derivatives at an internal point may be computed simply using Eq. (27). This scheme can avoid further integrations. However, it has the drawback of serious ‘boundary layer effect’, i.e. the accuracy of the result near the boundary is very sensitive to the proximity of the interior points near the boundary. To overcome this problem, Zhang et al. proposed an adaptive integration technique in the HBNM [15]. This technique is proven to be very accurate [15,17,18]. We calculate the internal values of the potential by the traditional BIEs, i.e. Z Z ~ ~ dG Q ðx; yÞuðyÞ dG uðxÞ ¼ U ðx; yÞqðyÞ G
¼
G
NG Z X i¼1
~ U ðx; yÞqðyÞ dG
NG Z X
Gi
i¼1
~ Q ðx; yÞuðyÞ dG
ð41Þ
Gi
where NG denotes the number of the segments which compose the whole boundary. Then the derivatives can be computed by differentiating u at the internal points, namely Z Z ~ ~ u;I ðxÞ ¼ U;I ðx; yÞqðyÞ dG Q;I ðx; yÞuðyÞ dG G
¼
G
NG Z X i¼1
Gi
~ U;I ðx; yÞqðyÞ dG
NG Z X i¼1
Gi
~ Q;I ðx; yÞuðyÞ dG;
I ¼ 1; 2 ð42Þ
Because every segment can be represented by a unit sector in the parametric space, the integrals on each segment in Eqs. (41) and (42) can be evaluated easily. Here, an adaptive technique is introduced to compute these integrals on the segment. In this technique, the unit sector is first divided into four equal quarters (see Fig. 2). Then for each quarter, we calculate the quarter length l
1277
and the distance between the evaluation point and the center of the quarter d. If l is smaller than d, this quarter is considered as a regular integration segment. Otherwise, it will be further divided into sub-quarters, 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. We emphasize 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 found that all integrations are calculated along the boundary only and that no boundary elements are used either for interpolation or for integration. Thus, the present HRBNM is a truly boundary-type meshless method.
4.3. Domains of integration and interpolation For carrying out numerical integration, the local integration domain Gis should be defined. It can be found that the size of Gis will affect the performance of the HRBNM. For the boundary node yi , the size of the corresponding local integration domain Gis is defined as ri ¼ adi
ð43Þ
where a is a scaling parameter, di is the minimum distance between the node yi and adjacent nodes. In this work, Gauss quadrature is performed to calculate the numerical integrations. For each Gauss quadrature point, point interpolation is employed to obtain the integrand. Hence, an interpolation domain for a Gauss quadrature point is required to do the point interpolation. The size of interpolation domain may be written as dmi ¼ bdi
ð44Þ
where b is a coefficient chosen. We will investigate the influences of a and b by numerical examples.
Fig. 2. Sub-dividing of the segment Gi at an evaluation point x.
Fig. 3. Arrangement of boundary nodes and test nodes (a) Example 1; (b) Example 2 (the star and the dot, respectively, represent boundary nodes and test nodes).
ARTICLE IN PRESS 1278
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
4.4. Weight functions The weight function wi should monotonically decrease with the distance and it must vanish along the boundary Lis. A variety of weight functions have been investigated in the past for meshless methods. In this research, both Gaussian and spline weight functions are considered. The Gaussian weight function corresponding to the boundary node yi can be expressed by 8 2 2 ^ > > < exp½ðd i =ci Þ exp½ðri =ci Þ ; 0rd^ i rri 2 1 exp½ðri =ci Þ wi ðxÞ ¼ ð45Þ > > : 0; d^ Zr i
i
segment boundary, so the problem of the discontinuity at corners can be avoided. However, if no node is arranged in the corner of the boundary, the node value of variable is needed extrapolation for points between the corner and its adjacent node, and this may produce a not very accurate computed value for variable at the corner. As q ¼ @u=@n cannot be uniquely defined at corners, a simple approach proposed here to avoid this problem is to assume there are two nodes very near each other but which belong to different segments. This approach is very easy to implement and is employed in the following numerical experiments.
5. Numerical experiments
Here, ri is given by Eq. (43); d^ i is the distance between a sampling point x and the node yi ; and ci is a constant controlling the shape of the weight function wi . It should be pointed out that ri can be utilized as the size of the support of the weight function wi , since wi must be zero on Lis . Besides, a spline weight function is defined as ( 1 ðd^ i =ri Þ2 þ 8ðd^ i =ri Þ3 3ðd^ i =ri Þ4 ; 0rd^ i rri ð46Þ wi ðxÞ ¼ 0; d^ i Zri
To demonstrate the validity of our method, we present some numerical experiments. In all examples, the following relative error is defined to evaluate the present method performance: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2ffi u Nt u1 X uðeÞ uðnÞ t i i ð47Þ eðuÞ ¼ Nt i¼1 uðeÞ
4.5. Treatment of corners
5.1. Parameters study
In the present HRBNM, the boundary variables are interpolated by the radial basis point interpolation on the independent smooth
In the HRBNM, there exist some parameters which can affect the accuracy of solutions. However, as in many meshless methods,
i
where Nt is the number of test nodes inside the domain; uðeÞ and i refer to the nodal values obtained by the exact solutions and uðnÞ i numerical method, respectively.
Fig. 4. Effect of the parameter Z of MQ radial function (a) Example 1; (b) Example 2.
Fig. 5. Effect of the parameter c0 of MQ radial function (a) Example 1; (b) Example 2.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
the absence of the knowledge how to select properly these free parameters is a drawback of this method. In this section, these parameters will be studied numerically. Example 1. A harmonic problem whose exact solution is uðx1 ; x2 Þ ¼ sinðx1 Þsinhðx2 Þ
in O ¼ fðx1 ; x2 Þ : 1rx1 ; x2 r1g
Example 2. A harmonic problem with exact solution uðr; yÞ ¼ 3 þ r 2 cosð2yÞ
in a peanut shape domain O: 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > > x ¼ x ð y Þ ¼ cosð2yÞ þ 1:1 sin2 ð2yÞcosy; 1 > < 1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > cosð2 ¼ x ð y Þ ¼ y Þ þ 1:1 sin2 ð2yÞsiny; x > 2 2 > > :
1279
0ryr2p
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ x21 þ x22 is the distance from the origin to any point ðx1 ; x2 Þ 2 O and y ¼ tan1 ðx2 =x1 Þ.
Fig. 6. Effect of the parameter b of EXP radial function (a) Example 1; (b) Example 2.
Fig. 7. Effect of the parameter a of the integration domain (a) Example 1; (b) Example 2.
Fig. 8. Effect of the parameter b of the interpolation domain (a) Example 1; (b) Example 2.
ARTICLE IN PRESS 1280
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
Both examples are subjected to the Dirichlet boundary conditions. The distribution of boundary nodes and test nodes are shown in Fig. 3.
5.1.1. RBFs A variety of RBFs have been investigated in the past for meshless methods. To obtain a better understanding of the proposed approach, a detailed numerical study of three RBFs has been implemented in this work. These are MQ, EXP and TPS radial functions, given by Eqs. (5)–(7), respectively. However, some parameters (C and Z for MQ, b for EXP) need to be determined for MQ radial function and EXP radial function. These parameters will affect the performance of the radial functions. We first investigate the effect of C and Z for the MQ radial function. It is found that nodal spacing affects the choice of the shape parameter C [2]. As a result, C is defined as C ¼ c 0 di
There is only one parameter b in the EXP radial function. Results of relative errors for different b are obtained and depicted in Fig. 6. From this figure, it can be observed that b 2 ½0:2; 0:8 leads better results. For the TPS radial function, we obtain eðuÞ ¼ 0:02057, eðu;x1 Þ ¼ 0:02058 and eðu;x2 Þ ¼ 0:02222 for Example 1, and eðuÞ ¼ 0:02414, eðu;x1 Þ ¼ 0:02418 and eðu;x2 Þ ¼ 0:02153 for Example 2. Comparing
ð48Þ
where c0 is a dimensionless shape parameter, di is the minimum distance between the node yi and adjacent nodes. It is also found that the parameter Z has great influence on the performance of the radial basis point interpolation than that of the parameter C. Therefore, Z is firstly studied with c0 fixed at 0.5, 1.0, 1.5 and 2.0. Results for various values of Z are shown in Fig. 4. It can be found that Z 2 ½0:5; 1Þ [ ð1; 2Þ yields a satisfied result. Errors for different values of c0 are plotted in Fig. 5. It can be observed that c0 2 ½1:0; 5:0 can lead to satisfactory results.
Fig. 11. Heat conduction in a quadrant of an annular.
Fig. 9. Effect of the parameter g of the Gauss weights (a) Example 1; (b) Example 2.
Fig. 10. Convergence in relative error by using the Spline and Gauss weights ðg ¼ 1:5Þ (a) Example 1; (b) Example 2.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
the results of errors between MQ, EXP and TPS, we can see that the computational accuracy is little influenced by RBFs. For the following problems, the MQ radial function with parameters Z ¼ 0:5 and c0 ¼ 2:0 will be used. 5.1.2. Integration domain parameter a To study how the local sub-domain size affects the solution accuracy of the present method, relative errors with respect to the parameter a is shown in Fig. 7. It should be mentioned here that the union of fGis gN i¼1 does not cover the whole bounding surface when ao0:50. Fig. 7 indicates a 2 ½0:80; 0:99 leads to a satisfactory result. For convenience, a ¼ 0:90 will be used in the future computations. 5.1.3. Interpolation domain parameter b The size of interpolation domain of a quadrature point is decided by the parameter b in Eq. (44). Fig. 8 displays the relationship between the parameter b and relative errors. The results show that b 2 ½2:5; 4:5 yields a satisfactory result. For the following calculations, b is chosen as 3.0.
where g is a dimensionless shape parameter. Relative errors for different values of g are plotted in Fig. 9. It can be found that the errors change with g and the results for g 2 ½1:0; 2:0 are very good. When the parameter is too small ðgo1:0Þ or too big ðg42:0Þ, the errors increase. Both the Spline weight function and the Gaussian weight function are used to study the convergence of the HRBNM. The convergence results are shown in Fig. 10, where three kinds of nodal arrangement of 16, 32 and 64 uniformly distributed boundary nodes are employed. The results show that the convergence rates of the Gaussian weights are higher than those of the Spline weights.
5.2. Heat conduction A heat conduction problem in Fig. 11 is studied here, u represents the temperature field. Because of symmetry, only onequarter of the structure is modeled. For this problem, the boundary conditions are uða; yÞ ¼ Ta cosð4yÞ;
5.1.4. Weight functions In the Gaussian weight function, the parameter ci will influence the performance of the HRBNM. This parameter is defined as ci ¼ gdi
ð49Þ
Fig. 12. u and q on BC of the plate ðr ¼ 2Þ.
1281
@u ¼0 @y
uðb; yÞ ¼ Tb cosð4yÞ
on y ¼ 0 and y ¼
p 2
The analytical solution for this problem is [32] 4 4 a4 b4 r a4 r b4 Tb 4 4 Ta 4 4 cosð4yÞ uðr; yÞ ¼ 8 8 b a a r b r
The parameters are chosen as a ¼ 1, b ¼ 2, Ta ¼ 10 and Tb ¼ 5 in computation. Besides, we applied 32 boundary nodes (four regularly distributed nodes on AB and CD, 16 regularly distributed nodes on BC, and eight regularly distributed nodes on AD). Fig. 12 displays a comparison between the present numerical results with the exact solutions for u and q along the curved boundary BC. The results of the HBNM are also presented in this figure for comparison. It can be found that excellent agreement between the analytical and numerical solutions is achieved. The numerical results for the temperature u and its derivatives along the arc r ¼ 1:5 obtained by the HRBNM, the HBNM and the BRPIM are shown in Fig. 13 together with the analytical solutions. It is clearly shown that the numerical solutions for all approaches are seen to capture the behavior of the exact solutions very well for points far away from the boundary, while at points close to the boundary, the BRPIM has a little lower accuracy than the HRBNM and the HBNM.
Fig. 13. Results of (a) u and (b) u;x1 on the arc r ¼ 1:5.
ARTICLE IN PRESS 1282
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
Fig. 14. Flow in a confined porous layer. Fig. 15. Potential u and its derivatives at the section x2 ¼ 2:5.
5.3. Fluid flow In this example, we solve the problem of fluid flow in a thin layer of a porous geological material (see Fig. 14). The thickness of the layer erL or b. The layer is contained between two impervious strata. The flow takes place in the plane of the porous layer and the hydraulic conductivity of the layer is isotropic and accounts for any non-uniformity in the fluid velocity distribution over its thickness. The boundary conditions for the hydraulic potential are uðx1 ; 0Þ ¼ 0;
uðx1 ; bÞ ¼ 0;
uð0; x2 Þ ¼ 0;
uðL; x2 Þ ¼ P0
where P0 is a constant pressure head imposed on x1 ¼ L. It can be noted that the potential u is composed of the datum potential ud and the pressure potential up . Because the porous layer is located on the plane x3 ¼ 0, we can set the datum at x3 ¼ 0 and u can be identified directly as the pressure potential. With these conditions, the exact solution of this problem can be given by [32] ( ) 1 2P0 X 1 þ ð1Þnþ1 sinhðnpx1 =bÞsinðnpx2 =bÞ uðx1 ; x2 Þ ¼ p n¼1 n sinhðnpL=bÞ The parameters that are used in our analysis are L ¼ 10, b ¼ 5 and P0 ¼ 10. The exact solutions and numerical solutions of the present HRBNM for u and its derivatives are depicted in Fig. 15. In this study, the boundary of the domain is discretized by 48 uniformly distributed nodes (16 nodes on AB and CD, and eight nodes on AD and BC). From this figure, we can see that the numerical solutions are in excellent agreement with the exact solutions. To show the convergence of the present method, regularly distributed 24, 48 and 96 nodes are used. The ratio of the number of nodes on AB and CD to that on AD and BC is kept constant and equal to 2. The results of convergence of @u=@x1 are shown in Fig. 16. In the figure, the convergences of the HBNM and the BRPIM are also shown for comparison. It can be found that the solution accuracy of the present approach is slightly higher than that of the HBNM and the BRPIM.
6. Conclusions A hybrid radial boundary node method is presented in this paper. It is based on a modified variational principle, which involves three independent variables, two of them defined on the boundary and the other inside the domain. The domain variables
Fig. 16. Convergence for the fluid flow.
are interpolated by the classical fundamental solutions and the boundary variables by the radial basis point interpolation. The shape functions constructed by this point interpolation possess the delta function property, thus boundary conditions can be applied directly. The HRBNM is a truly boundary-type meshless method, in which the interpolation does not use any element and the integral does not require any background mesh. Thus compared with the BNM, the BRPIM and the hybrid BRPIM, the HRBNM is truly meshless. Besides, in contrast with the LBIE method and the LRPIM, the HRBNM only requires a nodal data structure on the bounding surface of the domain to be solved. The differences among the HRBNM and the HBNM or the BNM are as follows. In the HRBNM, the basic unknown quantities are the real solutions of the nodal variables, but in the HBNM and the BNM, the basic unknown quantities are approximations of the nodal variables. In the HRBNM, boundary conditions can be directly and easily implemented, but in the HBNM and the BNM, boundary conditions are applied after they are transformed into their approximations on the boundary nodes with the MLS approximation. Hence, the HRBNM is expected to have higher computational precision.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
Some strongly singular integration in the HRBNM is circumvented with the help of a uniform potential field. Besides, an adaptive integration scheme is employed for alleviating the ‘boundary layer effect’. Numerical examples are presented for several potential problems. The numerical results have demonstrated the accuracy, convergence and effectiveness of the HRBNM. The effects of shape parameters of RBFs, sizes of interpolation domains and integration domains, and weight functions are investigated in detail. The present method can easily be implemented to other problems that can be solved by the BEM.
Acknowledgment The authors would like to thank the anonymous reviewers whose constructive comments improve the quality of this article.
References [1] Belytschko T, Krongauz Y, Organ D, Fleming M. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng 1996;139:3–47. [2] Liu GR, Gu YT. An introduction to meshfree methods and their programming. Dordrecht: Springer; 2005. [3] Lancaster P, Salkauskas K. Surface generated by moving least squares methods. Math Comput 1981;37:141–58. [4] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse element. Comput Mech 1992;10:307–18. [5] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [6] Atluri SN, Shen SP. The meshless local Petrov–Galerkin (MLPG) method: a simple & less-costly alternative to the finite element and boundary element methods. Comput Model Eng Sci 2002;3:11–51. [7] Atluri SN, Sladek J, Sladek V, Zhu T. The local boundary integral equation (LBIE) and it’s meshless implementation for linear elasticity. Comput Mech 2000;25:180–98. [8] Mukherjee YX, Mukherjee S. The boundary node method for potential problems. Int J Numer Methods Eng 1997;40:797–15. [9] Zhang J, Yao Z, Li H. A hybrid boundary node method. Int J Numer Methods Eng 2002;53:751–63. [10] Kitipornchai S, Liew KM, Cheng Y. A boundary element-free method (BEFM) for three-dimensional elasticity problems. Comput Mech 2005;36:13–20.
1283
[11] Peng M, Cheng Y. A boundary element-free method (BEFM) for twodimensional potential problems. Eng Anal Bound Elem 2009;33:77–82. [12] Li X, Zhu J. A Galerkin boundary node method and its convergence analysis. J Comput Appl Math 2009;230:314–28. [13] Li X, Zhu J. A Galerkin boundary node method for biharmonic problems. Eng Anal Bound Elem 2009;33:858–65. [14] Zhang J, Yao Z. Meshless regular hybrid boundary node method. Comput Model Eng Sci 2001;2:301–18. [15] Zhang J, Tanaka M, Matsumoto T. Meshless analysis of potential problems in three dimensions with the hybrid boundary node method. Int J Numer Methods Eng 2004;59:1147–68. [16] Zhang J, Yao Z. The regular hybrid boundary node method for threedimensional linear elasticity. Eng Anal Bound Elem 2004;28:525–34. [17] Miao Y, Wang Y. Development of hybrid boundary node method in two dimensional elasticity. Eng Anal Bound Elem 2005;29:703–12. [18] Yan F, Wang Y, Tham LG, Cheung YK. Dual reciprocity hybrid boundary node method for 2-D elasticity with body force. Eng Anal Bound Elem 2008;32: 713–75. [19] Yan F, Wang Y, Miao Y, Cheung YK. Dual reciprocity hybrid boundary node method for free vibration analysis. J Sound Vib 2009;321:1036–57. [20] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 2002;54:1623–48. [21] Liu GR, Gu YT. A local radial point interpolation method (LRPIM) for free vibration analysis of 2-D solids. J Sound Vib 2001;246:29–46. [22] Liu GR, Yan L, Wang JG, Gu YT. Point interpolation method based on local residual formulation using radial basis functions. Struct Eng Mech 2002;14:713–32. [23] Gu YT, Liu GR. A boundary radial point interpolation method (BRPIM) for 2-D structural analyses. Struct Eng Mech 2002;15:535–50. [24] Liu GR, Gu YT. Boundary meshfree methods based on the boundary point interpolation methods. Eng Anal Bound Elem 2004;28:475–87. [25] Wu YL, Liu GR. A meshfree formulation of local radial point interpolation method (LRPIM) for incompressible flow simulation. Comput Mech 2003;30:355–65. [26] Liu GR, Gu YT. Comparisons of two meshfree local point interpolation methods for structural analyses. Comput Mech 2002;29:107–21. [27] Dumont NA. The hybrid boundary element method. In: Brebbia CA, Wendland WL, Kuhn G, editors. Boundary elements IX, vol. 1. Berlin: Springer; 1988. [28] DeFiigueiredo TGB. A new boundary element formulation in engineering. Berlin: Springer; 1991. [29] Gu YT, Liu GR. Hybrid boundary point interpolation method and their coupling with the element free Galerkin method. Eng Anal Bound Elem 2003;27:905–17. [30] Madych WR, Nelson SA. Multivariate interpolation and conditionally positive definite functions II. Math Comput Math Comput 1990;54:211–30. [31] Powell MJD. The theory of radial basis function approximation in 1990. In: Light W, editor. Advances in numerical analysis, vol. II. Oxford: Oxford Science Publications; 1997. [32] Selvadurai APS. Partial differential equations in mechanics, vol. 1. Berlin: Springer; 2000.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A hybrid radial boundary node method based on radial basis point interpolation Xiaolin Li a,, Jialin Zhu a, Shougui Zhang b a b
College of Mathematics and Physics, Chongqing University, Chongqing 400044, PR China College of Mathematics and Computer Science, Chongqing Normal University, Chongqing 400047, PR China
a r t i c l e in fo
abstract
Article history: Received 14 January 2009 Accepted 4 June 2009 Available online 17 July 2009
The hybrid boundary node method (HBNM) retains the meshless attribute of the moving least squares (MLS) approximation and the reduced dimensionality advantages of the boundary element method. However, the HBNM inherits the deficiency of the MLS approximation, in which shape functions lack the delta function property. Thus in the HBNM, boundary conditions are implemented after they are transformed into their approximations on the boundary nodes with the MLS scheme. This paper combines the hybrid displacement variational formulation and the radial basis point interpolation to develop a direct boundary-type meshless method, the hybrid radial boundary node method (HRBNM) for two-dimensional potential problems. The HRBNM is truly meshless, i.e. absolutely no elements are required either for interpolation or for integration. The radial basis point interpolation is used to construct shape functions with delta function property. So unlike the HBNM, the HRBNM is a direct numerical method 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. Some selected numerical tests illustrate the efficiency of the method proposed. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Meshless Hybrid radial boundary node method Radial basis function Boundary integral equations
1. Introduction In recent years, the meshless or meshfree methods have generated considerable interest among researchers for solving boundary value problems [1,2]. Compared to the finite element method (FEM) and the boundary element method (BEM), the main feature of this type of method is the absence of an explicit mesh, and the approximate solutions are constructed entirely based on a cluster of scattered nodes. Meshless methods have been found to have special advantages on problems to which the conventional mesh-based methods are difficult to be applied. These include problems with complicated boundary, moving boundary, large deformations, mesh adaptively, and so on. The moving least squares (MLS) as an approximation method has been introduced by Lancaster and Salkauskas [3]. Since the numerical approximations of the MLS start from a cluster of scattered nodes instead of meshes, there have many meshless methods based on the MLS scheme for the numerical solution of differential equations. Typical of them are the diffuse element method [4], the element-free Galerkin (EFG) method [5] and the meshless local Petrov–Galerkin (MLPG) method [6]. The EFG method uses a global weak form. Although no mesh is required in
Corresponding author. Fax: +86 23 65120520.
E-mail address: [email protected] (X. Li). 0955-7997/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2009.06.003
this method for variable interpolation, background cells are necessary for ‘energy’ integration in general. The MLPG approach is formulated on a local weak form. The main advantage of the MLPG approach over the FEM and the EFG method is that it does not need a ‘finite element mesh’, either for purposes of interpolation of the solution variables or for the integration of the ‘energy’. The MLS approximation has also been used in boundary integral equations (BIEs). Among them are the meshless local boundary integral equation (LBIE) method [7], the boundary node method (BNM) [8], the hybrid boundary node method (HBNM) [9], the boundary element-free method (BEFM) [10,11] and the Galerkin boundary node method (GBNM) [12,13]. The first boundary-type meshless method was the BNM. This method takes the advantages of both the BIE in dimension reduction and the MLS approximation in elements elimination. Nevertheless, background cells are required for integration. To avoid this hindrance, Zhu et al. proposed the LBIE method. It is equivalent to a sort of MLPG approach, which needs absolutely no domain and boundary elements. The LBIE method, however, is not strictly a boundary method since it requires evaluation of integrals over certain surfaces (called Ls in [7]) that can be regarded as ‘closure surfaces’ of boundary elements. To obtain a truly boundary-type meshless method, Zhang et al. [9] proposed the HBNM that combines the MLS approximation with the hybrid displacement variational formulation. In this approach, the integration is limited to a fixed local region, thus no cells are needed either
ARTICLE IN PRESS 1274
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
for the interpolation purposes or for integration process. The HBNM is a very promising method, which has been used for problems in potential theory [9,14,15], linear elasticity [16–18] and elastodynamics [19]. These MLS-based meshless methods have gained remarkable progress in solving a broad class of boundary value problems. However, since the MLS approximations lack the delta function property of the usual FEM and BEM shape functions, boundary conditions in these meshless techniques are difficult to be enforced. This problem becomes even more serious in the boundary-type meshless methods, since a large number of boundary conditions need to be satisfied. The strategy proposed in the BNM [8] involves a new definition of the discrete norm used for the construction of the MLS approximations, making the BNM computationally much more expensive than the traditional BEM. 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. However, they are not the real nodal variables, and thus boundary conditions could not be directly applied. The BNM and the HBNM are the indirect boundary-type meshless methods. In addition to the abovementioned strategy, Kitipornchai et al. [10] presented an improved MLS approximation that uses weighted orthogonal polynomials as basis functions to restore the delta function property of the MLS. The BEFM is based on the improved MLS, thus it is a direct boundary-type meshless method. Recently, via using an equivalent variational form of BIEs, another effective technique is developed in the GBNM to impose boundary conditions [12]. The radial point interpolation method (RPIM) [2,20] is another meshless interpolation technique proposed by Liu et al. The RPIM approximates its field variable by letting the problem function pass through each scattered node in its local support domain and radial basis functions (RBFs) are used for point interpolation. Different from the MLS approximation, the shape functions so constructed have the delta function property. As a result, the construction of a RPIM shape function is more efficient than the MLS procedures. Some meshless methods using the RPIM shape functions have been developed. Typical of them are the RPIM [20], the local radial point interpolation method (LRPIM) [21,22] and the boundary radial point interpolation method (BRPIM) [23,24]. As the MLPG approach, the LRPIM is a truly domain-type meshless method. The BRPIM is a boundary-type method, but it requires background cells for numerical integration. In these RPIM-based methods, boundary conditions can be enforced as easily as in the conventional FEM and BEM, thus they are much more efficient than the methods using MLS shape functions. Up to date, these approaches have been successfully implemented for problems of solid [20–24], fluid [25] and elastodynamics [26]. An interesting question arises here—Is there possibly a boundary-type meshless method that only involves the real solution of nodal variables and does not need any element either for interpolation or for integration? Such a method is a direct numerical approach, in which boundary conditions can be applied directly and easily, compared with the HBNM and the BNM; and it is a truly meshless method, which require no mesh for either interpolation or integration purpose, compared with the BNM and the BRPIM. The answer to this question is positive. The new method is called the hybrid radial boundary node method (HRBNM) and is the subject of this paper. The HRBNM combines the RPIM with the hybrid displacement variational formulation. Hybrid boundary element formulations are based on generalized variational principles. In the late eighties, two hybrid models have been developed: the hybrid stress boundary element model [27] and the hybrid displacement boundary element model [28].
The first is based on the Hellinger–Reissner principle with the stresses in the domain and displacements on the boundary as independent functions, and the second on a modified functional using three independent variables, two of them defined on the boundary and the other inside the domain. Both models use the classical fundamental solution to approximate the domain variable and thus allowing for the transfer of the domain integrals to the boundary. Besides, the resulting coefficient matrices are symmetric. Recently, with combination of the RPIM and the hybrid displacement boundary element model, a hybrid BRPIM has been developed [29]. The hybrid BRPIM results in symmetric matrices, but an underlying cell structure is still used for integration. As in the HBNM, the aim of using the hybrid displacement variational formulation in the present HRBNM is to gain a truly boundary-type meshless method by localizing the integration domain to a regular sub-domain, other than to gain the symmetrical system of equations. The rest of this paper is organized as follows. Section 2 describes the radial basis point interpolation scheme. In Section 3, the HRBNM is developed. Then, a detailed numerical implementation is provided in the next section. Section 5 gives some numerical results. Finally, the conclusion is presented in Section 6.
2. Radial basis point interpolation Let O be a bounded domain of R2 with smooth boundary G. A generic point in R2 is denoted by x ¼ ðx1 ; x2 Þ or y ¼ ðy1 ; y2 Þ. Let n ¼ ðn1 ; n2 Þ denote the outward normal direction on G. It is assumed that G is the union of piecewise smooth segments Gi , i ¼ 1; 2; . . . ; NG . In contrast to the BRPIM, the point interpolation in the present method is independently constructed on each segment, respectively, so that the discontinuity at corners can be avoided. As in the traditional BEM formulation, the point interpolation for u and its normal derivative q ¼ @u=@n on each Gi can be constructed independently. In the following interpolation scheme, let the variable v denote u or q for simplicity. In an interpolation domain with Ns scattered nodes, the function v can be approximated as vðsÞ vh ðsÞ ¼
Ns X
aj Bj ðsÞ ¼ BT ðsÞa
ð1Þ
j¼1
Gi ; where s is a parametric coordinate on BT ðsÞ ¼ ½B1 ðsÞ; B2 ðsÞ; . . . ; BNs ðsÞ; Bj ðsÞ is a RBF; aT ¼ ½a1 ; a2 ; . . . ; aNs ; aj are normally obtained by collocation; i.e. via solving vðsk Þ ¼ BT ðsk Þa;
1rkrNs
ð2Þ
Notice that the theory of RBFs provides a solid theoretical foundation to approximate v by vh. If the condition vðsk Þ ¼ vh ðsk Þ, 1rkrNs , is imposed, then the system of Eq. (2) is uniquely solvable providing that the coefficient matrix BT0 ¼ ½BT ðs1 Þ; BT ðs2 Þ; . . . ; BT ðsNs Þ
ð3Þ
is non-singular. It is well known that positive definiteness of the matrix is sufficient to guarantee its invertibility. However, most of the globally defined RBFs are only conditionally positive definite [30]. To guarantee the unique solvability of the interpolation problem and achieve polynomial precision, a polynomial term has to be added to Eq. (1), that is vðsÞ vh ðsÞ ¼
Ns X j¼1
aj Bj ðsÞ þ
m X i¼1
bi pi ðsÞ ¼ BT ðsÞa þ pT ðsÞb
ð4Þ
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
where pi ðsÞ is a monomial; m is the number of monomials; T bi is an unknown coefficient; b ¼ ½b1 ; b2 ; . . . ; bm and pT ðsÞ ¼ ½p1 ðsÞ; p2 ðsÞ; . . . ; pm ðsÞ. The effectiveness and accuracy of the interpolation depends on the choice of the RBFs. A number of different types of RBFs have been proposed and may be used for this purpose. Three usually used RBFs are: (1) Multiquadric (MQ): Bj ðsÞ ¼ ½ðs sj Þ2 þ C 2 Z
ð5Þ
(2) Gaussian (EXP): Bj ðsÞ ¼ exp½bðs sj Þ2
ð6Þ
(3) Thin plate spline (TPS): Bj ðsÞ ¼ ðs sj Þ2 lnjs sj j
ð7Þ
where C, Z and b are shape parameters. For the MQ and TPS, the following linear polynomial basis is required: pT ðsÞ ¼ ½1; s;
m¼2
ð8Þ
But the EXP can be implemented without any additional polynomial terms as it is positive definite [31]. Let the approximation represented by Eq. (4) pass through all the Ns boundary nodes, namely v^ ¼ B0 a þ p0 b ð9Þ where v^ T ¼ ½vðs1 Þ; vðs2 Þ; . . . ; vðsNs Þ, B0 is given in Eq. (3) and pT0 ¼ ½pT ðs1 Þ; pT ðs2 Þ; . . . ; pT ðsNs Þ. In order to guarantee unique solutions, the following constraints are imposed:
the functional P given by Z Z 1 P¼ u;I u;J dO qu dG 2 Gq
1275
ð17Þ
O
with the boundary condition (15). The variational formulation (17) should also satisfy the subsidiary compatibility conditions ð18Þ u~ ¼ u on G where u~ denotes the boundary potential field and u is the potential in the domain but very close to the boundary. Introduce the compatibility condition (18) in the variational expression (17) in terms of a Lagrange multiplier l. It can be concluded that the Lagrange multiplier l represents the boundary normal flux q~ [28]. As a result, the modified variational principle can be written as Z Z Z 1 ~ uÞ ~ dG P¼ u;I u;J dO qðu q u~ dG ð19Þ 2 O G Gq By carrying out the first variation of Eq. (19), we obtain Z Z Z ~ du dG ðu uÞ ~ dq~ dG dP ¼ u;II du dO þ ðq qÞ O G G Z þ ðq~ qÞdu~ dG
ð20Þ
Gq
Then, one can get the following integral equations via setting dP to zero: Z Z ~ du dG u;II du dO ¼ 0 ðq qÞ ð21Þ G
Z
O
~ dq~ dG ¼ 0 ðu uÞ
ð22Þ
G
pT0 a ¼ 0
ð10Þ
Solving Eqs. (9) and (10), we obtain ^ a ¼ Sa v; b ¼ Sb v^
ð11Þ
ð12Þ
where the shape function vector is
UðsÞ ¼ BT ðsÞSa þ pT ðsÞSb ¼ ½j1 ðsÞ; j2 ðsÞ; . . . ; ji ðsÞ; . . . ; jNs ðsÞ
ðq~ qÞdu~ dG ¼ 0
ð23Þ
Gq
1 T 1 where Sb ¼ ½pT0 B1 and Sa ¼ B1 0 p 0 p 0 B0 0 ½I p0 Sb . Then substituting Eq. (11) into Eq. (4) yields
vh ðsÞ ¼ ½BT ðsÞSa þ pT ðsÞSb v^ ¼ UðsÞv^
Z
ð13Þ
Obviously, the shape function UðsÞ obtained through the above process possesses the property of a delta function.
Since Eq. (23) will be satisfied if the flux boundary condition, q~ ¼ q, is imposed, we will ignore it in what follows. We can show that Eqs. (21) and (22) should be satisfied in any sub-domain Os , which is bounded by Gs and Ls (see Fig. 1). As in the LBIE method and the LRPIM, Eqs. (21) and (22) can be replaced by the following weak formulations: Z Z ~ i dG ðq qÞw u;II wi dO ¼ 0 ð24Þ Gis þLis
Z Gis þLis
Ois
~ i dG ¼ 0 ðu uÞw
3. The hybrid radial boundary node method 3.1. The boundary variational principle We consider a mixed boundary value problem in potential theory
DuðxÞ ¼ 0 for x 2 O
ð14Þ
uðxÞ ¼ uðxÞ
ð15Þ
for x 2 Gu
qðxÞ ¼ u;I ðxÞnI ¼ qðxÞ
for x 2 Gq
ð16Þ
where ð Þ;I denotes @ð Þ=@xI ; u and q are the given potential and normal flux on the essential boundary Gu and on the flux boundary Gq ¼ G Gu, respectively. Based on the classical one-field variational principle, the solution of a potential problem is the function u which minimizes
Fig. 1. The sub-domain Ois for the boundary node yi .
ð25Þ
ARTICLE IN PRESS 1276
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
where wi is the weight function or the test function; Ois denotes a sub-domain, which is located entirely inside O and contains the boundary node yi ; Gis and Lis are the boundary of Ois . Notice that in the above equations, the shape and dimension of the sub-domains Ois could be arbitrary. Thus we can choose a simple regular shape for Ois. Since a circle is the simplest regularly shaped sub-domain in R2 , we choose the sub-domain Ois as the intersection of the domain O and a circle centered at the boundary node yi (see Fig. 1).
In Eqs. (24) and (25), u~ and q~ on Gis can be represented by the point interpolation, namely 8 Ns X > > ~ > uðsÞ jj ðsÞuj ; ¼ > > > < j¼1 xðsÞ 2 G ð26Þ Ns X > > > qðsÞ ~ j ðsÞq ; ¼ > j j > > : j¼1 where jj ðsÞ is the contribution from the boundary node yj to the evaluation point xðsÞ, which is given by Eq. (13) and is not equal to zero in the interpolation domain of the j-th node only. It is worth noting that in the HRBNM, since the shape functions given by Eq. (13) have the delta function property, uj and qj of Eq. (26) are equal to the nodal values uðyj Þ and qðyj Þ, respectively. In the HBNM where the MLS is used, the shape functions lack the delta function property, which may lead to difficulties in imposition of boundary conditions, so that uj and qj are approximate values of the nodal values uðyj Þ and qðyj Þ, respectively. This is the main difference between the HRBNM and the HBNM. Bear in mind that u~ and q~ on the local boundary Lis are not known a priori. To get rid of these unknowns in the integral over Lis , we should select the weight function wi such that no integration is needed along Lis . As the hybrid boundary element model [28], the domain variables u and q can be approximated as 8 N X > > > uðxÞ ¼ U ðx; yj Þsj ; > > > < j¼1 x2O ð27Þ N X > > > qðxÞ ¼ Q ðx; y Þ s ; > j j > > : j¼1 where sj is the unknown parameter; N is the total number of boundary nodes; U ðx; yÞ ¼ ð1=2pÞlnjx yj is the fundamental solution of the Laplace operator; and Q ¼ @U =@n. 3.3. System of equations As u is represented by Eq. (27), the last integral in Eq. (24) vanishes if we exclude the node yi from the sub-domain Ois at which the singularity occurs. Indeed, this integral term is only attributed to the principal diagonal of the matrix. And it will be taken into account when calculating the boundary integrals. Thus by substituting Eqs. (26) and (27) into Eqs. (24) and (25), we obtain N Z N Z X X Q ðx; yj Þsj wi ðxÞ dG ¼ jj ðxðsÞÞqj wi ðxÞ dG ð28Þ Gis
N Z X j¼1
Gis
Gis
j¼1
U ðx; yj Þsj wi ðxÞ dG ¼
Q r ¼ Hq
ð30Þ
Ur ¼ Hu
ð31Þ
where Z Q ðx; yj Þwi ðxÞ dG Qij ¼
ð32Þ
Gis
Uij ¼
3.2. Variables interpolation
j¼1
Applying the above equations to all boundary nodes yields
N Z X j¼1
Gis
jj ðxðsÞÞuj wi ðxÞ dG
ð29Þ
Hij ¼
Z Gis
Z Gis
U ðx; yj Þwi ðxÞ dG
ð33Þ
jj ðsÞwi ðxÞ dG
ð34Þ
rT ¼ ½s1 ; s2 ; . . . ; sN
ð35Þ
q T ¼ ½q1 ; q2 ; . . . ; qN
ð36Þ
u T ¼ ½u1 ; u2 ; . . . ; uN
ð37Þ
The matrix U is a non-singular square matrix, thus it is invertible and the unknown vector r can be obtained from Eq. (31) as
r ¼ U1 Hu
ð38Þ
then substituting Eq. (38) into Eq. (30) yields QU1 Hu ¼ Hq
ð39Þ
Since the point interpolation satisfies the delta function property, boundary conditions in the present HRBNM can be applied as easily as in the traditional BEM. For a well-defined problem, either u or q is known at each boundary node. Thus the system equation (39) can be solved in a standard way to get the unknown boundary values. Then the unknown vector r is obtained via Eq. (38).
4. Numerical aspects 4.1. Integration = Gis , and In Eqs. (32) and (33), the integrals are regular when yj2 these integrals can be evaluated by usual Gaussian quadrature. If yj 2 Gis , a special logarithmic weighted numerical integration formula can be employed to compute the integrals in Eq. (33). For Eq. (32), the integral is difficult to compute when yj 2 Gis . In the situation, if iaj, namely the region Gis includes more than one boundary node, it is very hard to obtain this 1=r type singular integrals. To avoid this drawback, we assume that the region Gis covers only one boundary node in this work. Consequently, only the main diagonal terms of the matrix Q involve a strongly singular integral. This integral can be obtained by applying a uniform potential. More specifically, from Eq. (39) we have QU1 HI ¼ 0
ð40Þ
where I is a unit vector. As a result, the value of the main diagonal terms of matrix Q can be easily evaluated by the off-diagonal terms using Eq. (40). It should be pointed out that the uniform potential can be utilized in this research, since the shape function obtained by the point interpolation possesses the delta function property.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
4.2. Computation of potentials and its derivatives Once the unknowns u and q are found, the values of potential u~ and its normal derivative q~ at any boundary point could be evaluated using Eq. (26). The potential and its derivatives at an internal point may be computed simply using Eq. (27). This scheme can avoid further integrations. However, it has the drawback of serious ‘boundary layer effect’, i.e. the accuracy of the result near the boundary is very sensitive to the proximity of the interior points near the boundary. To overcome this problem, Zhang et al. proposed an adaptive integration technique in the HBNM [15]. This technique is proven to be very accurate [15,17,18]. We calculate the internal values of the potential by the traditional BIEs, i.e. Z Z ~ ~ dG Q ðx; yÞuðyÞ dG uðxÞ ¼ U ðx; yÞqðyÞ G
¼
G
NG Z X i¼1
~ U ðx; yÞqðyÞ dG
NG Z X
Gi
i¼1
~ Q ðx; yÞuðyÞ dG
ð41Þ
Gi
where NG denotes the number of the segments which compose the whole boundary. Then the derivatives can be computed by differentiating u at the internal points, namely Z Z ~ ~ u;I ðxÞ ¼ U;I ðx; yÞqðyÞ dG Q;I ðx; yÞuðyÞ dG G
¼
G
NG Z X i¼1
Gi
~ U;I ðx; yÞqðyÞ dG
NG Z X i¼1
Gi
~ Q;I ðx; yÞuðyÞ dG;
I ¼ 1; 2 ð42Þ
Because every segment can be represented by a unit sector in the parametric space, the integrals on each segment in Eqs. (41) and (42) can be evaluated easily. Here, an adaptive technique is introduced to compute these integrals on the segment. In this technique, the unit sector is first divided into four equal quarters (see Fig. 2). Then for each quarter, we calculate the quarter length l
1277
and the distance between the evaluation point and the center of the quarter d. If l is smaller than d, this quarter is considered as a regular integration segment. Otherwise, it will be further divided into sub-quarters, 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. We emphasize 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 found that all integrations are calculated along the boundary only and that no boundary elements are used either for interpolation or for integration. Thus, the present HRBNM is a truly boundary-type meshless method.
4.3. Domains of integration and interpolation For carrying out numerical integration, the local integration domain Gis should be defined. It can be found that the size of Gis will affect the performance of the HRBNM. For the boundary node yi , the size of the corresponding local integration domain Gis is defined as ri ¼ adi
ð43Þ
where a is a scaling parameter, di is the minimum distance between the node yi and adjacent nodes. In this work, Gauss quadrature is performed to calculate the numerical integrations. For each Gauss quadrature point, point interpolation is employed to obtain the integrand. Hence, an interpolation domain for a Gauss quadrature point is required to do the point interpolation. The size of interpolation domain may be written as dmi ¼ bdi
ð44Þ
where b is a coefficient chosen. We will investigate the influences of a and b by numerical examples.
Fig. 2. Sub-dividing of the segment Gi at an evaluation point x.
Fig. 3. Arrangement of boundary nodes and test nodes (a) Example 1; (b) Example 2 (the star and the dot, respectively, represent boundary nodes and test nodes).
ARTICLE IN PRESS 1278
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
4.4. Weight functions The weight function wi should monotonically decrease with the distance and it must vanish along the boundary Lis. A variety of weight functions have been investigated in the past for meshless methods. In this research, both Gaussian and spline weight functions are considered. The Gaussian weight function corresponding to the boundary node yi can be expressed by 8 2 2 ^ > > < exp½ðd i =ci Þ exp½ðri =ci Þ ; 0rd^ i rri 2 1 exp½ðri =ci Þ wi ðxÞ ¼ ð45Þ > > : 0; d^ Zr i
i
segment boundary, so the problem of the discontinuity at corners can be avoided. However, if no node is arranged in the corner of the boundary, the node value of variable is needed extrapolation for points between the corner and its adjacent node, and this may produce a not very accurate computed value for variable at the corner. As q ¼ @u=@n cannot be uniquely defined at corners, a simple approach proposed here to avoid this problem is to assume there are two nodes very near each other but which belong to different segments. This approach is very easy to implement and is employed in the following numerical experiments.
5. Numerical experiments
Here, ri is given by Eq. (43); d^ i is the distance between a sampling point x and the node yi ; and ci is a constant controlling the shape of the weight function wi . It should be pointed out that ri can be utilized as the size of the support of the weight function wi , since wi must be zero on Lis . Besides, a spline weight function is defined as ( 1 ðd^ i =ri Þ2 þ 8ðd^ i =ri Þ3 3ðd^ i =ri Þ4 ; 0rd^ i rri ð46Þ wi ðxÞ ¼ 0; d^ i Zri
To demonstrate the validity of our method, we present some numerical experiments. In all examples, the following relative error is defined to evaluate the present method performance: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2ffi u Nt u1 X uðeÞ uðnÞ t i i ð47Þ eðuÞ ¼ Nt i¼1 uðeÞ
4.5. Treatment of corners
5.1. Parameters study
In the present HRBNM, the boundary variables are interpolated by the radial basis point interpolation on the independent smooth
In the HRBNM, there exist some parameters which can affect the accuracy of solutions. However, as in many meshless methods,
i
where Nt is the number of test nodes inside the domain; uðeÞ and i refer to the nodal values obtained by the exact solutions and uðnÞ i numerical method, respectively.
Fig. 4. Effect of the parameter Z of MQ radial function (a) Example 1; (b) Example 2.
Fig. 5. Effect of the parameter c0 of MQ radial function (a) Example 1; (b) Example 2.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
the absence of the knowledge how to select properly these free parameters is a drawback of this method. In this section, these parameters will be studied numerically. Example 1. A harmonic problem whose exact solution is uðx1 ; x2 Þ ¼ sinðx1 Þsinhðx2 Þ
in O ¼ fðx1 ; x2 Þ : 1rx1 ; x2 r1g
Example 2. A harmonic problem with exact solution uðr; yÞ ¼ 3 þ r 2 cosð2yÞ
in a peanut shape domain O: 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > > x ¼ x ð y Þ ¼ cosð2yÞ þ 1:1 sin2 ð2yÞcosy; 1 > < 1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > cosð2 ¼ x ð y Þ ¼ y Þ þ 1:1 sin2 ð2yÞsiny; x > 2 2 > > :
1279
0ryr2p
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ x21 þ x22 is the distance from the origin to any point ðx1 ; x2 Þ 2 O and y ¼ tan1 ðx2 =x1 Þ.
Fig. 6. Effect of the parameter b of EXP radial function (a) Example 1; (b) Example 2.
Fig. 7. Effect of the parameter a of the integration domain (a) Example 1; (b) Example 2.
Fig. 8. Effect of the parameter b of the interpolation domain (a) Example 1; (b) Example 2.
ARTICLE IN PRESS 1280
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
Both examples are subjected to the Dirichlet boundary conditions. The distribution of boundary nodes and test nodes are shown in Fig. 3.
5.1.1. RBFs A variety of RBFs have been investigated in the past for meshless methods. To obtain a better understanding of the proposed approach, a detailed numerical study of three RBFs has been implemented in this work. These are MQ, EXP and TPS radial functions, given by Eqs. (5)–(7), respectively. However, some parameters (C and Z for MQ, b for EXP) need to be determined for MQ radial function and EXP radial function. These parameters will affect the performance of the radial functions. We first investigate the effect of C and Z for the MQ radial function. It is found that nodal spacing affects the choice of the shape parameter C [2]. As a result, C is defined as C ¼ c 0 di
There is only one parameter b in the EXP radial function. Results of relative errors for different b are obtained and depicted in Fig. 6. From this figure, it can be observed that b 2 ½0:2; 0:8 leads better results. For the TPS radial function, we obtain eðuÞ ¼ 0:02057, eðu;x1 Þ ¼ 0:02058 and eðu;x2 Þ ¼ 0:02222 for Example 1, and eðuÞ ¼ 0:02414, eðu;x1 Þ ¼ 0:02418 and eðu;x2 Þ ¼ 0:02153 for Example 2. Comparing
ð48Þ
where c0 is a dimensionless shape parameter, di is the minimum distance between the node yi and adjacent nodes. It is also found that the parameter Z has great influence on the performance of the radial basis point interpolation than that of the parameter C. Therefore, Z is firstly studied with c0 fixed at 0.5, 1.0, 1.5 and 2.0. Results for various values of Z are shown in Fig. 4. It can be found that Z 2 ½0:5; 1Þ [ ð1; 2Þ yields a satisfied result. Errors for different values of c0 are plotted in Fig. 5. It can be observed that c0 2 ½1:0; 5:0 can lead to satisfactory results.
Fig. 11. Heat conduction in a quadrant of an annular.
Fig. 9. Effect of the parameter g of the Gauss weights (a) Example 1; (b) Example 2.
Fig. 10. Convergence in relative error by using the Spline and Gauss weights ðg ¼ 1:5Þ (a) Example 1; (b) Example 2.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
the results of errors between MQ, EXP and TPS, we can see that the computational accuracy is little influenced by RBFs. For the following problems, the MQ radial function with parameters Z ¼ 0:5 and c0 ¼ 2:0 will be used. 5.1.2. Integration domain parameter a To study how the local sub-domain size affects the solution accuracy of the present method, relative errors with respect to the parameter a is shown in Fig. 7. It should be mentioned here that the union of fGis gN i¼1 does not cover the whole bounding surface when ao0:50. Fig. 7 indicates a 2 ½0:80; 0:99 leads to a satisfactory result. For convenience, a ¼ 0:90 will be used in the future computations. 5.1.3. Interpolation domain parameter b The size of interpolation domain of a quadrature point is decided by the parameter b in Eq. (44). Fig. 8 displays the relationship between the parameter b and relative errors. The results show that b 2 ½2:5; 4:5 yields a satisfactory result. For the following calculations, b is chosen as 3.0.
where g is a dimensionless shape parameter. Relative errors for different values of g are plotted in Fig. 9. It can be found that the errors change with g and the results for g 2 ½1:0; 2:0 are very good. When the parameter is too small ðgo1:0Þ or too big ðg42:0Þ, the errors increase. Both the Spline weight function and the Gaussian weight function are used to study the convergence of the HRBNM. The convergence results are shown in Fig. 10, where three kinds of nodal arrangement of 16, 32 and 64 uniformly distributed boundary nodes are employed. The results show that the convergence rates of the Gaussian weights are higher than those of the Spline weights.
5.2. Heat conduction A heat conduction problem in Fig. 11 is studied here, u represents the temperature field. Because of symmetry, only onequarter of the structure is modeled. For this problem, the boundary conditions are uða; yÞ ¼ Ta cosð4yÞ;
5.1.4. Weight functions In the Gaussian weight function, the parameter ci will influence the performance of the HRBNM. This parameter is defined as ci ¼ gdi
ð49Þ
Fig. 12. u and q on BC of the plate ðr ¼ 2Þ.
1281
@u ¼0 @y
uðb; yÞ ¼ Tb cosð4yÞ
on y ¼ 0 and y ¼
p 2
The analytical solution for this problem is [32] 4 4 a4 b4 r a4 r b4 Tb 4 4 Ta 4 4 cosð4yÞ uðr; yÞ ¼ 8 8 b a a r b r
The parameters are chosen as a ¼ 1, b ¼ 2, Ta ¼ 10 and Tb ¼ 5 in computation. Besides, we applied 32 boundary nodes (four regularly distributed nodes on AB and CD, 16 regularly distributed nodes on BC, and eight regularly distributed nodes on AD). Fig. 12 displays a comparison between the present numerical results with the exact solutions for u and q along the curved boundary BC. The results of the HBNM are also presented in this figure for comparison. It can be found that excellent agreement between the analytical and numerical solutions is achieved. The numerical results for the temperature u and its derivatives along the arc r ¼ 1:5 obtained by the HRBNM, the HBNM and the BRPIM are shown in Fig. 13 together with the analytical solutions. It is clearly shown that the numerical solutions for all approaches are seen to capture the behavior of the exact solutions very well for points far away from the boundary, while at points close to the boundary, the BRPIM has a little lower accuracy than the HRBNM and the HBNM.
Fig. 13. Results of (a) u and (b) u;x1 on the arc r ¼ 1:5.
ARTICLE IN PRESS 1282
X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
Fig. 14. Flow in a confined porous layer. Fig. 15. Potential u and its derivatives at the section x2 ¼ 2:5.
5.3. Fluid flow In this example, we solve the problem of fluid flow in a thin layer of a porous geological material (see Fig. 14). The thickness of the layer erL or b. The layer is contained between two impervious strata. The flow takes place in the plane of the porous layer and the hydraulic conductivity of the layer is isotropic and accounts for any non-uniformity in the fluid velocity distribution over its thickness. The boundary conditions for the hydraulic potential are uðx1 ; 0Þ ¼ 0;
uðx1 ; bÞ ¼ 0;
uð0; x2 Þ ¼ 0;
uðL; x2 Þ ¼ P0
where P0 is a constant pressure head imposed on x1 ¼ L. It can be noted that the potential u is composed of the datum potential ud and the pressure potential up . Because the porous layer is located on the plane x3 ¼ 0, we can set the datum at x3 ¼ 0 and u can be identified directly as the pressure potential. With these conditions, the exact solution of this problem can be given by [32] ( ) 1 2P0 X 1 þ ð1Þnþ1 sinhðnpx1 =bÞsinðnpx2 =bÞ uðx1 ; x2 Þ ¼ p n¼1 n sinhðnpL=bÞ The parameters that are used in our analysis are L ¼ 10, b ¼ 5 and P0 ¼ 10. The exact solutions and numerical solutions of the present HRBNM for u and its derivatives are depicted in Fig. 15. In this study, the boundary of the domain is discretized by 48 uniformly distributed nodes (16 nodes on AB and CD, and eight nodes on AD and BC). From this figure, we can see that the numerical solutions are in excellent agreement with the exact solutions. To show the convergence of the present method, regularly distributed 24, 48 and 96 nodes are used. The ratio of the number of nodes on AB and CD to that on AD and BC is kept constant and equal to 2. The results of convergence of @u=@x1 are shown in Fig. 16. In the figure, the convergences of the HBNM and the BRPIM are also shown for comparison. It can be found that the solution accuracy of the present approach is slightly higher than that of the HBNM and the BRPIM.
6. Conclusions A hybrid radial boundary node method is presented in this paper. It is based on a modified variational principle, which involves three independent variables, two of them defined on the boundary and the other inside the domain. The domain variables
Fig. 16. Convergence for the fluid flow.
are interpolated by the classical fundamental solutions and the boundary variables by the radial basis point interpolation. The shape functions constructed by this point interpolation possess the delta function property, thus boundary conditions can be applied directly. The HRBNM is a truly boundary-type meshless method, in which the interpolation does not use any element and the integral does not require any background mesh. Thus compared with the BNM, the BRPIM and the hybrid BRPIM, the HRBNM is truly meshless. Besides, in contrast with the LBIE method and the LRPIM, the HRBNM only requires a nodal data structure on the bounding surface of the domain to be solved. The differences among the HRBNM and the HBNM or the BNM are as follows. In the HRBNM, the basic unknown quantities are the real solutions of the nodal variables, but in the HBNM and the BNM, the basic unknown quantities are approximations of the nodal variables. In the HRBNM, boundary conditions can be directly and easily implemented, but in the HBNM and the BNM, boundary conditions are applied after they are transformed into their approximations on the boundary nodes with the MLS approximation. Hence, the HRBNM is expected to have higher computational precision.
ARTICLE IN PRESS X. Li et al. / Engineering Analysis with Boundary Elements 33 (2009) 1273–1283
Some strongly singular integration in the HRBNM is circumvented with the help of a uniform potential field. Besides, an adaptive integration scheme is employed for alleviating the ‘boundary layer effect’. Numerical examples are presented for several potential problems. The numerical results have demonstrated the accuracy, convergence and effectiveness of the HRBNM. The effects of shape parameters of RBFs, sizes of interpolation domains and integration domains, and weight functions are investigated in detail. The present method can easily be implemented to other problems that can be solved by the BEM.
Acknowledgment The authors would like to thank the anonymous reviewers whose constructive comments improve the quality of this article.
References [1] Belytschko T, Krongauz Y, Organ D, Fleming M. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng 1996;139:3–47. [2] Liu GR, Gu YT. An introduction to meshfree methods and their programming. Dordrecht: Springer; 2005. [3] Lancaster P, Salkauskas K. Surface generated by moving least squares methods. Math Comput 1981;37:141–58. [4] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse element. Comput Mech 1992;10:307–18. [5] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [6] Atluri SN, Shen SP. The meshless local Petrov–Galerkin (MLPG) method: a simple & less-costly alternative to the finite element and boundary element methods. Comput Model Eng Sci 2002;3:11–51. [7] Atluri SN, Sladek J, Sladek V, Zhu T. The local boundary integral equation (LBIE) and it’s meshless implementation for linear elasticity. Comput Mech 2000;25:180–98. [8] Mukherjee YX, Mukherjee S. The boundary node method for potential problems. Int J Numer Methods Eng 1997;40:797–15. [9] Zhang J, Yao Z, Li H. A hybrid boundary node method. Int J Numer Methods Eng 2002;53:751–63. [10] Kitipornchai S, Liew KM, Cheng Y. A boundary element-free method (BEFM) for three-dimensional elasticity problems. Comput Mech 2005;36:13–20.
1283
[11] Peng M, Cheng Y. A boundary element-free method (BEFM) for twodimensional potential problems. Eng Anal Bound Elem 2009;33:77–82. [12] Li X, Zhu J. A Galerkin boundary node method and its convergence analysis. J Comput Appl Math 2009;230:314–28. [13] Li X, Zhu J. A Galerkin boundary node method for biharmonic problems. Eng Anal Bound Elem 2009;33:858–65. [14] Zhang J, Yao Z. Meshless regular hybrid boundary node method. Comput Model Eng Sci 2001;2:301–18. [15] Zhang J, Tanaka M, Matsumoto T. Meshless analysis of potential problems in three dimensions with the hybrid boundary node method. Int J Numer Methods Eng 2004;59:1147–68. [16] Zhang J, Yao Z. The regular hybrid boundary node method for threedimensional linear elasticity. Eng Anal Bound Elem 2004;28:525–34. [17] Miao Y, Wang Y. Development of hybrid boundary node method in two dimensional elasticity. Eng Anal Bound Elem 2005;29:703–12. [18] Yan F, Wang Y, Tham LG, Cheung YK. Dual reciprocity hybrid boundary node method for 2-D elasticity with body force. Eng Anal Bound Elem 2008;32: 713–75. [19] Yan F, Wang Y, Miao Y, Cheung YK. Dual reciprocity hybrid boundary node method for free vibration analysis. J Sound Vib 2009;321:1036–57. [20] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 2002;54:1623–48. [21] Liu GR, Gu YT. A local radial point interpolation method (LRPIM) for free vibration analysis of 2-D solids. J Sound Vib 2001;246:29–46. [22] Liu GR, Yan L, Wang JG, Gu YT. Point interpolation method based on local residual formulation using radial basis functions. Struct Eng Mech 2002;14:713–32. [23] Gu YT, Liu GR. A boundary radial point interpolation method (BRPIM) for 2-D structural analyses. Struct Eng Mech 2002;15:535–50. [24] Liu GR, Gu YT. Boundary meshfree methods based on the boundary point interpolation methods. Eng Anal Bound Elem 2004;28:475–87. [25] Wu YL, Liu GR. A meshfree formulation of local radial point interpolation method (LRPIM) for incompressible flow simulation. Comput Mech 2003;30:355–65. [26] Liu GR, Gu YT. Comparisons of two meshfree local point interpolation methods for structural analyses. Comput Mech 2002;29:107–21. [27] Dumont NA. The hybrid boundary element method. In: Brebbia CA, Wendland WL, Kuhn G, editors. Boundary elements IX, vol. 1. Berlin: Springer; 1988. [28] DeFiigueiredo TGB. A new boundary element formulation in engineering. Berlin: Springer; 1991. [29] Gu YT, Liu GR. Hybrid boundary point interpolation method and their coupling with the element free Galerkin method. Eng Anal Bound Elem 2003;27:905–17. [30] Madych WR, Nelson SA. Multivariate interpolation and conditionally positive definite functions II. Math Comput Math Comput 1990;54:211–30. [31] Powell MJD. The theory of radial basis function approximation in 1990. In: Light W, editor. Advances in numerical analysis, vol. II. Oxford: Oxford Science Publications; 1997. [32] Selvadurai APS. Partial differential equations in mechanics, vol. 1. Berlin: Springer; 2000.