- Email: [email protected]
A new method for numerical evaluation of nearly singular integrals over high-order geometry elements in 3D BEM
Accepted Manuscript A new method for numerical evaluation of nearly singular integrals over high-order geometry elements in 3D BEM Yaoming Zhang, Xiaochao Li, Vladimir Sladek, Jan Sladek, Xiaowei Gao PII: DOI: Reference:
S0377-0427(14)00384-7 http://dx.doi.org/10.1016/j.cam.2014.08.027 CAM 9781
To appear in:
Journal of Computational and Applied Mathematics
Received date: 11 October 2013 Revised date: 20 May 2014 Please cite this article as: Y. Zhang, X. Li, V. Sladek, J. Sladek, X. Gao, A new method for numerical evaluation of nearly singular integrals over high-order geometry elements in 3D BEM, Journal of Computational and Applied Mathematics (2014), http://dx.doi.org/10.1016/j.cam.2014.08.027 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A new method for numerical evaluation of nearly singular integrals over high-order geometry elements in 3D BEM Yaoming Zhang a,b,*, Xiaochao Li a , Vladimir Sladek c, Jan Sladek c, Xiaowei Gao b a
c
Institute of Applied Mathematics, Shandong University of Technology,Zibo255049,PR China b State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China Slovak Academy of Sciences, Institute of Construction and Architecture, 84503 Bratislava, Slovakia
Abstract: This work presents a new method for numerical computation of the two-dimensional nearly singular integrals with using the eight-node second-order quadrilateral surface elements in 3D BEM. A new indirect regularized boundary element formulation excluding the CPV (Cauchy Principal Value) and HFP (Hadamard-Finite-Part) integrals is proposed. Based on this, a new approximation formula of the distance from the fixed calculation point to a generic point of the aforementioned surface geometry elements is developed firstly, and then the exponential transformation, which has been widely employed in 2D BEM, is extended to 3D BEM to remove the near singularities of integrands for considered integrals. Several numerical examples are given to verify the high efficiency and the stability of the proposed scheme. Key Words: BEM; Nearly singular integrals; High-order geometry elements; Transformation; 3D potential problems.
1 Introduction The BEM as an important numerical technique has been widely used in many areas, but the accuracy of its calculation depends heavily on the accuracy and efficiency of the computation of integrals with singular and nearly singular kernels. These kernels become singular when the collocation point belongs to the integration elements, and many effective methods have been developed to handle them[1-11]. Another important issue is the integration of the kernels for the collocation points which are close to but not on the integration element. Then, the considered integrals, although regular in the sense of mathematics, are termed nearly singular integrals since they exhibit similar characteristics of the singularity and generally can not be calculated accurately by standard integration quadratures. This is so-called boundary layer effect in BEM. In this work, we focus on numerical computation of the two-dimensional nearly singular integrals that arise in the *
Corresponding author at :Institute of Applied Mathematics, Shandong University of Technology, Zibo255049, PR China. . E-mail address: [email protected] (Y. M. Zhang).
1
solution of three-dimensional(3D) boundary element method(BEM) using eight-node second-order quadrilateral surface elements. Such problems would be more crucial to some of engineering applications, such as crack problems when the crack tip is deformed to have a small opening displacement[12], contact problems when the contact area of the two contacting bodies is very small[13], inverse and sensitivity problems[14], as well as thin or shell-like structure problems[15-17], etc. In these situations, the number of the nearly singular integrals, which have to be evaluated for a system matrix, can be much larger than that of the singular integrals. Effective evaluation of nearly singular integrals has attracted a great amount of attention in recent years. Various numerical techniques have been developed to remove the near singularities, such as rigid-body displacement and other global regularization methods[17-23], analytical or semi-analytical formulas[24-27], adaptive element subdivision techniques [28-29], the use of polar coordinate transformations[30-32], and various nonlinear transformations[33-43]. Among these methods, the variable transformations technique seems be a promising method due to its wide suitability and high accuracy and stability. The methods developed so far include, but are not limited to, polynomial transformation[33], sigmoidal transformation[34], coordinate optimal transformation[35], degenerate mapping method[13,36],
sinh
transformation[37-38],
rational
transformation[39],
distance
transformation[40-43], and exponential transformation[44-48]. A comprehensive review was given by present author[44]. Although great progresses have been achieved for each of the above methods, a number of drawbacks are associated with these methods and these include the fact that some, which benefit from the strategies for calculating the singular integrals, are failed to fully eliminate near singularities, some apply only to the planar elements, some essentially use the first-order functions to depict the surface geometry boundary although the authors claim that the problems are considered over the curved surface elements, and others are tailored specifically to the form of the kernel functions and therefore lack wide applicability. In previous work, the present author proposed a nonlinear transformation, named as the exponential transformation, to remove or damp out the near singularities of integrands for considered integrals in 2D BEM[44]. This transformation, based on the idea of diminishing the difference of the orders of magnitude or the scale of change of addition factors in the denominator of the kernels, is more simple and applicable than the previous works of the above researches according to extensive numerical experiments and applications. It is also demonstrated that the exponential transformation offers two primary advantages over most of existing methodologies. One is the accuracy and ease of implementation, the other one is applicable to a broad range of integrals with various kernels without extra computational 2
effort. In present work, we would extend such transformation to 3D BEM. The usage of high-order geometry elements, most usually of the eight-node second-order quadrilateral surface elements, for efficient evaluation of the nearly singular integrals is necessary in many cases. Clearly, the advantage of using the high-order geometry elements in complex geometrical domains does not concern only its power to achieve higher calculation accuracy, but even with a small number of computational elements. More crucial, computational models of thin structure problems or the case when the calculation points are very close to the actual boundary demand a higher level of geometry approximation, and the usage of high-order geometry in computational model can meet this requirement. For example, if the boundary geometry is depicted by using planar elements, the planar elements of the outer surface will contact or even pass through the inner boundary if the thickness of the considered structure is very small or the case where some interior points very close to the actual boundary may be located outside the boundary elements, and obviously this calculation model is fallacious. Consequently, the actual geometry of considered domain cannot be described faithfully, and thus lower-order geometry approximation will fail to yield reliable results for such problems. In order to avoid this phenomenon, very fine meshes must be used in this situation, but this would yields too much preprocessing and CPU time. In addition, a great number of meshes will produce a lot of artificial corners that will lead to the discontinuity of the tangent derivative of the boundary unknowns. This is fatal to many engineering problems. In this paper, a general methodology is proposed to compute the nearly singular integrals that arise in the solution of 3D BEM using the eight-node second-order quadrilateral surface elements. A new indirect regularized boundary element formulation excluding the CPV and HFP integrals is presented. Based on this, a new approximation formula of the distance between the fixed calculation point and a generic point of curved surface elements is developed here firstly, and then the exponential transformation, which has been widely applied in 2D BEM, is developed to remove the near singularities of the considered integrals in 3D BEM before applying the standard Gaussian quadrature to numerical integration. Four numerical examples with exact benchmark solutions are presented to test the proposed scheme, yielding very promising results even when the internal point is very close to the boundary. The outline of the rest of this paper is as follows. The indirect regularized BIEs excluding the CPV and HFP integrals are presented in Section2. Then, in Section 3, the expression of the nearly singular integrals over the eight-node second-order quadrilateral surface elements is investigated, including the determination of the projection point and construction of the distance functions. Section 4 considers the regularization of nearly 3
singular integrals over high-order geometrical elements using the exponential transformation. In Section 5, the accuracy and stability of the proposed scheme are tested in three 3D potential examples with known benchmark solutions. Finally, the conclusions are provided in Section 6.
2 Regularized boundary integral equations(RBIEs) It is well known that having known all relevant boundary quantities, the solutions at interior points can be obtained by using the integral representation. Thus, the accuracy of boundary quantities directly affects the validity of the domain quantities, but not opposite. This is a great advantage of pure boundary integral formulations as compared with the domain discretization methods such as FEM, for instance. On the other hand, we have to deal with the singular boundary integrals when calculating the boundary unknowns. Therefore, a good choice is to use the regularized BIEs with indirect unknowns. In this paper, we always assume that Ω is a bounded domain in R 3 , Ω c is its open complement, and Γ denotes the common boundary. n( x ) is the unit outward normal vector on Γ to the domain Ω at the point x . With omitting the body sources in potential ˆ can be expressed as problems, the regularized BIEs with indirect unknowns on Ω u ( y ) = ∫ φ ( x )u * ( x, y )d Γ, y ∈ Γ Γ
∂u ( y ) ˆ ∂u * ( x , y ) = Sφ ( y ) + ∫ [φ ( x ) − φ ( y ) ] dΓ Γ ∂nˆ ( y ) ∂nˆ ( y ) ⎡ ∂u * ( x , y ) ∂u * ( x , y ) ⎤ +φ ( y ) ∫ ⎢ + ⎥ d Γ, y ∈ Γ Γ ∂nˆ ( x ) ⎦ ⎣ ∂nˆ ( y )
(1)
(2)
For the internal point y , the integral equations can be written as ˆ u ( y ) = ∫ φ ( x )u * ( x , y ) d Γ x , y ∈ Ω Γ
ˆ ∇ y u ( y ) = ∫ φ ( x )∇ y u * ( x , y )d Γ x , y ∈ Ω Γ
(3) (4)
In the Eqs.(1)-(4), φ ( x ) is the density function to be determined; u * ( x , y ) denotes the ˆ = Ω , Sˆ = 1 , nˆ ( x ) is the unit outward Kelvin fundamental solution. For interior problems, Ω
ˆ = Ωc , Sˆ = 0 , normal vectors on Γ to domain Ω at point x . For exterior problems, Ω nˆ ( x ) is the unit outward normal vectors on Γ to domain Ωc at point x .
For the discretized form of the Eqs. (3) and (4), when the field point y is far enough from the integration elements, a straightforward application of Gaussian quadrature procedure suffices to evaluate such integrals. However, when the field point y is very close to the 4
integration elements Γ e , the distance r between the field point y and the source point x is almost zero. Hence, the integrals in discretized Eqs. (3)-(4) are nearly singular and the numerical integrations by the standard Gaussian quadrature fail. These nearly singular integrals can be expressed as Ι=∫
Γe
f ( x, y ) dΓ rα
(5)
where r = x − y 2 , α > 0 is a real constant, and f ( x, y ) denotes a well-behaved function.
3 Nearly singular integrals on high-order geometry elements The quintessence of the BEM is to discretize the boundary into a finite number of segments, not necessarily uniform, which are called boundary elements. Two approximations are made over each of these elements. One is about the geometry of the boundary, while the other has to interpolate the spatial variation of the field variables or the boundary unknowns over the element. The plane shape geometry element is not an ideal one as it can not approximate the curved surface boundaries with sufficient accuracy. For this reason, it is recommended to use higher order elements, which approximate geometry and boundary quantities by using higher order interpolation polynomials-usually of the second order. In this paper, the geometry segment is modeled by a continuous paraboloidal element, which has eight knots, namely, the boundary geometry is approximated by the piecewise continuous eight-node second-order quadrilateral surface elements, while the distribution of the boundary quantities over each of these segments is approximated using discontinuous elements, eight nodes of which are located away from the edges of the element. Assume x j = ( x1j , x2j , x3 j ), j = 1,L,8 are the eight knots of the segment Γ j , then
Cartesian coordinates of the points on the element Γ j can be interpolated as[50] 8
xk (ξ1 , ξ 2 ) = ∑ N j (ξ1 , ξ 2 ) xkj , k = 1, 2,3
(6)
1 1 N1 (ξ1 , ξ 2 ) = (1 − ξ1 )(1 − ξ 2 )(−ξ1 − ξ 2 − 1), N 2 (ξ1 , ξ 2 ) = (1 + ξ1 )(1 − ξ 2 )(ξ1 − ξ 2 − 1)， 4 4 1 1 N 3 (ξ1 , ξ 2 ) = (1 + ξ1 )(1 + ξ 2 )(ξ1 + ξ 2 − 1)，N 4 (ξ1 , ξ 2 ) = (1 − ξ1 )(1 + ξ 2 )(−ξ1 + ξ 2 − 1)， 4 4 1 1 N 5 (ξ1 , ξ 2 ) = (1 − ξ12 )(1 − ξ 2 )，N 6 (ξ1 , ξ 2 ) = (1 + ξ1 )(1 − ξ 2 2 ), 2 2 1 1 N 7 (ξ1 , ξ 2 ) = (1 + ξ 2 )(1 − ξ12 ), N8 (ξ1 , ξ 2 ) = (1 − ξ1 )(1 − ξ 2 2 ), − 1 ≤ ξ1 ≤ 1, − 1 ≤ ξ 2 ≤ 1 2 2
(7)
j =1
where
5
where ξ1 , ξ 2 are the local coordinates. 3.1
Determination of the projection point
The minimum distance d from the field point y to the integration element Γ e is defined as the length y − x p , where x p is the projection point of y onto integration element Γ e . Letting (η1 , η 2 ) be the local coordinates of the projection point x p , i.e. x p = ( x1 (η1 ,η2 ), x2 (η1 ,η2 ), x3 (η1 ,η2 )) , then η1 , η2 are the real roots of the following equation ∂xi ⎧ ⎪[ xi (η1 ,η2 ) − yi ] ∂ξ = 0 ⎪ 1 , ⎨ ⎪[ xi (η1 ,η2 ) − yi ] ∂xi = 0 ⎪⎩ ∂ξ 2 in which the summation convention is used, and
i = 1, 2,3
∂xi ∂xi = ∂ξ k ∂ξ k
ξ1 =η1 ξ 2 =η2
(8)
, k = 1, 2 . These assumptions
will be applied also in what follows unless specified otherwise. If the field point y is sufficiently close to the boundary Γ , then x p is inside the integration element, and Eq. (8) has a pair of the unique real roots (η1 , η 2 ) ∈ [−1,1] × [−1,1] . The real roots η1 , η2 can be evaluated numerically by using the Newton’s method. Setting
f1 (η1 ,η2 ) = [ xi (η1 ,η2 ) − yi ]
∂xi ∂x , f 2 (η1 ,η2 ) = [ xi (η1 ,η2 ) − yi ] i ∂ξ1 ∂ξ 2
The formula of the Newton’method can be expressed as F ′(η ( k ) )Δη ( k ) = − F (η ( k ) )
(9)
where Δη ( k ) = η ( k +1) − η ( k ) ,η ( k ) = (η1( k ) ,η2 ( k ) )T , η ( k +1) = (η1( k +1) ,η2( k +1) )T , ⎡ ⎢ ⎡ f1 (η ,η2 ) ⎤ (k ) (k ) ′ , F (η ) = ⎢ F (η ) = ⎢ (k ) (k ) ⎥ ⎢ ⎢⎣ f 2 (η1 ,η 2 ) ⎥⎦ ⎢ ⎣ (k ) 1
(k )
∂f1 ∂η1 ∂f 2 ∂η1
∂f1 ⎤ ∂η 2 ⎥ ⎥ ∂f 2 ⎥ ∂η 2 ⎥⎦
η =η ( k )
here 3 ⎛ ∂x ∂xi ∂ 2 xi ⎞ = ∑⎜ i + [ xi (η1 ,η 2 ) − yi ] ⎟, j , m = 1, 2 . ∂ηm i =1 ⎜⎝ ∂ξ m ∂ξ j ∂ξ j ∂ξ m ⎟⎠
∂f j
It should be emphasized that the effective and accurate computation of the projection 6
point x p , whose local coordinates in the parameter plane are (η1 ,η 2 ) , is crucial to implementation of the proposed scheme. In order to demonstrate the validity and convergence of the above Newton’s iterative scheme，the following numerical example 2 is investigated here. The discretization of boundary surface is as in example 2, and we randomly take an eight-node quadrilateral surface element for which the coordinates of eight nodes are given in the Table 1. Ten interior points Ai , i = 1,L ,10 are considered close to this element, as shown in Fig. 1. Pi , i = 1,L ,10 denote the projection points of Ai , i = 1,L ,10 onto this element, respectively, which have been computed by aforementioned Newton’s iterative scheme. The same one initial approximation η0 = (0, 0) is chosen for all iterative procedures, and the iterative results, including projection points coordinates and number of iterations, are listed in Table 2 based on the stopping criterion η ( k +1) − η ( k ) 2 ≤ 1 × 10−11 . From these results, we can conclude that it doesn’t matter if the field point y is very close to integration elements, 3-5 iteration steps are sufficient to approximate the real roots η1 , η2 of Eq. (9), and the procedure doesn't require a choice of harsh initial approximation. Moreover, many applications show that the computational accuracy of (η1 ,η 2 ) has already been good enough based on the stopping criterion η ( k +1) − η ( k ) 2 ≤ 1 × 10−7 , and the iteration step may be much less.
3.2
Construction of new distance functions
Because of xk (ξ1 , ξ 2 ) being at most a second-order polynomial with respect to ξ1 or ξ 2 , applying the Taylor’s expansion of xk (ξ1 , ξ2 ) in the neighborhood of projection point x p in the parameter plane, we have ∂xk 1 ∂ 2 xk xk (ξ1 , ξ 2 ) = xk (η1 ,η 2 ) + (ξα − ηα ) + (ξα − ηα )(ξ β − η β ) ∂ξα 2 ∂ξα ∂ξ β
(10)
Using Eq. (10), the distance square r 2 between the field point y and the source point
x (ξ1 , ξ2 ) can be written as r 2 (ξ1 , ξ 2 ) = [ xk (ξ1 , ξ 2 ) − yk ][ xk (ξ1 , ξ 2 ) − yk ] = d 2 + (ξα − ηα ) g%α + (ξα − ηα )(ξ β − η β ) g%αβ +(ξα − ηα )(ξ β − η β )(ξγ − ηγ ) g%αβγ
(11)
+(ξα − ηα )(ξ β − η β )(ξγ − ηγ )(ξ μ − η μ ) g%αβγμ
where the summation rule is applied with respect to the Latin indices (taken from the range 1,2,3) and Greek indices (taken from the range 1,2) 7
d 2 = ( yk − xkp )( yk − xkp )
g%α = 2( yk − xkp ) xk ,α , xk ,α =
∂xk ∂ξα
ξ1 =η1 ξ 2 =η2
g%αβ = ( yk − xkp ) xk ,αβ + xk ,α xk ,β = g% βα g%αβγ = xk ,αβ xk ,γ = g% βαγ , g%αβγμ =
1 xk ,αβ xk ,γμ = g% βαγμ 4
(12)
Recall that g%α ≡ 0 , since ( yk − xkp ) is orthogonal to the element and xk ,α is tangential to the element at the projection point x p . Thus, Eq. (11) can be rewritten as r 2 (ξ1 , ξ 2 ) = d 2 + (ξα − ηα )(ξ β − η β ) gˆαβ
= d 2 + (ξ1 − η1 ) 2 g11 + (ξ 2 − η2 ) 2 g 22 + (ξ1 − η1 )(ξ 2 − η2 ) g12 ,
(13)
where
gˆαβ = g%αβ + (ξγ − ηγ ) g%αβγ + (ξγ − ηγ )(ξ μ − η μ ) g%αβγμ , g11 = gˆ11 , g 22 = gˆ 22 , 3.3
g12 = gˆ12 + gˆ 21 = 2 gˆ12
Nearly singular integrals on the second-order elements
By some simple deductions and based on the expression form (13) of the distance function r 2 , the nearly singular integrals in Eq. (5) would be reduced to the following form
Ι=∫
B
0
∫
f ( x, y)
A
0
⎡⎣ d 2 + x 2 g11 ( x, y ) + y 2 g 22 ( x, y ) + xyg12 ( x, y) ⎤⎦
α
dxdy
(14)
where A, B are two constants which are possibly different values in different integrals;
f ( ⋅ ) is a regular function that consists of shape functions, Jacobian and terms which arise from talking the derivative of the integral kernels.
4 Variable transformation The exponential transformation has been proved to be feasible in dealing with the nearly singular linear integrals over curved elements in 2D BEM[32-34]. In this paper, the exponential transformation is extended to treat the nearly singular surface integrals over paraboloidal surface elements in 3D BEM. The extended transformation can be expressed as follows
x = d (em1 +m2s − 1), y = d (en1 +n2t − 1), 8
− 1 ≤ s ≤ 1, − 1 ≤ t ≤ 1
(15)
1 2
where m1 = m2 = ln(1 +
A 1 B ), n1 = n2 = ln(1 + ) . d 2 d
Substituting (15) into Eq. (14), we obtain the following equation
Ι=
1 d 2α −2
1
∫ ∫
1
−1 −1
f ( s, t )m2 n2 em1 + n1 + m2 s + n2t dsdt F α ( s, t )
(16)
where
F ( s, t ) = 1 + (e m1 + m2 s − 1) 2 g11 ( x, y) + (e n1 + n2t − 1) 2 g 22 ( x, y ) + (e m1 + m2 s − 1)(e n1 + n2t − 1) g12 ( x, y ) with x and y being given by Eq. (15). By following the procedures described above, the near singularity of the boundary integrals has been fully regularized. The final integral formulations over curved boundary elements are obtained as shown in Eq. (16), which can be computed straightforward by using standard Gaussian quadrature.
5 Numerical examples In this section, four benchmark numerical examples of 3D potential problems are examined to verify the methodology developed above. The eight-node second-order quadrilateral surface elements are employed to depict the geometry boundary, and the same type of discontinuous interpolation functions to approximate the boundary functions. The boundary quantities are computed by solving the RBIEs with indirect unknowns(1)-(2), while the potentials u and its partial derivatives ∂u ∂x1 at interior points by using the integral representations(3)-(4). The sixteen-points standard Gaussian integration is used for calculation of the various element integrals unless specified otherwise. For the nearly singular integrals the standard Guassian quadrature is employed after application of the proposed transformation. The results obtained by using the present method as well as by the conventional algorithm (without any transformation) and the exact solutions are all presented for convenience of comparison, in order to demonstrate the usefulness of the proposed method. The numerical solution accuracy at single computed point is assessed by means of the relative error defined by RE =
I exa − I num I exa
9
where I num and I exa denote the numerical and exact value at the evaluation points, respectively. Furthermore, the average relative error (ARE) of the multiple computational results is defined by ARE =
1 N
N
∑ RE ( j ) j =1
where RE( j ) denotes the relative error at the jth evaluation point, N is the number of the interior evaluation points. In what follows, d denotes the distance between the evaluation point and the integration boundary element. Example 1 This example concerns the evaluation of nearly singular integral on a curved
surface segment to test the proposed method. The chosen surface segment, named as spherical surface element[49], is represented in parametric form with the usual spherical polar system (θ， ϕ ) . And the element’s geometric parameters are given as follows: θ ∈ [0, π 4] ,
ϕ ∈ [π 4, π 2] , the sphere radius r = 0.1 , and with center (0,0,0) . The projection point of the
evaluation point is located at the center of the element. This example is taken from Ref [43]. Using the conception of Ref.[43], the relative distance between the evaluation point and the element is defined as r0 a1 2 , where a stands for element’s area and r0 is the minimum distance from the evaluation point to element. In Ref [43], a new distance transformation, depending on the above relative distance,
r0 a1 2 , is developed to remove the near singularity based on two local systems ( ρ ,θ ) and (α , β ) . It should be noted that this method possesses two major drawbacks: one is that when
the four subtriangles, each of which is with projection point as one of its vertices, have unsuitable shapes depending on the position of projection point, an adaptive element subdivision technique[49] is necessary to improve the computational accuracy. Such a procedure is sometime very cumbersome. The other drawback is that the distance transformation is not very effective for evaluation of nearly hypersingular integrals due to not fully eliminating their singularities. On the contrary, this study presents a general methodology for numerical evaluation of the nearly singular 2D integrals over high-order surface elements arising in 3D BEM. First, the curved surface segment is modeled by an eight-node second-order quadrilateral geometry element, and then an exact formula of distance is proposed, in which the ‘exact’ means that this formula is equal to the real distance from the evaluation point to a generic point of this curved surface element. Finally, an extended exponential transformation is developed to fully remove the near singularities of integrals with various kernels. Compared with the distance 10
transformation, the present method always works no matter where the position point is located. Furthermore, the present method essentially no need to define the aforementioned relative distance, which is used here only for comparison purpose. Table 3 and 4 list the relative errors for numerical evaluation of integrals with the kernels u * and q* , respectively, with the relative distance d a1 2 changing, using both the present
method and the distance transformation. We can observe from the Table 3 and 4 that for these very simple test problems, two methods can achieve very similar accuracies, and meanwhile, the relative errors for integrals with kernel u * are both very small with the order less than 10−5 , whereas for integrals with kernel q* , with the relative errors reaching 10−4 .
Example 2 As shown in Fig. 2, this example concerns a three dimensional spherical
structure with radius 1.0. The prescribed temperature on the boundary is
u = x12 2 − x22 2 + 2 x1 + x3 The spherical surface is divided into one hundred surface elements of the second-order, and the same type of discontinuous interpolation shape functions is adopted to approximate the boundary functions. The numerical solutions for the potentials u and its derivatives ∂u ∂x1 (in the x1 direction) at internal points are listed in Tables 5 and 6, respectively, hence we can see that when the evaluation points are not too close to the boundary, both the methods with and without transformation of the integration variables are effective and can give acceptable results. As the evaluation point approaches the boundary element of integration, i.e., when the distance of the internal point from the integration element is equal to or less than 0.001, the results of the conventional method become less satisfactory. On the other hand, the results of the proposed method are still steady and satisfactory even when the distance of the evaluation point to the integration element reaches 1E − 10 . This can be seen from the relative errors with respect to the exact solutions which are also shown in Tables 5 and 6 and demonstrate the efficiency and the usefulness of the developed algorithm. Furthermore, on the inner spherical surface
S1 : x12 + x22 + x32 = 0.9992 or
⎧ x1 = 0.999sin θ cos ϕ ⎪ ⎨ x2 = 0.999sin θ sin ϕ , 0 ≤ θ ≤ π , 0 ≤ ϕ ≤ 2π , ⎪ x = 0.999 cos θ ⎩ 3
400 interior points, uniformly spaced according to θ and ϕ , are taken into account. Fig. 3(a) and (b) display the profiles of the analytical solutions for the potentials u and its derivatives ∂u ∂x1 on the inner spherical surface S1, respectively, and Fig. 4(a) and (b) show the surfaces 11
of the numerical solutions for the potentials u and its derivatives ∂u ∂x1 at these 400 interior points, respectively. Hence we can see through the comparison of Fig. 3(a)-(b) and Fig. 4(a)-(b) that the numerical results match the exact solution very well. Fig. 5(a) and (b) show the relative error surfaces of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these 400 interior points, where their AREs are 8.7705 ×10−4 and 2.2353 ×10−3 at these 400 interior points, respectively. Hence it can be seen that the proposed method is accurate.
Example 2 In the second example, we consider a mixed boundary value problem. As
shown in Fig. 6, a central hollow sphere with inner and outer radii are r = 1 and R = 2 , respectively. On the outer surface, the Dirichlet boundary condition is assumed u = u , while on the inner surface, it is Neumann boundary condition q = q , with
u = 2 x1 + x2 + x3 + x1 x2 + 3 x2 x3 + 2 x1 x3 + 5 q = (2 + x2 + 2 x3 )(− x1 ) + (1 + x1 + 3x3 )(− x2 ) + (1 + 3 x2 + 2 x1 )(− x3 ) There are 120 total second-order surface elements spread on the boundaries of the sphere shell, 48 elements on the inner boundary, and 72 elements on the outer boundary, and the same type of discontinuous interpolation shape functions is adopted to approximate the boundary functions. In this section, only domain quantities near the outer boundary are presented, and for those near the inner boundary, similar results can be obtained. Taking one interior point near each boundary element on outer surface, we consider total N = 80 interior points, whose distance d from the integration elements varies from 1.0E − 1 to 1.0E − 9 . The dependences of the AREs for the potentials u and its derivative
∂u ∂x1 on d are shown in Fig. 7 and 8, respectively. Fig. 7 shows that AREs of results obtained by the conventional method (without any transformation) are rather large to be accepted when the distance of the evaluation point to the out integration element is 0.001 or smaller. On the other hand, the results of the proposed method are excellently consistent with the exact solutions with the largest ARE less than 0.0006 even when the distance of the evaluation point to the out integration elements reaches 1.0E − 9 . In Fig. 8, the results computed by the conventional method are unacceptable when d ≤ 0.01 , in contrast to the very high accuracy and stability of the results by the proposed method even when the distance of the evaluation point to the integration element is as small as 1.0E − 9 . Furthermore, the convergence curves of the computed u and its derivatives ∂u ∂x1 at points (1.999999, 0, 0) are shown in Fig. 9 and 10, respectively, hence we can observe that the convergence rates remain monotonic and rapid even when the distance from the field point to 12
the boundary is as small as 1.0E − 7 .
Example 3 The last example concerns a problem in a torus centered at origin, as shown
in Fig. 11, with the exterior radius and interior radius being R = 3 and r = 1 , respectively. The parametric equation of the boundary surface is x1 = ( R + r cos θ ) cos ϕ , x2 = ( R + r cos θ ) sin ϕ , x3 = r sin θ , 0 ≤ θ ≤ 2π , 0 ≤ ϕ ≤ 2π The prescribed potential distribution along the boundary is u = x12 − x32 + x1 x2 + x2 x3 + 2 .
To solve this problem numerically the boundary is discretized by eighty second-order quadrilateral surface elements. Near each boundary element, one interior point is chosen, and thus a total of N = 80 interior points are taken into account. The average relative error curves of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these points are shown in Fig. 12 and 13, respectively. Hence it can be seen that when the evaluation points are not too close to the integration element ( d = 0.1 ), the conventional method and the proposed method are both efficient, but the conventional method fails as the evaluation points are closer to the boundary. On the other hand, the results obtained by the proposed method are stable and satisfactory even when the distance of the evaluation point to the integration element is equal to 1.0E − 9 or even smaller. Furthermore, on the inner torus S2 : x1 = (3 + 0.99 cos θ ) cos ϕ , x2 = (3 + 0.99 cos θ ) sin ϕ , x3 = 0.99sin θ , 0 ≤ θ ≤ 2π , 0 ≤ ϕ ≤ 2π , 320 interior points, uniformly spaced according to θ and ϕ , are taken into account. Fig. 14 (a) and (b) display the profiles of the analytical solutions for the potentials u and its derivatives
∂u ∂x1 on the inner torus S2 , respectively, and Fig. 15 (a) and (b) show the surfaces of the numerical solutions for the potentials u and its derivatives ∂u ∂x1 at these 320 interior points, respectively. Hence we can observe from the comparison of the Fig. 14 (a)-(b) and Fig. 15 (a)-(b) that the numerical results match the exact solutions very well. Fig. 16 (a) and (b) show the relative error surfaces of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these 320 interior points, respectively, where the AREs for the potentials u and its partial derivatives ∂u ∂x1 at these 320 interior points are 1.1436 ×10-3 and 2.6382 ×10-3 , respectively. On the other hand, when the boundary is divided into one hundred and sixty second-order quadrilateral surface elements, the relative error surfaces of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these 320 13
interior points are shown in Fig. 17(a) and (b), respectively, in which their AREs are 4.2060 × 10-4 and 6.7721 × 10-4 at these 320 interior points, respectively. Hence we can see through the comparison of Fig. 16(a)-(b) and Fig. 17(a)-(b) that the proposed method is accurate and convergent as the number of discretization boundary elements increases.
6 Conclusions In this paper, a general scheme is proposed in order to calculate the nearly singular integrals occurring on high-order geometrical elements, which arise in the regularized indirect BIEs formulation. The new distance function r 2 is constructed when the high-order surface elements are adopted for modeling the boundary geometry, and then a nonlinear transformation, based on the exponential functions, is developed to remove or damp out the near singularities of integrands for considered integrals before applying the standard Guassian quadrature to numerical integration. Three numerical examples with exact benchmark solutions are presented to test the proposed scheme, yielding very promising results even when the distance between the evaluation point and the integration element is as small as 1.0E − 9 . The results verify the feasibility and the effectiveness of the proposed scheme, with
eliminating the boundary layer effect in the regularized BEM formulation. Compared with existing approaches, the presented scheme makes the first attempt to evaluate the nearly singular integrals arising in the regularized indirect BEM formulations when the high-order quadrilateral surface elements are used. The proposed scheme also can be applied to other problems in BEM, such as sensitivity analyses, contact problems, and thin-body problems. Some work for thin structures is already underway and will be reported in a subsequent paper. Acknowledgements
The support from the Opening Fund of the State Key Laboratory of Structural Analysis for Industrial Equipment (GZ1017) and the National Natural Science Foundation of Shangdong Province of China (ZR2010AZ003) are gratefully acknowledged. This article has been produced partially with the financial assistance of the European Regional Development Fund (ERDF) under the Operational Programme Research and Development/Measure 4.1 Support of networks of excellence in research and development as the pillars of regional development and support to international cooperation in the Bratislava region/ Project No. 26240120020 Building the centre of excellence for research and 14
development of structural composite materials – 2nd stage.
REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
[17]
[18]
M. Tanaka, V. Sladek, J. Sladek, Regularization techniques applied to boundary element methods, Applied Mechanics Reviews 47 (1994) 457-499. Y.J. Liu, On the simple solution and non-singular nature of the BIE/BEM-a review and some new results, Engineering Analysis with Boundary Elements 24 (2000) 286-292. J.T. Chen, Hong H.K, Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series, Applied Mechanics Reviews 52 (1999) 17-33. Y.J. Liu, T.J. Rudolphi, New identities for fundamental solutions and their applications to non-singular boundary element formulations, Comput Mech 24 (1999) 286-292. M. Guiggiani, G. Krishnasamy, T.J. Rudolphi, et al, A general algorithm for the numerical solution of hyper-singular boundary integral equations, (ASME) J. Appl. Mech, 59 (1992) 604-614. M. Guiggiani, A. Gigante, A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method, Journal of Applied Mechanics (ASME) 57 (1990) 906-915. V. Sladek, J. Sladek, Singular integrals and boundary elements, Computer Methods in Applied Mechanics and Engineering 157 (1998) 251-266. J.T. Chen, A.C. Wu, Null-field integral equation approach for multi-inclusion problem under antiplane shear, ASME, Journal of Applied Mechanics 74 (2007) 469-487. X.W. Gao, An effective method for numerical evaluation of general 2D and 3D high order singular boundary elements, Comput Methods Appl. Mech, Engrg 199 (2010) 2856-2864. G. Krishnasamy, L.W. Schmerr, T.J. Rudolphi, F.J. Rizzo, Hypersingular boundary integral equations: some applications in acoustic and elastic wave scattering, J Appl Mech 57 (1990) 404-414. Y.M. Zhang, V. Sladek, J. Sladek, etc, A new boundary integral equation formulation for plane orthotropic elastic media, Applied Mathematical Modelling 36 (2012) 4862-4875. T. Dirgantara, M.H. Aliabadi, Crack growth analysis of plates loaded by bending and tension using dual boundary element method, International Journal of Fracture 105 (2000) 27-74. M.H. Aliabadi, D. Martin, Boundary element hyper-singular formulation for elastoplastic contact problems, International Journal for Numerical Methods in Engineering 48 (2000) 995-1014. D. Zhang, F.J. Rizzo, Y.J. Rudolphi, Stress intensity sensitivities via hypersingular boundary element integral equations, Computational Mechanics 23 (1999) 389-396. J.F. Luo, Y.J. Liu, E.J. Berger, Analysis of two-dimensional thin structures (from micro- to nano-scales) using the boundary element method, Computational Mechanics 22(5) (1998) 404-412. Y.M. Zhang, Y. Gu, J.T. Chen, Boundary element analysis of 2D thin walled structures with high-order geometry elements using transformation, Engineering Analysis with Boundary Elements 35 (2011) 581-586. Y.J. Liu, Analysis of shell-like structures by the boundary element method based on 3-D elasticity: formulation and verification, International Journal for Numerical Methods in Engineering 41 (1998) 541-558. H.B. Chen, P. Lu, E. Schnack, Regularized algorithms for the calculation of values on and near boundaries in 2D elastic BEM, Engineering Analysis with Boundary Elements 25(10) (2001) 851-876. 15
[19] V. Sladek, J. Sladek, M. Tanaka, Regularization of hypersingular and nearly singular integrals in the potential theory and elasticity, International Journal for Numerical Methods in Engineering 36(10) (1993) 1609-1628. [20] C.Z. Cheng, Z.R. Niu, N. Recho, Z.Y. Yang, R.Y. Ge, A natural stress boundary integral equation for calculating the near boundary stress field, Computers Structures 89 (2011) 1449-1455. [21] S. Mukerjee, M.K. Chati, X.L. Shi, Evaluation of nearly singular integrals in boundary element contour and node methods for three-dimensional linear elasticity, International Journal of Solids and Structures 37(51) (2000) 7633-7654. [22] J.J. Granados, G. Gallego, Regularization of nearly hypersingular integrals in the boundary element method, Engineering Analysis with Boundary Elements 25(3) (2001) 165-184. [23] Y.C. Shiah, Y.S. Shih, Regularization of nearly singular integrals in the boundary element analysis for interior anisotropic thermal field near the boundary, Engineering Analysis with Boundary Element 30(2) (2007) 219-230. [24] Z.R. Niu, W.L. Wendland, X.X. Wang, H.L. Zhou, A semi-analytic algorithm for the evaluation of the nearly singular integrals in three-dimensional boundary element methods, Comput Methods Appl Mech Eng 194 (2005) 1057-1074. [25] H.L. Zhou, Z.R. Niu, C.Z. Cheng, Z.W. Guan, Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for anisotropic potential problems, Comput Struct 86 (2008) 1656-1671. [26] S. Nintcheu Fata, Semi-analytic treatment of nearly-singular Galerkin surface integrals, Applied Numerical Mathematics 60 (2010) 974-993. [27] Y.M. Zhang, Y. Gu, J.T. Chen, Stress analysis for multilayered coating systems using semi-analytical BEM with geometric non-linearities, Computational Mechanics 47 (2011) 493-504. [28] X.W. Gao, K. Yang, J. Wang, An adaptive element subdivision technique for evaluation of various 2D singular boundary integrals, Engineering Analysis with Boundary Elements 32 (2008) 692-696. [29] Telles, JCF, A self-adaptive coordinate transformation for efficient numerical evaluation of general boundary element integral, International Journal for Numerical Methods in Engineering 24 (1987) 959-973. [30] Letizia Scuderi, A new smoothing strategy for computing nearly singular integrals in 3D Galerkin BEM, Journal of Computational and Applied Mathematics 225 (2009) 406-427. [31] G.Z. Xie, F.L. Zhou, J.M. Zhang, X.S. Zheng, etc, New variable transformations for evaluating nearly singular integrals in 3D boundary element method, Engineering Analysis with Boundary Elements 37 (2013) 1169-1178. [32] M. Barbara, Johnston, P. R. Johnston, David Elliott, A new method for the numerical evaluation of nearly singular integrals on triangular elements in the 3D boundary element method, Journal of Computational and Applied Mathematics 245 (2013) 148-161. [33] M. Cerrolaza, E. Alarcon, A bi-cubic transformation for the numerical evaluation of the Cauchy principal value integrals in boundary methods, International Journal for Numerical Methods in Engineering 28(5) (1989) 987-999. [34] P.R. Johnston, Application of Sigmoidal transformations to weakly singular and near singular boundary element integrals, International Journal for Numerical Methods in Engineering 45(10) (1999) 1333-1348. [35] V. Sladek, J. Sladek, M. Tanaka, Optimal transformations of the integration variables in computation of singular integrals in BEM, International Journal for Numerical Methods in Engineering 47(7) 16
[36] [37] [38]
[39] [40]
[41]
[42]
[43]
[44]
[45] [46] [47]
[48] [49] [50]
(2000) 1263-1283. W.J. Ye, A new transformation technique for evaluating nearly singular integrals. Comput Mech 42 (2008) 457-466. Peter R. Johnston, David Elliott, A sinh transformation for evaluating nearly singular boundary element integrals, Int. J. Numer. Meth. Engng 62 (2005) 564-578. M. Barbara, Johnston, Peter R. Johnston, David Elliott, A sinh transformation for evaluating two-dimensional nearly singular boundary element integrals, Int. J. Numer. Meth. Engng 69 (2007) 1460-1479. Q. Huang, T.A. Cruse, Some notes on singular integral techniques in boundary element analysis, International Journal for Numerical Methods in Engineering 36(15) (1993) 2643- 2659. H. Ma, N. Kamiya, Distance transformation for the numerical evaluation of near singular boundary integrals with various kernels in boundary element method, Engineering Analysis with Boundary Elements 26(4) (2002) 329-339. H. Ma, N. Kamiya, A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two-and three-dimensional elasticity, Comput Mech 29 (2002) 277-288. Y. Miao, W. Li, J.H. Lv, X.H. Long, Distance transformation for the numerical evaluation of nearly singular integrals on triangular elements, Engineering Analysis with Boundary Elements 37 (2013) 1311-1317. X.Y. Qin, J.M. Zhang, G.X. Xie, etc, A general algorithm for the numerical evaluation of nearly singular integrals on 3D boundary element, Journal of Computational and Applied Mathematics 235 (2011) 4174-4186. Y.M. Zhang, C.L. Sun, A general algorithm for the numerical evaluation of nearly singular boundary integrals in the equivalent non-singular BIEs with indirect unknowns, Journal of the Chinese Institute of Engineers 31(3) (2008) 437-447. Y.M. Zhang, Y. Gu, J.T. Chen, Boundary layer effect in BEM with high order geometry elements using transformation, Computer Modeling in Engineering & Sciences 45(3) (2009) 227-247. Y.M. Zhang, Y. Gu, J.T. Chen, Boundary element analysis of the thermal behaviour in thin-coated cutting tools, Engineering Analysis with Boundary Elements 34(9) (2010) 775-784. Y.M. Zhang, Y. Gu, J.T. Chen, Internal stress analysis for single and multilayered coating systems using the boundary element method, Engineering Analysis with Boundary Elements 35(4) (2011) 708-717. Y.M. Zhang, W.Z. Qu, J.T. Chen, BEM analysis of thin structures for thermoelastic problems, Engineering Analysis with Boundary Elements 37 (2013) 441-452. X.Y. Qin, J.M. Zhang, G.Y. Li, et al., An element implementation of the boundary face method for 3D potential problems, Eng. Anal. Bound. Elem. 34 (2010) 934–943. Brebbia, C. A., Tells, J. C. F., and Wrobel, L. C. Boundary Element Technique, Springer-Verlag, Berlin, Germany(1984).
17
Figure
Figures：
Fig. 1. Eight-noded quadrilateral surface element, interior points and projection points
Example 2: Fig. 2. Discretization of a unit sphere with 100 second-order surface elements
(a)
(b)
Example 2: Fig. 3. Profiles of analytical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 2: Fig. 4. Surfaces of numerical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 2: Fig. 5. Surfaces of REs for the potentials (a) and its derivatives (b)
Example 3: Fig. 6. Discretization for one-eight of a hollow sphere with 15 second-order surface elements
Example 3: Fig. 7. AREs for the potentials u at interior points approaching to the boundary
Example 3: Fig. 8. AREs for the derivatives u x1 at interior points approaching to the boundary
Example 3: Fig. 9. Convergence curve of the computed u
Example 3: Fig. 10. Convergence curve of the computed u x1
Example 4: Fig. 11. Discretization of the torus with 80 second-order surface elements
Example 4: Fig. 12. AREs of the potentials
u
at interior points approaching to the boundary
Example 4: Fig. 13. AREs of the fluxes
(a)
u x1
at interior points approaching to the boundary
(b)
Example 4: Fig. 14. Profiles of the analytical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 4: Fig. 15. Profiles of numerical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 4: Fig. 16. Surfaces of REs for the potentials (a) and its derivatives (b) with eighty disretization boundary elements
(a)
(b)
Example 4: Fig. 17. Surfaces of REs for the potentials (a) and its derivatives (b) with hundred and sixty disretization boundary elements
Table
Tables：
Table 1 Coordinates ( x1j , x3j , x3j ) of eight nodes x j . Nodes x j
x1j
x2j
x3j
x1
0.304322331872978
0.936607830800249
0.173648177666930
x
2
-0.304322331872978
0.936607830800249
0.173648177666930
x
3
-0.304322331872978
0.936607830800249
-0.173648177666930
x
4
0.304322331872978
0.936607830800249
-0.173648177666930
x
5
6.0300091581E-017
0.984807753012208
0.173648177666930
x
6
-0.309016994374947
0.951056516295154
6.1230317691E-017
x
7
6.0300091581E-017
0.984807753012208
-0.173648177666930
x8
0.309016994374947
0.951056516295154
6.12303176911E-017
Table 2 Iterative results of projection points (1 ,2 ) with criterion ( k 1) ( k )
2
1 1011 .
Projection points (1 ,2 )
Interior points ( y1, y2 , y3 )
Iterations
-0.2290578 -0.2022083
0.9608331 0.9600097
0.0573038 -0.1289337
0.749627609474804 0.666556026866918
-0.33316956351993 0.749130717587247
5 5
0.1508432
0.9737583
-0.0859558
-0.500012936834761
0.499128077463858
4
-0.1528704
0.9767622
-0.0343823
0.499957919680965
0.199811834063831
4
0.1018787
0.9830678
0.0429779
-0.333392866230069
-0.249608365402485
3
-0.1018034
0.9823725
0.0573038
0.333394259894939
-0.332871412077530
3
-0.1274570
0.9808801
-0.0143202
0.416671464829870
8.32713312822E-002
4
0.1015781
0.9814146
0.1002760
-0.332996240934481
-0.581755480500003
4
-0.0760863
0.9715013
0.0992632
0.252591808052607
-0.582748457267639
4
-0.0509764
0.9875403
-0.028657
0.166712963275079
0.166626480502265
3
Table 3 Relative errors of various integrals with kernel
u * on a spherical surface element.
r0 a1 2
101
102
103
104
105
106
Reference solution
0.0077988
0.0079031
0.0079136
0.0079147
0.0079148
0.0079148
Ref [43] : ( , )
6.42E-9
2.52E-7
1.99E-6
6.09E-6
7.44E-6
9.11E-6
Ref [43] : ( , )
7.16E-9
2.52E-7
1.99E-6
6.09E-6
7.44E-6
9.11E-6
Present method 10
7.2694E-9
2.5498E-7
1.8427E-6
5.0576E-6
4.8640E-6
4.5634E-6
Table 4 Relative errors of various integrals with kernel
q* on a spherical surface element.
r0 a1 2
101
102
103
104
105
106
Reference solution
0.2839962
0.2890126
0.2895177
0.2895683
0.2895733
0.2895783
Ref [43] : ( , )
1.81E-7
4.48E-6
1.59E-6
7.89E-5
1.77E-4
1.50E-4
Ref [43] : ( , )
1.76E-7
4.48E-6
1.59E-6
7.89E-5
1.77E-4
1.50E-4
Present method 10
1.6743E-7
4.1356E-6
1.5059E-6
7.7185E-5
2.6966E-4
2.1417E-4
Table 5 Potentials u at internal points increasingly close to the boundary. Present Distance d
Exact
No transform Numerical
Relative error
0.1 0.01
0.2203216E+01 0.2468210E+01
0.2203587E+01 0.2478618E+01
0.2201609E+01 0.2467192E+01
5.604263E-04 4.126834E-04
0.001
0.2495155E+01
0.2472128E+01
0.2494217E+01
3.761218E-04
0.0001
0.2497854E+01
0.2469416E+01
0.2496925E+01
3.720910E-04
0.00001
0.2498151E+01
0.2469117E+01
0.2497223E+01
3.716440E-04
0.000001
0.2498154E+01
0.2469086E+01
0.2497225E+01
3.716563E-04
0.0000001
0.2497769E+01
0.2469083E+01
0.2496801E+01
3.717048E-04
0.00000001
0.2498154E+01
0.2469083E+01
0.2497225E+01
3.718859E-04
0.000000001
0.2498154E+01
0.2497224E+01
0.2496801E+01
3.721940E-04
0.0000000001
0.2498154E+01
0.2469083E+01
0.2497223E+01
3.726790E-04
Table 6 Potential derivatives u x1 at internal points increasingly close to the boundary. Present Distance d
Exact
No transform Numerical
Relative error
0.1 0.01
0.2899256E+01 0.2989385E+01
0.2910363E+01 0.2507812E+01
0.2896660E+01 0.2997269E+01
8.955288E-04 2.637319E-03
0.001
0.2998385E+01
-0.2714809E+01
0.3008335E+01
3.318621E-03
0.0001
0.2999285E+01
-0.3296846E+01
0.3009453E+01
3.390136E-03
0.00001
0.2999375E+01
-0.3350076E+01
0.3009565E+01
3.397357E-03
0.000001
0.2999384E+01
-0.3355347E+01
0.3009575E+01
3.397986E-03
0.0000001
0.2999385E+01
-0.3355873E+01
0.3009568E+01
3.395174E-03
0.00000001
0.2999385E+01
-0.3355926E+01
0.3009454E+01
3.357243E-03
0.000000001
0.2999385E+01
-0.3355931E+01
0.3009376E+01
3.331104E-03
0.0000000001
0.2999385E+01
-0.3355932E+01
0.3009223E+01
3.280253E-03
Research highlights • We present a new indirect regularized boundary element formulation. •
A new distance formula for high-order geometry elements is constructed.
•
The exponential transformation is developed to remove the near singularities.
•
The efficiency and applicability of the scheme are verified by some examples.
S0377-0427(14)00384-7 http://dx.doi.org/10.1016/j.cam.2014.08.027 CAM 9781
To appear in:
Journal of Computational and Applied Mathematics
Received date: 11 October 2013 Revised date: 20 May 2014 Please cite this article as: Y. Zhang, X. Li, V. Sladek, J. Sladek, X. Gao, A new method for numerical evaluation of nearly singular integrals over high-order geometry elements in 3D BEM, Journal of Computational and Applied Mathematics (2014), http://dx.doi.org/10.1016/j.cam.2014.08.027 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A new method for numerical evaluation of nearly singular integrals over high-order geometry elements in 3D BEM Yaoming Zhang a,b,*, Xiaochao Li a , Vladimir Sladek c, Jan Sladek c, Xiaowei Gao b a
c
Institute of Applied Mathematics, Shandong University of Technology,Zibo255049,PR China b State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China Slovak Academy of Sciences, Institute of Construction and Architecture, 84503 Bratislava, Slovakia
Abstract: This work presents a new method for numerical computation of the two-dimensional nearly singular integrals with using the eight-node second-order quadrilateral surface elements in 3D BEM. A new indirect regularized boundary element formulation excluding the CPV (Cauchy Principal Value) and HFP (Hadamard-Finite-Part) integrals is proposed. Based on this, a new approximation formula of the distance from the fixed calculation point to a generic point of the aforementioned surface geometry elements is developed firstly, and then the exponential transformation, which has been widely employed in 2D BEM, is extended to 3D BEM to remove the near singularities of integrands for considered integrals. Several numerical examples are given to verify the high efficiency and the stability of the proposed scheme. Key Words: BEM; Nearly singular integrals; High-order geometry elements; Transformation; 3D potential problems.
1 Introduction The BEM as an important numerical technique has been widely used in many areas, but the accuracy of its calculation depends heavily on the accuracy and efficiency of the computation of integrals with singular and nearly singular kernels. These kernels become singular when the collocation point belongs to the integration elements, and many effective methods have been developed to handle them[1-11]. Another important issue is the integration of the kernels for the collocation points which are close to but not on the integration element. Then, the considered integrals, although regular in the sense of mathematics, are termed nearly singular integrals since they exhibit similar characteristics of the singularity and generally can not be calculated accurately by standard integration quadratures. This is so-called boundary layer effect in BEM. In this work, we focus on numerical computation of the two-dimensional nearly singular integrals that arise in the *
Corresponding author at :Institute of Applied Mathematics, Shandong University of Technology, Zibo255049, PR China. . E-mail address: [email protected] (Y. M. Zhang).
1
solution of three-dimensional(3D) boundary element method(BEM) using eight-node second-order quadrilateral surface elements. Such problems would be more crucial to some of engineering applications, such as crack problems when the crack tip is deformed to have a small opening displacement[12], contact problems when the contact area of the two contacting bodies is very small[13], inverse and sensitivity problems[14], as well as thin or shell-like structure problems[15-17], etc. In these situations, the number of the nearly singular integrals, which have to be evaluated for a system matrix, can be much larger than that of the singular integrals. Effective evaluation of nearly singular integrals has attracted a great amount of attention in recent years. Various numerical techniques have been developed to remove the near singularities, such as rigid-body displacement and other global regularization methods[17-23], analytical or semi-analytical formulas[24-27], adaptive element subdivision techniques [28-29], the use of polar coordinate transformations[30-32], and various nonlinear transformations[33-43]. Among these methods, the variable transformations technique seems be a promising method due to its wide suitability and high accuracy and stability. The methods developed so far include, but are not limited to, polynomial transformation[33], sigmoidal transformation[34], coordinate optimal transformation[35], degenerate mapping method[13,36],
sinh
transformation[37-38],
rational
transformation[39],
distance
transformation[40-43], and exponential transformation[44-48]. A comprehensive review was given by present author[44]. Although great progresses have been achieved for each of the above methods, a number of drawbacks are associated with these methods and these include the fact that some, which benefit from the strategies for calculating the singular integrals, are failed to fully eliminate near singularities, some apply only to the planar elements, some essentially use the first-order functions to depict the surface geometry boundary although the authors claim that the problems are considered over the curved surface elements, and others are tailored specifically to the form of the kernel functions and therefore lack wide applicability. In previous work, the present author proposed a nonlinear transformation, named as the exponential transformation, to remove or damp out the near singularities of integrands for considered integrals in 2D BEM[44]. This transformation, based on the idea of diminishing the difference of the orders of magnitude or the scale of change of addition factors in the denominator of the kernels, is more simple and applicable than the previous works of the above researches according to extensive numerical experiments and applications. It is also demonstrated that the exponential transformation offers two primary advantages over most of existing methodologies. One is the accuracy and ease of implementation, the other one is applicable to a broad range of integrals with various kernels without extra computational 2
effort. In present work, we would extend such transformation to 3D BEM. The usage of high-order geometry elements, most usually of the eight-node second-order quadrilateral surface elements, for efficient evaluation of the nearly singular integrals is necessary in many cases. Clearly, the advantage of using the high-order geometry elements in complex geometrical domains does not concern only its power to achieve higher calculation accuracy, but even with a small number of computational elements. More crucial, computational models of thin structure problems or the case when the calculation points are very close to the actual boundary demand a higher level of geometry approximation, and the usage of high-order geometry in computational model can meet this requirement. For example, if the boundary geometry is depicted by using planar elements, the planar elements of the outer surface will contact or even pass through the inner boundary if the thickness of the considered structure is very small or the case where some interior points very close to the actual boundary may be located outside the boundary elements, and obviously this calculation model is fallacious. Consequently, the actual geometry of considered domain cannot be described faithfully, and thus lower-order geometry approximation will fail to yield reliable results for such problems. In order to avoid this phenomenon, very fine meshes must be used in this situation, but this would yields too much preprocessing and CPU time. In addition, a great number of meshes will produce a lot of artificial corners that will lead to the discontinuity of the tangent derivative of the boundary unknowns. This is fatal to many engineering problems. In this paper, a general methodology is proposed to compute the nearly singular integrals that arise in the solution of 3D BEM using the eight-node second-order quadrilateral surface elements. A new indirect regularized boundary element formulation excluding the CPV and HFP integrals is presented. Based on this, a new approximation formula of the distance between the fixed calculation point and a generic point of curved surface elements is developed here firstly, and then the exponential transformation, which has been widely applied in 2D BEM, is developed to remove the near singularities of the considered integrals in 3D BEM before applying the standard Gaussian quadrature to numerical integration. Four numerical examples with exact benchmark solutions are presented to test the proposed scheme, yielding very promising results even when the internal point is very close to the boundary. The outline of the rest of this paper is as follows. The indirect regularized BIEs excluding the CPV and HFP integrals are presented in Section2. Then, in Section 3, the expression of the nearly singular integrals over the eight-node second-order quadrilateral surface elements is investigated, including the determination of the projection point and construction of the distance functions. Section 4 considers the regularization of nearly 3
singular integrals over high-order geometrical elements using the exponential transformation. In Section 5, the accuracy and stability of the proposed scheme are tested in three 3D potential examples with known benchmark solutions. Finally, the conclusions are provided in Section 6.
2 Regularized boundary integral equations(RBIEs) It is well known that having known all relevant boundary quantities, the solutions at interior points can be obtained by using the integral representation. Thus, the accuracy of boundary quantities directly affects the validity of the domain quantities, but not opposite. This is a great advantage of pure boundary integral formulations as compared with the domain discretization methods such as FEM, for instance. On the other hand, we have to deal with the singular boundary integrals when calculating the boundary unknowns. Therefore, a good choice is to use the regularized BIEs with indirect unknowns. In this paper, we always assume that Ω is a bounded domain in R 3 , Ω c is its open complement, and Γ denotes the common boundary. n( x ) is the unit outward normal vector on Γ to the domain Ω at the point x . With omitting the body sources in potential ˆ can be expressed as problems, the regularized BIEs with indirect unknowns on Ω u ( y ) = ∫ φ ( x )u * ( x, y )d Γ, y ∈ Γ Γ
∂u ( y ) ˆ ∂u * ( x , y ) = Sφ ( y ) + ∫ [φ ( x ) − φ ( y ) ] dΓ Γ ∂nˆ ( y ) ∂nˆ ( y ) ⎡ ∂u * ( x , y ) ∂u * ( x , y ) ⎤ +φ ( y ) ∫ ⎢ + ⎥ d Γ, y ∈ Γ Γ ∂nˆ ( x ) ⎦ ⎣ ∂nˆ ( y )
(1)
(2)
For the internal point y , the integral equations can be written as ˆ u ( y ) = ∫ φ ( x )u * ( x , y ) d Γ x , y ∈ Ω Γ
ˆ ∇ y u ( y ) = ∫ φ ( x )∇ y u * ( x , y )d Γ x , y ∈ Ω Γ
(3) (4)
In the Eqs.(1)-(4), φ ( x ) is the density function to be determined; u * ( x , y ) denotes the ˆ = Ω , Sˆ = 1 , nˆ ( x ) is the unit outward Kelvin fundamental solution. For interior problems, Ω
ˆ = Ωc , Sˆ = 0 , normal vectors on Γ to domain Ω at point x . For exterior problems, Ω nˆ ( x ) is the unit outward normal vectors on Γ to domain Ωc at point x .
For the discretized form of the Eqs. (3) and (4), when the field point y is far enough from the integration elements, a straightforward application of Gaussian quadrature procedure suffices to evaluate such integrals. However, when the field point y is very close to the 4
integration elements Γ e , the distance r between the field point y and the source point x is almost zero. Hence, the integrals in discretized Eqs. (3)-(4) are nearly singular and the numerical integrations by the standard Gaussian quadrature fail. These nearly singular integrals can be expressed as Ι=∫
Γe
f ( x, y ) dΓ rα
(5)
where r = x − y 2 , α > 0 is a real constant, and f ( x, y ) denotes a well-behaved function.
3 Nearly singular integrals on high-order geometry elements The quintessence of the BEM is to discretize the boundary into a finite number of segments, not necessarily uniform, which are called boundary elements. Two approximations are made over each of these elements. One is about the geometry of the boundary, while the other has to interpolate the spatial variation of the field variables or the boundary unknowns over the element. The plane shape geometry element is not an ideal one as it can not approximate the curved surface boundaries with sufficient accuracy. For this reason, it is recommended to use higher order elements, which approximate geometry and boundary quantities by using higher order interpolation polynomials-usually of the second order. In this paper, the geometry segment is modeled by a continuous paraboloidal element, which has eight knots, namely, the boundary geometry is approximated by the piecewise continuous eight-node second-order quadrilateral surface elements, while the distribution of the boundary quantities over each of these segments is approximated using discontinuous elements, eight nodes of which are located away from the edges of the element. Assume x j = ( x1j , x2j , x3 j ), j = 1,L,8 are the eight knots of the segment Γ j , then
Cartesian coordinates of the points on the element Γ j can be interpolated as[50] 8
xk (ξ1 , ξ 2 ) = ∑ N j (ξ1 , ξ 2 ) xkj , k = 1, 2,3
(6)
1 1 N1 (ξ1 , ξ 2 ) = (1 − ξ1 )(1 − ξ 2 )(−ξ1 − ξ 2 − 1), N 2 (ξ1 , ξ 2 ) = (1 + ξ1 )(1 − ξ 2 )(ξ1 − ξ 2 − 1)， 4 4 1 1 N 3 (ξ1 , ξ 2 ) = (1 + ξ1 )(1 + ξ 2 )(ξ1 + ξ 2 − 1)，N 4 (ξ1 , ξ 2 ) = (1 − ξ1 )(1 + ξ 2 )(−ξ1 + ξ 2 − 1)， 4 4 1 1 N 5 (ξ1 , ξ 2 ) = (1 − ξ12 )(1 − ξ 2 )，N 6 (ξ1 , ξ 2 ) = (1 + ξ1 )(1 − ξ 2 2 ), 2 2 1 1 N 7 (ξ1 , ξ 2 ) = (1 + ξ 2 )(1 − ξ12 ), N8 (ξ1 , ξ 2 ) = (1 − ξ1 )(1 − ξ 2 2 ), − 1 ≤ ξ1 ≤ 1, − 1 ≤ ξ 2 ≤ 1 2 2
(7)
j =1
where
5
where ξ1 , ξ 2 are the local coordinates. 3.1
Determination of the projection point
The minimum distance d from the field point y to the integration element Γ e is defined as the length y − x p , where x p is the projection point of y onto integration element Γ e . Letting (η1 , η 2 ) be the local coordinates of the projection point x p , i.e. x p = ( x1 (η1 ,η2 ), x2 (η1 ,η2 ), x3 (η1 ,η2 )) , then η1 , η2 are the real roots of the following equation ∂xi ⎧ ⎪[ xi (η1 ,η2 ) − yi ] ∂ξ = 0 ⎪ 1 , ⎨ ⎪[ xi (η1 ,η2 ) − yi ] ∂xi = 0 ⎪⎩ ∂ξ 2 in which the summation convention is used, and
i = 1, 2,3
∂xi ∂xi = ∂ξ k ∂ξ k
ξ1 =η1 ξ 2 =η2
(8)
, k = 1, 2 . These assumptions
will be applied also in what follows unless specified otherwise. If the field point y is sufficiently close to the boundary Γ , then x p is inside the integration element, and Eq. (8) has a pair of the unique real roots (η1 , η 2 ) ∈ [−1,1] × [−1,1] . The real roots η1 , η2 can be evaluated numerically by using the Newton’s method. Setting
f1 (η1 ,η2 ) = [ xi (η1 ,η2 ) − yi ]
∂xi ∂x , f 2 (η1 ,η2 ) = [ xi (η1 ,η2 ) − yi ] i ∂ξ1 ∂ξ 2
The formula of the Newton’method can be expressed as F ′(η ( k ) )Δη ( k ) = − F (η ( k ) )
(9)
where Δη ( k ) = η ( k +1) − η ( k ) ,η ( k ) = (η1( k ) ,η2 ( k ) )T , η ( k +1) = (η1( k +1) ,η2( k +1) )T , ⎡ ⎢ ⎡ f1 (η ,η2 ) ⎤ (k ) (k ) ′ , F (η ) = ⎢ F (η ) = ⎢ (k ) (k ) ⎥ ⎢ ⎢⎣ f 2 (η1 ,η 2 ) ⎥⎦ ⎢ ⎣ (k ) 1
(k )
∂f1 ∂η1 ∂f 2 ∂η1
∂f1 ⎤ ∂η 2 ⎥ ⎥ ∂f 2 ⎥ ∂η 2 ⎥⎦
η =η ( k )
here 3 ⎛ ∂x ∂xi ∂ 2 xi ⎞ = ∑⎜ i + [ xi (η1 ,η 2 ) − yi ] ⎟, j , m = 1, 2 . ∂ηm i =1 ⎜⎝ ∂ξ m ∂ξ j ∂ξ j ∂ξ m ⎟⎠
∂f j
It should be emphasized that the effective and accurate computation of the projection 6
point x p , whose local coordinates in the parameter plane are (η1 ,η 2 ) , is crucial to implementation of the proposed scheme. In order to demonstrate the validity and convergence of the above Newton’s iterative scheme，the following numerical example 2 is investigated here. The discretization of boundary surface is as in example 2, and we randomly take an eight-node quadrilateral surface element for which the coordinates of eight nodes are given in the Table 1. Ten interior points Ai , i = 1,L ,10 are considered close to this element, as shown in Fig. 1. Pi , i = 1,L ,10 denote the projection points of Ai , i = 1,L ,10 onto this element, respectively, which have been computed by aforementioned Newton’s iterative scheme. The same one initial approximation η0 = (0, 0) is chosen for all iterative procedures, and the iterative results, including projection points coordinates and number of iterations, are listed in Table 2 based on the stopping criterion η ( k +1) − η ( k ) 2 ≤ 1 × 10−11 . From these results, we can conclude that it doesn’t matter if the field point y is very close to integration elements, 3-5 iteration steps are sufficient to approximate the real roots η1 , η2 of Eq. (9), and the procedure doesn't require a choice of harsh initial approximation. Moreover, many applications show that the computational accuracy of (η1 ,η 2 ) has already been good enough based on the stopping criterion η ( k +1) − η ( k ) 2 ≤ 1 × 10−7 , and the iteration step may be much less.
3.2
Construction of new distance functions
Because of xk (ξ1 , ξ 2 ) being at most a second-order polynomial with respect to ξ1 or ξ 2 , applying the Taylor’s expansion of xk (ξ1 , ξ2 ) in the neighborhood of projection point x p in the parameter plane, we have ∂xk 1 ∂ 2 xk xk (ξ1 , ξ 2 ) = xk (η1 ,η 2 ) + (ξα − ηα ) + (ξα − ηα )(ξ β − η β ) ∂ξα 2 ∂ξα ∂ξ β
(10)
Using Eq. (10), the distance square r 2 between the field point y and the source point
x (ξ1 , ξ2 ) can be written as r 2 (ξ1 , ξ 2 ) = [ xk (ξ1 , ξ 2 ) − yk ][ xk (ξ1 , ξ 2 ) − yk ] = d 2 + (ξα − ηα ) g%α + (ξα − ηα )(ξ β − η β ) g%αβ +(ξα − ηα )(ξ β − η β )(ξγ − ηγ ) g%αβγ
(11)
+(ξα − ηα )(ξ β − η β )(ξγ − ηγ )(ξ μ − η μ ) g%αβγμ
where the summation rule is applied with respect to the Latin indices (taken from the range 1,2,3) and Greek indices (taken from the range 1,2) 7
d 2 = ( yk − xkp )( yk − xkp )
g%α = 2( yk − xkp ) xk ,α , xk ,α =
∂xk ∂ξα
ξ1 =η1 ξ 2 =η2
g%αβ = ( yk − xkp ) xk ,αβ + xk ,α xk ,β = g% βα g%αβγ = xk ,αβ xk ,γ = g% βαγ , g%αβγμ =
1 xk ,αβ xk ,γμ = g% βαγμ 4
(12)
Recall that g%α ≡ 0 , since ( yk − xkp ) is orthogonal to the element and xk ,α is tangential to the element at the projection point x p . Thus, Eq. (11) can be rewritten as r 2 (ξ1 , ξ 2 ) = d 2 + (ξα − ηα )(ξ β − η β ) gˆαβ
= d 2 + (ξ1 − η1 ) 2 g11 + (ξ 2 − η2 ) 2 g 22 + (ξ1 − η1 )(ξ 2 − η2 ) g12 ,
(13)
where
gˆαβ = g%αβ + (ξγ − ηγ ) g%αβγ + (ξγ − ηγ )(ξ μ − η μ ) g%αβγμ , g11 = gˆ11 , g 22 = gˆ 22 , 3.3
g12 = gˆ12 + gˆ 21 = 2 gˆ12
Nearly singular integrals on the second-order elements
By some simple deductions and based on the expression form (13) of the distance function r 2 , the nearly singular integrals in Eq. (5) would be reduced to the following form
Ι=∫
B
0
∫
f ( x, y)
A
0
⎡⎣ d 2 + x 2 g11 ( x, y ) + y 2 g 22 ( x, y ) + xyg12 ( x, y) ⎤⎦
α
dxdy
(14)
where A, B are two constants which are possibly different values in different integrals;
f ( ⋅ ) is a regular function that consists of shape functions, Jacobian and terms which arise from talking the derivative of the integral kernels.
4 Variable transformation The exponential transformation has been proved to be feasible in dealing with the nearly singular linear integrals over curved elements in 2D BEM[32-34]. In this paper, the exponential transformation is extended to treat the nearly singular surface integrals over paraboloidal surface elements in 3D BEM. The extended transformation can be expressed as follows
x = d (em1 +m2s − 1), y = d (en1 +n2t − 1), 8
− 1 ≤ s ≤ 1, − 1 ≤ t ≤ 1
(15)
1 2
where m1 = m2 = ln(1 +
A 1 B ), n1 = n2 = ln(1 + ) . d 2 d
Substituting (15) into Eq. (14), we obtain the following equation
Ι=
1 d 2α −2
1
∫ ∫
1
−1 −1
f ( s, t )m2 n2 em1 + n1 + m2 s + n2t dsdt F α ( s, t )
(16)
where
F ( s, t ) = 1 + (e m1 + m2 s − 1) 2 g11 ( x, y) + (e n1 + n2t − 1) 2 g 22 ( x, y ) + (e m1 + m2 s − 1)(e n1 + n2t − 1) g12 ( x, y ) with x and y being given by Eq. (15). By following the procedures described above, the near singularity of the boundary integrals has been fully regularized. The final integral formulations over curved boundary elements are obtained as shown in Eq. (16), which can be computed straightforward by using standard Gaussian quadrature.
5 Numerical examples In this section, four benchmark numerical examples of 3D potential problems are examined to verify the methodology developed above. The eight-node second-order quadrilateral surface elements are employed to depict the geometry boundary, and the same type of discontinuous interpolation functions to approximate the boundary functions. The boundary quantities are computed by solving the RBIEs with indirect unknowns(1)-(2), while the potentials u and its partial derivatives ∂u ∂x1 at interior points by using the integral representations(3)-(4). The sixteen-points standard Gaussian integration is used for calculation of the various element integrals unless specified otherwise. For the nearly singular integrals the standard Guassian quadrature is employed after application of the proposed transformation. The results obtained by using the present method as well as by the conventional algorithm (without any transformation) and the exact solutions are all presented for convenience of comparison, in order to demonstrate the usefulness of the proposed method. The numerical solution accuracy at single computed point is assessed by means of the relative error defined by RE =
I exa − I num I exa
9
where I num and I exa denote the numerical and exact value at the evaluation points, respectively. Furthermore, the average relative error (ARE) of the multiple computational results is defined by ARE =
1 N
N
∑ RE ( j ) j =1
where RE( j ) denotes the relative error at the jth evaluation point, N is the number of the interior evaluation points. In what follows, d denotes the distance between the evaluation point and the integration boundary element. Example 1 This example concerns the evaluation of nearly singular integral on a curved
surface segment to test the proposed method. The chosen surface segment, named as spherical surface element[49], is represented in parametric form with the usual spherical polar system (θ， ϕ ) . And the element’s geometric parameters are given as follows: θ ∈ [0, π 4] ,
ϕ ∈ [π 4, π 2] , the sphere radius r = 0.1 , and with center (0,0,0) . The projection point of the
evaluation point is located at the center of the element. This example is taken from Ref [43]. Using the conception of Ref.[43], the relative distance between the evaluation point and the element is defined as r0 a1 2 , where a stands for element’s area and r0 is the minimum distance from the evaluation point to element. In Ref [43], a new distance transformation, depending on the above relative distance,
r0 a1 2 , is developed to remove the near singularity based on two local systems ( ρ ,θ ) and (α , β ) . It should be noted that this method possesses two major drawbacks: one is that when
the four subtriangles, each of which is with projection point as one of its vertices, have unsuitable shapes depending on the position of projection point, an adaptive element subdivision technique[49] is necessary to improve the computational accuracy. Such a procedure is sometime very cumbersome. The other drawback is that the distance transformation is not very effective for evaluation of nearly hypersingular integrals due to not fully eliminating their singularities. On the contrary, this study presents a general methodology for numerical evaluation of the nearly singular 2D integrals over high-order surface elements arising in 3D BEM. First, the curved surface segment is modeled by an eight-node second-order quadrilateral geometry element, and then an exact formula of distance is proposed, in which the ‘exact’ means that this formula is equal to the real distance from the evaluation point to a generic point of this curved surface element. Finally, an extended exponential transformation is developed to fully remove the near singularities of integrals with various kernels. Compared with the distance 10
transformation, the present method always works no matter where the position point is located. Furthermore, the present method essentially no need to define the aforementioned relative distance, which is used here only for comparison purpose. Table 3 and 4 list the relative errors for numerical evaluation of integrals with the kernels u * and q* , respectively, with the relative distance d a1 2 changing, using both the present
method and the distance transformation. We can observe from the Table 3 and 4 that for these very simple test problems, two methods can achieve very similar accuracies, and meanwhile, the relative errors for integrals with kernel u * are both very small with the order less than 10−5 , whereas for integrals with kernel q* , with the relative errors reaching 10−4 .
Example 2 As shown in Fig. 2, this example concerns a three dimensional spherical
structure with radius 1.0. The prescribed temperature on the boundary is
u = x12 2 − x22 2 + 2 x1 + x3 The spherical surface is divided into one hundred surface elements of the second-order, and the same type of discontinuous interpolation shape functions is adopted to approximate the boundary functions. The numerical solutions for the potentials u and its derivatives ∂u ∂x1 (in the x1 direction) at internal points are listed in Tables 5 and 6, respectively, hence we can see that when the evaluation points are not too close to the boundary, both the methods with and without transformation of the integration variables are effective and can give acceptable results. As the evaluation point approaches the boundary element of integration, i.e., when the distance of the internal point from the integration element is equal to or less than 0.001, the results of the conventional method become less satisfactory. On the other hand, the results of the proposed method are still steady and satisfactory even when the distance of the evaluation point to the integration element reaches 1E − 10 . This can be seen from the relative errors with respect to the exact solutions which are also shown in Tables 5 and 6 and demonstrate the efficiency and the usefulness of the developed algorithm. Furthermore, on the inner spherical surface
S1 : x12 + x22 + x32 = 0.9992 or
⎧ x1 = 0.999sin θ cos ϕ ⎪ ⎨ x2 = 0.999sin θ sin ϕ , 0 ≤ θ ≤ π , 0 ≤ ϕ ≤ 2π , ⎪ x = 0.999 cos θ ⎩ 3
400 interior points, uniformly spaced according to θ and ϕ , are taken into account. Fig. 3(a) and (b) display the profiles of the analytical solutions for the potentials u and its derivatives ∂u ∂x1 on the inner spherical surface S1, respectively, and Fig. 4(a) and (b) show the surfaces 11
of the numerical solutions for the potentials u and its derivatives ∂u ∂x1 at these 400 interior points, respectively. Hence we can see through the comparison of Fig. 3(a)-(b) and Fig. 4(a)-(b) that the numerical results match the exact solution very well. Fig. 5(a) and (b) show the relative error surfaces of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these 400 interior points, where their AREs are 8.7705 ×10−4 and 2.2353 ×10−3 at these 400 interior points, respectively. Hence it can be seen that the proposed method is accurate.
Example 2 In the second example, we consider a mixed boundary value problem. As
shown in Fig. 6, a central hollow sphere with inner and outer radii are r = 1 and R = 2 , respectively. On the outer surface, the Dirichlet boundary condition is assumed u = u , while on the inner surface, it is Neumann boundary condition q = q , with
u = 2 x1 + x2 + x3 + x1 x2 + 3 x2 x3 + 2 x1 x3 + 5 q = (2 + x2 + 2 x3 )(− x1 ) + (1 + x1 + 3x3 )(− x2 ) + (1 + 3 x2 + 2 x1 )(− x3 ) There are 120 total second-order surface elements spread on the boundaries of the sphere shell, 48 elements on the inner boundary, and 72 elements on the outer boundary, and the same type of discontinuous interpolation shape functions is adopted to approximate the boundary functions. In this section, only domain quantities near the outer boundary are presented, and for those near the inner boundary, similar results can be obtained. Taking one interior point near each boundary element on outer surface, we consider total N = 80 interior points, whose distance d from the integration elements varies from 1.0E − 1 to 1.0E − 9 . The dependences of the AREs for the potentials u and its derivative
∂u ∂x1 on d are shown in Fig. 7 and 8, respectively. Fig. 7 shows that AREs of results obtained by the conventional method (without any transformation) are rather large to be accepted when the distance of the evaluation point to the out integration element is 0.001 or smaller. On the other hand, the results of the proposed method are excellently consistent with the exact solutions with the largest ARE less than 0.0006 even when the distance of the evaluation point to the out integration elements reaches 1.0E − 9 . In Fig. 8, the results computed by the conventional method are unacceptable when d ≤ 0.01 , in contrast to the very high accuracy and stability of the results by the proposed method even when the distance of the evaluation point to the integration element is as small as 1.0E − 9 . Furthermore, the convergence curves of the computed u and its derivatives ∂u ∂x1 at points (1.999999, 0, 0) are shown in Fig. 9 and 10, respectively, hence we can observe that the convergence rates remain monotonic and rapid even when the distance from the field point to 12
the boundary is as small as 1.0E − 7 .
Example 3 The last example concerns a problem in a torus centered at origin, as shown
in Fig. 11, with the exterior radius and interior radius being R = 3 and r = 1 , respectively. The parametric equation of the boundary surface is x1 = ( R + r cos θ ) cos ϕ , x2 = ( R + r cos θ ) sin ϕ , x3 = r sin θ , 0 ≤ θ ≤ 2π , 0 ≤ ϕ ≤ 2π The prescribed potential distribution along the boundary is u = x12 − x32 + x1 x2 + x2 x3 + 2 .
To solve this problem numerically the boundary is discretized by eighty second-order quadrilateral surface elements. Near each boundary element, one interior point is chosen, and thus a total of N = 80 interior points are taken into account. The average relative error curves of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these points are shown in Fig. 12 and 13, respectively. Hence it can be seen that when the evaluation points are not too close to the integration element ( d = 0.1 ), the conventional method and the proposed method are both efficient, but the conventional method fails as the evaluation points are closer to the boundary. On the other hand, the results obtained by the proposed method are stable and satisfactory even when the distance of the evaluation point to the integration element is equal to 1.0E − 9 or even smaller. Furthermore, on the inner torus S2 : x1 = (3 + 0.99 cos θ ) cos ϕ , x2 = (3 + 0.99 cos θ ) sin ϕ , x3 = 0.99sin θ , 0 ≤ θ ≤ 2π , 0 ≤ ϕ ≤ 2π , 320 interior points, uniformly spaced according to θ and ϕ , are taken into account. Fig. 14 (a) and (b) display the profiles of the analytical solutions for the potentials u and its derivatives
∂u ∂x1 on the inner torus S2 , respectively, and Fig. 15 (a) and (b) show the surfaces of the numerical solutions for the potentials u and its derivatives ∂u ∂x1 at these 320 interior points, respectively. Hence we can observe from the comparison of the Fig. 14 (a)-(b) and Fig. 15 (a)-(b) that the numerical results match the exact solutions very well. Fig. 16 (a) and (b) show the relative error surfaces of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these 320 interior points, respectively, where the AREs for the potentials u and its partial derivatives ∂u ∂x1 at these 320 interior points are 1.1436 ×10-3 and 2.6382 ×10-3 , respectively. On the other hand, when the boundary is divided into one hundred and sixty second-order quadrilateral surface elements, the relative error surfaces of the computational results for the potentials u and its partial derivatives ∂u ∂x1 at these 320 13
interior points are shown in Fig. 17(a) and (b), respectively, in which their AREs are 4.2060 × 10-4 and 6.7721 × 10-4 at these 320 interior points, respectively. Hence we can see through the comparison of Fig. 16(a)-(b) and Fig. 17(a)-(b) that the proposed method is accurate and convergent as the number of discretization boundary elements increases.
6 Conclusions In this paper, a general scheme is proposed in order to calculate the nearly singular integrals occurring on high-order geometrical elements, which arise in the regularized indirect BIEs formulation. The new distance function r 2 is constructed when the high-order surface elements are adopted for modeling the boundary geometry, and then a nonlinear transformation, based on the exponential functions, is developed to remove or damp out the near singularities of integrands for considered integrals before applying the standard Guassian quadrature to numerical integration. Three numerical examples with exact benchmark solutions are presented to test the proposed scheme, yielding very promising results even when the distance between the evaluation point and the integration element is as small as 1.0E − 9 . The results verify the feasibility and the effectiveness of the proposed scheme, with
eliminating the boundary layer effect in the regularized BEM formulation. Compared with existing approaches, the presented scheme makes the first attempt to evaluate the nearly singular integrals arising in the regularized indirect BEM formulations when the high-order quadrilateral surface elements are used. The proposed scheme also can be applied to other problems in BEM, such as sensitivity analyses, contact problems, and thin-body problems. Some work for thin structures is already underway and will be reported in a subsequent paper. Acknowledgements
The support from the Opening Fund of the State Key Laboratory of Structural Analysis for Industrial Equipment (GZ1017) and the National Natural Science Foundation of Shangdong Province of China (ZR2010AZ003) are gratefully acknowledged. This article has been produced partially with the financial assistance of the European Regional Development Fund (ERDF) under the Operational Programme Research and Development/Measure 4.1 Support of networks of excellence in research and development as the pillars of regional development and support to international cooperation in the Bratislava region/ Project No. 26240120020 Building the centre of excellence for research and 14
development of structural composite materials – 2nd stage.
REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
[17]
[18]
M. Tanaka, V. Sladek, J. Sladek, Regularization techniques applied to boundary element methods, Applied Mechanics Reviews 47 (1994) 457-499. Y.J. Liu, On the simple solution and non-singular nature of the BIE/BEM-a review and some new results, Engineering Analysis with Boundary Elements 24 (2000) 286-292. J.T. Chen, Hong H.K, Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series, Applied Mechanics Reviews 52 (1999) 17-33. Y.J. Liu, T.J. Rudolphi, New identities for fundamental solutions and their applications to non-singular boundary element formulations, Comput Mech 24 (1999) 286-292. M. Guiggiani, G. Krishnasamy, T.J. Rudolphi, et al, A general algorithm for the numerical solution of hyper-singular boundary integral equations, (ASME) J. Appl. Mech, 59 (1992) 604-614. M. Guiggiani, A. Gigante, A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method, Journal of Applied Mechanics (ASME) 57 (1990) 906-915. V. Sladek, J. Sladek, Singular integrals and boundary elements, Computer Methods in Applied Mechanics and Engineering 157 (1998) 251-266. J.T. Chen, A.C. Wu, Null-field integral equation approach for multi-inclusion problem under antiplane shear, ASME, Journal of Applied Mechanics 74 (2007) 469-487. X.W. Gao, An effective method for numerical evaluation of general 2D and 3D high order singular boundary elements, Comput Methods Appl. Mech, Engrg 199 (2010) 2856-2864. G. Krishnasamy, L.W. Schmerr, T.J. Rudolphi, F.J. Rizzo, Hypersingular boundary integral equations: some applications in acoustic and elastic wave scattering, J Appl Mech 57 (1990) 404-414. Y.M. Zhang, V. Sladek, J. Sladek, etc, A new boundary integral equation formulation for plane orthotropic elastic media, Applied Mathematical Modelling 36 (2012) 4862-4875. T. Dirgantara, M.H. Aliabadi, Crack growth analysis of plates loaded by bending and tension using dual boundary element method, International Journal of Fracture 105 (2000) 27-74. M.H. Aliabadi, D. Martin, Boundary element hyper-singular formulation for elastoplastic contact problems, International Journal for Numerical Methods in Engineering 48 (2000) 995-1014. D. Zhang, F.J. Rizzo, Y.J. Rudolphi, Stress intensity sensitivities via hypersingular boundary element integral equations, Computational Mechanics 23 (1999) 389-396. J.F. Luo, Y.J. Liu, E.J. Berger, Analysis of two-dimensional thin structures (from micro- to nano-scales) using the boundary element method, Computational Mechanics 22(5) (1998) 404-412. Y.M. Zhang, Y. Gu, J.T. Chen, Boundary element analysis of 2D thin walled structures with high-order geometry elements using transformation, Engineering Analysis with Boundary Elements 35 (2011) 581-586. Y.J. Liu, Analysis of shell-like structures by the boundary element method based on 3-D elasticity: formulation and verification, International Journal for Numerical Methods in Engineering 41 (1998) 541-558. H.B. Chen, P. Lu, E. Schnack, Regularized algorithms for the calculation of values on and near boundaries in 2D elastic BEM, Engineering Analysis with Boundary Elements 25(10) (2001) 851-876. 15
[19] V. Sladek, J. Sladek, M. Tanaka, Regularization of hypersingular and nearly singular integrals in the potential theory and elasticity, International Journal for Numerical Methods in Engineering 36(10) (1993) 1609-1628. [20] C.Z. Cheng, Z.R. Niu, N. Recho, Z.Y. Yang, R.Y. Ge, A natural stress boundary integral equation for calculating the near boundary stress field, Computers Structures 89 (2011) 1449-1455. [21] S. Mukerjee, M.K. Chati, X.L. Shi, Evaluation of nearly singular integrals in boundary element contour and node methods for three-dimensional linear elasticity, International Journal of Solids and Structures 37(51) (2000) 7633-7654. [22] J.J. Granados, G. Gallego, Regularization of nearly hypersingular integrals in the boundary element method, Engineering Analysis with Boundary Elements 25(3) (2001) 165-184. [23] Y.C. Shiah, Y.S. Shih, Regularization of nearly singular integrals in the boundary element analysis for interior anisotropic thermal field near the boundary, Engineering Analysis with Boundary Element 30(2) (2007) 219-230. [24] Z.R. Niu, W.L. Wendland, X.X. Wang, H.L. Zhou, A semi-analytic algorithm for the evaluation of the nearly singular integrals in three-dimensional boundary element methods, Comput Methods Appl Mech Eng 194 (2005) 1057-1074. [25] H.L. Zhou, Z.R. Niu, C.Z. Cheng, Z.W. Guan, Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for anisotropic potential problems, Comput Struct 86 (2008) 1656-1671. [26] S. Nintcheu Fata, Semi-analytic treatment of nearly-singular Galerkin surface integrals, Applied Numerical Mathematics 60 (2010) 974-993. [27] Y.M. Zhang, Y. Gu, J.T. Chen, Stress analysis for multilayered coating systems using semi-analytical BEM with geometric non-linearities, Computational Mechanics 47 (2011) 493-504. [28] X.W. Gao, K. Yang, J. Wang, An adaptive element subdivision technique for evaluation of various 2D singular boundary integrals, Engineering Analysis with Boundary Elements 32 (2008) 692-696. [29] Telles, JCF, A self-adaptive coordinate transformation for efficient numerical evaluation of general boundary element integral, International Journal for Numerical Methods in Engineering 24 (1987) 959-973. [30] Letizia Scuderi, A new smoothing strategy for computing nearly singular integrals in 3D Galerkin BEM, Journal of Computational and Applied Mathematics 225 (2009) 406-427. [31] G.Z. Xie, F.L. Zhou, J.M. Zhang, X.S. Zheng, etc, New variable transformations for evaluating nearly singular integrals in 3D boundary element method, Engineering Analysis with Boundary Elements 37 (2013) 1169-1178. [32] M. Barbara, Johnston, P. R. Johnston, David Elliott, A new method for the numerical evaluation of nearly singular integrals on triangular elements in the 3D boundary element method, Journal of Computational and Applied Mathematics 245 (2013) 148-161. [33] M. Cerrolaza, E. Alarcon, A bi-cubic transformation for the numerical evaluation of the Cauchy principal value integrals in boundary methods, International Journal for Numerical Methods in Engineering 28(5) (1989) 987-999. [34] P.R. Johnston, Application of Sigmoidal transformations to weakly singular and near singular boundary element integrals, International Journal for Numerical Methods in Engineering 45(10) (1999) 1333-1348. [35] V. Sladek, J. Sladek, M. Tanaka, Optimal transformations of the integration variables in computation of singular integrals in BEM, International Journal for Numerical Methods in Engineering 47(7) 16
[36] [37] [38]
[39] [40]
[41]
[42]
[43]
[44]
[45] [46] [47]
[48] [49] [50]
(2000) 1263-1283. W.J. Ye, A new transformation technique for evaluating nearly singular integrals. Comput Mech 42 (2008) 457-466. Peter R. Johnston, David Elliott, A sinh transformation for evaluating nearly singular boundary element integrals, Int. J. Numer. Meth. Engng 62 (2005) 564-578. M. Barbara, Johnston, Peter R. Johnston, David Elliott, A sinh transformation for evaluating two-dimensional nearly singular boundary element integrals, Int. J. Numer. Meth. Engng 69 (2007) 1460-1479. Q. Huang, T.A. Cruse, Some notes on singular integral techniques in boundary element analysis, International Journal for Numerical Methods in Engineering 36(15) (1993) 2643- 2659. H. Ma, N. Kamiya, Distance transformation for the numerical evaluation of near singular boundary integrals with various kernels in boundary element method, Engineering Analysis with Boundary Elements 26(4) (2002) 329-339. H. Ma, N. Kamiya, A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two-and three-dimensional elasticity, Comput Mech 29 (2002) 277-288. Y. Miao, W. Li, J.H. Lv, X.H. Long, Distance transformation for the numerical evaluation of nearly singular integrals on triangular elements, Engineering Analysis with Boundary Elements 37 (2013) 1311-1317. X.Y. Qin, J.M. Zhang, G.X. Xie, etc, A general algorithm for the numerical evaluation of nearly singular integrals on 3D boundary element, Journal of Computational and Applied Mathematics 235 (2011) 4174-4186. Y.M. Zhang, C.L. Sun, A general algorithm for the numerical evaluation of nearly singular boundary integrals in the equivalent non-singular BIEs with indirect unknowns, Journal of the Chinese Institute of Engineers 31(3) (2008) 437-447. Y.M. Zhang, Y. Gu, J.T. Chen, Boundary layer effect in BEM with high order geometry elements using transformation, Computer Modeling in Engineering & Sciences 45(3) (2009) 227-247. Y.M. Zhang, Y. Gu, J.T. Chen, Boundary element analysis of the thermal behaviour in thin-coated cutting tools, Engineering Analysis with Boundary Elements 34(9) (2010) 775-784. Y.M. Zhang, Y. Gu, J.T. Chen, Internal stress analysis for single and multilayered coating systems using the boundary element method, Engineering Analysis with Boundary Elements 35(4) (2011) 708-717. Y.M. Zhang, W.Z. Qu, J.T. Chen, BEM analysis of thin structures for thermoelastic problems, Engineering Analysis with Boundary Elements 37 (2013) 441-452. X.Y. Qin, J.M. Zhang, G.Y. Li, et al., An element implementation of the boundary face method for 3D potential problems, Eng. Anal. Bound. Elem. 34 (2010) 934–943. Brebbia, C. A., Tells, J. C. F., and Wrobel, L. C. Boundary Element Technique, Springer-Verlag, Berlin, Germany(1984).
17
Figure
Figures：
Fig. 1. Eight-noded quadrilateral surface element, interior points and projection points
Example 2: Fig. 2. Discretization of a unit sphere with 100 second-order surface elements
(a)
(b)
Example 2: Fig. 3. Profiles of analytical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 2: Fig. 4. Surfaces of numerical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 2: Fig. 5. Surfaces of REs for the potentials (a) and its derivatives (b)
Example 3: Fig. 6. Discretization for one-eight of a hollow sphere with 15 second-order surface elements
Example 3: Fig. 7. AREs for the potentials u at interior points approaching to the boundary
Example 3: Fig. 8. AREs for the derivatives u x1 at interior points approaching to the boundary
Example 3: Fig. 9. Convergence curve of the computed u
Example 3: Fig. 10. Convergence curve of the computed u x1
Example 4: Fig. 11. Discretization of the torus with 80 second-order surface elements
Example 4: Fig. 12. AREs of the potentials
u
at interior points approaching to the boundary
Example 4: Fig. 13. AREs of the fluxes
(a)
u x1
at interior points approaching to the boundary
(b)
Example 4: Fig. 14. Profiles of the analytical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 4: Fig. 15. Profiles of numerical solutions for the potentials (a) and its derivatives (b)
(a)
(b)
Example 4: Fig. 16. Surfaces of REs for the potentials (a) and its derivatives (b) with eighty disretization boundary elements
(a)
(b)
Example 4: Fig. 17. Surfaces of REs for the potentials (a) and its derivatives (b) with hundred and sixty disretization boundary elements
Table
Tables：
Table 1 Coordinates ( x1j , x3j , x3j ) of eight nodes x j . Nodes x j
x1j
x2j
x3j
x1
0.304322331872978
0.936607830800249
0.173648177666930
x
2
-0.304322331872978
0.936607830800249
0.173648177666930
x
3
-0.304322331872978
0.936607830800249
-0.173648177666930
x
4
0.304322331872978
0.936607830800249
-0.173648177666930
x
5
6.0300091581E-017
0.984807753012208
0.173648177666930
x
6
-0.309016994374947
0.951056516295154
6.1230317691E-017
x
7
6.0300091581E-017
0.984807753012208
-0.173648177666930
x8
0.309016994374947
0.951056516295154
6.12303176911E-017
Table 2 Iterative results of projection points (1 ,2 ) with criterion ( k 1) ( k )
2
1 1011 .
Projection points (1 ,2 )
Interior points ( y1, y2 , y3 )
Iterations
-0.2290578 -0.2022083
0.9608331 0.9600097
0.0573038 -0.1289337
0.749627609474804 0.666556026866918
-0.33316956351993 0.749130717587247
5 5
0.1508432
0.9737583
-0.0859558
-0.500012936834761
0.499128077463858
4
-0.1528704
0.9767622
-0.0343823
0.499957919680965
0.199811834063831
4
0.1018787
0.9830678
0.0429779
-0.333392866230069
-0.249608365402485
3
-0.1018034
0.9823725
0.0573038
0.333394259894939
-0.332871412077530
3
-0.1274570
0.9808801
-0.0143202
0.416671464829870
8.32713312822E-002
4
0.1015781
0.9814146
0.1002760
-0.332996240934481
-0.581755480500003
4
-0.0760863
0.9715013
0.0992632
0.252591808052607
-0.582748457267639
4
-0.0509764
0.9875403
-0.028657
0.166712963275079
0.166626480502265
3
Table 3 Relative errors of various integrals with kernel
u * on a spherical surface element.
r0 a1 2
101
102
103
104
105
106
Reference solution
0.0077988
0.0079031
0.0079136
0.0079147
0.0079148
0.0079148
Ref [43] : ( , )
6.42E-9
2.52E-7
1.99E-6
6.09E-6
7.44E-6
9.11E-6
Ref [43] : ( , )
7.16E-9
2.52E-7
1.99E-6
6.09E-6
7.44E-6
9.11E-6
Present method 10
7.2694E-9
2.5498E-7
1.8427E-6
5.0576E-6
4.8640E-6
4.5634E-6
Table 4 Relative errors of various integrals with kernel
q* on a spherical surface element.
r0 a1 2
101
102
103
104
105
106
Reference solution
0.2839962
0.2890126
0.2895177
0.2895683
0.2895733
0.2895783
Ref [43] : ( , )
1.81E-7
4.48E-6
1.59E-6
7.89E-5
1.77E-4
1.50E-4
Ref [43] : ( , )
1.76E-7
4.48E-6
1.59E-6
7.89E-5
1.77E-4
1.50E-4
Present method 10
1.6743E-7
4.1356E-6
1.5059E-6
7.7185E-5
2.6966E-4
2.1417E-4
Table 5 Potentials u at internal points increasingly close to the boundary. Present Distance d
Exact
No transform Numerical
Relative error
0.1 0.01
0.2203216E+01 0.2468210E+01
0.2203587E+01 0.2478618E+01
0.2201609E+01 0.2467192E+01
5.604263E-04 4.126834E-04
0.001
0.2495155E+01
0.2472128E+01
0.2494217E+01
3.761218E-04
0.0001
0.2497854E+01
0.2469416E+01
0.2496925E+01
3.720910E-04
0.00001
0.2498151E+01
0.2469117E+01
0.2497223E+01
3.716440E-04
0.000001
0.2498154E+01
0.2469086E+01
0.2497225E+01
3.716563E-04
0.0000001
0.2497769E+01
0.2469083E+01
0.2496801E+01
3.717048E-04
0.00000001
0.2498154E+01
0.2469083E+01
0.2497225E+01
3.718859E-04
0.000000001
0.2498154E+01
0.2497224E+01
0.2496801E+01
3.721940E-04
0.0000000001
0.2498154E+01
0.2469083E+01
0.2497223E+01
3.726790E-04
Table 6 Potential derivatives u x1 at internal points increasingly close to the boundary. Present Distance d
Exact
No transform Numerical
Relative error
0.1 0.01
0.2899256E+01 0.2989385E+01
0.2910363E+01 0.2507812E+01
0.2896660E+01 0.2997269E+01
8.955288E-04 2.637319E-03
0.001
0.2998385E+01
-0.2714809E+01
0.3008335E+01
3.318621E-03
0.0001
0.2999285E+01
-0.3296846E+01
0.3009453E+01
3.390136E-03
0.00001
0.2999375E+01
-0.3350076E+01
0.3009565E+01
3.397357E-03
0.000001
0.2999384E+01
-0.3355347E+01
0.3009575E+01
3.397986E-03
0.0000001
0.2999385E+01
-0.3355873E+01
0.3009568E+01
3.395174E-03
0.00000001
0.2999385E+01
-0.3355926E+01
0.3009454E+01
3.357243E-03
0.000000001
0.2999385E+01
-0.3355931E+01
0.3009376E+01
3.331104E-03
0.0000000001
0.2999385E+01
-0.3355932E+01
0.3009223E+01
3.280253E-03
Research highlights • We present a new indirect regularized boundary element formulation. •
A new distance formula for high-order geometry elements is constructed.
•
The exponential transformation is developed to remove the near singularities.
•
The efficiency and applicability of the scheme are verified by some examples.