# High Order Projection Plane Method for Evaluation of Supersingular Curved Boundary Integrals in BEM

Authors: Miao Cui; Wei-zhe Feng; Xiao-wei Gao; Kai Yang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Boundary element method (BEM) is a very promising approach for solving various engineering problems, in which accurate evaluation of boundary integrals is required. In the present work, the direct method for evaluating singular curved boundary integrals is developed by considering the third-order derivatives in the projection plane method when expanding the geometry quantities at the field point as Taylor series. New analytical formulas are derived for geometry quantities defined on the curved line/plane, and unified expressions are obtained for both two-dimensional and three-dimensional problems. For the two-dimensional boundary integrals, analytical expressions for the third-order derivatives are derived and are employed to verify the complex-variable-differentiation method (CVDM) which is used to evaluate the high order derivatives for three-dimensional problems. A few numerical examples are given to show the effectiveness and the accuracy of the present method.


---

## Body

## 1. Introduction

Boundary element method (BEM) [1] is a very promising approach for solving various engineering problems, such as heat transfer problems [2–5] in energy science and engineering, in which accurate evaluation of boundary integrals [6–10] is required. Weak, strong, or hypersingularities are involved in these boundary integrals, if the source point is located on the element under integration [11]. These issues have been focused on by researchers, and many techniques, such as Gauss logarithmic quadrature rule for weak singularities [12], rigid body displacement algorithm for strong singularities [12], and regularization for hypersingularities [13], have been proposed to remove the singularities. Each method has its advantage and disadvantage [11–19]. A comprehensive review has been given in [11] on the treatment of singularities. In general, the efficient methods for singular curved line and surface integrals are those based on the operation over the intrinsic coordinate system [14–16].Recently, Gao et al. [11] proposed another efficient direct method for numerically evaluating all kinds of singular curved boundary integrals for two-dimensional (2D) and three-dimensional (3D) BEM analyses, and numerical examples were given to demonstrate the stability of the proposed method. However, only quadratic terms were truncated when expanding the coordinate of the field point as Taylor series, which means the method will be invalid if the curved boundary is beyond the quadratic, such as the cubed curve.In the present work, the direct method in [11] for evaluating arbitrary high order singular curved boundary integrals is developed, by considering the third-order derivatives when expanding the geometry quantities at the field point as Taylor series. Corresponding analytical expressions are derived. The first innovation is that the present method will extend the application range of the method in [11]. That is, the singular integrals over higher order elements, such as cubed curve, can be directly evaluated. Another innovation is that the complex-variable-differentiation method [20] is introduced to determine the third-order derivatives. For 2D boundary integral, the third-order derivatives are analytically derived, which are used to verify the CVDM.The paper is organized as follows. In Section2, review and modification of the direct method in [11] are provided. In Section 3, numerical examples are given to show the effectiveness and the accuracy of the present method. Finally, conclusions are drawn from computational results and discussions.


## 2. Development of the Direct Method for Evaluating Arbitrary High Order Singular Curved Boundary Integrals

### 2.1. Review of Singular Integrals [11]

Singular integrals over an element can be classified into the following form:(1)
I
e
x
p
=
∫
Γ
e
f
-
x
p
,
x
r
β
x
p
,
x
d
Γ
xin which x
p and x represent the source and field points, respectively; Γ
e is the boundary element under integration, which is a curved line or a curved surface for 2D or 3D problems, respectively; and r is the distance between the source and field points.The termf
-
(
x
p
,
x
) in (1) is a regular function and β is the order of the singularity.In 2D problems, the following type of singular line integrals is also frequently encountered:(2)
J
e
x
p
=
∫
Γ
e
f
-
x
p
,
x
ln
⁡
r
x
p
,
x
r
β
x
p
,
x
d
Γ
x
.It is assumed that the integrals in (1) and (2) always exist: that is, each integration has a finite value.


### 2.2. Review of Expansions of Geometry Quantities on the Projection Line/Plane [11]

In [11], a projection line for 2D or a projection plane for 3D problems, which is the tangential line/plane of the element to the origin of the intrinsic coordinate system, was introduced. A local orthogonal coordinate system was established on the projection line/plane. The coordinate transformation between the local and global systems can be performed using the following relationships: (3)
x
i
′
=
L
i
j
x
j
-
x
j
o
,
(4)
x
i
=
x
i
o
+
L
j
i
x
j
′
,where L
i
j is the direction cosine of the local coordinate axes with respect to the global one. The repeated subscripts represent summation; x
i
o is the global coordinates of the origin of the local coordinate system.The original curved element can be projected onto the projection line or plane to form a straight or a flat projection element by (3), and then all geometry quantities can be expressed in terms of variables defined on the projection line/plane.


### 2.3. High Order Expansion of Geometry Quantities in terms of Variables on the Projection Line for 2D Problems

For 2D and 3D problems, derivations for the modification are different, which are separately described in detail.In 2D problems, the boundary element is a curved line as shown in Figure1. The local orthogonal coordinate system is denoted by (x
1
′
,
x
2
′) and the projection line is along the axis x
1
′.Figure 1
The curved line boundary and the projection line.There is only one independent variable for a line, and all other geometry quantities can be expressed byx
1
′. Similar to the treatment in [11], the coordinate x
2
′ at the field point over the curved line can be expressed by x
1
′ using Taylor series. Different from the treatment in [11], the third-order derivative is considered in the present work:(5)
x
2
′
=
∂
x
2
′
∂
x
1
′
x
1
′
+
1
2
∂
2
x
2
′
∂
x
1
′
2
x
1
′
2
+
1
6
∂
3
x
2
′
∂
x
1
′
3
x
1
′
3
=
a
1
x
2
′
x
1
′
+
a
2
x
2
′
x
1
′
2
+
a
3
x
2
′
x
1
′
3
,where(6)
a
1
x
2
′
=
∂
x
2
′
∂
x
1
′
=
∂
x
2
′
∂
ξ
1
∂
ξ
1
∂
x
1
′
,
a
2
x
2
′
=
1
2
∂
2
x
2
′
∂
x
1
′
2
=
1
2
∂
2
x
2
′
∂
ξ
2
∂
ξ
∂
x
1
′
∂
ξ
∂
x
1
′
+
∂
x
2
′
∂
ξ
∂
2
ξ
∂
x
1
′
2
,
a
3
x
2
′
=
1
6
∂
3
x
2
′
∂
x
1
′
3
=
1
6
∂
∂
x
1
′
∂
2
x
2
′
∂
x
1
′
2
=
1
6
∂
3
x
2
′
∂
ξ
3
∂
ξ
∂
x
1
′
∂
ξ
∂
x
1
′
∂
ξ
∂
x
1
′
+
3
∂
2
x
2
′
∂
ξ
2
∂
ξ
∂
x
1
′
∂
2
ξ
∂
x
1
′
2
+
∂
x
2
′
∂
ξ
∂
3
ξ
∂
x
1
′
3
.It can be seen that the expressions are more complex. For 2D problem, the corresponding derivatives are analytically derived, and the detailed derivation process is omitted here.In (5) and (6), (7a)
∂
x
2
′
∂
ξ
=
∑
α
=
1
N
node
∂
N
α
∂
ξ
x
2
′
α
,
(7b)
∂
2
x
2
′
∂
ξ
2
=
∑
α
=
1
N
node
∂
2
N
α
∂
ξ
2
x
2
′
α
,
(7c)
∂
3
x
2
′
∂
ξ
3
=
∑
α
=
1
N
node
∂
3
N
n
∂
ξ
3
x
2
′
α
,
(8a)
∂
ξ
∂
x
1
′
=
1
m
1
,
(8b)
∂
2
ξ
∂
x
1
′
2
=
-
1
m
1
4
∂
x
1
∂
ξ
∂
2
x
1
∂
ξ
2
+
∂
x
2
∂
ξ
∂
2
x
2
∂
ξ
2
,
(8c)
∂
3
ξ
∂
x
1
′
3
=
∂
∂
x
1
′
∂
2
ξ
∂
x
1
′
2
=
-
m
1
2
B
-
4
C
m
1
7
.In (7a)–(8c),(9)
m
1
=
∂
x
1
∂
ξ
2
+
∂
x
2
∂
ξ
2
,
(10a)
B
=
∂
2
x
1
∂
ξ
2
2
+
∂
2
x
2
∂
ξ
2
2
+
∂
x
1
∂
ξ
∂
3
x
1
∂
ξ
3
+
∂
x
2
∂
ξ
∂
3
x
2
∂
ξ
3
,
(10b)
C
=
∂
x
1
∂
ξ
2
∂
2
x
1
∂
ξ
2
2
+
2
∂
x
1
∂
ξ
∂
x
2
∂
ξ
∂
2
x
1
∂
ξ
2
∂
2
x
2
∂
ξ
2
+
∂
x
2
∂
ξ
2
∂
2
x
2
∂
ξ
2
2
.Define a local distance:(11)
ρ
=
x
1
′
-
x
1
′
p
,where x
1
′
p is the local coordinate of point P
′ which is the projection point of the source point P and x
1
′ is the local coordinate of point Q
′ which is the projection point of the field point Q.The derivative ofρ with respect to local coordinate x
1
′ is(12a)
ρ
,
1
=
∂
ρ
∂
x
1
′
=
x
1
′
-
x
1
′
p
ρ
.In (12a), ρ
,
1
=
1 or -
1.Then, we can obtain(12b)
x
1
′
=
x
1
′
p
+
ρ
,
1
ρ
.By substituting (12a) into (5), we can obtain(13)
x
2
′
=
p
x
2
′
+
q
x
2
′
,
ρ
ρ
,where(14)
p
x
2
′
=
a
1
x
1
′
p
+
a
2
x
2
′
x
1
′
p
x
1
′
p
+
a
3
x
2
′
x
1
′
p
x
1
′
p
x
1
′
p
=
x
2
′
p
,
(15)
q
x
2
′
,
ρ
=
q
1
x
2
′
+
q
2
x
2
′
ρ
+
q
3
x
2
′
ρ
2
.In (15), (16a)
q
1
x
2
′
=
a
1
x
2
′
+
2
a
2
x
2
′
x
1
′
p
+
3
a
3
x
2
′
x
1
′
p
2
ρ
,
1
,
(16b)
q
2
x
2
′
=
a
2
x
2
′
ρ
,
1
2
+
3
a
3
x
2
′
x
1
′
p
=
a
2
x
2
′
+
3
a
3
x
2
′
x
1
′
p
,
(16c)
q
3
x
2
′
,
ρ
=
a
3
x
2
′
ρ
,
1
.Similarly, we can express the intrinsic coordinate in terms of the local distance as(17)
ξ
=
p
ξ
+
q
ξ
,
ρ
ρ
,where p
(
ξ
) and q
(
ξ
,
ρ
) can be determined by replacing x
2
′ with ξ in (13)–(16c).Substituting (12b) and (13) into (4), we can obtain(18)
x
i
=
x
i
o
+
L
j
i
x
j
′
=
x
i
o
+
L
1
i
x
1
′
+
L
2
i
x
2
′
=
x
i
o
+
L
1
i
x
1
′
p
+
ρ
,
1
ρ
+
L
2
i
x
2
′
p
+
q
x
2
′
,
ρ
ρ
=
x
i
p
+
L
1
i
ρ
,
1
ρ
+
L
2
i
q
x
2
′
,
ρ
ρ
.Then, substituting (15) into (18), we can obtain(19)
x
i
=
x
i
p
+
b
i
+
c
i
ρ
+
d
i
ρ
2
ρ
,where(20a)
b
i
=
L
1
i
ρ
,
1
+
L
2
i
q
1
y
′
,
(20b)
c
i
=
L
2
i
q
2
y
′
,
(20c)
d
i
=
L
2
i
q
3
y
′
.Referring to Figure1, the global distance r can be easily expressed as(21)
r
=
x
i
-
x
i
p
x
i
-
x
i
p
=
g
ρ
ρ
,where(22)
g
ρ
=
b
i
+
c
i
ρ
+
d
i
ρ
2
b
i
+
c
i
ρ
+
d
i
ρ
2
=
G
0
+
G
1
ρ
+
G
2
ρ
2
+
G
3
ρ
3
+
G
4
ρ
4
.In (22), G
0
=
b
i
b
i, G
1
=
2
b
i
c
i, G
2
=
c
i
c
i
+
2
b
i
d
i, G
3
=
2
c
i
d
i, and G
4
=
d
i
d
i.It can be easily derived that(23)
r
,
i
=
∂
r
∂
x
i
=
x
i
-
x
i
p
r
=
b
i
+
c
i
ρ
+
d
i
ρ
2
g
ρ
.


### 2.4. The Complex-Variable-Differentiation Method (CVDM)

From the above derivations, it can be seen that the corresponding expressions have been changed, by considering the third-order derivatives when expanding the geometry quantities at the field point as Taylor series. The third-order derivatives should be first determined.In the present work, the complex-variable-differentiation method (CVDM) [20] is introduced to determine the high order derivatives. For 2D analysis, the third-order derivatives are determined explicitly by (7c), (8c), (9), (10a), and (10b), which are also calculated by the CVDM.In the CVDM, the variablex of a real function f
(
x
) is replaced by a complex one x
+
i
h, with the imaginary part h being very small (usually h = 10−20). The function f
(
x
+
i
h
) can be expanded into Taylor’s series as(24)
f
x
+
i
h
=
f
x
+
i
h
f
′
x
-
h
2
2
f
′
′
x
+
o
h
3
.Since h is small enough, (25) can be obtained: (25)
f
′
x
=
Im
⁡
f
x
+
i
h
h
,where Im denotes the imaginary part.The CVDM is a very promising method, since the derivatives only require direct evaluations [21]. The CVDM has been validated to be an efficient and accurate approach for determination of derivatives [22–24].For example, considering∂
3
x
2
′
/
∂
x
1
′
3
=
∂
/
∂
x
1
′
∂
2
x
2
′
/
∂
x
1
′
2, we can add the very small imaginary part h to x
1
′; then the real value of ∂
3
x
2
′
/
∂
x
1
′
3 yields Im
⁡
(
∂
2
x
2
′
/
∂
x
1
′
2
)
/
h. In addition, the value of ∂
2
x
2
′
/
∂
x
1
′
2 can be easily obtained using the analytical expressions in the previous work [11]. The CVDM and the derived expressions in (7c)–(8c) can be validated by each other.


### 2.5. Expansion of Geometry Quantities in terms of Variables on the Projection Plane for 3D Problems

In 3D problems, the boundary element is a curved surface marked asA
B
C
D in Figure 2. A local orthogonal coordinate system (x
1
′
,
x
2
′
,
x
3
′) is established on the projection plane with its origin being at the point (ξ
1
=
0, ξ
2
=
0), in which axes x
1
′ and x
2
′ are located within the plane. The axis x
1
′ is along ξ
1 direction and axis x
3
′ is along the outward normal direction to the element. By (4), we can project the original curved surface element onto the projection plane, to form a flat projection element marked as A
′
B
′
C
′
D
′ in Figure 2.Figure 2
The curved surface boundary and the projection plane.In a plane, there are only two independent variables. We choosex
1
′ and x
2
′ as the independent variables, and x
3
′ over the curved surface can be expressed in terms of x
1
′ and x
2
′. Compared with the treatment in [11], x
3
′ is expanded as Taylor series by x
1
′ and x
2
′, by considering the third-order derivatives:(26)
x
3
′
=
a
I
x
3
′
x
I
′
+
a
I
J
x
3
′
x
I
′
x
J
′
+
a
I
J
K
x
3
′
x
I
′
x
J
′
x
K
′
,where I, J, and K are integers varying from 1 to 2, and(27a)
a
I
x
3
′
=
∂
x
3
′
∂
x
I
′
=
∂
x
3
′
∂
ξ
K
∂
ξ
K
∂
x
I
′
,
(27b)
a
I
J
x
3
′
=
1
2
∂
2
x
3
′
∂
x
I
′
∂
x
J
′
=
1
2
∂
2
x
3
′
∂
ξ
K
∂
ξ
L
∂
ξ
K
∂
x
I
′
∂
ξ
L
∂
x
J
′
+
∂
x
3
′
∂
ξ
K
∂
2
ξ
K
∂
x
I
′
∂
x
J
′
,
(27c)
a
I
J
K
x
3
′
=
1
6
∂
3
x
3
′
∂
x
I
′
∂
x
J
′
∂
x
K
′
=
1
3
∂
a
I
J
x
3
′
∂
x
K
′
=
1
3
∂
a
I
J
x
3
′
∂
ξ
L
∂
ξ
L
∂
x
K
′
.The treatment is the same as that in [11], to determine a
I
(
x
3
′
) in (27a) and a
I
J
(
x
3
′
) in (27b). For a
I
J
K
(
x
3
′
) in (27c), we use the CVDM, to avoid rather complicated analytical derivations.The local distanceρ projected from the global distance r onto the projection plane and its derivative ρ
,
I are introduced as follows: (28)
ρ
=
x
1
′
-
x
1
′
p
2
+
x
2
′
-
x
2
′
p
2
,
ρ
,
I
=
∂
ρ
∂
x
I
′
=
x
I
′
-
x
I
′
p
ρ
.The local coordinates at the field point can be expressed as(29)
x
I
′
=
x
I
′
p
+
ρ
,
I
ρ
.Similar to the 2D analysis, substituting (29) into (26) yields(30)
x
3
′
=
p
x
3
′
+
q
x
3
′
,
ρ
ρ
,where(31)
p
x
3
′
=
a
I
x
3
′
x
I
′
p
+
a
I
J
x
3
′
x
I
′
p
x
J
′
p
+
a
I
J
K
x
3
′
x
I
′
p
x
J
′
p
x
K
′
p
=
x
3
′
p
,
(32)
q
x
3
′
,
ρ
=
q
1
x
3
′
+
q
2
x
3
′
ρ
+
q
3
x
3
′
ρ
2
.In (32),(33a)
q
1
x
3
′
=
a
I
x
3
′
ρ
,
I
+
a
I
J
x
3
′
x
I
′
p
ρ
,
J
+
a
I
J
x
3
′
x
J
′
p
ρ
,
I
+
a
I
J
K
x
3
′
x
I
′
p
x
K
′
p
ρ
,
J
+
a
I
J
K
x
3
′
x
J
′
p
x
K
′
p
ρ
,
I
+
a
I
J
K
x
3
′
x
I
′
p
x
J
′
p
ρ
,
K
,
(33b)
q
2
x
3
′
=
a
I
J
x
3
′
ρ
,
I
ρ
,
J
+
a
I
J
K
x
3
′
x
K
′
p
ρ
,
I
ρ
,
J
+
a
I
J
K
x
3
′
x
I
′
p
ρ
,
J
ρ
,
K
+
a
I
J
K
x
3
′
x
J
′
p
ρ
,
I
ρ
,
K
,
(33c)
q
3
x
3
′
=
a
I
J
K
x
3
′
ρ
,
I
ρ
,
J
ρ
,
K
.Similar to the 2D analysis, intrinsic coordinates can be expressed as(34)
ξ
K
=
p
ξ
K
+
q
ξ
K
,
ρ
ρ
,where p
(
ξ
K
) and q
(
ξ
K
,
ρ
) can be determined by replacing x
3
′ with ξ
K in (31)–(33c).Substituting (29) and (30) into (4) yields(35)
x
i
=
x
i
o
+
L
J
i
x
J
′
+
L
3
i
x
3
′
=
x
i
o
+
L
J
i
x
J
′
p
+
ρ
,
J
ρ
+
L
3
i
x
3
′
p
+
q
x
3
′
,
ρ
ρ
=
x
i
p
+
L
J
i
ρ
,
J
ρ
+
L
3
i
q
x
3
′
,
ρ
ρ
.Substituting (32) into (35), we can obtain (36)
x
i
=
x
i
p
+
b
i
+
c
i
ρ
+
d
i
ρ
2
ρ
,where(37a)
b
i
=
L
J
i
ρ
,
J
+
L
3
i
q
1
x
3
′
,
(37b)
c
i
=
L
3
i
q
2
x
3
′
,
(37c)
d
i
=
L
3
i
q
3
x
3
′
.The expression ofr is as follows:(38)
r
=
x
i
-
x
i
p
x
i
-
x
i
p
=
g
ρ
ρ
,where(39)
g
ρ
=
b
i
+
c
i
ρ
+
d
i
ρ
2
b
i
+
c
i
ρ
+
d
i
ρ
2
=
G
0
+
G
1
ρ
+
G
2
ρ
2
+
G
3
ρ
3
+
G
4
ρ
4
.In (39), G
0
=
b
i
b
i, G
1
=
2
b
i
c
i, G
2
=
c
i
c
i
+
2
b
i
d
i, G
3
=
2
c
i
d
i, and G
4
=
d
i
d
i.It can be easily derived thatr
,
i has the same expression as that in 2D analysis.Considering thatI
=
J
=
K
=
1 and ρ
,
1
=
1 or -
1 for 2D problems, it can be seen that unified expressions are obtained for both two-dimensional and three-dimensional problems.Then, we can handle the singular boundary integrals in (1) and (2), with the above modifications. The 2D problem is taken as an example to describe the procedure to determine the integral in (1). The differential relationship between the local distance and the real curved line can be written as (40)
d
ρ
=
d
Γ
cos
⁡
ϕ
=
d
Γ
n
i
0
n
i
,where n
i
0 and n
i are the outward normals to the tangential lines of passing the origin of the local coordinate system and the field point, respectively, and ϕ is the angle between the two outward normals.Substituting (40) and (38) into (1) leads to(41)
I
e
x
p
=
lim
ε
→
0
⁡
∫
ε
ρ
E
F
-
ρ
ρ
β
d
ρ
,where ρ
E is the distance from the projected source point P
′ to one of the ends of the projected line element and F
-
(
ρ
) is the regular part: that is, (42)
F
-
ρ
=
f
-
x
p
,
x
ρ
g
β
ρ
n
i
0
n
i
ρ
.In order to integrate (41), the nonsingular part F
- is expanded as a power series in ρ, such that(43)
F
-
ρ
=
∑
n
=
0
N
B
n
ρ
nin which N is the order of the power series, usually taking a value between 2 and 7 depending on the size of ρ
E, and B
(
n
) are constants which are determined by collocating N
+
1 points over the integration region (
0
,
ρ
E
). In the paper, N
+
1 equally spaced points are used: that is, (
0
,
ρ
1
,
…
,
ρ
N
). The coefficient for the first point n
=
0 is B
(
0
)
=
F
-
(
0
) and other coefficients can be solved using the following equation set:(44)
R
B
=
Y
,where R is a square matrix with the order of N:(45)
R
=
1
ρ
1
⋯
ρ
1
N
-
1
1
ρ
2
⋯
ρ
2
N
-
1
⋮
⋮
⋯
⋮
1
ρ
N
⋯
ρ
N
N
-
1
.
B and Y are vectors as follows:(46)
B
=
B
1
B
2
⋮
B
N
,
Y
=
F
-
ρ
1
-
B
0
ρ
1
F
-
ρ
2
-
B
0
ρ
2
⋮
F
-
ρ
N
-
B
0
ρ
N
.Solving (44) for coefficient vector B and then substituting (43) into (41) yield (47)
I
e
x
p
=
∑
n
=
0
N
B
n
lim
ε
→
0
⁡
∫
ε
ρ
E
ρ
n
-
β
d
ρ
=
∑
n
=
0
N
B
n
E
n
,where (48)
E
n
=
1
n
-
β
+
1
1
ρ
E
β
-
n
-
1
-
lim
ρ
ε
→
0
⁡
1
ρ
ε
β
-
n
-
1
n
≠
β
-
1
ln
⁡
ρ
E
-
lim
ρ
ε
→
0
⁡
ln
⁡
ρ
ε
n
=
β
-
1
.Then, singularities can be analytically removed.


## 2.1. Review of Singular Integrals [11]

Singular integrals over an element can be classified into the following form:(1)
I
e
x
p
=
∫
Γ
e
f
-
x
p
,
x
r
β
x
p
,
x
d
Γ
xin which x
p and x represent the source and field points, respectively; Γ
e is the boundary element under integration, which is a curved line or a curved surface for 2D or 3D problems, respectively; and r is the distance between the source and field points.The termf
-
(
x
p
,
x
) in (1) is a regular function and β is the order of the singularity.In 2D problems, the following type of singular line integrals is also frequently encountered:(2)
J
e
x
p
=
∫
Γ
e
f
-
x
p
,
x
ln
⁡
r
x
p
,
x
r
β
x
p
,
x
d
Γ
x
.It is assumed that the integrals in (1) and (2) always exist: that is, each integration has a finite value.


## 2.2. Review of Expansions of Geometry Quantities on the Projection Line/Plane [11]

In [11], a projection line for 2D or a projection plane for 3D problems, which is the tangential line/plane of the element to the origin of the intrinsic coordinate system, was introduced. A local orthogonal coordinate system was established on the projection line/plane. The coordinate transformation between the local and global systems can be performed using the following relationships: (3)
x
i
′
=
L
i
j
x
j
-
x
j
o
,
(4)
x
i
=
x
i
o
+
L
j
i
x
j
′
,where L
i
j is the direction cosine of the local coordinate axes with respect to the global one. The repeated subscripts represent summation; x
i
o is the global coordinates of the origin of the local coordinate system.The original curved element can be projected onto the projection line or plane to form a straight or a flat projection element by (3), and then all geometry quantities can be expressed in terms of variables defined on the projection line/plane.


## 2.3. High Order Expansion of Geometry Quantities in terms of Variables on the Projection Line for 2D Problems

For 2D and 3D problems, derivations for the modification are different, which are separately described in detail.In 2D problems, the boundary element is a curved line as shown in Figure1. The local orthogonal coordinate system is denoted by (x
1
′
,
x
2
′) and the projection line is along the axis x
1
′.Figure 1
The curved line boundary and the projection line.There is only one independent variable for a line, and all other geometry quantities can be expressed byx
1
′. Similar to the treatment in [11], the coordinate x
2
′ at the field point over the curved line can be expressed by x
1
′ using Taylor series. Different from the treatment in [11], the third-order derivative is considered in the present work:(5)
x
2
′
=
∂
x
2
′
∂
x
1
′
x
1
′
+
1
2
∂
2
x
2
′
∂
x
1
′
2
x
1
′
2
+
1
6
∂
3
x
2
′
∂
x
1
′
3
x
1
′
3
=
a
1
x
2
′
x
1
′
+
a
2
x
2
′
x
1
′
2
+
a
3
x
2
′
x
1
′
3
,where(6)
a
1
x
2
′
=
∂
x
2
′
∂
x
1
′
=
∂
x
2
′
∂
ξ
1
∂
ξ
1
∂
x
1
′
,
a
2
x
2
′
=
1
2
∂
2
x
2
′
∂
x
1
′
2
=
1
2
∂
2
x
2
′
∂
ξ
2
∂
ξ
∂
x
1
′
∂
ξ
∂
x
1
′
+
∂
x
2
′
∂
ξ
∂
2
ξ
∂
x
1
′
2
,
a
3
x
2
′
=
1
6
∂
3
x
2
′
∂
x
1
′
3
=
1
6
∂
∂
x
1
′
∂
2
x
2
′
∂
x
1
′
2
=
1
6
∂
3
x
2
′
∂
ξ
3
∂
ξ
∂
x
1
′
∂
ξ
∂
x
1
′
∂
ξ
∂
x
1
′
+
3
∂
2
x
2
′
∂
ξ
2
∂
ξ
∂
x
1
′
∂
2
ξ
∂
x
1
′
2
+
∂
x
2
′
∂
ξ
∂
3
ξ
∂
x
1
′
3
.It can be seen that the expressions are more complex. For 2D problem, the corresponding derivatives are analytically derived, and the detailed derivation process is omitted here.In (5) and (6), (7a)
∂
x
2
′
∂
ξ
=
∑
α
=
1
N
node
∂
N
α
∂
ξ
x
2
′
α
,
(7b)
∂
2
x
2
′
∂
ξ
2
=
∑
α
=
1
N
node
∂
2
N
α
∂
ξ
2
x
2
′
α
,
(7c)
∂
3
x
2
′
∂
ξ
3
=
∑
α
=
1
N
node
∂
3
N
n
∂
ξ
3
x
2
′
α
,
(8a)
∂
ξ
∂
x
1
′
=
1
m
1
,
(8b)
∂
2
ξ
∂
x
1
′
2
=
-
1
m
1
4
∂
x
1
∂
ξ
∂
2
x
1
∂
ξ
2
+
∂
x
2
∂
ξ
∂
2
x
2
∂
ξ
2
,
(8c)
∂
3
ξ
∂
x
1
′
3
=
∂
∂
x
1
′
∂
2
ξ
∂
x
1
′
2
=
-
m
1
2
B
-
4
C
m
1
7
.In (7a)–(8c),(9)
m
1
=
∂
x
1
∂
ξ
2
+
∂
x
2
∂
ξ
2
,
(10a)
B
=
∂
2
x
1
∂
ξ
2
2
+
∂
2
x
2
∂
ξ
2
2
+
∂
x
1
∂
ξ
∂
3
x
1
∂
ξ
3
+
∂
x
2
∂
ξ
∂
3
x
2
∂
ξ
3
,
(10b)
C
=
∂
x
1
∂
ξ
2
∂
2
x
1
∂
ξ
2
2
+
2
∂
x
1
∂
ξ
∂
x
2
∂
ξ
∂
2
x
1
∂
ξ
2
∂
2
x
2
∂
ξ
2
+
∂
x
2
∂
ξ
2
∂
2
x
2
∂
ξ
2
2
.Define a local distance:(11)
ρ
=
x
1
′
-
x
1
′
p
,where x
1
′
p is the local coordinate of point P
′ which is the projection point of the source point P and x
1
′ is the local coordinate of point Q
′ which is the projection point of the field point Q.The derivative ofρ with respect to local coordinate x
1
′ is(12a)
ρ
,
1
=
∂
ρ
∂
x
1
′
=
x
1
′
-
x
1
′
p
ρ
.In (12a), ρ
,
1
=
1 or -
1.Then, we can obtain(12b)
x
1
′
=
x
1
′
p
+
ρ
,
1
ρ
.By substituting (12a) into (5), we can obtain(13)
x
2
′
=
p
x
2
′
+
q
x
2
′
,
ρ
ρ
,where(14)
p
x
2
′
=
a
1
x
1
′
p
+
a
2
x
2
′
x
1
′
p
x
1
′
p
+
a
3
x
2
′
x
1
′
p
x
1
′
p
x
1
′
p
=
x
2
′
p
,
(15)
q
x
2
′
,
ρ
=
q
1
x
2
′
+
q
2
x
2
′
ρ
+
q
3
x
2
′
ρ
2
.In (15), (16a)
q
1
x
2
′
=
a
1
x
2
′
+
2
a
2
x
2
′
x
1
′
p
+
3
a
3
x
2
′
x
1
′
p
2
ρ
,
1
,
(16b)
q
2
x
2
′
=
a
2
x
2
′
ρ
,
1
2
+
3
a
3
x
2
′
x
1
′
p
=
a
2
x
2
′
+
3
a
3
x
2
′
x
1
′
p
,
(16c)
q
3
x
2
′
,
ρ
=
a
3
x
2
′
ρ
,
1
.Similarly, we can express the intrinsic coordinate in terms of the local distance as(17)
ξ
=
p
ξ
+
q
ξ
,
ρ
ρ
,where p
(
ξ
) and q
(
ξ
,
ρ
) can be determined by replacing x
2
′ with ξ in (13)–(16c).Substituting (12b) and (13) into (4), we can obtain(18)
x
i
=
x
i
o
+
L
j
i
x
j
′
=
x
i
o
+
L
1
i
x
1
′
+
L
2
i
x
2
′
=
x
i
o
+
L
1
i
x
1
′
p
+
ρ
,
1
ρ
+
L
2
i
x
2
′
p
+
q
x
2
′
,
ρ
ρ
=
x
i
p
+
L
1
i
ρ
,
1
ρ
+
L
2
i
q
x
2
′
,
ρ
ρ
.Then, substituting (15) into (18), we can obtain(19)
x
i
=
x
i
p
+
b
i
+
c
i
ρ
+
d
i
ρ
2
ρ
,where(20a)
b
i
=
L
1
i
ρ
,
1
+
L
2
i
q
1
y
′
,
(20b)
c
i
=
L
2
i
q
2
y
′
,
(20c)
d
i
=
L
2
i
q
3
y
′
.Referring to Figure1, the global distance r can be easily expressed as(21)
r
=
x
i
-
x
i
p
x
i
-
x
i
p
=
g
ρ
ρ
,where(22)
g
ρ
=
b
i
+
c
i
ρ
+
d
i
ρ
2
b
i
+
c
i
ρ
+
d
i
ρ
2
=
G
0
+
G
1
ρ
+
G
2
ρ
2
+
G
3
ρ
3
+
G
4
ρ
4
.In (22), G
0
=
b
i
b
i, G
1
=
2
b
i
c
i, G
2
=
c
i
c
i
+
2
b
i
d
i, G
3
=
2
c
i
d
i, and G
4
=
d
i
d
i.It can be easily derived that(23)
r
,
i
=
∂
r
∂
x
i
=
x
i
-
x
i
p
r
=
b
i
+
c
i
ρ
+
d
i
ρ
2
g
ρ
.


## 2.4. The Complex-Variable-Differentiation Method (CVDM)

From the above derivations, it can be seen that the corresponding expressions have been changed, by considering the third-order derivatives when expanding the geometry quantities at the field point as Taylor series. The third-order derivatives should be first determined.In the present work, the complex-variable-differentiation method (CVDM) [20] is introduced to determine the high order derivatives. For 2D analysis, the third-order derivatives are determined explicitly by (7c), (8c), (9), (10a), and (10b), which are also calculated by the CVDM.In the CVDM, the variablex of a real function f
(
x
) is replaced by a complex one x
+
i
h, with the imaginary part h being very small (usually h = 10−20). The function f
(
x
+
i
h
) can be expanded into Taylor’s series as(24)
f
x
+
i
h
=
f
x
+
i
h
f
′
x
-
h
2
2
f
′
′
x
+
o
h
3
.Since h is small enough, (25) can be obtained: (25)
f
′
x
=
Im
⁡
f
x
+
i
h
h
,where Im denotes the imaginary part.The CVDM is a very promising method, since the derivatives only require direct evaluations [21]. The CVDM has been validated to be an efficient and accurate approach for determination of derivatives [22–24].For example, considering∂
3
x
2
′
/
∂
x
1
′
3
=
∂
/
∂
x
1
′
∂
2
x
2
′
/
∂
x
1
′
2, we can add the very small imaginary part h to x
1
′; then the real value of ∂
3
x
2
′
/
∂
x
1
′
3 yields Im
⁡
(
∂
2
x
2
′
/
∂
x
1
′
2
)
/
h. In addition, the value of ∂
2
x
2
′
/
∂
x
1
′
2 can be easily obtained using the analytical expressions in the previous work [11]. The CVDM and the derived expressions in (7c)–(8c) can be validated by each other.


## 2.5. Expansion of Geometry Quantities in terms of Variables on the Projection Plane for 3D Problems

In 3D problems, the boundary element is a curved surface marked asA
B
C
D in Figure 2. A local orthogonal coordinate system (x
1
′
,
x
2
′
,
x
3
′) is established on the projection plane with its origin being at the point (ξ
1
=
0, ξ
2
=
0), in which axes x
1
′ and x
2
′ are located within the plane. The axis x
1
′ is along ξ
1 direction and axis x
3
′ is along the outward normal direction to the element. By (4), we can project the original curved surface element onto the projection plane, to form a flat projection element marked as A
′
B
′
C
′
D
′ in Figure 2.Figure 2
The curved surface boundary and the projection plane.In a plane, there are only two independent variables. We choosex
1
′ and x
2
′ as the independent variables, and x
3
′ over the curved surface can be expressed in terms of x
1
′ and x
2
′. Compared with the treatment in [11], x
3
′ is expanded as Taylor series by x
1
′ and x
2
′, by considering the third-order derivatives:(26)
x
3
′
=
a
I
x
3
′
x
I
′
+
a
I
J
x
3
′
x
I
′
x
J
′
+
a
I
J
K
x
3
′
x
I
′
x
J
′
x
K
′
,where I, J, and K are integers varying from 1 to 2, and(27a)
a
I
x
3
′
=
∂
x
3
′
∂
x
I
′
=
∂
x
3
′
∂
ξ
K
∂
ξ
K
∂
x
I
′
,
(27b)
a
I
J
x
3
′
=
1
2
∂
2
x
3
′
∂
x
I
′
∂
x
J
′
=
1
2
∂
2
x
3
′
∂
ξ
K
∂
ξ
L
∂
ξ
K
∂
x
I
′
∂
ξ
L
∂
x
J
′
+
∂
x
3
′
∂
ξ
K
∂
2
ξ
K
∂
x
I
′
∂
x
J
′
,
(27c)
a
I
J
K
x
3
′
=
1
6
∂
3
x
3
′
∂
x
I
′
∂
x
J
′
∂
x
K
′
=
1
3
∂
a
I
J
x
3
′
∂
x
K
′
=
1
3
∂
a
I
J
x
3
′
∂
ξ
L
∂
ξ
L
∂
x
K
′
.The treatment is the same as that in [11], to determine a
I
(
x
3
′
) in (27a) and a
I
J
(
x
3
′
) in (27b). For a
I
J
K
(
x
3
′
) in (27c), we use the CVDM, to avoid rather complicated analytical derivations.The local distanceρ projected from the global distance r onto the projection plane and its derivative ρ
,
I are introduced as follows: (28)
ρ
=
x
1
′
-
x
1
′
p
2
+
x
2
′
-
x
2
′
p
2
,
ρ
,
I
=
∂
ρ
∂
x
I
′
=
x
I
′
-
x
I
′
p
ρ
.The local coordinates at the field point can be expressed as(29)
x
I
′
=
x
I
′
p
+
ρ
,
I
ρ
.Similar to the 2D analysis, substituting (29) into (26) yields(30)
x
3
′
=
p
x
3
′
+
q
x
3
′
,
ρ
ρ
,where(31)
p
x
3
′
=
a
I
x
3
′
x
I
′
p
+
a
I
J
x
3
′
x
I
′
p
x
J
′
p
+
a
I
J
K
x
3
′
x
I
′
p
x
J
′
p
x
K
′
p
=
x
3
′
p
,
(32)
q
x
3
′
,
ρ
=
q
1
x
3
′
+
q
2
x
3
′
ρ
+
q
3
x
3
′
ρ
2
.In (32),(33a)
q
1
x
3
′
=
a
I
x
3
′
ρ
,
I
+
a
I
J
x
3
′
x
I
′
p
ρ
,
J
+
a
I
J
x
3
′
x
J
′
p
ρ
,
I
+
a
I
J
K
x
3
′
x
I
′
p
x
K
′
p
ρ
,
J
+
a
I
J
K
x
3
′
x
J
′
p
x
K
′
p
ρ
,
I
+
a
I
J
K
x
3
′
x
I
′
p
x
J
′
p
ρ
,
K
,
(33b)
q
2
x
3
′
=
a
I
J
x
3
′
ρ
,
I
ρ
,
J
+
a
I
J
K
x
3
′
x
K
′
p
ρ
,
I
ρ
,
J
+
a
I
J
K
x
3
′
x
I
′
p
ρ
,
J
ρ
,
K
+
a
I
J
K
x
3
′
x
J
′
p
ρ
,
I
ρ
,
K
,
(33c)
q
3
x
3
′
=
a
I
J
K
x
3
′
ρ
,
I
ρ
,
J
ρ
,
K
.Similar to the 2D analysis, intrinsic coordinates can be expressed as(34)
ξ
K
=
p
ξ
K
+
q
ξ
K
,
ρ
ρ
,where p
(
ξ
K
) and q
(
ξ
K
,
ρ
) can be determined by replacing x
3
′ with ξ
K in (31)–(33c).Substituting (29) and (30) into (4) yields(35)
x
i
=
x
i
o
+
L
J
i
x
J
′
+
L
3
i
x
3
′
=
x
i
o
+
L
J
i
x
J
′
p
+
ρ
,
J
ρ
+
L
3
i
x
3
′
p
+
q
x
3
′
,
ρ
ρ
=
x
i
p
+
L
J
i
ρ
,
J
ρ
+
L
3
i
q
x
3
′
,
ρ
ρ
.Substituting (32) into (35), we can obtain (36)
x
i
=
x
i
p
+
b
i
+
c
i
ρ
+
d
i
ρ
2
ρ
,where(37a)
b
i
=
L
J
i
ρ
,
J
+
L
3
i
q
1
x
3
′
,
(37b)
c
i
=
L
3
i
q
2
x
3
′
,
(37c)
d
i
=
L
3
i
q
3
x
3
′
.The expression ofr is as follows:(38)
r
=
x
i
-
x
i
p
x
i
-
x
i
p
=
g
ρ
ρ
,where(39)
g
ρ
=
b
i
+
c
i
ρ
+
d
i
ρ
2
b
i
+
c
i
ρ
+
d
i
ρ
2
=
G
0
+
G
1
ρ
+
G
2
ρ
2
+
G
3
ρ
3
+
G
4
ρ
4
.In (39), G
0
=
b
i
b
i, G
1
=
2
b
i
c
i, G
2
=
c
i
c
i
+
2
b
i
d
i, G
3
=
2
c
i
d
i, and G
4
=
d
i
d
i.It can be easily derived thatr
,
i has the same expression as that in 2D analysis.Considering thatI
=
J
=
K
=
1 and ρ
,
1
=
1 or -
1 for 2D problems, it can be seen that unified expressions are obtained for both two-dimensional and three-dimensional problems.Then, we can handle the singular boundary integrals in (1) and (2), with the above modifications. The 2D problem is taken as an example to describe the procedure to determine the integral in (1). The differential relationship between the local distance and the real curved line can be written as (40)
d
ρ
=
d
Γ
cos
⁡
ϕ
=
d
Γ
n
i
0
n
i
,where n
i
0 and n
i are the outward normals to the tangential lines of passing the origin of the local coordinate system and the field point, respectively, and ϕ is the angle between the two outward normals.Substituting (40) and (38) into (1) leads to(41)
I
e
x
p
=
lim
ε
→
0
⁡
∫
ε
ρ
E
F
-
ρ
ρ
β
d
ρ
,where ρ
E is the distance from the projected source point P
′ to one of the ends of the projected line element and F
-
(
ρ
) is the regular part: that is, (42)
F
-
ρ
=
f
-
x
p
,
x
ρ
g
β
ρ
n
i
0
n
i
ρ
.In order to integrate (41), the nonsingular part F
- is expanded as a power series in ρ, such that(43)
F
-
ρ
=
∑
n
=
0
N
B
n
ρ
nin which N is the order of the power series, usually taking a value between 2 and 7 depending on the size of ρ
E, and B
(
n
) are constants which are determined by collocating N
+
1 points over the integration region (
0
,
ρ
E
). In the paper, N
+
1 equally spaced points are used: that is, (
0
,
ρ
1
,
…
,
ρ
N
). The coefficient for the first point n
=
0 is B
(
0
)
=
F
-
(
0
) and other coefficients can be solved using the following equation set:(44)
R
B
=
Y
,where R is a square matrix with the order of N:(45)
R
=
1
ρ
1
⋯
ρ
1
N
-
1
1
ρ
2
⋯
ρ
2
N
-
1
⋮
⋮
⋯
⋮
1
ρ
N
⋯
ρ
N
N
-
1
.
B and Y are vectors as follows:(46)
B
=
B
1
B
2
⋮
B
N
,
Y
=
F
-
ρ
1
-
B
0
ρ
1
F
-
ρ
2
-
B
0
ρ
2
⋮
F
-
ρ
N
-
B
0
ρ
N
.Solving (44) for coefficient vector B and then substituting (43) into (41) yield (47)
I
e
x
p
=
∑
n
=
0
N
B
n
lim
ε
→
0
⁡
∫
ε
ρ
E
ρ
n
-
β
d
ρ
=
∑
n
=
0
N
B
n
E
n
,where (48)
E
n
=
1
n
-
β
+
1
1
ρ
E
β
-
n
-
1
-
lim
ρ
ε
→
0
⁡
1
ρ
ε
β
-
n
-
1
n
≠
β
-
1
ln
⁡
ρ
E
-
lim
ρ
ε
→
0
⁡
ln
⁡
ρ
ε
n
=
β
-
1
.Then, singularities can be analytically removed.


## 3. Numerical Examples

Based on the Fortran code for singular integrals in [11], a new Fortran code has been developed with the above modification presented in this paper. It should be emphasized that the functions of the old Fortran code are fully included in the new Fortran code. If all the third-order derivatives are set to be zero, the functions are the same with those of the method in [11]. Besides, the new code has its unique characteristics.For 2D analysis, the third-order derivatives are computed by both CVDM and analytical expression. For 3D analysis, the high order derivatives are computed by the CVDM, because analytical derivations are very difficult and complicated to be obtained.In this section, the first example is given to show the advantage of CVDM over a difference method. Then, two singular integrals are calculated, to show the effectiveness and the accuracy of the present method.


### 3.1. The Advantage of the CVDM

A singular function is selected:(49)
f
x
,
x
p
=
1
x
-
x
p
2
.Then, the exact derivative with respect to the variablex
p can be obtained: (50)
∂
f
∂
x
p
=
2
x
-
x
p
3
.For comparison, the central difference formulation is adopted: that is,(51)
∂
f
∂
x
p
=
f
x
,
x
p
+
h
-
f
x
,
x
p
-
h
2
h
.Whenx
=
1.001, x
p
=
1, the exact first-order derivative is 2 × 109. The dependence of the derivatives on h for the CVDM and the difference method is examined, and the results are listed in Table 1.Table 1
The first-order derivatives with differenth for f
(
x
,
x
p
)
=
1
/
(
x
-
x
p
)
2.
h
CVDM
Difference method
1
.
00
E
-
04
1.9603210
E
+
09
2.0403182
E
+
09
1.00
E
-
06
1.9997157
E
+
09
1.9997236
E
+
09
1.00
E
-
08
1.9997197
E
+
09
1.9997197
E
+
09
1.00
E
-
10
1.9997197
E
+
09
1.9997194
E
+
09
1.00
E
-
12
1.9997197
E
+
09
1.9997857
E
+
09
1.00
E
-
14
1.9997197
E
+
09
1.9981112
E
+
09
1.00
E
-
16
1.9997197
E
+
09
1.1100197
E
+
09
1.00
E
-
17
1.9997197
E
+
09
0.0000000
E
+
00
1.00
E
-
28
1.9997197
E
+
09
0.0000000
E
+
00It can be seen that the difference method can give satisfactory results only forh between 10−5 and 10−10. For smaller h, its result becomes unacceptable. However, the result of the CVDM tends to be independent of h as the step size becomes sufficiently small. This clearly verifies the advantage of the CVDM.


### 3.2. A Weakly Singular Line Integral over a Cubed Curve

In this example, a singular line integral for a 2D BEM problem is chosen:(52)
J
=
∫
Γ
ln
⁡
r
d
Γ
.The integration line is a cubed curve as shown in Figure3.Figure 3
The cubed curve for integration.The four typical points are 1(
-
1
,
-
1
), 2 (
1,1
), 3 (
-
1
/
3
,
-
1
/
27
), and 4 (
1
/
3,1
/
27
). The computational point (source point) is at 1. The curve is taken as one boundary element. It should be emphasized here that the direct method in [11] has been invalid, which cannot be used for the singular line integration in (52) over the cubed curve in Figure 4, with only one element, because only quadratic terms were truncated when expanding the geometry quantities at the field point as Taylor series.Figure 4
The cylindrical surface.Firstly, the accuracy of the CVDM for the proposed technique is validated. We use the derived analytical expressions in (7c)–(8c) to compute the third-order derivatives. Then, the CVDM is employed to determine the same derivatives.Using the derived analytical expressions in this paper, the third-order derivatives can be obtained:(53)
∂
3
x
2
′
∂
x
1
′
3
=
-
6.0000
,
∂
3
ξ
∂
x
1
′
3
=
-
6.7493
×
1
0
-
6
.The results of the CVDM are exactly the same with that in (53). It can be seen that the CVDM and the analytical expressions in (7c)–(8c) are validated by each other.Secondly, the singular integral is evaluated using the present method. The result is 0.229935, with the Gauss order being 6.This integral is also evaluated using the Gauss quadrature rule for logarithmic function [12]. The result is J
=
0.226920. It can be seen that the present method is effective and with high accuracy, for higher order curved boundary.


### 3.3. A Supersingular Integral over a Curved Surface

In this example, the accuracy of the present method in handling singular integrals over a highly curved surface is examined. The boundary element is a cylindrical surface with radius 1 and the rotational angle is9
0
°. The geometrical dimensions are shown in Figure 4.The boundary integral is as follows:(54)
I
x
p
=
∫
Γ
-
1
4
π
r
β
3
r
,
3
∂
r
∂
n
-
n
3
d
Γ
.Two source points,a and b, are located at positions with intrinsic coordinates of (
0,0
) and (
0.66,0
).Guiggiani et al. [14] computed this integral for the case of β
=
3. Table 2 lists the results using both methods. For the present method, the Gauss order is 4.Table 2
Results for two pointsa and b(β
=
3).
Pointa
Pointb
Method in [14]
−0.3438
−0.4971
Present method
−0.3458
−0.4981From Table2, it can be seen that he results using the present method agree well with the Guiggiani et al.’s results, which indicate the high accuracy of the present method for singular integrals over a highly curved surface.


## 3.1. The Advantage of the CVDM

A singular function is selected:(49)
f
x
,
x
p
=
1
x
-
x
p
2
.Then, the exact derivative with respect to the variablex
p can be obtained: (50)
∂
f
∂
x
p
=
2
x
-
x
p
3
.For comparison, the central difference formulation is adopted: that is,(51)
∂
f
∂
x
p
=
f
x
,
x
p
+
h
-
f
x
,
x
p
-
h
2
h
.Whenx
=
1.001, x
p
=
1, the exact first-order derivative is 2 × 109. The dependence of the derivatives on h for the CVDM and the difference method is examined, and the results are listed in Table 1.Table 1
The first-order derivatives with differenth for f
(
x
,
x
p
)
=
1
/
(
x
-
x
p
)
2.
h
CVDM
Difference method
1
.
00
E
-
04
1.9603210
E
+
09
2.0403182
E
+
09
1.00
E
-
06
1.9997157
E
+
09
1.9997236
E
+
09
1.00
E
-
08
1.9997197
E
+
09
1.9997197
E
+
09
1.00
E
-
10
1.9997197
E
+
09
1.9997194
E
+
09
1.00
E
-
12
1.9997197
E
+
09
1.9997857
E
+
09
1.00
E
-
14
1.9997197
E
+
09
1.9981112
E
+
09
1.00
E
-
16
1.9997197
E
+
09
1.1100197
E
+
09
1.00
E
-
17
1.9997197
E
+
09
0.0000000
E
+
00
1.00
E
-
28
1.9997197
E
+
09
0.0000000
E
+
00It can be seen that the difference method can give satisfactory results only forh between 10−5 and 10−10. For smaller h, its result becomes unacceptable. However, the result of the CVDM tends to be independent of h as the step size becomes sufficiently small. This clearly verifies the advantage of the CVDM.


## 3.2. A Weakly Singular Line Integral over a Cubed Curve

In this example, a singular line integral for a 2D BEM problem is chosen:(52)
J
=
∫
Γ
ln
⁡
r
d
Γ
.The integration line is a cubed curve as shown in Figure3.Figure 3
The cubed curve for integration.The four typical points are 1(
-
1
,
-
1
), 2 (
1,1
), 3 (
-
1
/
3
,
-
1
/
27
), and 4 (
1
/
3,1
/
27
). The computational point (source point) is at 1. The curve is taken as one boundary element. It should be emphasized here that the direct method in [11] has been invalid, which cannot be used for the singular line integration in (52) over the cubed curve in Figure 4, with only one element, because only quadratic terms were truncated when expanding the geometry quantities at the field point as Taylor series.Figure 4
The cylindrical surface.Firstly, the accuracy of the CVDM for the proposed technique is validated. We use the derived analytical expressions in (7c)–(8c) to compute the third-order derivatives. Then, the CVDM is employed to determine the same derivatives.Using the derived analytical expressions in this paper, the third-order derivatives can be obtained:(53)
∂
3
x
2
′
∂
x
1
′
3
=
-
6.0000
,
∂
3
ξ
∂
x
1
′
3
=
-
6.7493
×
1
0
-
6
.The results of the CVDM are exactly the same with that in (53). It can be seen that the CVDM and the analytical expressions in (7c)–(8c) are validated by each other.Secondly, the singular integral is evaluated using the present method. The result is 0.229935, with the Gauss order being 6.This integral is also evaluated using the Gauss quadrature rule for logarithmic function [12]. The result is J
=
0.226920. It can be seen that the present method is effective and with high accuracy, for higher order curved boundary.


## 3.3. A Supersingular Integral over a Curved Surface

In this example, the accuracy of the present method in handling singular integrals over a highly curved surface is examined. The boundary element is a cylindrical surface with radius 1 and the rotational angle is9
0
°. The geometrical dimensions are shown in Figure 4.The boundary integral is as follows:(54)
I
x
p
=
∫
Γ
-
1
4
π
r
β
3
r
,
3
∂
r
∂
n
-
n
3
d
Γ
.Two source points,a and b, are located at positions with intrinsic coordinates of (
0,0
) and (
0.66,0
).Guiggiani et al. [14] computed this integral for the case of β
=
3. Table 2 lists the results using both methods. For the present method, the Gauss order is 4.Table 2
Results for two pointsa and b(β
=
3).
Pointa
Pointb
Method in [14]
−0.3438
−0.4971
Present method
−0.3458
−0.4981From Table2, it can be seen that he results using the present method agree well with the Guiggiani et al.’s results, which indicate the high accuracy of the present method for singular integrals over a highly curved surface.


## 4. Conclusion and Future Work

In the present work, high order projection plane method for evaluation of supersingular curved boundary integrals in BEM has been developed, by considering the third-order derivatives when expanding the geometry quantities at the field point as Taylor series. Corresponding analytical expressions are derived, and unified expressions are obtained for both two-dimensional and three-dimensional problems.The present method extends the application range of the method in [11]. That is, the singular integrals over higher order elements, such as cubed curve, can be directly evaluated. The complex-variable-differentiation method is introduced to determine the third-order derivatives, which effectively avoids rather complicated analytical derivations for the three-dimensional analysis. For the two-dimensional boundary integral, the third-order derivatives are analytically derived and calculated, which are used to verify the complex-variable-differentiation method.The results show that the present method is effective and with high accuracy for high order singular curved boundary integrals.Future work is to apply the mathematical derivations and equations to evaluate the singularities involved in radial integration boundary element method in solving nonhomogeneous thermoelastic problems.



---

*Source: 9523405-2016-03-08.xml*