An Approach Based on Generalized Functions to Regularize Divergent Integrals

Engineering Analysis with Boundary Elements 40 (2014) 162–180

Contents lists available at ScienceDirect

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

An Approach Based on Generalized Functions to Regularize Divergent Integrals V.V. Zozulya n Centro de Investigacion Cientifica de Yucatan A.C., Calle 43, No 130, Colonia: Chuburná de Hidalgo, C.P. 97200, Mérida, Yucatán, México

art ic l e i nf o

a b s t r a c t

Article history: Received 13 February 2013 Accepted 2 November 2013 Available online 20 January 2014

This article addresses weakly singular, hypersingular integrals, which arise when the boundary integral equation (BIE) methods are used for 3-D potential theory problem solutions. An approach based on the theory of distributions and the application of the second Green theorem has been explored for the calculation of such divergent integrals. The divergent integrals have been transformed to a form that allows easy and uniform calculation of weakly singular and hypersingular integrals. For flat boundary elements (BE), piecewise constants and piecewise linear approximations, only regular integrals over the contour of the BE have to be evaluated. Furthermore, all calculations can be done analytically, so no numerical integration is required. In the case of 3-D, rectangular and triangular BE have been considered. The behavior of divergent integrals has been studied in the context that the collocation point moves to the contour of the BE. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Weakly singular Singular Hypersingular integrals Regularization Boundary integral equations

1. Introduction In recent years, many publications have been devoted to discussing BIE methods and their application to science and engineering because the BIE and the discrete analogy boundary element method (BEM) are very powerful tools for solving the mathematical problems that arise in science and engineering [1,2,7,17]. A historical account of the main stages of BIE and BEM development can be found in [9]. One of the main problems of the BIE solution is the calculation of divergent integrals as the numerical methods developed for regular integrals cannot be used and special methods have to be applied. Many methods have been developed using a classical approach for calculating the divergent integrals of different types. Usually, divergent integrals of a different type require different methods for their mathematical interpretation and numerical calculation. An analysis of most known methods used for the treatment of the different divergent integrals can be found in several books [1,7,17,20,25,27,33] and review articles [8,18,22,26,34,40]. From these publications, it follows that a significant impact on the topic under consideration has been made by Cruse [12–15], Mukherjee [28–31], Atluri [10,11,32,35] and their coauthors. There are some additional references that may also be of interest [3,4]. This publication is not intended to be a review; therefore, we do not present a long list of publications on the topic and will not critically analyze the problem. Instead, we concentrate on problems directly related to the topic of this publication. We begin by making one historical remark: the method of regularization for singular and hypersingular integrals,

n

Tel.: þ 52 9999428330. E-mail address: [email protected]

0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.11.001

which was developed by Guiggiani and coauthors (see [16] and references there), is very popular in the BEM community. Roughly speaking, the method consists of the following: A singularity is extracted from the divergent integral and divided into several parts. One part is relatively simple and contains singularities, while the other parts are regular. The regular parts are calculated using established methods, and the singular part is usually known and can be calculated. Such an approach is not new; it has long been used for divergent series calculations [21] and for calculating divergent integrals. For example, Kantorovich in [24] calculated 1-D divergent integrals with different types of singularities, including hypersingular integrals, and then extracted the divergent parts and calculated them analytically. Michlin extended this approach for the n-dimensional case and the theoretical foundation of that method [26,27]. Since that time, such methods have been used for divergent integral treatment. For example, we used it to derive a solution for the elastodynamic contact problems for cracked bodies in [36]. The classic approach for treating divergent integrals has one significant disadvantage; divergent integrals with a different singularity need different definitions, different theoretical justifications and different methods to be calculated. For example, weakly singular integrals are considered improper integrals, singular integrals are considered in the sense of the Cauchy principal value (PV), and hypersingular integrals have been considered by Hadamard to be finite part integrals (FP). In modern mathematics, divergent integrals have a strong theoretical foundation based on the theory of generalized functions (distributions), which permits us to apply the same approach for divergent integrals as for different types of singularities. According to this theory, divergent integrals with any type of singularity can be considered to be functionals (generalized functions) defined in

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

163

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ x21 þ x22 þ x23 is the distance in Euclidian space, x1 and x2 are the coordinates of integration and function fðxÞis regular. As soon as the point that corresponds to the origin of the coordinated system belongs to Ω, the integral (1) is divergent. The type of singularity depend on the value of the indices i, j and k. Because of the aforementioned singularity integrals, (1) cannot be transformed using the regular Green’s theorem. Therefore, a generalized function approach will be used here. In fact, the proposed method can be applied toward the regularization of an even wider class of divergent integrals. To illustrate this concept and to demonstrate the power of the proposed methods, we consider BIE and corresponding divergent integrals that appear in the general elliptic boundary-value problem. Let us consider a homogeneous region in which 3-D Euclidean space R3 occupies volume V with a smooth boundary ∂V. The region V is an open bounded subset of the Euclidean space with a C 0;1 Lipschitzian regular boundary ∂V . In the region V, we consider vector functions uðxÞand bðxÞthat are subject to a system of second order elliptic partial differential equations expressed in general form as

special functional spaces and on specially defined test functions [6]. We have shown in several publications that the approach based on the theory of distributions has not only theoretical meaning but is also an important application for practical calculations of divergent integrals. In our previous publications [17–19,37–44], the approach based on the theory of distributions has been developed for the regularization of divergent integrals with different singularities that arise in BEM applications. We explore the approach presented in [5], which interprets definite integrals as distributions and applies them to the solution of problems of fracture mechanics in [36]. Then, this approach was further developed for the regularization of 2-D hypersingular integrals, which appear in static and dynamic problems of fracture mechanics in [43] and [44], respectively. Regularized formulae for different types of divergent integrals have been reported in [38,39,40–41]. Additional applications of the regularization method can be found in review articles [18,19,42] and in a book [17]. Further development of that approach and the application of the second Green’s theorems in the sense of the theory of distributions have been described in [37,38]. The case of piecewise linear approximation has been considered in [39,40], while regularized formulae obtained in [38,40,41] permit us to transform weakly divergent singular, singular and hypersingular integrals over any polygonal area into regular integrals over the area contour. The developed approach can be applied to the regularization of a wide class of divergent integral regularization. In addition to hypersingular integrals, it is also suitable for the regularization of a variety of divergent integrals and any polynomial approximation. In this paper, the approach based on the theory of distributions and Green’s second theorem is developed and applied to the regularization of the divergent integral that appears in 2-D and 3-D potential theory problems solved by BIE methods. Generalized Gauss-Ostrogradski and Green theorems, which are applicable for the case of singular functions, have been obtained using methods presented in [5]. Then, the generalized second Green theorem was used for the development of the regularized formulae for divergent integral calculation. A special case of 2-D divergent integral regularization has been assessed for weakly and hypersingular integrals, and regular formulae for their calculation have been obtained. The weakly singular and hypersingular integrals for piecewise constant approximation have been considered for arbitrary convex polygons and for piecewise linear approximations that have been considered for rectangular and triangular BE. Using regularized formulae obtained here, divergent integrals were calculated for circular, quadratic and triangular areas. It is important to mention that in all presented equations, all calculations can be done analytically and no numerical integration is needed. The behavior of the divergent integrals was studied in the context of the collocation point being situated inside and outside of the BE and also when moving to the contour of the BE; the results are illustrated by corresponding diagrams created with Mathematica software.

where uand b are vector-functions and Lis a matrix differential operator; they have the form        L11 ⋯ L1n   u1  f1        ⋮    ⋯ ⋮ ; u ¼  ⋮ ; f ¼  ⋮  L¼ ð3Þ        Ln1 ⋯ Lnn   un  fn 

2. Divergent integrals and boundary integral equations

‖uðxÞ‖ ¼ Oðr  1 Þ;

Divergent integrals occur in various mathematical and engineering applications and can be of various types of singularity that exhibit different behavior in the vicinity of singular points. In this study, we concentrate our efforts on calculating the divergent integrals that appear in BIE when the corresponding system of integral equations is solved numerically using BEM. Usually, changing coordinates corresponding to divergent integrals can be presented as J i;j k

Z

xi1 xj2 fðxÞ ¼ dS rk S

ð1Þ

LUu ¼ b

ð2Þ

The coefficients of the matrix differential operator have the form Llk ¼

∂ ∂ ∂ c þ blki þ alk ∂xj lkji ∂xi ∂xi

ð4Þ

Coefficients clkji , blki and alk can be constants or can depend on coordinates. If the region Vis finite, it is necessary to establish boundary conditions. We consider the mixed boundary conditions in the form uðxÞ ¼ /ðxÞ;

8 x A ∂V u ;

pðxÞ ¼ P UuðxÞ ¼ ψðxÞ;

8 x A ∂V p

ð5Þ

The boundary contains two parts: ∂V u and ∂V p such that ∂V u \ ∂V p ¼ ∅ and ∂V u [ ∂V p ¼ ∂V. On the part ∂Vu is prescribed an unknown function uðxÞ, and on the part ∂V p is prescribed its generalized normal derivative pðxÞ. The generalized normal derivative is defined by the matrix differential operator with coefficients P lk ¼ nj clkji

∂ ∂xi

ð6Þ

here, ni are components of the outward normal vector to the surface ∂V p . If the region Vis infinite, then the solution of eq. (2) instead of the boundary conditions must satisfy additional conditions at the infinity in the form ‖P UuðxÞ‖ ¼ Oðr  2 Þ

for

r-1

ð7Þ

where ris the distance in the Euclidian space. According to the generalized second Green’s theorem Z Z ðun U L U u u U Ln Uun ÞdV ¼ ðun UP U u  u UPn Uun ÞdS

ð8Þ

we can obtain the following integral identity Z Z Z u U Ln Uun dV ¼ ðu UPn Uun  un U P UuÞdS  un U bdV

ð9Þ

∂V

V

V

∂V

V

where Ln is the operator adjoined to L. In the case of the scalar Poisson equation [5,7,39] and the system of the Lame linear equations of elasticity [1,17,40],

164

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

boundary, respectively, the BIE will have the form

we have ∂ L ¼ Ln ¼ Δ; P ¼ P n ¼ ni ∂xi

and

Llk ¼ Lnlk ¼ cljki ∂j ∂i ; P ij ¼ P nlk ¼ nk cikjl ∂l

ð10Þ respectively. Here, cikjl are elastic moduli. Eq. (9) is usually used to obtain integral representations for the solution of the boundary-value problems (2)–(5). To find those integral representations, we have to find the fundamental solution for the adjoined operator Ln Ln UU ¼ δ

ð11Þ

where U is a matrix of fundamental solution and δis the matrix delta function. Fundamental solutions for the system of partial differential equations with constant coefficients can be calculated analytically. If the coefficients are not constant, then only in special circumstances can fundamental solutions be derived. Analytical expressions for different types of differential operators and systems operators with applications in science and engineering have been described in numerous publications. In the case of the 3-D scalar Poisson equation, the fundamental solutions are expressed as 1 ; ð12Þ 4πr qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ ðx1  y1 Þ2 þ ðx2  y2 Þ2 þ ðx3 y3 Þ2 , which is the distance between points x and y in 3-D Euclidean space. Substituting a fundamental solution into (9), we obtain the integral representation for the solution of the general system of second order elliptic partial differential eq. (2) and its generalized normal derivative Uðx  yÞ ¼

uðyÞ ¼ Uðp; x; ∂VÞ  Wðu; x; ∂VÞ þ Uðf; x; VÞ pðyÞ ¼ Kðp; x; ∂V Þ  Fðu; x; ∂V Þ þ Kðf; x; VÞ

ð13Þ

For the case of the 3-D scalar Poisson equation, the fundamental solutions have the form n ðxÞðxi  yi Þ n ðyÞðxi  yi Þ Wðx; yÞ ¼  i ; Kðx; yÞ ¼ i ; 4πr 2 4πr 2   n ðxÞn ðyÞðx  y Þðx y Þ 1 i j i j i j 3  n ðxÞn ðyÞ Fðx; yÞ ¼ i i 4πr 3 r2

ð14Þ

and ni ðxÞ; ni ðyÞ are normal to the boundary ∂V vectors at points xand y, respectively. The corresponding fundamental solutions for the system of linear equations of elasticity have been described [1,2,17,40], etc. To be succinct, we introduce here the following notations for integral operators Z Z Uðf; x; V Þ ¼ fðxÞ U Uðx  yÞdV; Kðf; x; VÞ ¼ fðxÞ UKðx  yÞdV; V V Z Z pðxÞ U Uðx  yÞdS; Wðu; x; ∂VÞ ¼ uðxÞWðx; yÞdS; Uðp; x; ∂V Þ ¼ Z ∂V Z∂V pðxÞ UKðx; yÞdS; Fðu; x; ∂VÞ ¼ uðxÞFðx; yÞdS; Kðp; x; ∂VÞ ¼ ∂V

∂V

ð15Þ It is well known that boundary integral operators map between the following functional spaces U:

H  1=2 ð∂V Þ-H  1=2 ð∂VÞ;

W:

K:

H  1=2 ð∂VÞ-H  1=2 ð∂VÞ;

F:

H  1=2 ð∂VÞ-H  1=2 ð∂VÞ; H  1=2 ð∂VÞ-H  1=2 ð∂V Þ

ð16Þ

Together with boundary conditions (5), integral representations (13) are used to compose the BIE for the general boundary value problem (2), (5). Using (13) and boundary properties of the integral operators (15), we can construct various BIEs for the BV problem (2), (5). For example, if the first and second boundary integral representations are used on the parts ∂V u and ∂V p of the

1 2 uðxÞ Uðp; x; ∂V u Þ þ Wðu; x; ∂V p Þ ¼ ΦðxÞ; 1 2 pðxÞ  Kðp; x; ∂V u Þ þ Fðu; x; ∂V p ÞÞ

¼ ΨðxÞ;

8 x A ∂V u 8 x A ∂V p

ð17Þ

where ΦðxÞ ¼ Uðf; x; VÞ þUðψ; x; ∂V p Þ  Wð/; x; ∂V u Þ; ΨðxÞ ¼ Kðf; x; VÞ þKðψ; x; ∂V p Þ  Fð/; x; ∂V u Þ

ð18Þ

More possibilities for the creation of different types of the BIE have been considered in numerous books devoted to the BIE and BEM and their applications [1,2,7.17]. The kernels in (15) contain different types of singularities, so the corresponding integrals are divergent. Here, we will investigate these singularities and develop methods for the calculation of divergent integrals.

3. Classical definition of the divergent integrals Simple analysis of the fundamental solutions (12), (14) and corresponding fundamental solutions for elastostatics and general elliptic system (2) show that when x-y, U ij ðx  yÞ-r 1 ; W ij ðx; yÞ-r 2 ; K ij ðx; yÞ-r 2 ; F ij ðx; yÞ-r  3 ð19Þ Following Michlin, a definition and classification of the integrals with various singularities can be presented as follows. Definition 3.1. Let us consider two points with coordinates x; y A Rn (where n ¼ 3 or n¼2) and region V with smooth boundary ∂V of the classC 0;1 . The boundary integrals of type Z Gðx; yÞ fðxÞdS; α 40 ð20Þ rα ∂V where Gðx; yÞ is a finite function in Rn  ∂Vand fðxÞ is a finite function in ∂V are weakly singular for α o n  1, strongly singular for α ¼ n 1 and hypersingular for α 4 n 1. The integrals with singularities cannot be considered in the usual (Riemann or Lebegue) sense. For such integrals to make sense, special considerations must be made. Let us first consider the traditional approach to defining divergent integrals. Following [20,25,27], we define weakly singular integrals as improper, strongly singular in terms of Cauchy principal values and hypersingular in terms of Hadamard’s finite parts. Definition 3. 2. An integral with kernel U ij ðx  yÞ is weakly singular and must be considered as improper Z Z pi ðxÞU ij ðx  yÞdS ð21Þ W:S: pi ðxÞU ij ðx yÞdS ¼ lim ∂V

ε-0 ∂V n∂V ε

Here, ∂V ε is a part of the boundary, a projection of which on the tangential plane is contained in the circle C ε ðxÞ of the radio ε with its center at point x. Definition 3.3. Integrals with kernels W ij ðx; yÞ and K ij ðx; yÞ are singular and must be considered in the sense of the Cauchy principal values as Z Z ui ðxÞW ij ðx; yÞdS; P:V: ui ðxÞW ij ðx; yÞdS ¼ lim ∂V

Z P:V:

∂V

ε-0 ∂Vn∂Vðr o εÞ

Z pi ðxÞK ij ðx; yÞdS ¼ lim

ε-0 ∂Vn∂Vðr o εÞ

pi ðxÞK ij ðx; yÞdS

ð22Þ

Here, ∂Vðr o εÞ is a part of the boundary, a projection of which on the tangential plane is the circle C ε ðxÞ of the radio ε with its center at point x.

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

Definition 3.4. Integrals with kernels F ij ðx; yÞ are hypersingular and must be considered in sense of the Hadamard’s finite part as Z

Z

F:P: ∂V

ui ðxÞF ij ðx  yÞdS ¼ lim

ε-0

∂Vn∂V ðr o εÞ

ui ðxÞW ij ðx  yÞdS þ 2ui ðxÞ

f i ðxÞ ∂Vðr o εÞ



ð23Þ Here, functions f j ðxÞ are chosen from the condition of the limit existence. From these definitions, it follows that classical approaches require special consideration and special techniques for calculating each type of divergent integral. 4. Divergent integrals and generalized functions The methods of divergent integral regularization developed in our previous publications are based on the method of integration by parts and the theorems of Gauss-Ostrogradski and Green. In practice, their mathematically correct application implies some smoothness properties of the integrands. Unfortunately, divergent integrals do not possess such properties. Therefore, this approach is not suitable for divergent integral regularization in the frame of classical analysis using traditional approaches. Therefore, all of our publications drew attention to this situation and used generalized function theory. Generalized functions are widely used in modern analysis and applications, especially in theoretical and numerical studies of differential equations. Well-known generalized derivatives include Sobolev’s spaces, Dirac’s delta function, and others. The interpretation of divergent integrals as generalized functions is an uncommon approach, especially in the BEM community. Therefore, in this section, we consider the interpretation of divergent integrals as distributions and the application of the integral theorems of classical analysis for their regularization in the frame of a generalized function approach. We first introduce the basic concept of generalized functions or distributions and explain how they can be used for the regularization of divergent integrals. Instead of rigorous deductions and proofs for the corresponding theorems, we will concentrate on the central ideas and physical interpretation that have been used for the derivation of the formulae for the regularization of divergent integrals, directing the reader to special references for further details. To make the situation more clear, let us start from the 1-D case.

165

singularities and manipulation with concentrated loads. There are several mathematical definitions of generalized functions (see [5,6] for details). Here, we consider generalized functions to be linear continuous functionals defined on so-called space of test functions. In order to define correctly generalized function let us following [5] some additional mathematical constructions. Definition 4.1. Linear functional space K consists of real-value functions fðxÞ A C 1 , which go to zero together with all its derivatives faster than any power of the function 1=jxj. Function fðxÞ is called the test function, and C 1 is a functional space continuous with all its derivative functions. Definition 4.2. A generalized function f ðxÞ is defined as any linear functional on functional space Kof the test functions fðxÞ A K, such as the following Z ðf ; fÞ ¼ f ðxÞfðxÞdx ð26Þ R1

where R1 ¼ ð  1; 1Þis one dimensional Euclidian space. This definition means the following: 1. with any fðxÞ A K, functional (26) is associated with a real number, 2. for any α1 and α2 and f1 ðxÞ and f2 ðxÞ A K, there exists the relation      f ðxÞ; α1 f1 ðxÞ þ α2 f2 ðxÞ ¼ α1 f ðxÞ; f1 ðxÞ þ α2 f ðxÞ; f2 ðxÞ 3. from the convergence lim fk ðxÞ ¼ fðxÞ follows the convergence  k-1 lim f ðxÞ; fk ðxÞ ¼ ðf ðxÞ; fðxÞÞ. k-1

It is important to mention that one cannot speak of the value of a generalized function at a given point because they are defined as functionals on the whole space. For instance, one cannot, say that “a generalized function is equal to zero at point x0 ”. However, the statement that “a generalized function f ðxÞ is equal to zero in a neighborhood of x0 ” can be given with the following meaning: for every fðxÞ A K that vanish in a neighborhood of x0 , we have ðf ðxÞ; fðxÞÞ ¼ 0. The next important step is a definition of the derivative of a generalized function.

4.1. 1-D divergent integrals Generalized functions or distributions were introduced in order to extend and validate operation of differentiation to nonsmooth and even discontinuous functions. Let us consider function of one variable f ðxÞ ¼ x log ðxÞ x

ð24Þ

Its derivatives df ðxÞ ¼ log ðxÞ; dx

2

d f ðxÞ 1 ¼ ; x dx2

Definition 4.3. The derivative of the generalized function f ðxÞ is defined by the equation     df ðxÞ df ; fðxÞ ¼ f ðxÞ; ðxÞ ; 8 fðxÞ A K ð27Þ Dx f ðxÞ ¼ dx dx

3

n

d f ðxÞ 1 d f ðxÞ 1 ¼  2 ; ::: ; ¼ ð  1Þn n  1 ; dxn dx3 x x

n ¼ 2; 3; :::; 1:

ð25Þ have isolated singularities at point x ¼ 0. At this singular point, the differentiation and interpretation of these expressions cannot be achieved in the classical sense. However, such functions often appear in analytical science and engineering applications. The classical approach consists of considering the function and its derivatives on two intervals ð  1; 0Þ and ð0; 1Þ separately with some limit transitions for x- 7 0. In this way, the mechanics appear concentrated on loads such as forces and moments of different order. Generalized functions were invented to achieve accurate mathematical calculation of the derivatives of the functions with

The derivative of order kof the generalized function f ðxÞ is defined by ! ! k k d f ðxÞ d fðxÞ k k Dx f ðxÞ ¼ ; fðxÞ ¼ ð  1Þ f ðxÞ; ; 8 fðxÞ A K ð28Þ dxk dxk For a regular function, the ordinary derivative is equal to its generalized derivative. For our purposes, we do not need to study a generalized function in a general sense. The situation is greatly simplified if we consider generalized functions that can be presented in the form r

f ðxÞ ¼

d gðxÞ dxr

ð29Þ

where gðxÞis any continuous – or merely piecewise continuous – function; this is referred as a generating function.

166

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

Formula (28) for the derivative of order k of the generalized function f ðxÞ can be presented in the form r þk

Dkx f ðxÞ ¼

d

gðxÞ dxr þ k

ð30Þ

An antiderivative or indefinite integral of the generalized function f ðxÞ of the form (29) can be defined as Z r1 d gðxÞ ð31Þ f ðxÞdx ¼ Dx 1 f ðxÞ ¼ dxr  1 In the case r ¼ 1, Z f ðxÞdx ¼ Dx 1 f ðxÞ ¼ gðxÞ

ð32Þ

function, is regular and possesses all necessary derivatives. Therefore, function f ðxÞ is a generalized function. Let as consider a definite integral of the function f ðxÞ over finite interval Ω Z a I0 ¼ f ðxÞ dx ð38Þ a

What does the I 0 symbol means for such a singular function? The classical approach cannot be used to answer this question generally because the answer exists only for a special type of singularity and each type of singularity requires specific consideration. We introduce a finite test function fðxÞ A C 1 ðRÞ, such that fðxÞ ¼ 1 and 8 x A Ω, extended arbitrarily to the region Ω0 . Obviously, the derivatives for the defined function fðxÞ are equal to zero in the region Ω, including end points ∂Ω ¼ f a; ag. k

For our purposes, it is important to consider generalized functions of type (30), which are ordinary functions everywhere except in a subdomain Ωε of a larger domain Ω ¼ ½a; b such that outside of Ωε , the generating function gðxÞ possesses continuous derivatives. Frequently, the subdomain Ωε consists of isolated points. For example, if function f ðxÞ has piece-wise continuity, its generalized derivative is equal to Dx f ðxÞ ¼

m df ðxÞ þ ∑ ½½f ðxÞi δðx  xi Þ dx i¼1

ð33Þ

where ½½f ðxÞi is a jump of the function at point xi and δðx xi Þ is a Dirac’s delta function. Let us consider weakly singular, singular, and hypersingular integrals, which usually appear in the numerical solution of the boundary integral equations Z 1 Z 1 Z 1 1 1 dx; log j 1=xj dx; dx ð34Þ 2 1 1 x 1 x According to the traditional approach, each type of these divergent integrals has to be considered separately: weakly singular as improper (21) ! Z Z Z 1

W:S: 1

1

log j 1=xj dx ¼ lim ε-0

ε

log j 1=xj dx 

ε

log j 1=xj dx

¼ 2;

1

ð35Þ singular in the sense of Cauchy’s principal value (22) ! Z ε Z 1 Z 1 1 1 1 dx ¼ lim dx  log dx ¼ 0; P:V : x ε-0 1 x ε x 1

d fðxÞ ¼ 0; dxk

x A Ω ¼ ½  a; a

ð39Þ

Now, let us consider the scalar product, which is functional, and define function f ðxÞ in the sense of distributions Z Z k d gðxÞ ðf ; fÞ ¼ f ðxÞfðxÞ dx ¼ fðxÞ dx ð40Þ dxk R R Integrating by parts, taking into account properties (39) and the finiteness of the test function, we obtain Z Z Z k k k d gðxÞ d fðxÞ d fðxÞ k k fðxÞ dx ¼ ð  1Þ gðxÞ dx ¼ ð  1Þ gðxÞ dx k k 0 0 dx dx dxk R Ω[Ω Ω ð41Þ The last equality occurs because fðxÞ ¼ 1 and (39) in the region Ω. By integrating the last integral by the path again in inverse order and by taking into account that fð aÞ ¼ fðaÞ ¼ 1, we obtain x ¼ a Z Z k i1 k d fðxÞ d gðxÞ d gðxÞ k gðxÞ dx ¼ þ ð 1Þ fðxÞ dx ð42Þ  dxi  1  dxk dxk Ω0 Ω0 x ¼ a

Integral (38) in the sense of distribution can be written as Z Z f ðxÞfðxÞ dx ð43Þ I 0 ¼ f ðxÞfðxÞ dx  Ω0

R

The first term here is ð36Þ

and hypersingular in the sense of Hadamard’s finite part (23) ! Z 1 Z ε Z 1 1 1 1 F:P: dx ¼ lim dx  log dx ¼  2: ð37Þ 2 2 ε-0 x2 1 x ε x 1 For further details, see our previous publications. All of the above divergent integrals and those belonging to a much wider class of divergent integrals can be considered as generalized functions or, more specifically, as a linear functional defined on the corresponding test functions. To do this, we must introduce an accurate definition of the finite integral of the generalized function. As shown above, the indefinite integral of the generalized function can be represented in the form of (31) or (32). Now following [5,6], we consider definite integrals of the generalized function, which are very important for the derivation of the regular representation of the divergent integrals. Let f ðxÞ be a function of one variable defined in the region x A Ω ¼ ½  a; a. All singularities of the function f ðxÞ are concentrated in sub region Ωε ¼ ½ ε; ε  Ω. The region Ω n Ωε , including the boundary

x ¼ a Z k gðxÞ d gðxÞ f ðxÞfðxÞdx ¼ þ fðxÞ dx  i1  0 dx dxk R Ω

Z

i1

d

ð44Þ

x ¼ a

and the second term is Z Z k d gðxÞ f ðxÞfðxÞ dx ¼ fðxÞ dx dxk Ω0 Ω0

ð45Þ

As result of (43)–(45), we find the finite part of the divergent integral in the form: x ¼ a Z a k1 d gðxÞ I 0 ¼ F:P: f ðxÞ dx ¼ ð46Þ  dxk  1  a x ¼ a

This is the definition of the divergent integrals in the sense of distributions. We can use this equation to calculate weakly singular, singular and hypersingular integrals in the same way. Obviously, for r ¼ 1, we have Z a x ¼ a F:P: f ðxÞ dx ¼ gðxÞx ¼  a ð47Þ a

For regular functions this is an usual formula from integral calculus for calculation of the definite integrals, it usually refer as Newton-Leibnitz formula.

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

Thus, formulae (46) and (47), the generalized Newton-Leibnitz formula for the above-defined singular functions, provide a powerful tool for their calculation. Using this formula, one can easily calculate the 1-D divergent integrals, which can also be calculated using standard techniques, and the results can be compared. We start with a hypersingular function. Evidently, function f ðxÞ ¼ 1=ðy  xÞ2 can be represented in the form f ðxÞ ¼ ðd=dxÞð1=ðy  xÞÞ and gðxÞ ¼  ð1=ðy xÞÞ; then, from (47), it follows that  Z a Z a dx d 1 1 x ¼ a dx ¼  F:P: ¼  F:P: 2 x  yx ¼  a  a ðx yÞ  a dxðy  xÞ 1 1  ¼ ða þyÞ ða  yÞ

ð  a o y oaÞ

ð48Þ

  If we consider the functions f ðxÞ ¼ ðy xÞ  1 and gðxÞ ¼ lny  x, the singular integral, which is considered in the sense of Cauchy PV, follows from (47),   Z a Z a a x dy d lnj y xj  ð  a o x o aÞ ¼ P:V: dx ¼ ln P:V : dx a þx  a ðy xÞ a ð49Þ   If we consider the functions f ðxÞ ¼ lny  x and gðxÞ ¼  ðy  xÞ  lny  x  ðy xÞ, the WS integral, which is considered as an improper integral, follows from (47), Z a W:S: lnj y  xj dy ¼ ða  xÞlnj a xj þ ða þ xÞlnj aþ xj ; ð  a o xo aÞ a

Integrating the integral over domain Ω by the path in inverse order, we obtain Z Z dfðxÞ dfðxÞ a dx ¼ fðxÞfðxÞgðxÞj xx ¼ dx gðxÞfðxÞ fðxÞgðxÞ ¼ a  dx dx Ω0 Ω0 Z dgðxÞ dx ð56Þ  fðxÞfðxÞ 0 dx Ω Substituting into (43) Z Z dfðxÞ a dx f ðxÞfðxÞ dx ¼ fðxÞfðxÞgðxÞj xx ¼ fðxÞgðxÞ ¼ a  dx R Ω Z dgðxÞ dx  fðxÞfðxÞ 0 dx Ω and Z Z dgðxÞ dx f ðxÞfðxÞ dx ¼ fðxÞfðxÞ dx Ω0 Ω0

Z

1

dx ¼ 2 2 1 x

ð51Þ

Obviously, this result is consistent with the one obtained using the traditional approach (37). The presented approach can be applied and functions of type k

f ðxÞ ¼ fðxÞ

d gðxÞ dxk

ð52Þ

can be regularized using a distribution based method. Here, fðxÞ A C k ðΩÞ is a regular k-differentiable function. First, we consider simply the case of k ¼ 1 and functions of type f ðxÞ ¼ fðxÞ

dgðxÞ dx

x ¼ a

Z

Integrating by parts, we obtain Z

dgðxÞ dx ¼ fðxÞfðxÞ dx R

Z Ω [ Ω0

¼ 

Z

fðxÞfðxÞ

dgðxÞ dx ¼  dx

dfðxÞ dx  gðxÞfðxÞ dx Ω0

Z

Z

gðxÞ Ω [ Ω0

dfðxÞfðxÞ dx dx

gðxÞfðxÞ Ω [ Ω0

dfðxÞ dx dx ð55Þ

k

a

d fðxÞ gðxÞ dx dxk a

ð60Þ

With this formula, a divergent integral with any type of singularity can be represented as a sum of the boundary term and a regular integral. For BEM applications, it is very important to consider functions with singularities of type f ðxÞ ¼

1 rm

ð61Þ

This case has already been analyzed using the generalized function approach in [37]. The regularization formula for this specific case of singularities has the form x ¼ a Z a i k1i k1 fðxÞ Pi d fðxÞ iþ1 d F:P: dx ¼ ∑ ð  1Þ  m dxi r m  k dxk  1  i  a r i¼0 x ¼ a

ð53Þ

where fðxÞ A CðΩÞ is a regular differentiable function. All other assumptions are the same as in the previous case. In the same way as in (40), we consider the functional in the sense of distributions Z Z dgðxÞ ðf ; fÞ ¼ f ðxÞfðxÞ dx ¼ fðxÞfðxÞ dx ð54Þ dx R R

ð58Þ

To consider the regularization of the divergent integrals with a singular function of type (52), we can apply integration by parts k times or directly use the formula x ¼ a Z a i k1i k1 fðxÞ i þ 1 d gðxÞd I 0 ¼ F:P: f ðxÞ dx ¼ ∑ ð  1Þ  dxi dxk  1  i  a i¼0 þ ð 1Þk

F:P:

ð57Þ

we obtain the extension for the case of singular function formula for integration by parts in the form Z a Z dfðxÞ ¼a dx ð59Þ I 0 ¼ F:P: f ðxÞ dx ¼ fðxÞgðxÞj xx ¼ gðxÞ a  dx Ω a

ð50Þ A famous example of a Hadamard positive function with a negative definite integral comes from (48)

167

Z þ ð  1Þk

a a

k

P k d fðxÞ r m  k dxk

ð62Þ

1 1 where gðxÞ ¼ rmPk k and P k ¼ ð  1Þk ∏ki ¼ 0 ðm þ iÞ for k; m 4 1. It is important to mention that eq. (62) is very important for the regularization of a wide class of divergent integrals that appear in the numerical solution of the boundary integral equations. It can be used for regularization divergent integrals with any singularity of type (61) for plat and curvilinear boundary elements. Indeed, function fðxÞ can contain information about the interpolation polynomial and shape function of the curvilinear element.

4.2. n-D divergent integrals Now, let us extend the above consideration to a general n-D case. Definitions 1 and 2 are extended for the n-D case, so we have to replace R1 withRn and the function of one real variable f ðxÞ with the function of n variablesf ðxÞ. In the n-D case, we have to consider

168

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

partial derivatives of type Drx ¼

r 1 þ r 2 þ :::r n

∂ ∂xr11 ∂xr22 :::∂xrnn

ð63Þ

Rn

where r ¼ r 1 þ r 2 þ ::: þ r n is the order of the partial derivative. Definition 4.4. The partial derivative of order r of the generalized function f ðxÞ is defined by equation

Drx f ðxÞ ¼

∂r1 þ r2 þ :::rn f ðxÞ ; fðxÞ ∂xr11 ∂xr22 :::∂xrnn

! ¼ ð  1Þr f ðxÞ;

! ∂r1 þ r2 þ :::rn fðxÞ ; ∂xr11 ∂xr22 :::∂xrnn

Let us consider the scalar product, which is a definition of the singular function f ðxÞ in the sense of distributions, Z ðf ; fÞ ¼ f ðxÞfðxÞ dx ð70Þ

8 fðxÞ A K

ð64Þ

Taking into account that the test function is finite, Z Z fðxÞΔk gðxÞ dx ¼ fðxÞΔk gðxÞ dΩ Rn

Because the derivatives of the test function fðxÞ equal zero on Ω, the integration by parts using the second Green theorem in the sense of distribution gives Z

For generating functions of type (29) in the n-D case, there are more possibilities. In the general case, they can be presented in the form f ðxÞ ¼

∂r1 þ r2 þ :::rn gðxÞ ∂xr11 ∂xr22 :::∂xrnn

ð65Þ

f ðxÞ ¼ ∇gðxÞ

and

f ðxÞ ¼ ΔgðxÞ

Ω

As mentioned before, the classical approach cannot provide an accurate meaning for the integral I 0 . To consider this integral in the sense of the distribution, we introduce the function gðxÞ ð68Þ

where Δk ¼ Δ U Δ UΔ U ::: U Δ Δk is called the k – dimensional Laplace |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} operator. k This representation of the function f ðxÞ can be considered in the classical sense in the region Ω0 , but in the region Ω, it has to be considered in the sense of distribution. In the region Ω0 ¼ Rn =Ω, the function is regular and smooth up to the boundary, meaning that f ðxÞ A C k ðΩ0 Þ. The boundary ∂Ω satisfies the usual conditions for smoothness, which are discussed in the course of any standard analysis. We also introduce test function fðxÞ A C 1 ðRn Þ, such that fðxÞ ¼ 1 and 8 x A Ω. Function fðxÞ is finite and extends smoothly to the region Ω0 . In this case, its derivatives are equal to zero at the region Ω, including its boundary ∂Ω, i.e., Δk fðxÞ ¼ 0;

xAΩ

Ω [ Ω0

fðxÞΔk gðxÞ dx ¼ ð  1Þk

ð69Þ

Z Ω [ Ω0

gðxÞΔk fðxÞ dΩ ¼ ð  1Þk

Ω0

gðxÞΔk fðxÞ dV

ð72Þ The integration by parts in reverse order for the last integrals above leads to the following result

Ω

k1

0

fðxÞΔk gðxÞ dΩ ¼ ∑ ð  1Þi þ 1 i¼0

Z ∂Ω0

½fðxÞ∂n Δk  i  1 gðxÞ

 gðxÞ∂n Δk  i  1 fðxÞ dS þ ð 1Þk

ð66Þ

where ∇ ¼ ð∂=∂x1 Þ þ ð∂=∂x2 Þ þ⋯ þ ð∂=∂xn Þ and Δ ¼ ð∂2 =∂x21 Þ þ ð∂2 = ∂x22 Þ þ ⋯ þ ð∂2 =∂x2n Þ are Hamilton’s and Laplace’s operators. In the same way as in the 1-D case, in the n-D case, weakly singular, singular and hypersingular integrals have to be considered separately following definitions (21)–(23). However, such divergent integrals and many wide classes of divergent integrals can be considered as generalized functions and, more specifically, as linear functionals defined on the corresponding test functions. To do this, we have to introduce a definition for the finite integral of the generalized function in the n-D case. Following [5,6], we consider definite integrals of the generalized function in the n-D case. Let us consider function f ðxÞ, which is defined in the finite region with Ω  Rn such that all its singularities are concentrated in the sub region Ωε  Ω. In the region Ω n Ωε , the included boundary function is regular and possesses all necessary derivatives. Let us consider a definite integral over a finite region Z I 0 ¼ f ðxÞ dx ð67Þ

f ðxÞ ¼ Δk gðxÞ

Z

Z

For our purpose, a narrower class of the generating functions is preferable. For example:

ð71Þ

Ω [ Ω0

Z Ω0

gðxÞΔk fðxÞ dΩ ð73Þ

Here, ∂n ¼ ni ∂i is the normal derivative on the surface with respect to x and ni ðxÞ is a unit normal to the surface. Taking into account that Z Z Z f ðxÞfðxÞ dΩ ¼ f ðxÞfðxÞ dΩ f ðxÞfðxÞ dΩ ð74Þ Rn

Ω [ Ω0

Ω0

and considering (72) and (73), we find the formula for calculation of the divergent integral for functions of type (68) with any singularity: Z Z Z F:P: f ðxÞ dV ¼ F:P: Δk gðxÞ dV ¼ ∂n Δk  1 gðxÞ dS ð75Þ Ω

Ω

∂Ω

For example, weakly singular, singular and hypersingular integrals of any dimension can be calculated in the same way using formula (75). In the specific case of k ¼ 1, we get the generalization of the well-known Gauss-Ostrogradski theorem for the case of functions of type (68) with any singularity. Z Z Z F:P: f ðxÞ dΩ ¼ F:P: ΔgðxÞ dΩ ∂n gðxÞ dS ð76Þ Ω

Ω

∂Ω

Thus, formulae (75) and (76), the generalized Gauss-Ostrogradski theorem for the singular functions defined above, are powerful tools for their calculation. Using this formula, one can easily calculate the 2-D divergent integrals, which can also be calculated using standard techniques to compare the results. Let us consider integrals of type Z 1 Ik ¼ dΩ; k 4 0; k a 2 ð77Þ k Ω r The integrals of this type may be regularized using the OstrogradskiGauss theorem (76). In this case, f ðxÞ ¼ r1k , gðxÞ ¼ ðk  2Þ12 rk  2 and Z Z Z dS 1 1 1 rn ¼ ∂ dS ¼  dS ð78Þ I k ¼ F:P: n k ðk 2Þ ∂Ω r k ðk  2Þ2 ∂Ω r k  2 Ωr Here, r n ¼ ðxα yα Þnα and α ¼ 1; 2.

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

This integral can be calculated analytically over a circular area. Introducing polar coordinates yields Z Z 2π 1 rn 1 r 1 2π Ik ¼  dS ¼  r df ¼  ð79Þ k ðk  2Þ ∂Ω r ðk  2Þ 0 r k ðk  2Þ r k  2 For specific cases of k ¼ 1 and k ¼ 3, we get the well-known results I 1 ¼ 2πrand I 3 ¼ 2πr  1 , respectively. In another example, we present integrals of type Z x1 x2 I 1;1 ¼ dΩ; k 4 3; 5 ð80Þ k k Ω r x1 x2 1 . The regularizaIn this case f ðxÞ ¼ x1rkx2 and ΔgðxÞ ¼ ðk  2Þðk  6Þ r k  2 tion formula in this case has the form  Z  Z x1 x2 1 rn x1 x2 r n I 1;1 dS ð81Þ ¼ dΩ ¼  k k ðk 2Þðk  6Þ ∂Ω r k  2 rk Ω r

where r n ¼ x1 n2 þ x2 n1 . Due to symmetry, these integrals are equal to zero over a circular area, as was previously shown [39,41]. We can extend the presented approach and regularize functions of type k

f ðxÞ ¼ fðxÞΔ gðxÞ

169

theorem k times. In this case, the following formula is obtained Z Z k1 I 0 ¼ F:P: f ðxÞ dΩ ¼ ∑ ð  1Þi þ 1 ðfðxÞ∂n gðxÞ Ω

i¼0

∂Ω

Z

 gðxÞ∂n fðxÞÞ dS þ ð  1Þk

Ω

gðxÞΔk fðxÞ dΩ

ð90Þ

With this formula, a divergent integral with any type of singularity can be represented as a sum of the boundary term and a regular integral. For BEM applications, it is very important to consider functions with singularities of type f ðxÞ ¼

1 rm

ð91Þ

This case has already been analyzed using the generalized function approach [37]. The regularization formula for this specific case of singularity has the form Z Z k1 fðxÞ P I 0 ¼ F:P: dV ¼ ∑ ð  1Þi þ 1 Δk  i  1 fðxÞ∂n m i 2i m r V r ∂V i¼0 

ð82Þ

Pi

r

∂ Δk  i  1 fðxÞ m  2i n

i

Z

dS þð  1Þk

1 Vr

m  2k

½Δk fðxÞ dV ð92Þ

First, we consider simply the case k ¼ 1 and functions of type f ðxÞ ¼ fðxÞΔgðxÞ

ð83Þ

where fðxÞ A CðΩÞ is a regular differentiable function. All other assumptions are the same as in the previous case. In the same way as in (70), we consider the functional in the sense of distributions Z Z ðf ; fÞ ¼ fðxÞf ðxÞfðxÞ dx ¼ fðxÞfðxÞΔgðxÞ dx ð84Þ Rn

Rn

Using the second Green theorem, we obtain

Z

Z

n

R

fðxÞfðxÞΔgðxÞ dx ¼

Z

Ω[Ω

0

fðxÞfðxÞΔgðxÞ dx ¼ 

Ω [ Ω0

∇gðxÞ∇ðfðxÞfðxÞÞ cdx

1 1 where gðxÞ ¼ rmPk 2k and P k ¼ ð  1Þk ∏ki ¼ 0 ðm þ 2iÞ2 for k; m 4 1. We have to mention that eq. (92) is very important for the regularization of a wide class of divergent integrals that appear in the numerical solution of the boundary integral equations. It can be used for the regularization of divergent integrals with any singularity of type (91) for plat and curvilinear boundary elements. Indeed, function fðxÞ can contain information about the interpolation polynomial and shape function of the curvilinear element. In the following sections, we will demonstrate an application of the developed approach to the regularization and calculation of the divergent integrals that appear in BEM.

ð85Þ 5. The BEM equations and approximation

The last integral can be represented in the following form Z

Z

Ω[Ω

0

∇gðxÞ∇ðfðxÞfðxÞÞ dx ¼

Z

Ω

0

fðxÞ∇gðxÞ∇fðxÞdx þ

Ω [ Ω0

fðxÞ∇gðxÞ∇fðxÞ dx

ð86Þ Integrating the integral over domain Ω0 by the path in the inverse order using the first Green theorem, we obtain Z Z Z fðxÞfðxÞΔgðxÞ dx ¼ fðxÞfðxÞ∂n gðxÞ dS  fðxÞ∇fðxÞ∇gðxÞ dx Ω0

∂Ω

Ω0

Z

þ Ω0

fðxÞ∇fðxÞ∇gðxÞ dx

ð87Þ

Substituting (85) and (87) into (43), we obtain an extension of the first Green theorem for the case of singular functions in the form Z Z Z fðxÞ∂n gðxÞ dS  ∇gðxÞ∇fðxÞ dΩ I 0 ¼ F:P: fðxÞΔgðxÞ dΩ ¼ Ω

∂Ω

Ω

ð88Þ In the same way as in classical analysis, one can find an extension of the second Green theorem for the case of singular functions in the form Z Z Z I 0 ¼ F:P: f ðxÞ dΩ ¼ ðfðxÞ∂n gðxÞ  gðxÞ∂n fðxÞÞ dS  gðxÞΔfðxÞ dΩ Ω

∂Ω

Ω

ð89Þ To consider the regularization of the divergent integrals with singular function of type (82), we can apply the second Green

The BIE are usually solved numerically, transforming them into a discrete system of finite dimensional equations called the BEM. To transform the BIE into finite dimensional BEM equations, we have to split the boundary ∂V into finite elements N

S ¼ [ Sn ; n¼1

Sn \ Sk ¼ ∅

for n a k:

ð93Þ

Because Sis the boundary of body, these elements are known as boundary elements (BE). On each BE, we shall choose Q nodes of interpolation and shape functions fq ðxÞ. The vectors of displacement and traction on the BE Sn will then be represented approximately in the form Q

ui ðxÞ  ∑ uqi fq ðxÞ; q¼1

Q

pi ðxÞ  ∑ pqi fq ðxÞ; q¼1

x A Sn

ð94Þ

and on the hold boundary ∂V in the form N

ui ðxÞ  ∑

Q

∑ uqi fq ðxÞ;

n¼1q¼1

N

pi ðxÞ  ∑

Q

∑ pqi fq ðxÞ;

n¼1q¼1

xAS

ð95Þ

Here and below, if the unit q belongs to several BEs, it is considered in these sums once. In these and subsequent equations, if the node of interpolation belongs to several boundary elements, it is accounted for only once. Substituting expressions (95) into (13) gives us the finitedimensional representations for the vectors of displacements

170

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

and traction on the boundary in the form 1 2uðyr Þ ¼



Z Z Q N ∑ ∑ pq ðxm Þ fq ðxÞUðyr ; xÞ dS  uq ðxm Þ fq ðxÞWðyr ; xÞ dS þ Uðf; y; V n Þ

n ¼ 1q ¼ 1

Sn

Sn



Z Z Q N 1 pðyr Þ ¼ ∑ ∑ pq ðxm Þ fq ðxÞKðyr ; xÞ dS  uq ðxm Þ fq ðxÞFðyr ; xÞ dS þ Kðf; y; V n Þ 2 Sn Sn n ¼ 1q ¼ 1

ð96Þ The volume potentials Uðf; y; V n Þ and Kðf; y; V n Þ depend on the discretization of the V domain. More detailed information about the transition from BIE to BEM equations can be found in [1,2,7]. Let us consider how to construct a BE model of the boundary S ¼ ∂V  R3 of the area V. We fix in the boundary S a finite number of points xq ðg¼ 1; :::; Q Þ; these points are referred to as global node points SðqÞ ¼ xq A V : g ¼ 1; :::; Q . We shall divide the area S into a finite number of sub-areas Sn ðn ¼ 1; :::; NÞ, which are BEs that must satisfy the following conditions Sn [ Sm ¼ ∅;

m a n;

m; n ¼ 1; 2; :::; N;

N

S ¼ [ Sn : n¼1

ð97Þ

We introduce a local coordinate system ξ on each BE. We designate the nodal points xq A Sn in the local system of coordinates as ξq . They are coordinates of nodal points in the local coordinate system. Local and global coordinates are related according to the following: N

xq ¼ ∑ Λn ξqn :

ð98Þ

n¼1

Functions Λn depend on the position of the nodal points in the BE. They join individual BEs together in a finite dimensional BE model. The borders of the BEs and the positions of the nodal points should be such that, after joining together, separate elements form a discrete model of the area S. Having constructed a BE model of the area S, we shall consider an approximation of the function f ðxÞthat belongs to some functional space. The BE model of the area S is the domain of the function that should be approximated. We denote function f ðxÞ on n the BE Sn by f ðxÞ. Then, N

n

f ðxÞ ¼ ∑ f ðxÞ:

ð99Þ

The differential elements along the coordinate axes are related by dxi ¼ ð∂xi =∂ξj Þ dξj ;

ð103Þ dξi ¼ ð∂ξi =∂xj Þ dxj ; dξ ¼ J  1 ðxÞ dx: The differential element of the length of a contour in ℜ2 is defined by the expression qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dl ¼ ðdx1 ðξÞ=dξÞ2 þ ðdx2 ðξÞ=dξÞ2 Þ dξ ð104Þ The differential element of area in ℜ2 is transformed via the formula   ∂xα  dA ¼ dx1 dx2 ¼ Det  dξ1 dξ2 ; α; β ¼ 1; 2 ð105Þ ∂ξβ The differential element of the surface located in ℜ3 is defined by the expression dS ¼ ðn21 þ n22 þ n23 Þ1=2 dξ1 dξ2 ;

n

n

Q

n

f ðxÞ  ∑ f ðxq Þfnq ðξÞ;

ð100Þ

q¼1

ð106Þ

where the components of the unit normal to the surface vector are ∂x2 ∂x3 ∂x2 ∂x3  ; ∂ξ1 ∂ξ2 ∂ξ2 ∂ξ1 ∂x3 ∂x1 ∂x1 ∂x3 n2 ¼  ; ∂ξ1 ∂ξ2 ∂ξ1 ∂ξ2 ∂x1 ∂x2 ∂x2 ∂x1  : n3 ¼ ∂ξ1 ∂ξ2 ∂ξ1 ∂ξ2

n1 ¼

ð107Þ

On the boundary element, we consider Q nodes with global coordinates xqi of the nodal points and fq ðξÞ to be shape or interpolation functions. Then, the global coordinates can be expressed as Q

xi ðξÞ ¼ ∑ xqi fq ðξÞ; q¼1

Q

yi ðξÞ ¼ ∑ yqi fq ðξÞ

ð108Þ

q¼1

In general, on the BE with Q nodes of interpolation, the shape function is an approximation polynomial of Q degree. The interpolation polynomial of Q degree has the form  fqi ðξi Þ ¼

n¼1

On each BE, the local functions f ðxÞ may be represented in the form

dx ¼ JðξÞ dξ;

ðξi  ξ1i Þ…ðξi  ξqi  1 Þðξi  ξqi þ 1 Þ…ðξi  ξQi Þ

ð109Þ

ðξqi  ξ1i Þ…ðξqi  ξqi  1 Þðξqi  ξqi þ 1 Þ…ðξqi  ξQi Þ

We can write the shape function, which depends on two variables fq ðξ1 ; ξ2 Þ, in terms of the product of two interpolation polynomials of one variable, fq1 ðξ1 Þ and fq2 ðξ2 Þ, as fq ðξ1 ; ξ2 Þ ¼ fq1 ðξ1 Þfq2 ðξ2 Þ

ð110Þ q

where fnq ðξÞ are interpolation polynomials or nodal functions of the BE with number n. At the nodal point with coordinates xq , they are equal to 1, and at other nodal points, they are equal to zero. Taking into account (99) and (100), the global approximation of the functionf ðxÞ looks like N

f ðxÞ  ∑

Q

n

∑ f ðxq Þfnq ðξÞ:

ð101Þ

n¼1q¼1

If the nodal point q belongs to several BEs, it is considered in these sums only once. We introduce a local system of coordinates on each BE as ξ ¼ ðξ1 ; ξ3 ; ξ3 Þ, such that the local coordinates ξi are functions of the global ðξi ðx1 ; x2 ; x3 ÞÞ ones. Conversely, the global ones xi are functions of the localðxi ðξ1 ; ξ2 ; ξ3 ÞÞ ones. For these maps to be one-to-one, it is necessary and sufficient that the Jacobians of transformations be nonzero     ∂x  ∂ξ  J ¼ Det  i  a 0; J  1 ¼ Det  i  a 0 ð102Þ ∂ξj ∂xj

r

The distance between points with coordinates x and y in local coordinates is rðξ; ζÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ððxq1 fq ðξÞ  yr1 fr ðζÞÞ2 þ ðxq2 fq ðξÞ  yr2 fr ðζÞÞ2 þ ðxq3 fq ðξÞ  yr3 fr ðζÞÞ2

ð111Þ Substituting the expression for global coordinates from (108) into (104)–(107), the derivatives Q ∂fq ðξÞ ∂xi ðξÞ ¼ ∑ xqi ∂ξj ∂ξj q¼1

ð112Þ

the coordinates of the unit normal vectors Q

n1 ¼ ∑ xq2 q¼1 Q

n2 ¼ ∑ xq3 q¼1 Q

n3 ¼ ∑ xq1 q¼1

Q ∂fq ðξÞ Q q ∂fq ðξÞ ∂fq ðξÞ Q q ∂fq ðξÞ  ∑ xq2 ∑ x ∑ x ∂ξ1 q ¼ 1 3 ∂ξ2 ∂ξ2 q ¼ 1 3 ∂ξ1 q¼1 Q ∂fq ðξÞ Q q ∂fq ðξÞ ∂fq ðξÞ Q q ∂fq ðξÞ  ∑ xq1 ∑ x ∑ x ∂ξ1 q ¼ 1 1 ∂ξ2 ∂ξ1 q ¼ 1 3 ∂ξ2 q¼1 Q ∂fq ðξÞ Q q ∂fq ðξÞ ∂fq ðξÞ Q q ∂fq ðξÞ  ∑ xq2 ∑ x2 ∑ x ∂ξ1 q ¼ 1 ∂ξ2 ∂ξ1 q ¼ 1 1 ∂ξ2 q¼1

ð113Þ

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

and the Jacobian    Q ∂fq ðξÞ  J ¼ Det  ∑ xqi  q ¼ 1 ∂ξj 

171

ð114Þ

All of these parameter are used when integrals in (96) are calculated in local coordinates.

6. Piecewise constant approximation The piecewise constant approximation is the simplest one. Interpolation functions in this case do not depend on the BE form and the dimension of the domain. They have the form ( 1 8 x A Sn ; fq ðxÞ ¼ ð115Þ 0 8 x2 = Sn : For the purposes of simplicity, we transform the global system of coordinates such that the origin is located at the nodal point, where y0 ¼ 0, the coordinate axes x1 and x2 are located in the plane of the element, and the axis x3 is perpendicular to that plane. In this case, x3 ¼ 0 and n1 ¼ 0, n2 ¼ 0, n3 ¼ 1. In this specific case, the fundamental solutions in (13) contain the following divergent integrals Z i j x1 x2 J i;j ¼ fq ðxÞ dS for i; j ¼ 0; 1; 3 and k ¼ 1; 3; 5: ð116Þ k k S r Regular representations for integrals with these kernels can easily be obtained from the general formula (92). They can be presented in the following form Weakly singular R R rn 0;0 J 1 ¼ W:S: Sn dS r ¼ ∂Sn r dl;  Z  2 Z 2 x1 x1 r n 2r n 2x1 n1 1  dl; dS ¼ þ J 2;0 3 ¼ W:S: 3 3 ∂Sn r 3 r r Sn r  Z  2 Z 2 x2 x2 r n 2r n 2x2 n2 1 J 0;2  dl dS ¼ þ 3 ¼ W:S: 3 3 ∂Sn r 3 r r Sn r Z Z x1 x2 1 x1 x2 r n r þ  J 1;1 dl ð117Þ dS ¼  3 ¼ W:S: 3 3 r r3 Sn r ∂Sn Singular  R R ¼ P:V: Sn xr31 dS ¼ ∂Sn x1r3rn  nr1 dl; Z Z x2 x2 r n n2  J 0;1 dl dS ¼  3 ¼ P:V: 3 r r r3 Sn ∂Sn

though the nodal points and unit vector normal to the contour xi ðkÞ ¼

xki þ 1 þ xki ; 2

n^ 1 ðkÞ ¼

xk2 þ 1  xk2 ; 2Δk

n^ 2 ðkÞ ¼

xk1 þ 1  xk1 2Δk

ð120Þ

The coordinates of an arbitrary point on the contour ∂V n may be represented in the form x1 ðξÞ ¼ x1 ðkÞ þ ξΔk n^ 2 ðkÞ and x2 ðξÞ ¼ x2 ðkÞ þ ξΔk n^ 1 ðkÞ

ð121Þ

where x1 ðkÞ and x2 ðkÞ are the coordinates of the k-th side of the ^ n^ 1 ; n^ 2 Þ is a unit vector normal to the contour, ξ A ½  1; 1 contour, nð is a parameter of integration along the k-th side and 2Δk is the length of a k-th side. The coordinates can be calculate though the nodal points and unit vector normal to the contour qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xk þ 1 þ xki ; 2Δk ¼ ðxk1 þ 1  xk1 Þ2 þ ðxk2 þ 1  xk2 Þ2 ; ð122Þ xi ðkÞ ¼ i 2 Some useful notations are the following

J 1;0 3

Hypersingular R dS R rn J 0;0 3 ¼ F:P: Sn r 3 ¼  ∂Sn r 3 dl;  Z 2 Z  2 x1 x1 r n 2r n 2x1 n1 J 2;0 ¼ F:P: dS ¼   dl; 5 5 r5 3r 3 3r 3 Sn r ∂Sn  Z 2 Z  2 x2 x2 r n 2r n 2x2 n2 dl; dS ¼  3 J 0;2 5 ¼ F:P: 5 5 3 r 3r 3r Sn r ∂Sn Z Z x1 x2 x1 x2 r n r þ  J 1;1 dS ¼  3 dl; 5 ¼ F:P: 5 r5 3r Sn r ∂Sn

Fig. 1. Convex polygon.

rðξÞ ¼

ð118Þ

ð119Þ

where r n ¼ xα nα , r þ ¼ x1 n2 þ x2 n1 . Let us consider the contour ∂V n as a polygon with K vertexes, as shown in Fig. 1. The approach developed here will be used to calculate the divergent integrals (117)-(119). All the calculations will be performed using the local rectangular coordinate system with its origin located at point yq , the x1 and x2 axes located in the plane of the polygon and the x3 axis perpendicular to this plane. The global coordinates of the vertexes are ðxk1 ; xk2 Þ. They can be calculated

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Δ2k ξ2 þ 2ξΔk r þ ðkÞ þr 2 ðkÞ;

r n ðξÞ ¼ r n þ 2Δ1 n^ 1 ðkÞn^ 2 ðkÞξ

rðkÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x21 ðkÞþ x22 ðkÞ;

r n ðkÞ ¼ xα ðkÞn^ α ðkÞ;

ð123Þ

Analytical formulae for the calculation of the divergent integrals (117)-(119) are presented in [40]. It is important to mention that the obtained formulae can be applied to the calculation of regular integrals, which appear when the collocation point yq is situated outside of the area of integration Sn . Those formulae are valid for any collocation point situated inside or outside of the BE. Special considerations are only required for points situated at the vertexes of the boundary element; this will be addressed in the following sections. Tables 1 and 2 present the results of the divergent and regular integral calculations for the unit side square with coordinates of the vertexes 1  f0:5;  0:5g; 2  f0:5; 0:5g; 3  f  0:5; 0:5g; 4  f  0:5;  0:5g and for the unit side triangle with coordinates of the vertexes 1  f0:5; 0:289g; 2  f0; 0:577g; 3 f  0:5;  0:289g, respectively. Calculations have been done for collocation points with coordinate y01 ¼ 0:0 and different coordinates y02 , which correspond to points situated both inside and outside of the integration area. To validate the regularized eqs. (117)-(119), we compare our calculations of the hypersingular integrals with the results reported by Ioakimidis in [23] and, for weakly singular and regular integrals,

172

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

Table 1 Divergent integrals calculated for the unit square at different collocation points. Points

J 0;0 1

J 2;0 3

J 0;2 3

J 1;1 3

J 1;0 3

J 0;1 3

J 0;0 3

J 2;0 5

J 0;2 5

J 1;1 5

0.0 0.1 0.25 2.5 5.0

3.525 3.496 3.335 0.402 0.200

1.762 1.748 1.659 0.005 0.001

1.762 1.748 1.675 0.397 0.199

0.0 0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0

0.0  0.5782  1.636  0.163  0.0401

 11.31  11.61  13.75 0.067 0.008

 5.656  5.727  6.289 0.0 0.0

 5.656  5.882  7.462 0.067 0.008

0.0 0.0 0.0 0.0 0.0

Table 2 Divergent integrals calculated for the unit triangle at different collocation points. Points

J 0;0 1

J 2;0 3

J 0;2 3

J 1;1 3

J 1;0 3

J 0;1 3

J 0;0 3

J 2;0 5

J 0;2 5

J 1;1 5

0.0 0.10 0.25 2.50 5.00

2.281 2.239 2.039 0.173 0.086

1.141 1.048 0.827 0.0010 0.0001

1.141 1.190 1.212 0.172 0.086

0.0 0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0

0.0  0.808  1.840  0.0702  0.0173

 18.0  19.081  25.020 0.0287 0.0034

 9.0  10.263  14.452 0.0 0.0

 9.0  8.817  10.567 0.0285 0.0034

0.0 0.0 0.0 0.0 0.0

0,0

2,0

J1

J3

3.5 1.5

3.0 2.5

1.0

2.0

0.5

1.5 −1.0

−0.5

0.0

0.5

1.0

y02

−1.0

−0.5

0,2

0.5

1.0

y02

1,1

J3

J3

2.0

0.10

1.8

0.05

1.6 1.4 −1.0

1.2

−0.5

0.0

0.5

1.0

y02

−0.05

1.0 −1.0

−0.5

0.5

1.0

y02

−0.10

Fig. 2. Calculations of weakly singular integrals for the unit square versus y02 .

with the results obtained using regular 2-D numerical Gauss quadrature. Our calculations show that the results obtained using the regularized equations presented here agree with the results obtained by other methods. It is also important to mention that there are two possibilities calculations of the integrals in regular representations (117)-(119). The first one calculates the corresponding regularized integrals using numerical integration, and the second calculates the corresponding integrals using the analytical formulae presented in [40]. Our calculations show that with the analytical formulae, the results are more accurate and the calculation time is 5-7 times faster than numerical integration of the integrals and 8-12 times faster than 2-D numerical calculation using Gaussian quadrature. As previously mentioned, regularized eqs. (117)-(119) are valid for collocation points situated inside or outside of the integration area Sn but are not valid in the vicinity of the boundary of the integration area, including the boundary ∂Sn . Figs. 2–7 present the results from the calculation of the divergent integrals for the collocation point situated on the vertical line that passes through

the point y0 ¼ 0 for the unit square and triangle, respectively. Our calculations show that regularized formula (117) for weakly singular integrals is valid everywhere in the area Sn , including the boundary ∂Sn . For singular and hypersingular integrals, regularized formulae (118), (119) give good results for the collocation point situated inside or outside of the integration area Sn , but in the vicinity of boundary ∂Sn , the corresponding integrals are divergent. For regular integrals, formulae (117)-(119) have very good coincidence with the calculation using Gaussian quadrature in both cases, weakly singular and hypersingular integrals. It is apparent that Gaussian quadrature cannot be used to calculate divergent integrals. It is important to mention that we can calculate the direct value of the hypersingular integrals when the collocation point belongs to the contour of integration y0 A ∂Sn . If the collocation point is situated inside of side k of the polygon, we have to exclude that side from summation in the corresponding formulae with number k and calculate the rest. For side k, we have a 1-D hypersingular

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180 1,0

0,1

J3

−1.0

J3

10

10

5

5

−0.5

173

0.5

1.0

y02

−1.0

−0.5

−5

−5

−10

−10

0.5

1.0

0.5

1.0

0.5

1.0

y02

Fig. 3. Calculations of singular integrals for the unit square versus y02 .

0,0

2,0

J3

J5

15

40

10

20

5 0

−1.0

−0.5

0.5

1.0

y2

0

−1.0

−0.5

−5

−20

−10

−40

−15 −20 0,2

1,1

J5 30

J5 0.02

20

0.01

10 −1.0

−0.5

y2

−10

0.5

1.0

y02

0

−1.0

−0.5

y2

−0.01

−20

−0.02

−30

Fig. 4. Calculations of hypersingular integrals for the unit square versus y02 .

0,0

2,0

J1

J3 1.0

2.0

0.8

1.5

0.6 0.4

1.0

0.2 -1.0

-0.5

0.0

0.5

1.0

y02

-1.0

J30,2

0.10

0.8

0.05

0.6 0.0

y02

0.15

1.0

-0.5

1.0

J31,1

1.2

-1.0

0.5

-0.5

-1.0

0.5

1.0

y02

-0.5

-0.05

0.5

1.0

y02

-0.10

Fig. 5. Calculations of weakly singular integrals for the unit triangle versus y02 .

integral of type (51), which must be accounted for in the final result. Similarly, integrals can be calculated when the collocation point belongs to the vertex for the polygon (see [40] for details).

The final results for a rectangle with the collocation point situated at the middle of the down side y01 ¼ 0; y02 ¼  0:5 are J 1 ¼ 2:406, J 3 ¼  6:471, for a triangle with the collocation point situated at

174

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180 1,0

0,1

J3

J3

6

5

4 2 −1.0

−0.5

−2

0.5

1.0

y02

−1.0

−0.5

0.5

1.0

y02

0.5

1.0

y02

0.5

1.0

y02

−5

−4 −6

Fig. 6. Calculations of singular integrals for the unit triangle versus y02 .

2,0

J30,0

J5 20

40

10

20 −1.0

−0.5

−20

0.5

1.0

y02

−1.0

−0.5

−40

−20

−60

−30

−80

−40

J50,2

J51,1

30

0.3

20

0.2 0.1

10 −1.0

−0.5

−10

−10

0.5

1.0

y02

−1.0

−0.5

−0.1

−20

−0.2

−30

−0.3

−40

−0.4

Fig. 7. Calculations of hypersingular integrals for the unit triangle versus y02 .

the middle of the down side y01 ¼ 0; y02 ¼ 0:289 are J 1 ¼ 1:616, J 3 ¼  8:309 and for a vertex y01 ¼ 0; y02 ¼ 0:577 are J 1 ¼ 0:951, J 3 ¼  3:154. It is important to mention that all calculations shown here can be done analytically and that no numerical integration is needed.

Of course, for collocation points that do not belong to the vertexes of the rectangular or triangular BE, we can use slightly modified formulae developed for piecewise constant approximation. However, for piecewise linear approximation, it is very important when the collocation point does not belong to the vertexes. Therefore, we developed equations that are suitable for that case. Regular representations for the integrals with kernels (116) in the case of piecewise linear approximation can easily be obtained from general formula (92). They have the form Weakly singular i R j ðξÞ R h ¼ W:S: Sn qr dS ¼ ∂Sn jq ðξÞrrn  r∂n jq ðξÞ dl

 2   2   Z Z  x2 x r n 2r n 2x1 n1 x 1 þ 1  2r ∂n jq ðξÞ dl J 2;0 jq ðξÞ 31 dS ¼ jq ðξÞ 13 þ  q;3 ¼ W:S: r r r 3 r r Sn ∂Sn   2   2  Z Z  x2 x r n 2r n 2x2 n2 x 1 jq ðξÞ 32 dS ¼ jq ðξÞ 23 þ J 0;2  þ 2  2r ∂n jq ðξÞ dl q;3 ¼ W:S: r r r 3 r r Sn ∂Sn Z Z

x x r  x1 x2 1 r þ  x1 x2 1 2 n þ J 1;1 jq ðξÞ 3 dS ¼ jq ðξÞ  ∂n jq ðξÞ dl: q;3 ¼ W:S: 3 3 r r r r Sn ∂Sn

Z

Z

x r  n1  x1 1 n þ ∂n jq ðξÞ dl;  3 r r r ∂Sn Z Z

 jq ðξÞx2 x2 r n n2  x2 0;1 þ ∂n jq ðξÞ dl: J q;3 ¼ P:V: dS ¼  3 3 r r r r Sn ∂Sn Hypersingular J 1;0 q;3 ¼ P:V:

jq ðξÞx1 dS ¼ r3 Sn

ð125Þ

 R jq ðξÞ R rn 1 J 0;0 q;3 ¼ F:P: Sn r 3 dS ¼  ∂Sn jq ðξÞr 3 þ r ∂n jq ðξÞ dl;

7. Piecewise linear approximation

J 0;0 q;1

Singular

ð124Þ

 2  Z  x2 x r n 2r n 2x1 n1 jq ðξÞ 51 dS ¼ jq ðξÞ 15  3  3 r r 3r 3r Sn ∂Sn  2   x1 2 ∂n jq ðξÞ dl; þ  3r 3 3r Z

J 2;0 q;5 ¼ F:P:

  2   2  Z Z  x2 x r n 2r n 2x2 n2 x2 2 ∂n jq ðξÞ dl; þ J 0;2 jq ðξÞ 52 dS ¼ jq ðξÞ 25  3   q;5 ¼ F:P: 3 3 3r r r 3r 3r 3r Sn ∂Sn Z Z

x x r  x1 x2 r þ  x1 x2 1 2 n jq ðξÞ 5 dS ¼ jq ðξÞ  3 þ 3 ∂n jq ðξÞ dl: J 1;1 q;5 ¼ F:P: r r5 3r 3r Sn ∂Sn

ð126Þ

For details, one can refer to [39,40]. Analysis of these expressions shows that we have to calculate the sum of integrals of the following type Z J l;m q;p ¼

∂Sn

fq ðξÞ

K xl1 ðξÞxm 2 ðξÞ dlðξÞ ¼ ∑ r p ðξÞ k¼1

Z lk

fq ðξÞ

xl1 ðξÞxm 2 ðξÞ dlðξÞ r p ðξÞ

ð127Þ

This formula will be used for calculation over the rectangular and triangular BE of the divergent integrals of type (124)-(126) that

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

occur when the boundary-value problems (2), (5) are solved by the BEM.

The quadrilateral BE is defined by its angular nodes, and its shape functions are f1 ¼ 1=4ð1  ξ1 Þð1  ξ2 Þ;

f2 ¼ 1=4ð1 þ ξ1 Þð1 ξ2 Þ;

f3 ¼ 1=4ð1 þ ξ1 Þð1 þ ξ2 Þ;

f4 ¼ 1=4ð1 ξ1 Þð1 þ ξ2 Þ;

for

7.1. Rectangular boundary elements To simplify the situation, we transform the global coordinate system such that the origins of the global and local systems of coordinates coincide. The coordinate axes x1 and x2 are located in the plane of the element and coincide with the local ones ξ1 and ξ2 , while the axis x3 is perpendicular to that plane. In this case, x3 ¼ 0 and n1 ¼ 0, n2 ¼ 0, n3 ¼ 1. Let us consider the rectangle BE that is shown in Fig. 8.

175

ξ1 A ½0; 1;

ξ2 A ½0; 1:

ð128Þ

Then, the global coordinates can be expressed as functions of the local ones in the form 4

xi ðξ1 ; ξ2 Þ ¼ ∑ xqi fq ðξ1 ; ξ2 Þ; q¼1

4

yi ðξ1 ; ξ2 Þ ¼ ∑ yqi fq ðξ1 ; ξ2 Þ

ð129Þ

q¼1

The derivatives of the shape functions are ∂f1 ðξÞ ∂f1 ðξÞ ∂f2 ðξÞ ¼  1=4ð1 ξ2 Þ; ¼  1=4ð1 ξ1 Þ; ¼ 1=4ð1 ξ2 Þ; ∂ξ1 ∂ξ2 ∂ξ1

∂f2 ðξÞ ¼  1=4ð1 þ ξ1 Þ ∂ξ2 ∂f3 ðξÞ ∂f3 ðξÞ ∂f4 ðξÞ ¼ 1=4ð1 þ ξ2 Þ; ¼ 1=4ð1 þ ξ1 Þ; ¼  1=4ð1 ξ2 Þ; ∂ξ1 ∂ξ2 ∂ξ1

∂f4 ðξÞ ¼ 1=4ð1  ξ1 Þ: ∂ξ2

ð130Þ

The coordinates of the nodal points are 1  fx11 ¼  Δ1 ; x12 ¼  Δ2 g; 3  fx31

¼ Δ1 ; x32

¼ Δ2 g;

2  fx21 ¼ Δ1 ; x22 ¼  Δ2 g;

4  fx41 ¼  Δ1 ; x42 ¼ Δ2 g

ð131Þ

We also introduce some additional useful notations below rðξ; yq Þ ¼

Fig. 8. Rectangular BE.

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx1  yq1 Þ2 þ ðx2  yq2 Þ2 ¼ ðΔ1 ð1 þ ξ1 Þ  yq1 Þ2 þ ðΔ2 ð1 þ ξ2 Þ  yq2 Þ2 ;

Table 3 Divergent integrals calculated for the unit square at different collocation points. Points

J 0;0 1;1

J 2;0 1;3

J 0;2 1;3

J 1;1 1;3

0.0 0.10 0.25 2.50 5.00

0.881 0.787 0.633 0.093 0.048

0.440 0.367 0.252 0.0011 0.0

0.440 0.420 0.380 0.092 0.048

Points

J 0;0 2;1

J 2;0 2;3

J 0;2 2;3

J 1;1 2;3

0.0 0.10 0.25 2.50 5.00

0.881 0.787 0.633 0.093 0.048

0.440 0.367 0.252 0.0011 0.0

0.440 0.420 0.380 0.092 0.048

 0.058  0.072  0.081  0.005  0.001

Points

J 0;0 3;1

J 2;0 3;3

J 0;2 3;3

J 1;1 3;3

0.0 0.10 0.25 2.50 5.00

0.881 0.961 1.034 0.107 0.051

0.440 0.506 0.577 0.0016 0.0

0.440 0.454 0.457 0.105 0.051

0.058 0.037  0.009  0.007  0.001

Points

J 0;0 4;1

J 2;0 4;3

J 0;2 4;3

J 1;1 4;3

J 1;0 4;3

0.0 0.10 0.25 2.50 5.00

0.881 0.961 1.034 0.107 0.051

0.440 0.506 0.577 0.0016 0.0

0.440 0.454 0.457 0.105 0.051

 0.058  0.037 0.009 0.007 0.001

 0.881  1.013  1.154  0.003 0.0

Points

J 0;0 1

J 2;0 3

J 0;2 3

J 1;1 3

J 1;0 3

0.0 0.10 0.25 2.50 5.00

3.525 3.449 3.335 0.402 0.200

1.762 1.748 1.659 0.005 0.0

1.762 1.748 1.675 0.397 0.199

0.058 0.072 0.081 0.005 0.001

0.0 0.0 0.0 0.0 0.0

J 1;0 1;3

J 0;1 1;3

 0.881  0.734  0.505  0.002 0.0 J 1;0 2;3

 0.881  0.989  1.042  0.035  0.009 J 0;1 2;3

J 0;1 3;3

0.881 1.013 1.154 0.003

0.0 0.0 0.0 0.0 0.0

 2.828  2.032  0.900 0.013 0.0018 J 0;0 2;3

 0.881  0.989  1.042  0.035  0.009

0.881 0.734 0.505 0.002 0.0 J 1;0 3;3

J 0;0 1;3

 2.828  2.032  0.900 0.013 0.0018 J 0;0 3;3

0.881 0.700 0.224  0.046  0.010 J 0;1 4;3

 2.828  3.772  5.975 0.020 0.0022 J 0;0 4;3

0.881 0.700 0.224  0.046  0.010

 2.828  3.772  5.975 0.020 0.0022

J 2;0 1;5

J 0;2 1;5

 1.414  1.095  0.623 0.0 0.0

 1.414  0.936  0.273 0.013 0.018

J 2;0 2;5

J 0;2 2;5

J 1;1 2;5

 1.414  1.095  0.623 0.0 0.0

 1.414  0.936  0.273 0.013 0.018

 0.471  0.503  0.484 0.0

J 2;0 3;5

J 0;2 3;5

J 1;1 3;5

 1.414  1.767  2.525 0.0 0.0

 1.414  2.004  3.457 0.020 0.022

J 2;0 4;5

J 0;2 4;5

J 1;1 4;5

 1.414  1.767  2.525 0.0 0.0

 1.414  2.004  3.457 0.020 0.022

 0.471  0.404  0.187 0.0 0.0 J 1;1 5

J 0;1 3

J 0;0 3

J 2;0 5

J 0;2 5

0.0  0.578  1636  0.163  0.040

 11.31  11.61  13.75 0.067 0.0081

 5.656  5.727  6.297 0.0 0.0

 5.656  5.882  7.462 0.066 0.081

J 1;1 1;5 0.471 0.503 0.484 0.0

0.471 0.404 0.187 0.0 0.0

0.0 0.0 0.0 0.0 0.0

176

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

r n ¼ ðxα  yqα Þn^ α ; n^ 2 ðkÞ ¼

dl ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Δ21 dξ21 þ Δ22 dξ22 ;

n^ 1 ðkÞ ¼

methods. As in the case of piecewise constant approximation, our calculations show that with the analytical formulae, the results are more accurate and that the calculation time is 5-7 times faster than numerical integration of the corresponding integrals and 812 times faster than 2-D numerical integration. In Figs. 9–11, the results of the calculation of the integrals (127) for collocation points situated on the vertical line that passes through the point y0 ¼ 0 for the unit square are presented. The

xk2 þ 1  xk2 ; Δk

xk1 þ 1  xk1 : Δk

ð132Þ

Calculations will be done side-by-side using formula (127). The results of the calculations are presented side-by-side in [40]. Table 3 presents the results of the divergent and regular integrals of type J q;p and also the integrals of type J p ¼ ∑4q ¼ 1 J q;p and the calculations for the unit side square with coordinates of the vertexes 1  f0:5; 0:5g; 2  f0:5; 0:5g; 3 f  0:5; 0:5g; 4  f 0:5;  0:5g. Calculations have been done for collocation points with coordinate y01 ¼ 0:0 and different coordinates y02 that correspond to points situated inside and outside of the integration area. To validate the regularized equations obtained here, we compare our calculations with the results presented in the Table 3 and with results obtained using regular 2-D Gaussian numerical calculation for regular integrals. Our calculations show that the results of the calculations obtained using the regularized equations presented here agree with the results obtained by other

Fig. 12. Triangular BE in the global frame.

0,0

0,0

J1

Jq,1 1.0

2.0

0.8 0.6

1.5

0.4 1.0 −1.0

−0.5

0.0

0.5

1.0

y02

0.2 −1.0

−0.5

0.0

0.5

1.0

y02

Fig. 9. Calculations of weakly singular integrals for the unit square versus y02 .

0,1

0,1

J3

Jq,3

6

2

8 4

1

2

−1.0

−0.5

0.5

−2

1.0

y02

−1.0

−0.5

0.5

1.0

y02

−1 −2

−4

Fig. 10. Calculations of singular integrals for the unit square versus y02 .

0,0 Jq,3

0,0

J3

10

20 −1.0

5 0.5

−0.5

1.0

y02

−1.0

0.5

−0.5

−20

−5

−40

−10

Fig. 11. Calculations of weakly singular integrals for the unit square versus y02 .

1.0

y02

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

177

The quadrilateral BE is defined by its angular nodes, and its shape functions are

following diagrams present results that correspond to each shape function on the right and to their sum on the left. Our calculations show that the regularized formulae for weakly singular integrals are valid everywhere, including the boundary ∂Sn . For singular and hypersingular integrals, the regularized formulae yield good results for collocation points situated inside or outside of the integration area Sn ; however, in the vicinity of the boundary ∂Sn , the corresponding integrals are divergent. It is important to note that all calculations shown here can be done analytically and that no numerical integration is needed.

f1 ðξ1 ; ξ2 Þ ¼ ξ1 ;

f2 ðξ1 ; ξ2 Þ ¼ ξ2 ;

f3 ðξ1 ; ξ2 Þ ¼ ð1  ξ1  ξ2 Þ;

ξ1 A ½0; 1;

ξ2 A ½0; 1

ð133Þ Then, the global coordinates can be expressed as functions of the local ones in the form 3

xi ðξ1 ; ξ2 Þ ¼ ∑ xqi fq ðξ1 ; ξ2 Þ; q¼1

3

yi ðξ1 ; ξ2 Þ ¼ ∑ yqi fq ðξ1 ; ξ2 Þ

ð134Þ

q¼1

or x1 ¼ x31 þ Δx31 ξ1  Δx21 ξ2 ;

7.2. Triangular boundary elements

Δx3i

x2 ¼ x32 þ Δx32 ξ1  Δx22 ξ2

∂f1 ðξÞ ¼ 1; ∂ξ1

∂f1 ðξÞ ¼ 0; ∂ξ2

Δx2i

ð135Þ

where and ¼ The derivatives of the shape functions are

To simplify the situation, we transform the global system of coordinates such that the origins of the global and local systems of coordinates coincide. The coordinate axes x1 and x2 are located in the plane of the element and coincide with the local ones ξ1 and ξ2 , while the axis x3 is perpendicular to that plane. In this case x3 ¼ 0 and n1 ¼ 0, n2 ¼ 0, n3 ¼ 1. Let us consider the triangluar BE that is shown in Fig. 12.

¼ ðx1i  x3i Þ

∂f2 ðξÞ ¼ 0; ∂ξ1

ðx2i  x3i Þ. ∂f2 ðξÞ ¼ 1; ∂ξ2

∂f3 ðξÞ ¼  1; ∂ξ1

∂f3 ðξÞ ¼ 1 ∂ξ2

ð136Þ Taking into account that the coordinates ξ1 and ξ2 are oblique, normal derivatives have to be calculated using the following

Table 4 Divergent integrals calculated for the unit triangle at different collocation points Points

J 0;0 1;1

0.289 0.389 0.538 2.789 5.289

0.760 0.680 0.530 0.056 0.028

Points

J 0;0 2;1

0.289 0.389 0.538 2.789 5.289

0.760 0.878 0.977 0.061 0.029

Points

J 0;0 3;1

J 2;0 3;3

J 0;2 3;3

0.289 0.389 0.538 2.789 5.289

0.760 0.680 0.530 0.056 0.028

0.397 0319 0.195 3.8  10  4 0.0

0.362 0360 0.335 0.055 0.028

Points

J 0;0 1

J 2;0 3

J 0;2 3

0.289 0.389 0.538 2.789 5.289

2.281 2.239 2.039 0.173 0.086

J 2;0 1;3

J 0;2 1;3 0.397 0.319 0.195 3.8  10  4 0.0

J 2;0 2;3

J 1;1 1;3  0.030  0.053  0.067  0.002 0.0

0.362 0.360 0.335 0.055 0.028 J 0;2 2;3

0.345 0.409 0.438 2.3  10  4 0.0

J 1;1 2;3 0.415 0.469 0.540 0.061 0.029

1.140 1.048 0.827 0.001 0.0

J 1;0 1;3

0.030 0.053 0.067 0.002 0.0 J 1;1 3 0.0 0.0 0.0 0.0 0.0

J 2;0 1;5

 6.0  4.792  3.666 0.008 0.001

J 0;2 1;5

 2.880  2.601  2.397 0.0 0.0

J 0;0 2;3

1.316 1.012 0.255  0.026  0.006

 6.0  9.496  17.68 0.011 0.001

J 1;0 3;3

J 0;1 3;3

J 0;0 3;3

 1.140  1.048  0.827  0.001 0.0

 0.658  0.910  1.047  0.021  0.005

J 1;0 3

J 0;1 3

J 0;0 3

J 2;0 5

J 0;2 5

0.0  0.808  1.840  0.070  0.017

 18.0  19.08  25.02 0.028 0.003

 9.0  10.26  14.45 0.0 0.0

 9.0  8.817  10.56 0.028 0.003

0.0 0.0 0.0 0.0 0.0

J 2;0 2;5

J 0;2 2;5

 3.239  5.060  9.658 0.0 0.0 J 2;0 3;5

 6.0  4.792  3.666 0.008 0.001

 2.880  2.601  2.397 0.0 0.0

J q,1 1.0

2.0

0.8 0.6

1.5

0.4

1.0

0.2 0.0

0.5

1.0

y02

−0.5

0.0

0.5

1.0

Fig. 13. Calculations of weakly singular integrals for the unit triangle versus y02 .

 2.760  4.436  8.029 0.011 0.001 J 0;2 3;5

0,0

0,0

J1

−0.5

 3.119  2.190  1.269 0.008 0.001

J 0;1 2;3

0.0 0.0 0.0 0.0 0.0

J 1;1 3;3

J 0;0 1;3

 0.658  0.910  1.047  0.021  0.005

1.140 1.048 0.827 0.001 0.0 J 1;0 2;3

0.0 0.0 0.0 0.0 0.0

1.140 1.190 1.212 0.172 0.086

J 0;1 1;3

y02

 3.119  2.190  1.269 0.008 0.001

J 1;1 1;5  0.207  0.392  0.578 0.0 0.0 J 1;1 2;5 0.0 0.0 0.0 0.0 0.0 J 1;1 3;5 0.207 0.392 0.578 0.0 0.0 J 1;1 5 0.0 0.0 0.0 0.0 0.0

178

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180 0,0

0,0

J1

J q,1 3

5

2 0.5

−0.5

1

y02

1.0

−0.5

−5

0.5

−1

1.0

y02

1.0

y02

−2 −10

−3 Fig. 14. Calculations of singular integrals for the unit triangle versus y02 . 0,0

0,0

Jq ,3

Jq ,3

10

50

5 0.5

y02

1.0

−1.0

0.5

−0.5 −5

−50

−10 −15

−100

Fig. 15. Calculations of weakly hypersingular integrals for the unit triangle versus y02 .

formula ∂n ¼ n^ 1

  ∂ ∂ ∂ n^ 1 Δ2 n^ 2 Δ3 n^ α ð2Þ þ n^ α ð3Þ þ n^ 2 ¼ ∂x1 ∂x2 ∂ξα Δ Δ

ð137Þ

where Δ ¼ Δx31 Δx22  Δx21 Δx32 . The normal derivatives of the shape functions are ∂n f1 ðkÞ ¼

n^ 1 ðkÞn^ 1 ð2ÞΔ2 n^ 2 ðkÞn^ 1 ð3ÞΔ3 þ ; Δ Δ

∂n f3 ðkÞ ¼ 

∂n f2 ðkÞ ¼

n^ 1 ðkÞn^ 2 ð2ÞΔ2 n^ 2 ðkÞn^ 2 ð3ÞΔ3 þ ; Δ Δ

n^ 1 ðkÞðn^ 1 ð2Þ þ n^ 2 ð2ÞÞΔ2 þ n^ 2 ðkÞðn^ 1 ð3Þ þ n^ 2 ð3ÞÞΔ3 Δ

ð138Þ

The coordinates of the nodal points are point 1-ðx11 ; x12 Þ, point 2-ðx21 ; x22 Þ and point 3-ðx31 ; x32 Þ. The length of the triangle sides and radius are Δk ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx1k þ 1  xk1 Þ2 þðx2k þ 1  xk2 Þ2 ;

rðξ; yq Þ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx1  yq1 Þ2 þ ðx2  yq2 Þ2

Fig. 16. Solution for 80 rectangular and 160 triangular BEs.

ð139Þ Calculations will be performed side-by-side using formula (127). The results from the calculations are presented side-by-side in [40]. Table 4 presents the results of the calculations of the divergent and regular integrals of type J q;p and also the integrals of type J p ¼ ∑3q ¼ 1 J q;p for the unit side triangle with coordinates of the vertexes 1  f0:0;  0:0g; 2  f0:0; 1:0g; 3  f0:5; 0:867g. Calculations have been done for collocation points with coordinate y01 ¼ 0:5 and different coordinates y02 that correspond to points situated inside and outside of the integration area. Fig. 13–15 present the results of the integral calculations (127) for collocation points situated on the vertical line that passes through the middle point of the unit triangle. On the right diagrams, we present results that correspond to each shape function and, on the left, their sum. Our calculations show that the regularized formulae for weakly singular integrals are valid everywhere, including the boundary ∂Sn . For hypersingular integrals, the regularized formulae give good results for collocation points situated inside or outside of the integration area Sn , but in the vicinity of the boundary ∂Sn , the corresponding integrals are divergent. It is important to note that all calculations can be done analytically and that no numerical integration is needed.

8. Numerical solution of the hypersingular integral equation Let us consider in a circular area located in the plane  R2 ¼ fx : x3 ¼ 0g defined by coordinates Ω ¼ x21 þ x22 r R; x3 ¼ 0 a hypersingular integral equation of type Z uðxÞ dS ¼ f ðyÞ; 8 x; y A S ð140Þ 3 S r qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ ðx1  y1 Þ2 þ ðx2  y2 Þ2 is the distance between points x and y in circle S of the unit radius, f ðyÞ is a known function smooth enough in S, and uðxÞ is an unknown function defined in S. In the case of piecewise constant approximation, the BE equations have the form N

f ðyr Þ ¼  ∑ Fðyr ; xm Þuðxm Þ

ð141Þ

n¼1

where Z Fðyr ; xm Þ ¼

dS

2 2 3 m r m r Sm ððx1  y1 þ ξ1 Þ þ ðx2  y2 þ ξ2 Þ Þ

Here, ξ1 and ξ2 are local coordinates.

ð142Þ

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

Fig. 17. Solution for 156 rectangular and 312 triangular BEs.

In the event that yr and xm belong to different BEs, the integrals in (142) are regular and can be calculated using regular quadrature formulae or the regularized formulae presented here (127). In the event that yr ¼ xm , the integrals in (142) are hypersingular and must be calculated using the regularization methods presented here and using regularized formulae (127). In the case of piecewise constant approximation, the BEM equations have the form N

f ðyr Þ ¼  ∑

k

∑ F q ðyr ; xm Þfq ðξÞuðxm Þ

ð143Þ

n¼1q¼1

where k ¼ 3 for triangular BEs and k ¼ 4for quadratic BEs Z fq ðξÞ dS F q ðyr ; xm Þ ¼ 2 2 3 m r m r Ωm ððx1  y1 þ ξ1 Þ þ ðx2  y2 þ ξ2 Þ Þ

ð144Þ

In the event that yr and xm belong to different BEs, the integrals in (144) are regular and can be calculated using regular quadrature formulae or regularized formulae. In the event that yr and xm belong to the same BEs, the integrals in (144) are hypersingular and have to be calculated using the regularization methods presented here. The results of the numerical calculations of the normalized value uðxÞ for piecewise constant triangular and piecewise linear triangular and rectangular approximation are shown in Fig. 16 for 80 rectangular and 160 triangular BEs and in Fig. 17 for 156 rectangular and 312 triangular BEs, respectively. From these diagrams, it follows that the method for divergent integral regularization developed here gives good accuracy for piecewise linear and piecewise constant approximation. Our study shows that numerical solutions have good convergence, such that there is a less than 1% difference between the analytical and numerical solution for 100 rectangular and 200 triangular boundary elements.

9. Conclusions One of the main goals of this work is to demonstrate to the engineering and BEM communities that a generalized function based approach to the regularization of divergent integrals is not modern speculation but is instead a theoretically based and efficient computational tool. We clearly demonstrate that the theorems of Gauss-Ostrogradski and Green can be correctly applied to singular functions only in the context of generalized function theory. Through this and our other publications, it has been shown that divergent integral regularization, which is based on the theory of distribution, provides both a good mathematical foundation and an efficient tool for divergent integral calculations. A distribution based approach considers divergent integrals with

179

different singularities as a functional defined in special functional spaces. Using methods developed in [5], we obtain generalized Gauss-Ostrogradski and Green theorems, which are applicable in the case of singular functions. Using these theorems, divergent integrals can be transformed over any regular area into regular integrals over the boundary of the area. Thus, with the help of the second Green theorem, we have obtained a formula for the regularization of divergent integrals with any order of singularity over a regular area of any dimension. Using that formula in the same way, we can calculate weakly singular, singular and hypersingular integrals and also integrals with even higher singularity. In relatively simple cases, regularized formulae present only regular integrals over the contour of the area of integration, while in more complicated cases, regular integrals over the area may also be presented. Using the obtained formulae in the same way, we have considered 2-D weakly singular and hypersingular integrals, and we have obtained regular formulae for their calculation. The weakly singular and hypersingular integrals piecewise constant approximation has been considered for an arbitrary convex polygon and for rectangular and triangular BEs for piecewise linear approximation. It is important to mention that all calculations can be done analytically and that no numerical integration is needed for all of the equations developed here. The obtained regularization formulae can also be used for regular integral calculations. Such integrals appear when the collocation point moves to the BEs that do not belong to the area of integration. Mathematica software was used to verify the developed equations and to investigate their behavior for collocation points inside and outside the BE and near its contour. We compared our calculations with others reported in the literature and with numerical calculations of the regular integrals using Gaussian quadrature. In all considered cases, good agreement was observed. Our study shows that calculation of the regular integrals with our regularized formulae requires much less time than with Gaussian quadrature. Therefore, the regularized formulae presented here, and in other our publications, are recommended for calculating both divergent and regular integrals. We wish to note that the regularization formulae described here are analytical and can therefore be used to carry out asymptotic analysis and to analyze the behavior of the corresponding integrals when the collocation point moves to the boundary of the domain of integration.

Acknowledgments This work was supported by the Comity of Science and Technology of Mexico (CONACYT) by the Research Grants (Project no. 000000000101415), which is gratefully acknowledged. References [1] Aliabadi MH. The boundary element method. Application in solids and structures, Vol. 2. John Wiley & Sons, LTD; 2002. [2] Banerjee PK. Boundary element method in engineering science. New York, London: McGraw Hill; 1994. [3] Carley M. Numerical quadratures for singular and hypersingular integrals in boundary element methods. SIAM Journal on Scientific Computing 2007;29 (3):1207–16. [4] Carley M. Numerical quadratures for near-singular and near-hypersingular integrals in boundary element methods. Mathematical Proceedings of the Royal Irish Academy 2009;109A(1) (49(60)). [5] Courant R, Hilbert D. Methods of mathematical physics, vol. II. New York: Jonh Wiley&Sons; 1968. [6] Gel0 fand IM, Shilov GE. Generalized functions, Vol. 1. New York: Academic Press; 1964. [7] Chen G, Zhou J. Boundary Element Method. London: Academic Press; 1992; 646.

180

V.V. Zozulya / Engineering Analysis with Boundary Elements 40 (2014) 162–180

[8] Chen JT, Hong H-K. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. Applied Mechanics Review 1999;52(1):17–33. [9] Cheng AH-D, Cheng DT. Heritage and early history of the boundary element method. Engineering Analysis with Boundary Elements 2005;29:268–302. [10] Chien CC, Rajiyah H, Atluri SN. An effective method for solving the hypersingular integral equations in 3-D acoustics. Journal of the Acoustical Society of America 1990;88(2):918–37. [11] Chien CC, Rajiyah H, Atluri SN. On the evaluation of hyper-singular integrals arising in the boundary element method for linear elasticity. Computational Mechanics 1991;8:57–70. [12] Cruse TA. Numerical solution in three dimensional elastostatics. International Journal of Solids and Structures 1969;5:1259–74. [13] Cruse TA. An improved boundary-integral equation method for three dimensional elastic stress analysis. Computer & Structures 1974;4:741–54. [14] Cruse TA, Aithal R. Non-singular boundary integral equation implementation. International Journal for Numerical Methods in Engineering 1993;36:237–54. [15] Cruse TA, Suwito W. On the Somigliana stress identity in elasticity. Computational Mechanics 1993;11:1–10. [16] Guiggiani M. Formulation and numerical treatment of boundary integral equations with hypersingular kernels. In: Sladek V, Sladek J, editors. Singular Integrals in Boundary Element Methods. Southampton: WIT Press; 1998. [17] Guz AN, Zozulya VV. Brittle fracture of constructive materials under dynamic loading. Kiev: Naukova Dumka; 1993 (in Russian). [18] Guz AN, Zozulya VV. Fracture dynamics with allowance for a crack edges contact interaction. International Journal of Nonlinear Sciences and Numerical Simulation 2001;2(3):173–233. [19] Guz AN, Zozulya VV. Elastodynamic unilateral contact problems with friction for bodies with cracks. International Applied Mechanics 2002;38(8):3–45. [20] Hadamard J. Lectures on Cauchy0 s problem in linear hyperbolic differential equations. New York: Dover; 1953. [21] Hardy GH. Divergent series. Oxford: The Clarendon Press; 1949. [22] Hong H-K, Chen JT. Derivations of Integral Equations of Elasticity. Journal of Engineering Mechanics, ASCE 1988;114(6):1028–44. [23] Ioakimidis NI. Application of finite-part integrals to the singular integral equations of crack problems in plane and three-dimensional elasticity. Acta Mechanica 1982;45:31–47. [24] Kantorovich LV. On approximate calculation of some type of definite integrals and other applications of the method of singularities extraction. Matematicheskii Sbornik 1934;41(2):235–45 (in Russian). [25] Lifanov IK, Poltavskii LN, Vainikko GM. Hypersingular integral equations and their applications. London: Chapman & Hall/CRC; 2004. [26] Michlin SG. Singular integral equations. Uspehi Maematiceskih Nauk 1948;3 (25):29–112 (in Russian). [27] Michlin SG. Multidimensional singular integrals and integral equations. Oxford: Pergamon Press; 1965.

[28] Mukherjee S. CPV and HFP integrals and their applications in the boundary element method. International Journal of Solids and Structures 2000;37: 6623–34. [29] Mukherjee S. Finite parts of singular and hypersingular integrals with irregular boundary source points. Engineering Analysis with Boundary Elements 2000;24:767–76. [30] Mukherjeea S, Mukherjee YX, Wenjing Y. Cauchy principal values and finite parts of boundary integrals—revisited. Engineering Analysis with Boundary Elements 2005;29:844–9. [31] Mukherjee YX, Mukherjee S. Error analysis and adaptivity in threedimensional linear elasticity by the usual and hypersingular boundary contour method. International Journal of Solids and Structures 2001;38:161–78. [32] Qian ZY, Han ZD, Atluri SN. Directly Derived Non-Hyper-Singular Boundary Integral Equations for Acoustic Problems, and Their Solution through PetrovGalerkin Schemes. CMES: Computer Modeling in Engineering& Sciences 2004;5(6):541–62. [33] Sladek V, Sladek J, editors. Singular Integrals in Boundary Element Methods. Southampton: WIT Press; 1998. [34] Tanaka M, Sladek V, Sladek J. Regularization techniques applied to boundary element methods. Applied Mechanics Reviews 1994;47(10):457–99. [35] Zhang JD, Atluri SN. A boundary / interior element method for quasistatic and transient response analysis of shallow shells. Computers & Structures 1986;24:213–24. [36] Zozulya VV. Integrals of Hadamard type in dynamic problem of the crack theory. Doklady Academii Nauk. UkrSSR, Ser. A. Physical Mathematical & Technical Sciences 1991;2:19–22 (in Russian). [37] Zozulya VV. Regularization of the divergent integrals. I. General consideration. Electronic Journal of Boundary Elements 2006;4(2):49–57. [38] Zozulya VV. Regularization of the divergent integrals. II. Application in Fracture Mechanics. Electronic Journal of Boundary Elements 2006;4(2) (58–56). [39] Zozulya VV. Regularization of hypersingular integrals in 3-D fracture mechanics: Triangular BE, and piecewise-constant and piecewise-linear approximations. Engineering Analysis with Boundary Elements 2010;34(2):105–13. [40] Zozulya VV. Divergent Integrals in Elastostatics: Regularization in 3-D Case. CMES: Computer Modeling in Engineering & Sciences 2010;70(3):253–349. [41] Zozulya VV, Gonzalez-Chi PI. Weakly singular, singular and hypersingular integrals in elasticity and fracture mechanics. Journal of the Chinese Institute of Engineers 1999;22(6):763–75. [42] Zozulya VV, Gonzalez-Chi PI. Dynamic fracture mechanics with crack edges contact interaction. Engineering Analysis with Boundary Elements 2000;24 (9):643–59. [43] Zozulya VV, Lukin AN. Solution of three-dimensional problems of fracture mechanics by the method of integral boundary equations. International Applied Mechanics 1998;34(6):544–51. [44] Zozulya VV, Men’shikov VA. Solution of tree dimensional problems of the dynamic theory of elasticity for bodies with cracks using hypersingular integrals. International Applied Mechanics 2000;36(1):74–81.