- Email: [email protected]
Elementary analytical integrals required in subtraction of singularity method for evaluation of weakly singular boundary integrals
ARTICLE IN PRESS
Engineering Analysis with Boundary Elements 31 (2007) 241–247 www.elsevier.com/locate/enganabound
Elementary analytical integrals required in subtraction of singularity method for evaluation of weakly singular boundary integrals Krishna M. Singha, Masataka Tanakab, a
Department of Mechanical and Industrial Engineering, Indian Institute of Technology, Roorkee 247667, India Department of Mechanical Systems Engineering, Faculty of Engineering, Shinshu University, 4-17-1 Wakasato, Nagano 380-8553, Japan
b
Received 19 December 2005; accepted 13 May 2006 Available online 30 January 2007
Abstract This paper presents closed form analytical formulae for the elementary integrals required in the numerical evaluation of weakly singular boundary integrals involving logarithmic singularity using the subtraction of singularity method. Results are limited to curved quadratic elements in view of their practical usefulness in modeling complex boundaries. Continuous as well as discontinuous elements have been considered. Sample numerical results obtained using these formulae together with a hybrid subtraction of singularity–nonlinear transformation method have been included to demonstrate the usefulness of these formulae. r 2006 Elsevier Ltd. All rights reserved. Keywords: Boundary element method; Weakly singular integrals; Subtraction of singularity method; Analytical formulae; Elementary integrals
1. Introduction Accurate evaluation of singular boundary integrals is very important in the boundary element analysis, since these contribute to the dominant (diagonal and near diagonal) terms of the global boundary element matrices. Various techniques have been suggested for evaluation of these integrals. These include analytical evaluation methods (mostly limited to straight line or planar elements) [1], specialized quadrature formulae [2], recursive subdivision [3], subtraction of singularity [4] and nonlinear coordinate transformations [5]. A brief review and related bibliography of these methods can be found in Singh and Tanaka [6]. Of the preceding techniques, the subtraction of singularity method (SSM) is very general and can be used in conjunction with nonlinear transformation methods for near-exact evaluation of weakly singular integrals [6]. This composite scheme based on the SSM and nonlinear transformation would be referred to as the singularity subtraction plus nonlinear transformation (SSNT) method herein. The SSM requires analytical evaluation of a few Corresponding author. Tel.: +81 26 269 5120; fax: +81 26 269 5124.
E-mail address: [email protected] (M. Tanaka). 0955-7997/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.05.003
simple singular integrals. The aim of the present paper is to provide a collection of analytical formulae for these elementary integrals required in the numerical evaluation of boundary integrals involving logarithmic singularity. We focus on curved quadratic elements in view of their widespread use in modeling complex boundaries. Results are presented for continuous isoparametric elements as well as discontinuous elements—the latter being especially useful in modeling geometric (corners) and physical (change in boundary conditions) singularities. We have not considered linear or straight-line elements since analytical formulae are already available in the literature, which would be preferred over the SSM. All the formulae have been obtained using Maple V [7] and double-checked using hand derivations. We have further checked the results obtained with these formulae by comparing them with high-precision numerical integration in Maple V. The analytical formulae presented here are rather simple, and can be easily derived. The sole objective of putting these together is to make them easily available to the users/ researchers, so that they would not have to spend time rederiving these formulae. We also hope that availability of these formulae for various types of quadratic elements would encourage the use of near-exact integration of the
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
242
weakly singular integrals using the composite SSNT scheme. We start with a brief outline of SSM (and the related hybrid scheme) for singular boundary integrals and description of the elements of different type and their shape functions. Elementary integrals required in the SSM are presented for each type of element for various possible combinations of location of the singular point and nodal shape functions. We conclude the paper with sample numerical results obtained using these formulae in conjunction with the hybrid method for weakly singular integrals for potential problem. 2. SSM for weakly singular integrals Weakly singular integrals appearing in the boundary element analysis can be expressed as Z Iðx0 ; kÞ ¼ fk ðxÞ u ðx0 ; xÞ dGðxÞ, (1) Ge
where Ge represents the integration element, u is the fundamental solution of the problem, x0 is the collocation point and fk denotes the shape function associated with the kth local node of the element. Owing to the nature of the fundamental solution u , the preceding integral becomes singular when x0 2 Ge . For 2-D problems, u normally involves a logarithmic singularity. Hence, the integral I is weakly singular. This weakly singular integral can be re-written in terms of the normalized local coordinates as Z 1 0 Iðx ; kÞ ¼ f ðZ; Zs Þ dZ; f ðZ; Zs Þ ¼ fk ðZÞ u ðx0 ; xÞJðZÞ, 1
(2)
To obtain a smoother integrand for more efficient evaluation of I r , we can use a nonlinear transformation ZðgÞ: g 2 ½1; 1 ! Z 2 ½1; 1, and thus evaluate it as Z 1 I r ðx0 ; kÞ ¼ f r ðZðgÞ; Zs ÞJ g ðgÞ dg 1
¼
NG X
f r ðZðgi Þ; Zs ÞJ g ðgi Þwi ,
ð5Þ
i¼1
where NG denotes the order of Gaussian quadrature, and gi and wi denote the ith quadrature point and associated weight, respectively. This combined scheme is termed as the SSNT method. For 2-D problems involving logarithmic singularity, us ðx0 ; xÞ ¼ ln r. To obtain a simple form for f s ðZ; Zs Þ, we can use the following Taylor series expansion of r around Zs : r ¼ r~ þ OððZ Zs Þ2 Þ;
r~ ¼ JðZs ÞjZ Zs j.
(6)
Now, the following two simple choices are possible for f s ðZ; Zs Þ: f s ðZ; Zs Þ ¼ fk ðZs ÞJðZs Þ ln r~ and f s ðZ; Zs Þ ¼ fk ðZÞJðZs Þ ln r~.
(i) LOGA: (ii) LOGB:
Since fk ðZs Þ and JðZs Þ would be constant for a particular singular point, therefore, the only elementary integral involved with the choice LOGA would be Z 1 ln jZ Zs j dZ 1 8 2 lnð2Þ 2 if jZs j ¼ 1; > > < ð7Þ ¼ ð1 Zs Þ lnð1 Zs Þ > > : þðZ þ 1Þ lnðZ þ 1Þ 2 if Z 2 ð1; 1Þ: s
s
s
0
where Zs ¼ Zðx Þ and JðZÞ is the Jacobian of the transformation xðZÞ. The integrand f ðZ; Zs Þ is singular at the point Zs 2 ½1; 1. Let us suppose that the kernel u ðx0 ; xÞ exhibits the singularity typical of the singular function us ðx0 ; xÞ as x ! x0 . For Zs 2 ½1; 1, let us define f s ðZ; Zs Þ ¼ pðZ; Zs Þ us ðx0 ; xÞ, in which pðZ; Zs Þ is a regular function, and f r ðZ; Zs Þ ¼ f ðZ; Zs Þ f s ðZ; Zs Þ. Further, let us assume that f s ðZ; Zs Þ is such that limZ!Zs f r ðZ; Zs Þ ¼ 0 (or C), where C is a bounded constant. Then, the function f r ðZ; Zs Þ is regular. We can now write the integral I as Iðx0 ; kÞ ¼ I r ðx0 ; kÞ þ I s ðx0 ; kÞ,
(3)
I r ðx0 ; kÞ ¼
1
f r ðZ; Zs Þ dZ 1
and
I s ðx0 ; kÞ ¼
Z
f s ðZ; Zs Þ ¼ JðZs ÞflnðJðZs ÞÞfk ðZÞ þ fk ðZÞ ln jZ Zs jg.
(8)
0
Hence, the singular integral I s ðx ; kÞ is given by I s ðx0 ; kÞ ¼ JðZs Þ½lnðJðZs ÞÞBk þ Dk ðZs Þ, where Bk and Dk are elementary integrals given by Z 1 Bk ¼ fk ðZÞ dZ,
(9)
ð10Þ
1 Z 1
Dk ðZs Þ ¼
fk ðZÞ ln jZ Zs j dZ.
ð11Þ
1
The expressions for these elementary integrals for quadratic elements are presented in the sequel.
where Z
With the choice LOGB,
1
f s ðZ; Zs Þ dZ. 1
(4) We assume that I s can be evaluated analytically. Function f r ðZ; Zs Þ being regular, I r can be evaluated using Gaussian quadrature. The preceding approach summarized by Eq. (3) is the SSM for evaluation of the singular integral Iðx0 ; kÞ.
3. Elementary integrals for quadratic elements Quadratic elements are very often used to model problems with curvilinear geometry. For these elements, both the geometry and the boundary quantities are approximated by second degree polynomials in the local or intrinsic coordinate Z. Thus, variation of a function f
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
along a quadratic element can be expressed as f ðZÞ ¼
3 X
f i fi ðZÞ,
(12)
i¼1
where f i is the value of the function at the local node i and fi ’s are the so-called shape functions. Geometry of an element is also defined using the nodal coordinates xi and the shape functions as xðZÞ ¼
3 X
We would use superscript C to denote the elementary integrals for a continuous element. Integrals BC k are given by 1 BC 1 ¼ 3;
4 BC 2 ¼ 3;
xi fi ðZÞ.
(13)
Generally the three geometric nodes of an element are used as collocation points. Such an element is referred to as a continuous quadratic element. To model the corner points or the points where there is a change in the boundary conditions, one sided discontinuous elements (commonly referred to as the partially discontinuous elements) are used. These elements can be of two types depending on which extreme collocation point has been moved inside the element to model the geometric or physical singularity. If it is the third collocation node which has been moved inside, then the resulting element is referred to as the partially discontinuous element of Type 1. Otherwise, if the first node is moved inside, then the element is called the partially discontinuous element of Type 2. One may opt for moving both extreme nodes inside which results in a discontinuous element. In the graphical representation of these elements, we have used a for a geometric node and a for a collocation node.
(15)
16 DC 2 ð0Þ ¼ 9 ;
1 DC 3 ð0Þ ¼ 9.
(16)
Elementary integrals DC k for Zs ¼ 1 and Zs ¼ þ1 are C 1 17 DC 1 ð1Þ ¼ 3 lnð2Þ 18 ¼ D3 ðþ1Þ, C 4 10 DC 2 ð1Þ ¼ 3 lnð2Þ 9 ¼ D2 ðþ1Þ, C 1 1 DC 3 ð1Þ ¼ 3 lnð2Þ þ 18 ¼ D1 ðþ1Þ.
ð17Þ
3.2. Partially discontinuous element of Type 1 Fig. 2 provides the graphical representation of a partially discontinuous quadratic element of Type 1, in which Z1 represents the shifted local coordinate of third collocation point. Commonly used values of Z1 are 23 or 12. Geometry of the element can still be modelled using the coordinates of the geometric nodes and shape functions for the continuous elements given by Eq. (14). However, the variation of a problem variable is modeled using the following shape functions: ZðZ Z1 Þ , 1 þ Z1 ðZ þ 1ÞðZ Z1 Þ P f2 1 ðZÞ ¼ , Z1 ZðZ þ 1Þ P , f3 1 ðZÞ ¼ Z1 ð1 þ Z1 Þ P
f1 1 ðZÞ ¼
3.1. Continuous quadratic elements Fig. 1 gives the description of a continuous quadratic element. The shape functions for this element are f1 ðZÞ ¼ 12ZðZ 1Þ, 2
f2 ðZÞ ¼ 1 Z , f3 ðZÞ ¼ 12ZðZ þ 1Þ.
ð14Þ
ð18Þ
where superscript P1 has been used to denote the partially P discontinuous element of Type 1. Elementary integrals Bk 1
3 X
y
3 y
1 BC 3 ¼ 3.
Elementary integrals DC k have very simple form for any of the three collocation points. For Zs ¼ 0, these are 1 DC 1 ð0Þ ¼ 9;
i¼1
243
X 2
2
1
X x
x 1
2
3
-1
0
1
1 X η
Fig. 1. Continuous (isoparametric) quadratic element in global (x–y) and local (Z) coordinate systems.
-1
2 X 0
3 X η1
1
1
η
Fig. 2. Partially discontinuous quadratic element of Type 1 in global (x–y) and local (Z) coordinate systems: local node 3 is located at Z ¼ Z1 .
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
244
are given by P
Table 2 Elementary integrals for partially continuous element of Type 1 for Z1 ¼ 23
2 , 3 ð1 þ Z1 Þ 2 ¼2 , 3 Z1 2 . ¼ 3 Z1 ð1 þ Z1 Þ
B1 1 ¼ P
B2 1 P B3 1
P
P
D1 1
ð19Þ P Dk 1 ,
P D2 1 P D3 1
y
For elementary integrals the three possible values of Zs are 1, 0 and Z1 corresponding to the three collocation P points. Elementary integrals Dk 1 for Zs ¼ 1 are 6 lnð2Þ 9Z1 8 , ¼ 9ð1 þ Z1 Þ 2 1 P D2 1 ð1Þ ¼ 2 1, lnð2Þ 3Z1 9Z1 6 lnð2Þ þ 1 P . D3 1 ð1Þ ¼ 9Z1 ð1 þ Z1 Þ
P D1 1 ð1Þ
P
P
B1 1 ¼ 25 Zs ¼ 1
B2 1 ¼ 1 Zs ¼ 0
B3 1 ¼ 35 Zs ¼ 23
2 5 lnð2Þ
14 15
2 15
7 6
53 15
10 2 2 27 lnð5Þ 5 lnð3Þ 45 125 8 108 lnð5Þ lnð3Þ 9 5 3 16 36 lnð5Þ 5 lnð3Þ 15
lnð2Þ 3 5 lnð2Þ
1 þ 10
3 X
2 X
X
1
ð20Þ
P
Elementary integrals Dk 1 for Zs ¼ 0 are
x
2 , 9 ð1 þ Z1 Þ 2 P D2 1 ð0Þ ¼ 2 þ , 9 Z1 2 P . D3 1 ð0Þ ¼ 9 Z1 ð1 þ Z1 Þ P
1 X -η1
D1 1 ð0Þ ¼
-1
ð21Þ P
General relation for Dk 1 for Zs ¼ Z1 ; Z1 2 ð0; 1Þ can be written as P Dk 1 ðZ1 Þ
¼
P ak 1
þ
P bk 1
lnð1 Z1 Þ þ P
P ck 1
lnð1 þ Z1 Þ.
P
(22)
P
3.3. Partially discontinuous element of Type 2 Fig. 3 depicts a quadratic element of Type 2, in which Z1 represents the shifted local coordinate of the first
collocation point. Commonly used values of Z1 are The shape functions for this element are ZðZ 1Þ , Z1 ð1 þ Z1 Þ ð1 ZÞðZ þ Z1 Þ P f2 2 ðZÞ ¼ , Z1 ZðZ þ Z1 Þ P f3 2 ðZÞ ¼ , 1 þ Z1
2 3
or 12.
f1 2 ðZÞ ¼
Table 1 P P P Coefficients ak 1 , bk 1 and ck 1 in Eq. (22) for partially continuous element of Type 1 P
P
k
ak 1
bk 1
ck 1
1
3Z21 2 9ð1 þ Z1 Þ ð3Z2 þ 9Z1 2Þ 1 9Z1 ð6Z21 þ 9Z1 þ 2Þ 9Z1 ð1 þ Z1 Þ
ð2 þ Z1 Þð1 Z1 Þ2 6ð1 þ Z1 Þ ðZ1 þ 5Þð1 Z1 Þ2 6Z1 ð2Z21 þ 5Z1 þ 5Þð1 Z1 Þ 6Z1 ð1 þ Z1 Þ
ð1 þ Z1 Þð2 Z1 Þ 6 ð1 þ Z1 Þ3 6Z1 ð1 þ Z1 Þð2Z1 1Þ 6Z1
ð23Þ
where superscript P2 has been used to denote the partially discontinuous element of Type 2. P Elementary integrals Bk 2 are given by P
2 , 3 Z1 ð1 þ Z1 Þ 2 ¼2 , 3 Z1 2 . ¼ 3 ð1 þ Z1 Þ
B1 2 ¼ P
3
η
Fig. 3. Partially discontinuous quadratic element of Type 2 in global (x–y) and local (Z) coordinate systems: local node 1 is located at Z ¼ Z1 .
B2 2
2
3 X 1
P
Values of the coefficients ak 1 , bk 1 and ck 1 are listed in Table 1. Based on Table 1 and preceding equations, simplified values of the elementary integrals for the most commonly used value of Z1 ¼ 23 are listed in Table 2.
P
2 X 0
P
B3 2
ð24Þ P
For elementary integrals Dk 2 , the three possible values of Zs are Z1 , 0 and þ1, corresponding to the three P collocation points. The elementary integrals Dk 2 for Zs ¼ Z1 can be written as P
P
P
P
Dk 2 ðZ1 Þ ¼ ak 2 þ bk 2 lnð1 Z1 Þ þ ck 2 lnð1 þ Z1 Þ. P
P
P
(25)
Values of the coefficients ak 2 , bk 2 and ck 2 are listed in Table 3.
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
P
P
ak 1
bk 1
ck 1
1
ð6Z21 þ 9Z1 þ 2Þ 9Z1 ð1 þ Z1 Þ ð3Z21 þ 9Z1 2Þ 9Z1 3Z21 2 9ð1 þ Z1 Þ
ð2Z21 þ 5Z1 þ 5Þð1 Z1 Þ 6Z1 ð1 þ Z1 Þ ðZ1 þ 5Þð1 Z1 Þ2 6Z1 ð2 þ Z1 Þð1 Z1 Þ2 6ð1 þ Z1 Þ
ð1 þ Z1 Þð2Z1 1Þ 6Z1
2 3
P
P
P D2 2 P D3 2
2
X
1
ð1 þ Z1 Þ3 6Z1 ð1 þ Z1 Þð2 Z1 Þ 6
Table 4 Elementary integrals for partially continuous element of Type 2 for Z1 ¼ 2=3
D1 2
X
P
k
3 X
y
Table 3 P P P Coefficients ak 2 , bk 2 and ck 2 in Eq. (25) for partially continuous element of Type 2
245
P
P
B1 2 ¼ 35 Zs ¼ 23
B2 2 ¼ 1 Zs ¼ 0
B3 2 ¼ 25 Zs ¼ 1
5 3 16 36 lnð5Þ 5 lnð3Þ 15 125 8 lnð5Þ lnð3Þ 108 9 10 2 2 27 lnð5Þ 5 lnð3Þ 45
15
3 5 lnð2Þ
53
lnð2Þ 76
2 15
2 5 lnð2Þ
1 þ 10
14 15
For the middle collocation point (Zs ¼ 0), elementary P integrals Dk 2 are 2 , 9 Z1 ð1 þ Z1 Þ 2 P D2 2 ð0Þ ¼ 2 þ , 9 Z1 2 P . D3 2 ð0Þ ¼ 9 ð1 þ Z1 Þ P
x
-1
1 X -η1
2 X 0
3 X η2
1
η
Fig. 4. Discontinuous quadratic element in global (x–y) and local (Z) coordinate systems: local node 1 is located at Z ¼ Z1 and node 3 at Z ¼ Z2 .
discontinuous elements. The shape functions are ZðZ Z2 Þ , Z1 ðZ1 þ Z2 Þ ðZ þ Z1 ÞðZ Z2 Þ fD , 2 ðZÞ ¼ Z1 Z2 ZðZ þ Z1 Þ . fD 3 ðZÞ ¼ Z2 ðZ1 þ Z2 Þ fD 1 ðZÞ ¼
ð28Þ
D1 2 ð0Þ ¼
Elementary integrals BD k are given by
ð26Þ P
The elementary integrals Dk 2 for Zs ¼ þ1 are 6 lnð2Þ þ 1 , 9Z1 ð1 þ Z1 Þ 2 1 P D2 2 ðþ1Þ ¼ 2 1, lnð2Þ 3Z1 9Z1 6 lnð2Þ 9Z1 8 P . D3 2 ðþ1Þ ¼ 9ð1 þ Z1 Þ P
D1 2 ðþ1Þ ¼
ð27Þ
Based on the preceding, elementary integrals for the partially discontinuous element of Type 2 with the discontinuous first node at Z ¼ 23 (i.e. Z1 ¼ 23) are listed in Table 4.
2 , 3 Z1 ðZ1 þ Z2 Þ 2 BD , 2 ¼2 3 Z1 Z2 2 BD . 3 ¼ 3 Z2 ðZ1 þ Z2 Þ BD 1 ¼
ð29Þ
For elementary integrals DD k , the three possible values of Zs are Z1 , 0 and þZ2 , corresponding to the three collocation points. The elementary integrals DD k for Zs ¼ Z1 can be written as D
D
D
3.4. Discontinuous quadratic elements Though the continuous and partially discontinuous quadratic elements are the most widely used, some one might opt to use fully discontinuous elements. Formulae for these elements are included here for sake of completeness. In these elements, both the extreme collocation points are shifted inside as depicted in Fig. 4. We would use superscript D to denote the quantities related to
D
1 1 1 DD k ðZ1 Þ ¼ ak þ bk lnð1 Z1 Þ þ ck lnð1 þ Z1 Þ,
D
(30)
D
where coefficients ak 1 , bk 1 and ck 1 are D
ð6Z21 þ 9Z1 Z2 þ 2Þ , 9Z1 ðZ1 þ Z2 Þ ð2Z21 þ 3Z1 Z2 þ 2Z1 þ 3Z2 þ 2Þð1 Z1 Þ , ¼ 6Z1 ðZ1 þ Z2 Þ ð2Z21 þ 3Z1 Z2 2Z1 3Z2 þ 2Þð1 þ Z1 Þ , ¼ 6Z1 ðZ1 þ Z2 Þ
a1 1 ¼ D
b1 1 D
c1 1
ð31Þ
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
246 D
a2 1 ¼ D
b2 1 ¼ D
c2 1 ¼
D
a3 1 ¼ D
b3 1 ¼ D
c3 1 ¼
ð3Z21 þ 9Z1 Z2 2Þ , 9Z1 Z2
D
c3 2 ¼
ðZ1 þ 3Z2 þ 2Þð1 Z1 Þ2 , 6Z1 Z2
ð1 þ Z1 Þ2 ðZ1 þ 3Z2 2Þ , 6Z1 Z2
ð32Þ
3Z21 2 , 9Z2 ðZ1 þ Z2 Þ
ð1 þ Z1 Þ2 ð2 Z1 Þ . 6Z2 ðZ1 þ Z2 Þ
We may note that the relations for the discontinuous element represent the general formulae from which the relations for the continuous and partially discontinuous elements can be obtained by taking care of the proper limits.
ð33Þ
Let us consider weakly singular integral arising in the boundary element analysis of Laplace equation in 2-D is given by Z I L ðx0 ; kÞ ¼ fk ðxÞ log jx x0 j dGðxÞ Ge 1
Z
For the middle collocation point (Zs ¼ 0), elementary integrals DD k are 2 , 9 Z1 ðZ1 þ Z2 Þ 2 DD , 2 ð0Þ ¼ 2 þ 9 Z 1 Z2 2 DD . 3 ð0Þ ¼ 9 Z2 ðZ1 þ Z2 Þ
DD 1 ð0Þ ¼
¼
þ
where coefficients D
a1 2 ¼ D
b1 2 ¼ D
c1 2 ¼
D
D
(35)
are
ð36Þ
ð3Z22 þ 9Z1 Z2 2Þ , 9Z1 Z2 ðZ2 þ 3Z1 þ 2Þð1 Z2 Þ2 , 6Z1 Z2
ð1 þ Z2 Þ2 ðZ2 þ 3Z1 2Þ , 6Z1 Z2
ð6Z22 þ 9Z1 Z2 þ 2Þ , 9Z2 ðZ1 þ Z2 Þ ð2Z22 þ 3Z1 Z2 þ 2Z2 þ 3Z1 þ 2Þð1 Z2 Þ , ¼ 6Z2 ðZ1 þ Z2 Þ
a3 2 ¼ b3 2
and
D ck 2
lnð1 þ Z2 Þ,
ð1 þ Z2 Þ2 ð2 Z2 Þ , 6Z1 ðZ1 þ Z2 Þ
b2 2 ¼
D
D bk 2
D ck 2
ð2 þ Z2 Þð1 Z2 Þ2 , 6Z1 ðZ1 þ Z2 Þ
D
D
lnð1 Z2 Þ þ
D ak 2 ,
fk ðZÞ log jx x0 jJðZÞ dZ.
ð39Þ
1
ð34Þ
3Z22 2 , 9Z1 ðZ1 þ Z2 Þ
a2 2 ¼
c2 2 ¼
D bk 2
¼
We take a curved quadratic element with nodes x1 ¼ ð1; 1Þ; x2 ¼ ð2:5; 3Þ and x3 ¼ ð2; 5Þ, in which x2 is the middle node. The higher order transformation of Sato et al. [5] with k ¼ 3 has been employed in the SSNT method.
Finally, the elementary integrals DD k for Zs ¼ þZ2 can be written as D ak 2
ð38Þ
4. Numerical examples
ð2 þ Z1 Þð1 Z1 Þ2 , 6Z2 ðZ1 þ Z2 Þ
DD k ðþZ2 Þ
ð2Z22 þ 3Z1 Z2 2Z2 3Z1 þ 2Þð1 þ Z2 Þ . 6Z2 ðZ1 þ Z2 Þ
ð37Þ
Table 5 Relative error (in percent) in numerical integration of weakly singular Iðx0 ; kÞ with choice LOGA and SSNT method k¼1
k¼2
k¼3
1:1E 07 4:2E 12 5:2E 13 2:6E 12 2:0E 14 2:0E 14 4:6E 08 4:5E 12 3:1E 14
3:3E 07 2:0E 11 2:5E 13 2:4E 12 1:4E 14 2:8E 14 3:9E 08 1:0E 11 2:2E 13
3:0E 07 1:7E 11 7:6E 14 2:1E 12 0:0E þ 00 0:0E þ 00 5:4E 08 1:0E 11 4:1E 13
Partially discontinuous element of Type 1 1:00 16 2:5E 07 24 2:3E 13 32 7:6E 13 0.00 16 2:5E 12 24 0:0E þ 00 32 1:8E 14 0.67 16 1:3E 09 24 1:1E 12 32 1:9E 14
1:1E 06 6:3E 11 5:9E 13 2:3E 12 1:2E 14 0:0E þ 00 4:6E 09 5:7E 12 2:2E 13
3:0E 07 1:7E 11 1:0E 13 2:1E 12 0:0E þ 00 1:1E 14 1:3E 09 3:5E 12 5:8E 14
Partially discontinuous element of Type 2 0:67 16 1:1E 08 24 4:6E 12 32 1:6E 13 0.00 16 2:5E 12 24 0:0E þ 00 32 0:0E þ 00 1.00 16 4:6E 08 24 4:5E 12 32 1:7E 14
2:9E 08 7:5E 12 2:7E 13 2:4E 12 2:1E 14 0:0E þ 00 3:2E 07 5:8E 11 1:1E 12
2:4E 10 2:2E 12 2:4E 14 2:2E 12 0:0E þ 00 2:9E 14 5:7E 08 1:6E 11 4:4E 13
Zs
NG
Continuous elements 1:00 16 24 32 0.00 16 24 32 1.00 16 24 32
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247 Table 6 Relative error (in percent) in numerical integration of weakly singular Iðx0 ; kÞ with choice LOGB and SSNT method
247
These results clearly demonstrate the usefulness of the formulae presented herein and the SSNT method.
k¼1
k¼2
k¼3
5. Conclusions
1:2E 07 7:4E 12 1:0E 13 2:5E 12 2:0E 14 0:0E þ 00 4:6E 08 2:5E 12 1:5E 14
3:3E 07 2:6E 11 2:0E 14 2:4E 12 1:4E 14 1:4E 14 4:0E 08 4:0E 12 2:4E 14
3:0E 07 2:0E 11 0:0E þ 00 2:0E 12 2:1E 14 2:1E 14 5:3E 08 2:1E 12 6:5E 14
Partially discontinuous element of Type 1 1:00 16 2:5E 07 24 1:6E 11 32 1:4E 13 0.00 16 2:4E 12 24 1:8E 14 32 1:8E 14 0.67 16 1:4E 09 24 0:0E þ 00 32 0:0E þ 00
We have presented closed form analytical formulae for the elementary integrals required in the numerical evaluation of weakly singular boundary integrals involving logarithmic singularity using the subtraction of singularity method. Results have been presented for continuous, partially discontinuous and discontinuous quadratic elements. Sample numerical results obtained with a hybrid subtraction of singularity–nonlinear transformation method clearly demonstrate the usefulness of these formulae.
1:1E 06 8:2E 11 1:0E 13 2:3E 12 1:2E 14 2:4E 14 5:3E 09 0:0E þ 00 0:0E þ 00
3:0E 07 2:0E 11 2:1E 14 2:0E 12 3:4E 14 2:3E 14 1:6E 09 5:0E 13 0:0E þ 00
Partially discontinuous element of Type 2 0:67 16 1:1E 08 24 8:6E 13 32 3:9E 14 0.00 16 2:5E 12 24 0:0E þ 00 32 2:2E 14 1.00 16 4:6E 08 24 2:5E 12 32 1:7E 14
3:0E 08 2:1E 13 1:6E 14 2:4E 12 2:1E 14 2:1E 14 3:2E 07 2:5E 11 5:1E 14
4:6E 10 3:7E 13 2:4E 14 2:1E 12 1:5E 14 0:0E þ 00 5:6E 08 3:9E 12 1:1E 13
Zs
NG
Continuous elements 1:00 16 24 32 0.00 16 24 32 1.00 16 24 32
The ‘‘exact values’’ of the integrals correspond to the numerical values obtained with Maple V [7] using precision of 25 digits, and truncated to 16 digits in double precision computations. Errors in the numerical evaluation of the weakly singular integral I L ðxi ; kÞ obtained with the choices LOGA and LOGB are presented in Tables 5 and 6, respectively, for continuous as well as partial continuous elements. Accuracy of results is evident. Near-exact results are obtained with 32 Gauss points for any element type.
Acknowledgments Part of this work was financially supported by the Ministry of Education, Science, Sports and Culture, Japan, Grant-in-Aid for Scientific Research (C), No. 16560066, 2005, granted to the author (M.T.). This financial support is gratefully acknowledged. References [1] Kuwahara T, Takeda T. On a formula of boundary integral for the higher-order element in boundary element method and a consideration of error of solution. Trans Inst Electr Eng Jpn 1985;105-A:1–8 [in Japanese]. [2] Aimi A, Diligenti M, Monegato G. New numerical integration schemes for applications of Galerkin BEM to 2-D problems. Int J Numer Meth Eng 1997;40:1977–99. [3] Doblare´ M. Computational aspects of the boundary element method. Topics in boundary element research, vol. 3. Berlin: Springer; 1987 (chapter 4). [4] Katsikadelis JT. Boundary Elements: Theory and Applications. Oxford: Elsevier Science Ltd; 2002. [5] Sato M, Yoshiyoka S, Tsukui K, Yuuki R. Accurate numerical integration of singular kernels in the two-dimensional boundary element method. In: Brebbia CA, editor. Boundary elements X, vol. 1. Berlin: Springer; 1988. p. 279–96. [6] Singh KM, Tanaka M. On non-linear transformations for accurate numerical evaluation of weakly singular boundary integrals. Int J Numer Meth Eng 2001;50:2007–30. [7] Heal KM, Hansen ML, Rickard KM. Maple V: Learning Guide. New York: Springer; 1996.
Engineering Analysis with Boundary Elements 31 (2007) 241–247 www.elsevier.com/locate/enganabound
Elementary analytical integrals required in subtraction of singularity method for evaluation of weakly singular boundary integrals Krishna M. Singha, Masataka Tanakab, a
Department of Mechanical and Industrial Engineering, Indian Institute of Technology, Roorkee 247667, India Department of Mechanical Systems Engineering, Faculty of Engineering, Shinshu University, 4-17-1 Wakasato, Nagano 380-8553, Japan
b
Received 19 December 2005; accepted 13 May 2006 Available online 30 January 2007
Abstract This paper presents closed form analytical formulae for the elementary integrals required in the numerical evaluation of weakly singular boundary integrals involving logarithmic singularity using the subtraction of singularity method. Results are limited to curved quadratic elements in view of their practical usefulness in modeling complex boundaries. Continuous as well as discontinuous elements have been considered. Sample numerical results obtained using these formulae together with a hybrid subtraction of singularity–nonlinear transformation method have been included to demonstrate the usefulness of these formulae. r 2006 Elsevier Ltd. All rights reserved. Keywords: Boundary element method; Weakly singular integrals; Subtraction of singularity method; Analytical formulae; Elementary integrals
1. Introduction Accurate evaluation of singular boundary integrals is very important in the boundary element analysis, since these contribute to the dominant (diagonal and near diagonal) terms of the global boundary element matrices. Various techniques have been suggested for evaluation of these integrals. These include analytical evaluation methods (mostly limited to straight line or planar elements) [1], specialized quadrature formulae [2], recursive subdivision [3], subtraction of singularity [4] and nonlinear coordinate transformations [5]. A brief review and related bibliography of these methods can be found in Singh and Tanaka [6]. Of the preceding techniques, the subtraction of singularity method (SSM) is very general and can be used in conjunction with nonlinear transformation methods for near-exact evaluation of weakly singular integrals [6]. This composite scheme based on the SSM and nonlinear transformation would be referred to as the singularity subtraction plus nonlinear transformation (SSNT) method herein. The SSM requires analytical evaluation of a few Corresponding author. Tel.: +81 26 269 5120; fax: +81 26 269 5124.
E-mail address: [email protected] (M. Tanaka). 0955-7997/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.05.003
simple singular integrals. The aim of the present paper is to provide a collection of analytical formulae for these elementary integrals required in the numerical evaluation of boundary integrals involving logarithmic singularity. We focus on curved quadratic elements in view of their widespread use in modeling complex boundaries. Results are presented for continuous isoparametric elements as well as discontinuous elements—the latter being especially useful in modeling geometric (corners) and physical (change in boundary conditions) singularities. We have not considered linear or straight-line elements since analytical formulae are already available in the literature, which would be preferred over the SSM. All the formulae have been obtained using Maple V [7] and double-checked using hand derivations. We have further checked the results obtained with these formulae by comparing them with high-precision numerical integration in Maple V. The analytical formulae presented here are rather simple, and can be easily derived. The sole objective of putting these together is to make them easily available to the users/ researchers, so that they would not have to spend time rederiving these formulae. We also hope that availability of these formulae for various types of quadratic elements would encourage the use of near-exact integration of the
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
242
weakly singular integrals using the composite SSNT scheme. We start with a brief outline of SSM (and the related hybrid scheme) for singular boundary integrals and description of the elements of different type and their shape functions. Elementary integrals required in the SSM are presented for each type of element for various possible combinations of location of the singular point and nodal shape functions. We conclude the paper with sample numerical results obtained using these formulae in conjunction with the hybrid method for weakly singular integrals for potential problem. 2. SSM for weakly singular integrals Weakly singular integrals appearing in the boundary element analysis can be expressed as Z Iðx0 ; kÞ ¼ fk ðxÞ u ðx0 ; xÞ dGðxÞ, (1) Ge
where Ge represents the integration element, u is the fundamental solution of the problem, x0 is the collocation point and fk denotes the shape function associated with the kth local node of the element. Owing to the nature of the fundamental solution u , the preceding integral becomes singular when x0 2 Ge . For 2-D problems, u normally involves a logarithmic singularity. Hence, the integral I is weakly singular. This weakly singular integral can be re-written in terms of the normalized local coordinates as Z 1 0 Iðx ; kÞ ¼ f ðZ; Zs Þ dZ; f ðZ; Zs Þ ¼ fk ðZÞ u ðx0 ; xÞJðZÞ, 1
(2)
To obtain a smoother integrand for more efficient evaluation of I r , we can use a nonlinear transformation ZðgÞ: g 2 ½1; 1 ! Z 2 ½1; 1, and thus evaluate it as Z 1 I r ðx0 ; kÞ ¼ f r ðZðgÞ; Zs ÞJ g ðgÞ dg 1
¼
NG X
f r ðZðgi Þ; Zs ÞJ g ðgi Þwi ,
ð5Þ
i¼1
where NG denotes the order of Gaussian quadrature, and gi and wi denote the ith quadrature point and associated weight, respectively. This combined scheme is termed as the SSNT method. For 2-D problems involving logarithmic singularity, us ðx0 ; xÞ ¼ ln r. To obtain a simple form for f s ðZ; Zs Þ, we can use the following Taylor series expansion of r around Zs : r ¼ r~ þ OððZ Zs Þ2 Þ;
r~ ¼ JðZs ÞjZ Zs j.
(6)
Now, the following two simple choices are possible for f s ðZ; Zs Þ: f s ðZ; Zs Þ ¼ fk ðZs ÞJðZs Þ ln r~ and f s ðZ; Zs Þ ¼ fk ðZÞJðZs Þ ln r~.
(i) LOGA: (ii) LOGB:
Since fk ðZs Þ and JðZs Þ would be constant for a particular singular point, therefore, the only elementary integral involved with the choice LOGA would be Z 1 ln jZ Zs j dZ 1 8 2 lnð2Þ 2 if jZs j ¼ 1; > > < ð7Þ ¼ ð1 Zs Þ lnð1 Zs Þ > > : þðZ þ 1Þ lnðZ þ 1Þ 2 if Z 2 ð1; 1Þ: s
s
s
0
where Zs ¼ Zðx Þ and JðZÞ is the Jacobian of the transformation xðZÞ. The integrand f ðZ; Zs Þ is singular at the point Zs 2 ½1; 1. Let us suppose that the kernel u ðx0 ; xÞ exhibits the singularity typical of the singular function us ðx0 ; xÞ as x ! x0 . For Zs 2 ½1; 1, let us define f s ðZ; Zs Þ ¼ pðZ; Zs Þ us ðx0 ; xÞ, in which pðZ; Zs Þ is a regular function, and f r ðZ; Zs Þ ¼ f ðZ; Zs Þ f s ðZ; Zs Þ. Further, let us assume that f s ðZ; Zs Þ is such that limZ!Zs f r ðZ; Zs Þ ¼ 0 (or C), where C is a bounded constant. Then, the function f r ðZ; Zs Þ is regular. We can now write the integral I as Iðx0 ; kÞ ¼ I r ðx0 ; kÞ þ I s ðx0 ; kÞ,
(3)
I r ðx0 ; kÞ ¼
1
f r ðZ; Zs Þ dZ 1
and
I s ðx0 ; kÞ ¼
Z
f s ðZ; Zs Þ ¼ JðZs ÞflnðJðZs ÞÞfk ðZÞ þ fk ðZÞ ln jZ Zs jg.
(8)
0
Hence, the singular integral I s ðx ; kÞ is given by I s ðx0 ; kÞ ¼ JðZs Þ½lnðJðZs ÞÞBk þ Dk ðZs Þ, where Bk and Dk are elementary integrals given by Z 1 Bk ¼ fk ðZÞ dZ,
(9)
ð10Þ
1 Z 1
Dk ðZs Þ ¼
fk ðZÞ ln jZ Zs j dZ.
ð11Þ
1
The expressions for these elementary integrals for quadratic elements are presented in the sequel.
where Z
With the choice LOGB,
1
f s ðZ; Zs Þ dZ. 1
(4) We assume that I s can be evaluated analytically. Function f r ðZ; Zs Þ being regular, I r can be evaluated using Gaussian quadrature. The preceding approach summarized by Eq. (3) is the SSM for evaluation of the singular integral Iðx0 ; kÞ.
3. Elementary integrals for quadratic elements Quadratic elements are very often used to model problems with curvilinear geometry. For these elements, both the geometry and the boundary quantities are approximated by second degree polynomials in the local or intrinsic coordinate Z. Thus, variation of a function f
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
along a quadratic element can be expressed as f ðZÞ ¼
3 X
f i fi ðZÞ,
(12)
i¼1
where f i is the value of the function at the local node i and fi ’s are the so-called shape functions. Geometry of an element is also defined using the nodal coordinates xi and the shape functions as xðZÞ ¼
3 X
We would use superscript C to denote the elementary integrals for a continuous element. Integrals BC k are given by 1 BC 1 ¼ 3;
4 BC 2 ¼ 3;
xi fi ðZÞ.
(13)
Generally the three geometric nodes of an element are used as collocation points. Such an element is referred to as a continuous quadratic element. To model the corner points or the points where there is a change in the boundary conditions, one sided discontinuous elements (commonly referred to as the partially discontinuous elements) are used. These elements can be of two types depending on which extreme collocation point has been moved inside the element to model the geometric or physical singularity. If it is the third collocation node which has been moved inside, then the resulting element is referred to as the partially discontinuous element of Type 1. Otherwise, if the first node is moved inside, then the element is called the partially discontinuous element of Type 2. One may opt for moving both extreme nodes inside which results in a discontinuous element. In the graphical representation of these elements, we have used a for a geometric node and a for a collocation node.
(15)
16 DC 2 ð0Þ ¼ 9 ;
1 DC 3 ð0Þ ¼ 9.
(16)
Elementary integrals DC k for Zs ¼ 1 and Zs ¼ þ1 are C 1 17 DC 1 ð1Þ ¼ 3 lnð2Þ 18 ¼ D3 ðþ1Þ, C 4 10 DC 2 ð1Þ ¼ 3 lnð2Þ 9 ¼ D2 ðþ1Þ, C 1 1 DC 3 ð1Þ ¼ 3 lnð2Þ þ 18 ¼ D1 ðþ1Þ.
ð17Þ
3.2. Partially discontinuous element of Type 1 Fig. 2 provides the graphical representation of a partially discontinuous quadratic element of Type 1, in which Z1 represents the shifted local coordinate of third collocation point. Commonly used values of Z1 are 23 or 12. Geometry of the element can still be modelled using the coordinates of the geometric nodes and shape functions for the continuous elements given by Eq. (14). However, the variation of a problem variable is modeled using the following shape functions: ZðZ Z1 Þ , 1 þ Z1 ðZ þ 1ÞðZ Z1 Þ P f2 1 ðZÞ ¼ , Z1 ZðZ þ 1Þ P , f3 1 ðZÞ ¼ Z1 ð1 þ Z1 Þ P
f1 1 ðZÞ ¼
3.1. Continuous quadratic elements Fig. 1 gives the description of a continuous quadratic element. The shape functions for this element are f1 ðZÞ ¼ 12ZðZ 1Þ, 2
f2 ðZÞ ¼ 1 Z , f3 ðZÞ ¼ 12ZðZ þ 1Þ.
ð14Þ
ð18Þ
where superscript P1 has been used to denote the partially P discontinuous element of Type 1. Elementary integrals Bk 1
3 X
y
3 y
1 BC 3 ¼ 3.
Elementary integrals DC k have very simple form for any of the three collocation points. For Zs ¼ 0, these are 1 DC 1 ð0Þ ¼ 9;
i¼1
243
X 2
2
1
X x
x 1
2
3
-1
0
1
1 X η
Fig. 1. Continuous (isoparametric) quadratic element in global (x–y) and local (Z) coordinate systems.
-1
2 X 0
3 X η1
1
1
η
Fig. 2. Partially discontinuous quadratic element of Type 1 in global (x–y) and local (Z) coordinate systems: local node 3 is located at Z ¼ Z1 .
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
244
are given by P
Table 2 Elementary integrals for partially continuous element of Type 1 for Z1 ¼ 23
2 , 3 ð1 þ Z1 Þ 2 ¼2 , 3 Z1 2 . ¼ 3 Z1 ð1 þ Z1 Þ
B1 1 ¼ P
B2 1 P B3 1
P
P
D1 1
ð19Þ P Dk 1 ,
P D2 1 P D3 1
y
For elementary integrals the three possible values of Zs are 1, 0 and Z1 corresponding to the three collocation P points. Elementary integrals Dk 1 for Zs ¼ 1 are 6 lnð2Þ 9Z1 8 , ¼ 9ð1 þ Z1 Þ 2 1 P D2 1 ð1Þ ¼ 2 1, lnð2Þ 3Z1 9Z1 6 lnð2Þ þ 1 P . D3 1 ð1Þ ¼ 9Z1 ð1 þ Z1 Þ
P D1 1 ð1Þ
P
P
B1 1 ¼ 25 Zs ¼ 1
B2 1 ¼ 1 Zs ¼ 0
B3 1 ¼ 35 Zs ¼ 23
2 5 lnð2Þ
14 15
2 15
7 6
53 15
10 2 2 27 lnð5Þ 5 lnð3Þ 45 125 8 108 lnð5Þ lnð3Þ 9 5 3 16 36 lnð5Þ 5 lnð3Þ 15
lnð2Þ 3 5 lnð2Þ
1 þ 10
3 X
2 X
X
1
ð20Þ
P
Elementary integrals Dk 1 for Zs ¼ 0 are
x
2 , 9 ð1 þ Z1 Þ 2 P D2 1 ð0Þ ¼ 2 þ , 9 Z1 2 P . D3 1 ð0Þ ¼ 9 Z1 ð1 þ Z1 Þ P
1 X -η1
D1 1 ð0Þ ¼
-1
ð21Þ P
General relation for Dk 1 for Zs ¼ Z1 ; Z1 2 ð0; 1Þ can be written as P Dk 1 ðZ1 Þ
¼
P ak 1
þ
P bk 1
lnð1 Z1 Þ þ P
P ck 1
lnð1 þ Z1 Þ.
P
(22)
P
3.3. Partially discontinuous element of Type 2 Fig. 3 depicts a quadratic element of Type 2, in which Z1 represents the shifted local coordinate of the first
collocation point. Commonly used values of Z1 are The shape functions for this element are ZðZ 1Þ , Z1 ð1 þ Z1 Þ ð1 ZÞðZ þ Z1 Þ P f2 2 ðZÞ ¼ , Z1 ZðZ þ Z1 Þ P f3 2 ðZÞ ¼ , 1 þ Z1
2 3
or 12.
f1 2 ðZÞ ¼
Table 1 P P P Coefficients ak 1 , bk 1 and ck 1 in Eq. (22) for partially continuous element of Type 1 P
P
k
ak 1
bk 1
ck 1
1
3Z21 2 9ð1 þ Z1 Þ ð3Z2 þ 9Z1 2Þ 1 9Z1 ð6Z21 þ 9Z1 þ 2Þ 9Z1 ð1 þ Z1 Þ
ð2 þ Z1 Þð1 Z1 Þ2 6ð1 þ Z1 Þ ðZ1 þ 5Þð1 Z1 Þ2 6Z1 ð2Z21 þ 5Z1 þ 5Þð1 Z1 Þ 6Z1 ð1 þ Z1 Þ
ð1 þ Z1 Þð2 Z1 Þ 6 ð1 þ Z1 Þ3 6Z1 ð1 þ Z1 Þð2Z1 1Þ 6Z1
ð23Þ
where superscript P2 has been used to denote the partially discontinuous element of Type 2. P Elementary integrals Bk 2 are given by P
2 , 3 Z1 ð1 þ Z1 Þ 2 ¼2 , 3 Z1 2 . ¼ 3 ð1 þ Z1 Þ
B1 2 ¼ P
3
η
Fig. 3. Partially discontinuous quadratic element of Type 2 in global (x–y) and local (Z) coordinate systems: local node 1 is located at Z ¼ Z1 .
B2 2
2
3 X 1
P
Values of the coefficients ak 1 , bk 1 and ck 1 are listed in Table 1. Based on Table 1 and preceding equations, simplified values of the elementary integrals for the most commonly used value of Z1 ¼ 23 are listed in Table 2.
P
2 X 0
P
B3 2
ð24Þ P
For elementary integrals Dk 2 , the three possible values of Zs are Z1 , 0 and þ1, corresponding to the three P collocation points. The elementary integrals Dk 2 for Zs ¼ Z1 can be written as P
P
P
P
Dk 2 ðZ1 Þ ¼ ak 2 þ bk 2 lnð1 Z1 Þ þ ck 2 lnð1 þ Z1 Þ. P
P
P
(25)
Values of the coefficients ak 2 , bk 2 and ck 2 are listed in Table 3.
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
P
P
ak 1
bk 1
ck 1
1
ð6Z21 þ 9Z1 þ 2Þ 9Z1 ð1 þ Z1 Þ ð3Z21 þ 9Z1 2Þ 9Z1 3Z21 2 9ð1 þ Z1 Þ
ð2Z21 þ 5Z1 þ 5Þð1 Z1 Þ 6Z1 ð1 þ Z1 Þ ðZ1 þ 5Þð1 Z1 Þ2 6Z1 ð2 þ Z1 Þð1 Z1 Þ2 6ð1 þ Z1 Þ
ð1 þ Z1 Þð2Z1 1Þ 6Z1
2 3
P
P
P D2 2 P D3 2
2
X
1
ð1 þ Z1 Þ3 6Z1 ð1 þ Z1 Þð2 Z1 Þ 6
Table 4 Elementary integrals for partially continuous element of Type 2 for Z1 ¼ 2=3
D1 2
X
P
k
3 X
y
Table 3 P P P Coefficients ak 2 , bk 2 and ck 2 in Eq. (25) for partially continuous element of Type 2
245
P
P
B1 2 ¼ 35 Zs ¼ 23
B2 2 ¼ 1 Zs ¼ 0
B3 2 ¼ 25 Zs ¼ 1
5 3 16 36 lnð5Þ 5 lnð3Þ 15 125 8 lnð5Þ lnð3Þ 108 9 10 2 2 27 lnð5Þ 5 lnð3Þ 45
15
3 5 lnð2Þ
53
lnð2Þ 76
2 15
2 5 lnð2Þ
1 þ 10
14 15
For the middle collocation point (Zs ¼ 0), elementary P integrals Dk 2 are 2 , 9 Z1 ð1 þ Z1 Þ 2 P D2 2 ð0Þ ¼ 2 þ , 9 Z1 2 P . D3 2 ð0Þ ¼ 9 ð1 þ Z1 Þ P
x
-1
1 X -η1
2 X 0
3 X η2
1
η
Fig. 4. Discontinuous quadratic element in global (x–y) and local (Z) coordinate systems: local node 1 is located at Z ¼ Z1 and node 3 at Z ¼ Z2 .
discontinuous elements. The shape functions are ZðZ Z2 Þ , Z1 ðZ1 þ Z2 Þ ðZ þ Z1 ÞðZ Z2 Þ fD , 2 ðZÞ ¼ Z1 Z2 ZðZ þ Z1 Þ . fD 3 ðZÞ ¼ Z2 ðZ1 þ Z2 Þ fD 1 ðZÞ ¼
ð28Þ
D1 2 ð0Þ ¼
Elementary integrals BD k are given by
ð26Þ P
The elementary integrals Dk 2 for Zs ¼ þ1 are 6 lnð2Þ þ 1 , 9Z1 ð1 þ Z1 Þ 2 1 P D2 2 ðþ1Þ ¼ 2 1, lnð2Þ 3Z1 9Z1 6 lnð2Þ 9Z1 8 P . D3 2 ðþ1Þ ¼ 9ð1 þ Z1 Þ P
D1 2 ðþ1Þ ¼
ð27Þ
Based on the preceding, elementary integrals for the partially discontinuous element of Type 2 with the discontinuous first node at Z ¼ 23 (i.e. Z1 ¼ 23) are listed in Table 4.
2 , 3 Z1 ðZ1 þ Z2 Þ 2 BD , 2 ¼2 3 Z1 Z2 2 BD . 3 ¼ 3 Z2 ðZ1 þ Z2 Þ BD 1 ¼
ð29Þ
For elementary integrals DD k , the three possible values of Zs are Z1 , 0 and þZ2 , corresponding to the three collocation points. The elementary integrals DD k for Zs ¼ Z1 can be written as D
D
D
3.4. Discontinuous quadratic elements Though the continuous and partially discontinuous quadratic elements are the most widely used, some one might opt to use fully discontinuous elements. Formulae for these elements are included here for sake of completeness. In these elements, both the extreme collocation points are shifted inside as depicted in Fig. 4. We would use superscript D to denote the quantities related to
D
1 1 1 DD k ðZ1 Þ ¼ ak þ bk lnð1 Z1 Þ þ ck lnð1 þ Z1 Þ,
D
(30)
D
where coefficients ak 1 , bk 1 and ck 1 are D
ð6Z21 þ 9Z1 Z2 þ 2Þ , 9Z1 ðZ1 þ Z2 Þ ð2Z21 þ 3Z1 Z2 þ 2Z1 þ 3Z2 þ 2Þð1 Z1 Þ , ¼ 6Z1 ðZ1 þ Z2 Þ ð2Z21 þ 3Z1 Z2 2Z1 3Z2 þ 2Þð1 þ Z1 Þ , ¼ 6Z1 ðZ1 þ Z2 Þ
a1 1 ¼ D
b1 1 D
c1 1
ð31Þ
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247
246 D
a2 1 ¼ D
b2 1 ¼ D
c2 1 ¼
D
a3 1 ¼ D
b3 1 ¼ D
c3 1 ¼
ð3Z21 þ 9Z1 Z2 2Þ , 9Z1 Z2
D
c3 2 ¼
ðZ1 þ 3Z2 þ 2Þð1 Z1 Þ2 , 6Z1 Z2
ð1 þ Z1 Þ2 ðZ1 þ 3Z2 2Þ , 6Z1 Z2
ð32Þ
3Z21 2 , 9Z2 ðZ1 þ Z2 Þ
ð1 þ Z1 Þ2 ð2 Z1 Þ . 6Z2 ðZ1 þ Z2 Þ
We may note that the relations for the discontinuous element represent the general formulae from which the relations for the continuous and partially discontinuous elements can be obtained by taking care of the proper limits.
ð33Þ
Let us consider weakly singular integral arising in the boundary element analysis of Laplace equation in 2-D is given by Z I L ðx0 ; kÞ ¼ fk ðxÞ log jx x0 j dGðxÞ Ge 1
Z
For the middle collocation point (Zs ¼ 0), elementary integrals DD k are 2 , 9 Z1 ðZ1 þ Z2 Þ 2 DD , 2 ð0Þ ¼ 2 þ 9 Z 1 Z2 2 DD . 3 ð0Þ ¼ 9 Z2 ðZ1 þ Z2 Þ
DD 1 ð0Þ ¼
¼
þ
where coefficients D
a1 2 ¼ D
b1 2 ¼ D
c1 2 ¼
D
D
(35)
are
ð36Þ
ð3Z22 þ 9Z1 Z2 2Þ , 9Z1 Z2 ðZ2 þ 3Z1 þ 2Þð1 Z2 Þ2 , 6Z1 Z2
ð1 þ Z2 Þ2 ðZ2 þ 3Z1 2Þ , 6Z1 Z2
ð6Z22 þ 9Z1 Z2 þ 2Þ , 9Z2 ðZ1 þ Z2 Þ ð2Z22 þ 3Z1 Z2 þ 2Z2 þ 3Z1 þ 2Þð1 Z2 Þ , ¼ 6Z2 ðZ1 þ Z2 Þ
a3 2 ¼ b3 2
and
D ck 2
lnð1 þ Z2 Þ,
ð1 þ Z2 Þ2 ð2 Z2 Þ , 6Z1 ðZ1 þ Z2 Þ
b2 2 ¼
D
D bk 2
D ck 2
ð2 þ Z2 Þð1 Z2 Þ2 , 6Z1 ðZ1 þ Z2 Þ
D
D
lnð1 Z2 Þ þ
D ak 2 ,
fk ðZÞ log jx x0 jJðZÞ dZ.
ð39Þ
1
ð34Þ
3Z22 2 , 9Z1 ðZ1 þ Z2 Þ
a2 2 ¼
c2 2 ¼
D bk 2
¼
We take a curved quadratic element with nodes x1 ¼ ð1; 1Þ; x2 ¼ ð2:5; 3Þ and x3 ¼ ð2; 5Þ, in which x2 is the middle node. The higher order transformation of Sato et al. [5] with k ¼ 3 has been employed in the SSNT method.
Finally, the elementary integrals DD k for Zs ¼ þZ2 can be written as D ak 2
ð38Þ
4. Numerical examples
ð2 þ Z1 Þð1 Z1 Þ2 , 6Z2 ðZ1 þ Z2 Þ
DD k ðþZ2 Þ
ð2Z22 þ 3Z1 Z2 2Z2 3Z1 þ 2Þð1 þ Z2 Þ . 6Z2 ðZ1 þ Z2 Þ
ð37Þ
Table 5 Relative error (in percent) in numerical integration of weakly singular Iðx0 ; kÞ with choice LOGA and SSNT method k¼1
k¼2
k¼3
1:1E 07 4:2E 12 5:2E 13 2:6E 12 2:0E 14 2:0E 14 4:6E 08 4:5E 12 3:1E 14
3:3E 07 2:0E 11 2:5E 13 2:4E 12 1:4E 14 2:8E 14 3:9E 08 1:0E 11 2:2E 13
3:0E 07 1:7E 11 7:6E 14 2:1E 12 0:0E þ 00 0:0E þ 00 5:4E 08 1:0E 11 4:1E 13
Partially discontinuous element of Type 1 1:00 16 2:5E 07 24 2:3E 13 32 7:6E 13 0.00 16 2:5E 12 24 0:0E þ 00 32 1:8E 14 0.67 16 1:3E 09 24 1:1E 12 32 1:9E 14
1:1E 06 6:3E 11 5:9E 13 2:3E 12 1:2E 14 0:0E þ 00 4:6E 09 5:7E 12 2:2E 13
3:0E 07 1:7E 11 1:0E 13 2:1E 12 0:0E þ 00 1:1E 14 1:3E 09 3:5E 12 5:8E 14
Partially discontinuous element of Type 2 0:67 16 1:1E 08 24 4:6E 12 32 1:6E 13 0.00 16 2:5E 12 24 0:0E þ 00 32 0:0E þ 00 1.00 16 4:6E 08 24 4:5E 12 32 1:7E 14
2:9E 08 7:5E 12 2:7E 13 2:4E 12 2:1E 14 0:0E þ 00 3:2E 07 5:8E 11 1:1E 12
2:4E 10 2:2E 12 2:4E 14 2:2E 12 0:0E þ 00 2:9E 14 5:7E 08 1:6E 11 4:4E 13
Zs
NG
Continuous elements 1:00 16 24 32 0.00 16 24 32 1.00 16 24 32
ARTICLE IN PRESS K.M. Singh, M. Tanaka / Engineering Analysis with Boundary Elements 31 (2007) 241–247 Table 6 Relative error (in percent) in numerical integration of weakly singular Iðx0 ; kÞ with choice LOGB and SSNT method
247
These results clearly demonstrate the usefulness of the formulae presented herein and the SSNT method.
k¼1
k¼2
k¼3
5. Conclusions
1:2E 07 7:4E 12 1:0E 13 2:5E 12 2:0E 14 0:0E þ 00 4:6E 08 2:5E 12 1:5E 14
3:3E 07 2:6E 11 2:0E 14 2:4E 12 1:4E 14 1:4E 14 4:0E 08 4:0E 12 2:4E 14
3:0E 07 2:0E 11 0:0E þ 00 2:0E 12 2:1E 14 2:1E 14 5:3E 08 2:1E 12 6:5E 14
Partially discontinuous element of Type 1 1:00 16 2:5E 07 24 1:6E 11 32 1:4E 13 0.00 16 2:4E 12 24 1:8E 14 32 1:8E 14 0.67 16 1:4E 09 24 0:0E þ 00 32 0:0E þ 00
We have presented closed form analytical formulae for the elementary integrals required in the numerical evaluation of weakly singular boundary integrals involving logarithmic singularity using the subtraction of singularity method. Results have been presented for continuous, partially discontinuous and discontinuous quadratic elements. Sample numerical results obtained with a hybrid subtraction of singularity–nonlinear transformation method clearly demonstrate the usefulness of these formulae.
1:1E 06 8:2E 11 1:0E 13 2:3E 12 1:2E 14 2:4E 14 5:3E 09 0:0E þ 00 0:0E þ 00
3:0E 07 2:0E 11 2:1E 14 2:0E 12 3:4E 14 2:3E 14 1:6E 09 5:0E 13 0:0E þ 00
Partially discontinuous element of Type 2 0:67 16 1:1E 08 24 8:6E 13 32 3:9E 14 0.00 16 2:5E 12 24 0:0E þ 00 32 2:2E 14 1.00 16 4:6E 08 24 2:5E 12 32 1:7E 14
3:0E 08 2:1E 13 1:6E 14 2:4E 12 2:1E 14 2:1E 14 3:2E 07 2:5E 11 5:1E 14
4:6E 10 3:7E 13 2:4E 14 2:1E 12 1:5E 14 0:0E þ 00 5:6E 08 3:9E 12 1:1E 13
Zs
NG
Continuous elements 1:00 16 24 32 0.00 16 24 32 1.00 16 24 32
The ‘‘exact values’’ of the integrals correspond to the numerical values obtained with Maple V [7] using precision of 25 digits, and truncated to 16 digits in double precision computations. Errors in the numerical evaluation of the weakly singular integral I L ðxi ; kÞ obtained with the choices LOGA and LOGB are presented in Tables 5 and 6, respectively, for continuous as well as partial continuous elements. Accuracy of results is evident. Near-exact results are obtained with 32 Gauss points for any element type.
Acknowledgments Part of this work was financially supported by the Ministry of Education, Science, Sports and Culture, Japan, Grant-in-Aid for Scientific Research (C), No. 16560066, 2005, granted to the author (M.T.). This financial support is gratefully acknowledged. References [1] Kuwahara T, Takeda T. On a formula of boundary integral for the higher-order element in boundary element method and a consideration of error of solution. Trans Inst Electr Eng Jpn 1985;105-A:1–8 [in Japanese]. [2] Aimi A, Diligenti M, Monegato G. New numerical integration schemes for applications of Galerkin BEM to 2-D problems. Int J Numer Meth Eng 1997;40:1977–99. [3] Doblare´ M. Computational aspects of the boundary element method. Topics in boundary element research, vol. 3. Berlin: Springer; 1987 (chapter 4). [4] Katsikadelis JT. Boundary Elements: Theory and Applications. Oxford: Elsevier Science Ltd; 2002. [5] Sato M, Yoshiyoka S, Tsukui K, Yuuki R. Accurate numerical integration of singular kernels in the two-dimensional boundary element method. In: Brebbia CA, editor. Boundary elements X, vol. 1. Berlin: Springer; 1988. p. 279–96. [6] Singh KM, Tanaka M. On non-linear transformations for accurate numerical evaluation of weakly singular boundary integrals. Int J Numer Meth Eng 2001;50:2007–30. [7] Heal KM, Hansen ML, Rickard KM. Maple V: Learning Guide. New York: Springer; 1996.