Abstract

This article considers weakly singular, singular, and hypersingular integrals, which arise when the boundary integral equation methods are used to solve problems in elastostatics. The main equations related to formulation of the boundary integral equation and the boundary element methods in 2D and 3D elastostatics are discussed in details. For their regularization, an approach based on the theory of distribution and the application of the Green theorem has been used. The expressions, which allow an easy calculation of the weakly singular, singular, and hypersingular integrals, have been constructed.

1. Introduction

A huge amount of publications is devoted to the boundary integral equation methods (BIM) and its application in science and engineering [15]. One of the main problems arising in the numerical solution of the BIE by the boundary element method (BEM) is a calculation of the divergent integrals. In mathematics divergent integrals have established theoretical basis. For example, the weakly singular integrals are considered as improper integrals; the singular integrals are considered in the sense of Cauchy principal value (PV); the hypersingular integrals had been considered by Hadamard as finite part integrals (FP). Different methods have been developed for calculation of the divergent integrals. Usually different divergent integrals need different methods for their calculation. Analysis of the most known methods used for treatment of the different divergent integrals has been done in [3, 4, 68]. The correct mathematical interpretation of the divergent integrals with different singularities has been done in terms of the theory of distributions (generalized functions). The theory of distributions provides a unified approach for the study of the divergent integrals and integral operators with kernels containing different kind of singularities.

In our previous publications [919], approach based on the theory of distributions has been developed for regularization of the divergent integrals with different singularities. The above mentioned approach for regularization of the hypersingular integrals has been used for the first time in [11]. Then it was further developed for static and dynamic problems of fracture mechanics in [15, 18, 19], respectively, and also in [9, 10, 13]. The equations presented in [12, 14] can be applied for a wide class of divergent integral regularizations for example for 2D and 3D elasticity [5, 14].

In the present paper, the approach for the divergent integral regularization based on the theory of distribution and Green’s theorem is further developed. The equations that enable easy calculation of the weakly singular, singular, and hypersingular integrals for different shape functions are presented here.

2. Main Equations of Elastostatics

Lets consider a homogeneous, lineally elastic body, which in three-dimensional Euclidean space 3 occupies volume 𝑉 with smooth boundary 𝜕𝑉. The region 𝑉 is an open bounded subset of the three-dimensional Euclidean space 3 with a 𝐶0,1 Lipschitzian regular boundary 𝜕𝑉. The boundary contains two parts 𝜕𝑉𝑢 and 𝜕𝑉𝑝 such that 𝜕𝑉𝑢𝜕𝑉𝑝= and 𝜕𝑉𝑢𝜕𝑉𝑝=𝜕𝑉. On the part 𝜕𝑉𝑢 are prescribed displacement 𝑢𝑖(𝐱)of the body points, and on the part 𝜕𝑉𝑝 are prescribed traction 𝑝𝑖(𝐱), respectively. The body may be affected by volume forces 𝑏𝑖(𝐱). We assume that displacements of the body points and their gradients are small, so its stress-strain state is described by small strain deformation tensor 𝜀𝑖𝑗(𝐱). The strain tensor and the displacement vector are connected by Cauchy relations𝜀𝑖𝑗=12𝜕𝑖𝑢𝑗+𝜕𝑗𝑢𝑖,(2.1) where 𝜕/𝜕𝑥𝑖 is a derivative with respect to space coordinates 𝑥𝑖. The components of the strain tensor must also satisfy the Saint-Venant’s relations𝜕2𝑘𝑙𝜀𝑖𝑗𝜕2𝑖𝑙𝜀𝑘𝑗=𝜕2𝑘𝑗𝜀𝑖𝑙𝜕2𝑖𝑗𝜀𝑘𝑙.(2.2)

From the balance of impulse and the moment of impulse lows follow that the stress tensor is symmetric one and satisfy the equations of equilibrium𝜕𝑗𝜎𝑖𝑗+𝑏𝑖=0𝑖,𝐱𝑉.(2.3) Here and throughout the paper the summation convention applies to repeated indices.

The tensor of deformation 𝜀𝑖𝑗(𝐱) and stress 𝜎𝑖𝑗(𝐱) are related by Hook’s law𝜎𝑖𝑗=𝑐𝑖𝑗𝑘𝑙𝜀𝑖𝑗.(2.4) Here 𝑐𝑖𝑗𝑘𝑙 are elastic modules. In the case of homogeneous anisotropic medium they are symmetric 𝑐𝑖𝑗𝑘𝑙=𝑐𝑗𝑖𝑘𝑙=𝑐𝑘𝑙𝑖𝑗.(2.5) and satisfy condition of ellipticity𝑐𝑖𝑗𝑘𝑙𝜀𝑖𝑗𝜀𝑘𝑙𝛼1𝜀𝑖𝑗𝜀𝑖𝑗,𝜀𝑖𝑗,𝛼1>0.(2.6)

In the case of homogeneous isotropic medium, the elastic modules have the form 𝑐𝑖𝑗𝑘𝑙=𝜆𝛿𝑖𝑗𝛿𝑘𝑙𝛿+𝜇𝑖𝑘𝛿𝑗𝑙+𝛿𝑖𝑙𝛿𝑗𝑘,(2.7) where 𝜆 and 𝜇 are Lame constants, 𝜇>0 and 𝜆>𝜇, 𝛿𝑖𝑗 is a Kronecker’s symbol.

Substituting stress tensor in (2.3) and using Hook’s law (2.4) and Cauchy relations (2.1), we obtain the differential equations of equilibrium in the form of displacements which may be presented in the form𝐴𝑖𝑗𝑢𝑗+𝑏𝑖=0,𝐱𝑉.(2.8)

The differential operator 𝐴𝑖𝑗 for homogeneous anisotropic medium has the form𝐴𝑖𝑗=𝑐𝑖𝑘𝑗𝑙𝜕𝑘𝜕𝑙,(2.9) and for homogeneous isotropic medium has the form𝐴𝑖𝑗=𝜇𝛿𝑖𝑗𝜕𝑘𝜕𝑘+(𝜆+𝜇)𝜕𝑖𝜕𝑗.(2.10)

If the problem is defined in an infinite region, then solution of (2.8) must satisfy additional conditions at the infinity in the form 𝑢𝑗𝑟(𝐱)=𝑂1,𝜎𝑖𝑗𝑟(𝐱)=𝑂2for𝑟,(2.11) where 𝑟=𝑥21+𝑥22+𝑥23 is the distance in the three-dimensional Euclidian space.

If the body occupied a finite region 𝑉with the boundary 𝜕𝑉, it is necessary to establish boundary conditions. We consider the mixed boundary conditions in the form 𝑢𝑖(𝐱)=𝜑𝑖(𝐱),𝐱𝜕𝑉𝑢,𝑝𝑖(𝐱)=𝜎𝑖𝑗(𝐱)𝑛𝑗(𝐱)=𝑃𝑖𝑗𝑢𝑗(𝐱)=𝜓𝑖(𝐱),𝐱𝜕𝑉𝑝.(2.12)

The differential operator 𝑃𝑖𝑗𝑢𝑗𝑝𝑖 is called stress operator. It transforms the displacements into the tractions. For homogeneous anisotropic and isotropic medium, they have the forms 𝑃𝑖𝑗=𝑐𝑖𝑘𝑗𝑙𝑛𝑘𝜕𝑙,𝑃𝑖𝑗=𝜆𝑛𝑖𝜕𝑘𝛿+𝜇𝑖𝑗𝜕𝑛+𝑛𝑘𝜕𝑖,(2.13) respectively. Here 𝑛𝑖 are components of the outward normal vector, 𝜕𝑛=𝑛𝑖𝜕𝑖 is a derivative in direction of the vector 𝐧(𝐱) normal to the surface 𝜕𝑉𝑝.

3. Integral Representations for Displacements and Traction

In order to establish integral representations for the displacements and tractions, let us consider bilinear form which depends on two fields of the strain tensor that correspond to two fields of the displacements 𝐮 and 𝐮𝑎𝐮,𝐮=𝑐𝑖𝑗𝑘𝑙𝜀𝑖𝑗(𝐮)𝜀𝑘𝑙𝐮.(3.1) Obviously that 𝑎𝐮,𝐮𝐮=𝑎,𝐮,𝑎(𝐮,𝐮)=𝜎𝑖𝑗(𝐮)𝜀𝑖𝑗(𝐮)𝛼1𝜀𝑖𝑗(𝐮).(3.2)

Integrating the equality (3.1) over the volume 𝑉 and applying the Gauss-Ostrogradskii formula, we will obtain𝑉𝑎𝐮,𝐮𝑑𝑉=𝑉𝜎𝑖𝑗(𝐮)𝜀𝑘𝑙𝐮𝑑𝑉=𝜕𝑉𝜎𝑖𝑗𝑛𝑗𝑢𝑖𝑑𝑆𝑉𝑢𝑖𝜕𝑗𝜎𝑖𝑗𝑑𝑉.(3.3) Taking into account that 𝐴𝑖𝑗𝑢𝑗=𝜕𝑗𝜎𝑖𝑗 and 𝑝𝑖=𝜎𝑖𝑗𝑛𝑗=𝑃𝑖𝑗[𝑢𝑗], we will find the first Betti’s theorem in the form𝑉𝑢𝑖𝐴𝑖𝑗𝑢𝑗𝑑𝑉=𝑉𝑎𝐮,𝐮𝑑𝑉𝜕𝑉𝑢𝑖𝑃𝑖𝑗𝑢𝑗𝑑𝑆.(3.4) We will replace 𝑢𝑖 and 𝑢𝑖 in (3.4) and subtract resulting equation from (3.4). Because the form (3.1) is symmetrical one, we will obtain the second Betti’s theorem in the form𝑉𝑢𝑖𝐴𝑖𝑗𝑢𝑗𝑢𝑖𝐴𝑖𝑗𝑢𝑗𝑑𝑉=𝜕𝑉𝑢𝑖𝑃𝑖𝑗𝑢𝑗𝑢𝑖𝑃𝑖𝑗𝑢𝑗𝑑𝑆.(3.5) Taking into the account the definition of differential operator 𝐴𝑖𝑗 given in (2.8) and the definition of differential operator 𝑃𝑖𝑗 given in (2.12), we obtain relation 𝑉𝑏𝑖𝑢𝑖𝑏𝑖𝑢𝑖𝑑𝑉=𝜕𝑉𝑝𝑖𝑢𝑖𝑝𝑖𝑢𝑖𝑑𝑆,(3.6) which is called the Betti’s reciprocal theorem.

This theorem is usually used for obtaining integral representations for the displacements and traction vectors. To do that, we consider solution of the elliptic partial differential equation (2.8) in an infinite space for the body force 𝑏𝑖(𝐱)=𝛿𝑖𝑗𝛿(𝐱𝐲)𝐴𝑖𝑗𝑈𝑘𝑗(𝐱𝐲)+𝛿𝑘𝑖(𝐱𝐲)=0,𝐱,𝐲3.(3.7) Now considering that 𝑢𝑖(𝐱)=𝑈𝑖𝑗(𝐱𝐲),𝑝𝑖(𝐱)=𝑃𝑖𝑗𝑢𝑗(𝐱)=𝑊𝑖𝑗(𝐱,𝐲),(3.8)

from (3.6) we obtain the integral representation for the displacements vector𝑢𝑖(𝐲)=𝜕𝑉𝑝𝑖(𝐱)𝑈𝑗𝑖(𝐱𝐲)𝑢𝑗(𝐱)𝑊𝑗𝑖(𝐱,𝐲)𝑑𝑆+𝑉𝑏𝑖(𝐱)𝑈𝑗𝑖(𝐱𝐲)𝑑𝑉,(3.9) which is called Somigliana’s identity. The kernels 𝑈𝑗𝑖(𝐱𝐲) and 𝑊𝑗𝑖(𝐱,𝐲) are called fundamental solutions for elastostatics.

Applying to (3.9) the differential operator 𝑃𝑖𝑗, we will find integral representation for the traction in the form𝑝𝑖(𝐲)=𝜕𝑉𝑝𝑖(𝐱)𝐾𝑗𝑖(𝐱,𝐲)𝑢𝑗(𝐱)𝐹𝑗𝑖(𝐱,𝐲)𝑑𝑆+𝑉𝑏𝑖(𝐱)𝐾𝑗𝑖(𝐱,𝐲)𝑑𝑉.(3.10) The kernels 𝐾𝑗𝑖(𝐱,𝐲) and 𝐹𝑗𝑖(𝐱,𝐲) may be obtained applying the differential operator 𝑃𝑖𝑗 to the kernels 𝑈𝑗𝑖(𝐱𝐲) and 𝑊𝑗𝑖(𝐱,𝐲), respectively.

The integral representations (3.9) and (3.10) are usually used for direct formulation of the boundary integral equations in elastostatics.

4. Fundamental Solutions

In order to find the fundamental solutions 𝑈𝑖𝑗(𝐱𝐲) for the differential operator 𝐴𝑖𝑗, we consider the differential equations of elastostatics in the form displacements (3.7). Solutions of these equations are called the fundamental solutions.

In 3D elastostatics, they have the form𝑈𝑖𝑗1(𝐱𝐲)=16𝜋𝜇(1𝜐)𝑟(34𝜐)𝛿𝑖𝑗+𝜕𝑖𝑟𝜕𝑗𝑟.(4.1)

In 2D elastostatics they have the form𝑈𝑖𝑗1(𝐱𝐲)=8𝜋𝜇(1𝜐)(34𝜐)𝛿𝑖𝑗1ln𝑟+𝜕𝑖𝑟𝜕𝑗𝑟.(4.2) Here 𝑟=(𝑥1𝑦1)2+(𝑥2𝑦2)2+(𝑥3𝑦3)2 and 𝑟=(𝑥1𝑦1)2+(𝑥2𝑦2)2 for 3D and 2D case, respectively, 𝜕𝑖𝑟=𝜕𝑟/𝜕𝑥𝑖=𝜕𝑟/𝜕𝑦𝑖=(𝑥𝑖𝑦𝑖)/𝑟.

The kernels 𝑊𝑖𝑗(𝐱,𝐲) from (3.9) may be obtained by applying to 𝑈𝑖𝑗(𝐱𝐲) differential operator 𝑃𝑖𝑘[],(𝐱)=𝜆𝑛𝑖(𝐱)𝜕𝑘[]𝛿+𝜇𝑖𝑘𝑛𝑗(𝐱)𝜕𝑗[]+𝑛𝑘(𝐱)𝜕𝑖[],(4.3) as it is shown here𝑊𝑖𝑗(𝐱,𝐲)=𝜆𝑛𝑖(𝐱)𝜕𝑘𝑈𝑘𝑗(𝐱,𝐲)+𝜇𝑛𝑘𝜕(𝐱)𝑘𝑈𝑖𝑗(𝐱,𝐲)+𝜕𝑖𝑈𝑘𝑗(𝐱,𝐲).(4.4) Then after some transformations and simplifications, the expression for the kernels 𝑊𝑖𝑗(𝐱,𝐲) will have the following form𝑊𝑖𝑗(𝐱,𝐲)=14𝜋(1𝜐)𝑟𝛼𝑛𝑘(𝐱)𝜕𝑘𝑟(12𝜐)𝛿𝑖𝑗+𝛽𝜕𝑖𝑟𝜕𝑗𝑟𝑛+(12𝜐)𝑖(𝐱)𝜕𝑗𝑟𝑛𝑗(𝐱)𝜕𝑖𝑟.(4.5)

The kernels 𝐾𝑖𝑗(𝐱,𝐲) and 𝐹𝑖𝑗(𝐱,𝐲) from (3.10) may be obtained by applying differential operator 𝑃𝑖𝑘,(𝐲)=𝜆𝑛𝑖(𝐲)𝜕𝑘𝛿+𝜇𝑖𝑘𝑛𝑗(𝐲)𝜕𝑗+𝑛𝑘(𝐲)𝜕𝑖(4.6) to 𝑈𝑖𝑗(𝐱𝐲) and 𝑊𝑖𝑗(𝐱,𝐲), respectively𝐾𝑖𝑗(𝐱,𝐲)=𝜆𝑛𝑖(𝐲)𝜕𝑘𝑈𝑗𝑘(𝐱,𝐲)+𝜇𝑛𝑘𝜕(𝐲)𝑘𝑈𝑗𝑖(𝐱,𝐲)+𝜕𝑖𝑈𝑗𝑘,𝐹(𝐱,𝐲)𝑖𝑗(𝐱,𝐲)=𝜆𝑛𝑖(𝐲)𝜕𝑘𝑊𝑗𝑘(𝐱,𝐲)+𝜇𝑛𝑘𝜕(𝐲)𝑘𝑊𝑗𝑖(𝐱,𝐲)+𝜕𝑖𝑊𝑗𝑘.(𝐱,𝐲)(4.7) Then after some transformations and simplifications the expression for the kernels 𝐾𝑖𝑗(𝐱,𝐲) will have the form𝐾𝑖𝑗1(𝐱,𝐲)=4𝜋(1𝜐)𝑟𝛼𝑛𝑘(𝐲)𝜕𝑘𝑟(12𝜐)𝛿𝑖𝑗+𝛽𝜕𝑖𝑟𝜕𝑗𝑟𝑛+(12𝜐)𝑖(𝐲)𝜕𝑗𝑟𝑛𝑗(𝐲)𝜕𝑖𝑟,(4.8) and for the kernels, 𝐹𝑖𝑗(𝐱,𝐲) will have the form𝐹𝑖𝑗=𝜇(𝐱,𝐲)2𝛼𝜋(1𝜐)𝑟𝛽𝛽𝑛𝑘(𝐱)𝜕𝑘𝑟(12𝜐)𝑛𝑖(𝐲)𝜕𝑗𝛿𝑟+𝜐𝑖𝑗𝑛𝑘(𝐲)𝜕𝑘𝑟+𝑛𝑗(𝐲)𝜕𝑖𝑟𝛾𝑛𝑘(𝐲)𝜕𝑘𝑟𝜕𝑖𝑟𝜕𝑗𝑟𝑛+𝛽𝜈𝑖(𝐱)𝑛𝑘(𝐲)𝜕𝑘𝑟𝜕𝑗𝑟+𝑛𝑘(𝐱)𝑛𝑘(𝐲)𝜕𝑖𝑟𝜕𝑗𝑟+(12𝜐)𝛽𝑛𝑗(𝐱)𝑛𝑘(𝐲)𝜕𝑘𝑟𝜕𝑖𝑟+𝑛𝑘(𝐱)𝑛𝑘(𝐲)𝛿𝑖𝑗+𝑛𝑖(𝐱)𝑛𝑗(𝐲)(14𝜐)𝑛𝑗(𝐱)𝑛𝑖.(𝐲)(4.9) In (4.5)–(4.9), 𝛼=2,1, 𝛽=3,2 and 𝛾=5,4 in 3D and 2D cases, respectively.

The kernels (4.1)–(4.9) contain different kind singularities; therefore, corresponding integrals are divergent. Here we will investigate there singularities and develop methods of divergent integrals calculation.

5. Singularities, Boundary Properties, and Boundary Integral Equations

Simple observation shows that kernels in the integral representations (3.9) and (3.10) tend to infinity when 𝑟0. More detailed analysis of (4.1)–(4.9) gives us the following results.

In the 3D case with 𝐱𝐲,𝑈𝑖𝑗(𝐱𝐲)𝑟1,𝑊𝑖𝑗(𝐱,𝐲)𝑟2,𝐾𝑖𝑗(𝐱,𝐲)𝑟2,𝐹𝑖𝑗(𝐱,𝐲)𝑟3.(5.1)

In the 2D case with 𝐱𝐲,𝑈𝑖𝑗𝑟(𝐱𝐲)ln1,𝑊𝑖𝑗(𝐱,𝐲)𝑟1,𝐾𝑖𝑗(𝐱,𝐲)𝑟1,𝐹𝑖𝑗(𝐱,𝐲)𝑟2.(5.2) In order to investigate these functions and integrals with singular kernels, definition and classification of the integrals with various singularities will be presented here.

Definition 5.1. Let one considers two points with coordinated 𝐱,𝐲𝑛 (where 𝑛=3 or 𝑛=2) and region 𝑉 with smooth boundary 𝜕𝑉 of the class 𝐶0,1. The boundary integrals of the types 𝜕𝑉𝐺(𝐱,𝐲)𝑟𝛼𝜑(𝐱)𝑑𝑆,𝛼>0,(5.3) where 𝐺(𝐱,𝐲) is a finite function in 𝑛×𝜕𝑉 and 𝜑(𝐱) is a finite function in 𝜕𝑉, are weakly singular for 𝛼<𝑛1, singular for 𝛼=𝑛1, and strongly singular of hypersingular for 𝛼>𝑛1.

Definition 5.2. Let we consider two points with coordinated 𝐱,𝐲𝑛 (where 𝑛=3 or 𝑛=2) and region 𝑉 with smooth boundary 𝜕𝑉 of the class 𝐶0,1. The boundary integrals of the types 𝜕𝑉𝐺(𝐱,𝐲)ln(𝑟)𝜑(𝐲)𝑑𝐲,(5.4) where 𝐺(𝐱,𝐲) is a finite function in 𝑛×𝜕𝑉 and 𝜑(𝐱) is a finite function in 𝜕𝑉, are weakly singular.
The integrals with singularities can not be considered in usual (Riemann or Lebegue) sense. In order to such integrals have sense, it is necessary special consideration of them. We will apply this definition of the integrals from (3.9) and (3.10).

Definition 5.3. Integrals in (3.9) with kernels 𝑈𝑖𝑗(𝐱𝐲) are weakly singular and must be considered as improper 𝑊.𝑆.𝜕𝑉𝑝𝑖(𝐱)𝑈𝑖𝑗(𝐱𝐲)𝑑𝑆=lim𝜀0𝜕𝑉𝜕𝑉𝜀𝑝𝑖(𝐱)𝑈𝑖𝑗(𝐱𝐲)𝑑𝑆.(5.5) Here 𝜕𝑉𝜀 is a part of the boundary, projection of which on tangential plane is contained in the circle 𝐶𝜀(𝐱) of the radio 𝜀 with center in the point 𝐱.

Definition 5.4. Integrals in (3.9) and (3.10) with kernels 𝑊𝑖𝑗(𝐱,𝐲) and 𝐾𝑖𝑗(𝐱,𝐲) are singular and must be considered in sense of the Cauchy principal values as 𝑃.𝑉.𝜕𝑉𝑢𝑖(𝐱)𝑊𝑖𝑗(𝐱,𝐲)𝑑𝑆=lim𝜀0𝜕𝑉𝜕𝑉(𝑟<𝜀)𝑢𝑖(𝐱)𝑊𝑖𝑗(𝐱,𝐲)𝑑𝑆,𝑃.𝑉.𝜕𝑉𝑝𝑖(𝐱)𝐾𝑖𝑗(𝐱,𝐲)𝑑𝑆=lim𝜀0𝜕𝑉𝜕𝑉(𝑟<𝜀)𝑝𝑖(𝐱)𝐾𝑖𝑗(𝐱,𝐲)𝑑𝑆.(5.6) Here 𝜕𝑉(𝑟<𝜀) is a part of the boundary, projection of which on tangential plane is the circle 𝐶𝜀(𝐱) of the radio 𝜀 with center in the point 𝐱.

Definition 5.5. Integrals in (3.10) with kernels 𝐹𝑖𝑗(𝐱,𝐲) are hypersingular and must be considered in sense of the Hadamard finite part as 𝐹.𝑃.𝜕𝑉𝑢𝑖(𝐱)𝐹𝑗𝑖(𝐱𝐲)𝑑𝑆=lim𝜀0𝜕𝑉𝜕𝑉(𝑟<𝜀)𝑢𝑖(𝐱)𝑊𝑗𝑖(𝐱𝐲)𝑑𝑆+2𝑢𝑗𝑓(𝐱)𝑗(𝐱).𝜕𝑉(𝑟<𝜀)(5.7) Here functions 𝑓𝑗(𝐱) are chosen from the condition of the limit existence.
Singular character of the channels in (3.9) and (3.10) determine boundary properties of the corresponding potentials. Analysis of these formulae show that the boundary potentials, with the kernels 𝑈𝑖𝑗(𝐱𝐲), are weakly singular, and, therefore, they are continuous everywhere in the 𝑛 and, therefore, may be continuously extended on the boundary 𝜕𝑉. The potentials with the kernels 𝑊𝑖𝑗(𝐱,𝐲) and 𝐾𝑖𝑗(𝐱,𝐲) contain singular kernels, and they jump when crossing the boundary 𝜕𝑉. The potential with the kernels 𝐹𝑖𝑗(𝐱,𝐲) contain hypersingular kernels. They continuously cross the boundary 𝜕𝑉.
Boundary properties of these potentials have been studied extensively in [1, 2, 5, 8, 20] and so forth. Summary of these studies may be expressed by the equations 𝜕𝑉𝑝𝑖(𝐱)𝑈𝑗𝑖(𝐱𝐲)𝑑𝑆±=𝜕𝑉𝑝𝑖(𝐱)𝑈𝑗𝑖(𝐱𝐲)𝑑𝑆0,𝜕𝑉𝑢𝑖(𝐱)𝑊𝑗𝑖(𝐱,𝐲)𝑑𝑆±1=2𝑢𝑗(𝐲)+𝜕𝑉𝑢𝑖(𝐱)𝑊𝑗𝑖(𝐱,𝐲)𝑑𝑆0,𝜕𝑉𝑝𝑖(𝐱)𝐾𝑗𝑖(𝐱,𝐲)𝑑𝑆±1=±2𝑝𝑗(𝐲)+𝜕𝑉𝑝𝑖(𝐱)𝐾𝑗𝑖(𝐱,𝐲)𝑑𝑆0,𝜕𝑉𝑢𝑖(𝐱)𝐹𝑗𝑖(𝐱𝐲)𝑑𝑆±=𝜕𝑉𝑢𝑖(𝐱)𝐹𝑗𝑖(𝐱𝐲)𝑑𝑆0.(5.8) The symbols “±” and “” denote that two equalities, one with the top sign and the other with the bottom sign, are considered. The up index “0” points out that the direct value of the corresponding potentials on the surface  𝜕𝑉 should be taken.
Now using integral representations for displacements and traction and boundary properties of the potentials, we can get boundary integral equations for elastostatics. Tending 𝐲 in (3.9) and (3.10) to the boundary 𝜕𝑉 and taking into consideration boundary properties of the potentials (5.8), we obtain representation of the displacements and traction vectors on the boundary surface 𝜕𝑉. On the smooth parts of the boundary, they have the following form ±12𝑢𝑖(𝐲)=𝜕𝑉𝑝𝑗(𝐱)𝑈𝑖𝑗(𝐱𝐲)𝑢𝑗(𝐱)𝑊𝑖𝑗(𝐱,𝐲)𝑑𝑆+𝑉𝑝𝑗(𝐱)𝑈𝑖𝑗(1𝐱𝐲)𝑑𝑉,2𝑝𝑖(𝐲)=𝜕𝑉𝑝𝑗(𝐱)𝐾𝑖𝑗(𝐱,𝐲)𝑢𝑗(𝐱)𝐹𝑖𝑗(𝐱,𝐲)𝑑𝑆+𝑉𝑝𝑗(𝐱)𝐾𝑖𝑗(𝐱,𝐲)𝑑𝑉.(5.9) The plus and minus signs in these equations are used for the interior and exterior problems, respectively. Together with boundary conditions, they are used for compositing the BIE for the problems of elastostatics. The BIE are usually solved numerical transforming them into discrete system of finite dimensional equations. Such method is called the boundary element method (BEM).

6. Projection Method and the BEM Equations

The main idea of the BEM consists in approximation of the BIE and further solution of that approximated finite dimensional BE system of equations. The mathematical essence of this approach is the so-called projection method. Let us outline some results from mathematical theory of the projection methods related to the approximation of the BIE. For more information, one can refer to [4, 21].

We consider two Banach spaces 𝐗 and 𝐘 and functional equation in those spaces𝐀𝐮=𝐟,𝐮𝐷(𝐀)𝐗,𝐟𝑅(𝐀)𝐘.(6.1) Here 𝐀𝐗𝐘 is the linear operator mapping from Banach space 𝐗 in Banach space 𝐘, 𝐷(𝐀) is a domain, and 𝑅(𝐀) is a range of the operator 𝐀. The equation (6.1) is named the exact equation, and its solution is the exact solution. We denote 𝐿(𝐗,𝐘) Banach space of the linear operators mapping from 𝐗 in 𝐘.

Let in 𝐗 and 𝐘 act sequences of projection operators 𝐏 and 𝐏 such that 𝐏2=𝐏,𝐏𝐗=𝐗,𝐗𝐏𝐗,2=𝐏,𝐏𝐘=𝐘,𝐘𝐘,(6.2) where 𝐗 and 𝐘 are finite dimension subspaces of the Banach spaces 𝐗 and 𝐘; 𝑅1 is a parameter of discretization.

Now we consider operator 𝐀𝐿(𝐗,𝐘) mapping in finite dimensional subspaces 𝐗 and 𝐘 and an approximate equation𝐀𝐮=𝐟,𝐀=𝐏𝐀𝐏,𝐮=𝐏𝐮,𝐟=𝐏𝐟.(6.3) Solution 𝐮 of (6.3) is the approximate solutions of (6.1). The general scheme of the approached equations construction (6.3) is illustrated by the following diagram.726402.fig.001(6.4)

Now let us consider operator 𝐀𝐿(𝐗,𝐘) mapping in finite-dimensional subspaces 𝐗 and 𝐘 and an approximate equation.

Existence of the exact solution, convergence of the approximate solution to the exact one, and stability of the approximations are the main problems which arrived in application of the projection methods. In order to solve these problems, we have to formulate them mathematically.

We assume that projection operators 𝐏 and 𝐏 converge to identity operators in 𝐗 and 𝐘, respectively. It means thatlim0𝐏𝐮𝐮𝐗=0𝐮𝐗,lim0𝐏𝐟𝐟𝐘=0𝐟𝐘.(6.5)

Definition 6.1. Let conditions (6.5) be satisfied and stating from some =0>0 for any 𝐟𝐘, (6.3) has unique solution 𝐮. In this case if lim0𝐀𝐮𝐀𝐮𝐘=0,(6.6) then the solution of the approximate problem (6.3) converges to the exact solution (6.1). It means that the projective method presented on diagram (6.4) is applicable to the initial problem (6.1).

Definition 6.2. Let for some sequence of the operators {𝐀} mapping from 𝐗 into 𝐘 there is a constant 𝛾>0 such that stating from some =0>0𝐀𝐮𝐀𝐮𝐘𝐮𝛾𝐗𝐮𝐗(6.7) then for sequence of the operators {𝐀}, the condition of stability of the approximate solution is satisfied.
Conditions (6.5)–(6.7) are very important for formulation conditions of existence, convergence and stability of the approximate solution. These conditions contain the following theorem [21, 22].

Theorem 6.3. Let the following conditions be satisfied: (1)the projection operators 𝐏 and 𝐏 converge to identity operators in the Banach spaces 𝐗 and 𝐘, respectively, as it stated in (6.5);(2)the sequence of approximate operators {𝐀} converges to 𝐀 on each exact solution;(3)the condition of stability (6.7) is satisfied for the sequence of operators {𝐀}.
Then the following consequences take place: (1)the exact solution exists and it is unique;(2)for all enough small h exists a unique solution 𝐮𝐗 of the approximate equation (6.3);(3)the sequence of the approximate equation {𝐮} converges to the exact one and takes place in the estimation 𝐮𝐏𝐮𝐗𝛾1𝐏𝐀𝐮𝐀𝐏𝐮𝐘.(6.8)

Thus, using a projective method instead of the exact solution of (6.1) in functional space 𝐗, we can find the solution of the approximate equation (6.3) in finite dimensional space 𝐗. The functional spaces {𝐗,𝐗} and {𝐘,𝐘} are related by means of the projection operators 𝐏𝐿(𝐗,𝐗) and 𝐏𝐿(𝐘,𝐘), respectively. It is also important to construct inverse operator 𝐏1𝐿(𝐗,𝐗) which maps the finite dimensional space 𝐗 into the initial functional space 𝐗. Such operator refers to as the operator of interpolation. Because of 𝐗𝐗, the interpolation operator is not unique; moreover, for any two functional spaces 𝐗 and 𝐗, it is possible to construct infinite set of interpolation operators.

Let us apply the projection method to the BIE of elastostatics and construct corresponding finite dimensional BE equations. It is known [21, 23] that integral operators in (5.9) map between two functional spaces 𝐗(𝜕𝑉) and 𝐻1/2(𝜕𝑉) that are trace of displacements and traction on the boundary of the region in the following way: 𝐔𝑖𝑗𝐩𝑗=𝜕𝑉𝑈𝑖𝑗(𝐱,𝐲)𝑝𝑗(𝐖𝐱)𝑑𝑆𝐘(𝜕𝑉)𝐗(𝜕𝑉),𝑖𝑗𝐮𝑗=𝜕𝑉𝑊𝑖𝑗(𝐱,𝐲)𝑢𝑗𝐊(𝐱)𝑑𝑆𝐗(𝜕𝑉)𝐘(𝜕𝑉),𝑖𝑗𝐩𝑗=𝜕𝑉𝐾𝑖𝑗(𝐱,𝐲)𝑝𝑗𝐅(𝐱)𝑑𝑆𝐘(𝜕𝑉)𝐗(𝜕𝑉),𝑖𝑗𝐮𝑗=𝜕𝑉𝐹𝑖𝑗(𝐱,𝐲)𝑢𝑗(𝐱)𝑑𝑆𝐗(𝜕𝑉)𝐘(𝜕𝑉).(6.9)

We have to construct finite dimensional functional spaces that correspond to 𝐗(𝜕𝑉) and 𝐘(𝜕𝑉) and the corresponding projection operators. To construct finite dimensional functional spaces, we shall apply approximation by finite functions and splitting 𝜕𝑉 into finite elements 𝜕𝑉=𝑁𝑛=1𝜕𝑉𝑛,𝜕𝑉𝑛𝜕𝑉𝑘=,if𝑛𝑘.(6.10)

Because 𝜕𝑉 is the boundary of the region, these elements are called boundary elements. On each boundary element, we shall choose 𝑄 nodes of interpolation. Local projection operators act from functional 𝐗(𝜕𝑉𝑛) and 𝐘(𝜕𝑉𝑛) to the finite directional ones 𝐗𝑞(𝜕𝑉𝑛) and 𝐘𝑞𝜕𝑉𝑛)𝐏𝑢𝑞𝐗𝜕𝑉𝑛𝐗𝑞𝜕𝑉𝑛𝐱𝜕𝑉𝑛,𝐏𝑝𝑞𝐘𝜕𝑉𝑛𝐘𝑞𝜕𝑉𝑛𝐱𝜕𝑉𝑛.(6.11) Global projection operators are defined as the sum of the local projection operators 𝐏𝑢𝑛𝑞=𝑁𝑛=1𝐏𝑢𝑞,𝐏𝑝𝑛𝑞=𝑁𝑛=1𝐏𝑝𝑞.(6.12) They map 𝐗(𝜕𝑉) and 𝐘(𝜕𝑉) to finite dimensional interpolations spaces 𝐏𝑢𝑛𝑞𝐗(𝜕𝑉)𝐗𝑞𝑁𝑛=1𝜕𝑉𝑛𝐏𝐱𝜕𝑉,𝑝𝑛𝑞𝐘(𝜕𝑉)𝐘𝑞𝑁𝑛=1𝜕𝑉𝑛𝐱𝜕𝑉.(6.13) The local projection operators 𝐏𝑢𝑛 and 𝐏𝑝𝑛 establish correspondence between vectors of displacements and traction and their value on the nodes of interpolation of the boundary elements 𝜕𝑉𝑛𝐏𝑢𝑞𝐮𝑖𝐮(𝐱)=𝑛𝑖𝐱𝑞,𝑞=1,,𝑄𝐱𝜕𝑉𝑛,𝐏𝑝𝑞𝐩𝑖𝐩(𝐱)=𝑛𝑖𝐱𝑞,𝑞=1,,𝑄𝐱𝜕𝑉𝑛.(6.14) Similarly for the global operators, we have 𝐏𝑢𝑛𝑞𝐮𝑖𝐮(𝐱)=𝑛𝑖𝐱𝑞𝐏,𝑞=1,,𝑄;𝑛=1,,𝑁𝐱𝜕𝑉,𝑝𝑛𝑞𝐩𝑖𝐩(𝐱)=𝑛𝑖𝐱𝑞,𝑞=1,,𝑄;𝑛=1,,𝑁𝐱𝜕𝑉.(6.15) Let us construct local interpolation operators (𝐏𝑢𝑞)1 and (𝐏𝑝𝑞)1. For this purpose, we will introduce systems of shape functions 𝜑𝑛𝑞(𝐱) and 𝜓𝑛𝑞(𝐱) in the finite dimensional functional spaces 𝐗𝑞(𝜕𝑉𝑛) and 𝐘𝑞(𝜕𝑉𝑛). Then the vectors of displacements and traction on the boundary element 𝜕𝑉𝑛 will be represented approximately in the form 𝑢𝑖(𝐱)𝑄𝑞=1𝑢𝑛𝑖𝐱𝑞𝜑𝑛𝑞(𝐱),𝐱𝜕𝑉𝑛,𝑝𝑖(𝐱)𝑄𝑞=1𝑝𝑛𝑖𝐱𝑞𝜑𝑛𝑞(𝐱),𝐱𝜕𝑉𝑛,(6.16) and on the whole crack surface 𝜕𝑉 in the form 𝑢𝑖(𝐱)𝑁𝑄𝑛=1𝑞=1𝑢𝑛𝑖𝐱𝑞𝜑𝑛𝑞(𝐱),𝐱𝑁𝑛=1𝜕𝑉𝑛,𝑝𝑖(𝐱)𝑁𝑄𝑛=1𝑞=1𝑝𝑛𝑖𝐱𝑞𝜑𝑛𝑞(𝐱),𝐱𝑁𝑛=1𝜕𝑉𝑛.(6.17)

Finite-dimensional analogies for the integral operators (6.9) are operators which map the finite-dimensional functional spaces 𝐗𝑞(𝑁𝑛=1𝜕𝑉𝑛) and 𝐘𝑞(𝑁𝑛=1𝜕𝑉𝑛), from one to another 𝐔𝑛𝑞𝑖𝑗=𝐏𝑝𝑛𝑞𝐔𝑖𝑗𝐏𝑢𝑛𝑞𝐘𝑞𝜕𝑉𝑛𝐗𝑞𝜕𝑉𝑛,𝐖𝑛𝑞𝑖𝑗=𝐏𝑝𝑛𝑞𝐖𝑖𝑗𝐏𝑢𝑛𝑞𝐗𝑞𝜕𝑉𝑛𝐘𝑞𝜕𝑉𝑛,𝐊𝑛𝑞𝑖𝑗=𝐏𝑝𝑛𝑞𝐊𝑖𝑗𝐏𝑢𝑛𝑞𝐘𝑞𝜕𝑉𝑛𝐗𝑞𝜕𝑉𝑛,𝐅𝑛𝑞𝑖𝑗=𝐏𝑝𝑛𝑞𝐅𝑖𝑗𝐏𝑢𝑛𝑞𝐗𝑞𝜕𝑉𝑛𝐘𝑞𝜕𝑉𝑛.(6.18)

Note that in contrast to differential operators, the integral operators are global, and they are defined in the entire space, that is, at every boundary element.

Substitution of the expressions (6.17) in (2.8) gives us the finite-dimensional representations for the vectors of displacements and traction on the boundary in the form 12𝑢𝑚𝑖𝐲𝑟=𝑁𝑄𝑛=1𝑞=1𝑈𝑛𝑗𝑖𝐲𝑟,𝐱𝑞𝑝𝑛𝑗𝐱𝑞𝑊𝑛𝑗𝑖𝐱𝑟,𝐱𝑞𝑢𝑛𝑗𝐱𝑞+𝑈𝑖𝐟,𝐲,𝑉𝑛,12𝑝𝑚𝑖𝐲𝑟=𝑁𝑄𝑛=1𝑞=1𝐾𝑛𝑗𝑖𝐲𝑟,𝐱𝑞𝑝𝑛𝑗𝐲𝑞𝐹𝑛𝑗𝑖𝐲𝑟,𝐱𝑞𝑢𝑛𝑗𝐱𝑞+𝐾𝑖𝐟,𝐲,𝑉𝑛,(6.19) where 𝑈𝑛𝑗𝑖𝐲𝑟,𝐱𝑞=𝜕𝑉𝑛𝑈𝑗𝑖𝐲𝑟𝜓,𝐱𝑛𝑞(𝐱)𝑑𝑆,𝑊𝑛𝑗𝑖𝐲𝑟,𝐱𝑞=𝜕𝑉𝑛𝑊𝑗𝑖𝐲𝑟𝜑,𝐱𝑛𝑞(𝐾𝐱)𝑑𝑆,𝑛𝑗𝑖𝐲𝑟,𝐱𝑞=𝜕𝑉𝑛𝐾𝑗𝑖𝐲𝑟𝜓,𝐱𝑛𝑞(𝐱)𝑑𝑆,𝐹𝑛𝑗𝑖𝐲𝑟,𝐱𝑞=𝜕𝑉𝑛𝐹𝑗𝑖𝐲𝑟𝜑,𝐱𝑛𝑞(𝐱)𝑑𝑆.(6.20) The volume potentials 𝑈𝑖(𝐟,𝐲,𝑉𝑛) and 𝐾𝑖(𝐟,𝐲,𝑉𝑛) depend on discretization of the 𝑉 domain. More detailed information on transition from the BIE to the BEM equations can be found in [1, 2, 4, 5, 20].

7. Boundary Elements and Approximation

The BEM can be treated as the approximate method for the BIE solution, which includes approximation of the functions that belong to some functional space by discrete finite model. This model comprises finite number of values of the considered functions which are used for approximation of these functions by the shape functions determined on small subdomains called boundary elements. In this senses the BEM is closely related to a finite element method where the functions also belong to corresponding functional spaces and are approximated by finite model. Below we shall speak about finite element approximations and finite elements (FE), keeping in mind that boundary elements are their specific case.

It is important to mention that the local approximation of the considered function on one FE can be done independently from other FEs. It means that it is possible to approximate function on an FE by means of its values on the nodes independently of the place occupied considered the FE in the FE model and how behave the function on other FEs. Hence, it is possible to create the catalogue of various FE or BE with arbitrary node values interpolation function. Then from this catalogue can be chosen FEs which are necessary for approximation of the function and domain of its definition. The same FE can be used for discrete models of various functions or physical fields by determination of the necessary position of nodes in the model and further definition of the node values of the function or physical field. Thus, finite models of an area and its boundary do not not depend on functions and physical fields for which this area can be a domain of definition.

Let us consider how to construct an FE model of an area 𝑉𝑅𝑛 and a BE model of its boundary 𝜕𝑉𝑅𝑛1 (𝑛=2,3). We fix in the area 𝑉 finite number of points 𝐱𝑞(𝑔=1,,𝑄); these points refer to as global node points 𝑉(𝑞)={𝐱𝑞𝑉𝑔=1,,𝑄}. We shall divide the area 𝑉 into finite number of subareas 𝑉𝑛(𝑛=1,,𝑁) which are FEs. They have to satisfy the following conditions: 𝑉𝑛𝑉𝑚=,𝑚𝑛,𝑚,𝑛=1,2,,𝑁,𝑉=𝑁𝑛=1𝑉𝑛.(7.1) On each FE we introduce a local coordinate system 𝜉. The nodal points 𝐱𝑞𝑉𝑛 in the local system of coordinates we designate by 𝜉𝑞. They are coordinates of nodal points in the local coordinate system. Local and global coordinate are related in the following way:𝐱𝑞=𝑁𝑛=1𝚲𝑛𝜉𝑞𝑛.(7.2) Functions Λ𝑛 depend on position of the nodal points in the FE and BE. They join individual FE together in a FE model. Borders of the FEs and position of the nodal points should be such that, after joining together, separate elements form discrete model of the area 𝑉.

Having constructed FE model of the area 𝑉, we shall consider approximation of the function 𝑓(𝐱) that belong to some functional space. The FE model of the area 𝑉 is the domain of function which should be approximated. We denote function 𝑓(𝐱) on the FE 𝑉𝑛 by 𝑓𝑛(𝐱). Then𝑓(𝐱)=𝑁𝑛=1𝑓𝑛(𝐱).(7.3) On each FE, the local functions 𝑓𝑛(𝐱) may be represented in the form 𝑓𝑛(𝐱)𝑄𝑞=1𝑓𝑛(𝐱𝑞)𝜑𝑛𝑞(𝜉),(7.4) where 𝜑𝑛𝑞(𝜉) are interpolation polynomials or nodal functions of the FE with number 𝑛. In the nodal point with coordinates 𝐱𝑞, they are equal to 1, and in other nodal points are equal to zero. Taking into account (7.3) and (7.4), the global approximation of the function 𝑓(𝐱) looks like𝑓(𝐱)𝑁𝑄𝑛=1𝑞=1𝑓𝑛𝐱𝑞𝜑𝑛𝑞(𝜉).(7.5) If the nodal point 𝑞 belongs to several Fes, it is considered in these sums only once.

The FE and BE elements can be of various form and sizes, their surfaces can be curvilinear. The curvilinear FEs are very important in BEM because the boundary surface is usually curvilinear. But it is more convenient to use standard FE, whose surfaces coincide with coordinate planes of the local coordinates system. Mathematically it means that it is necessary to establish relation between local coordinates 𝜉𝑖, in which the element has a simple appearance, and global 𝑥𝑖 where the FE represents more complex figure. Local coordinates 𝜉𝑖 should be functions of global (𝜉𝑖(𝑥1,𝑥2,𝑥3)) ones, and on the contrary, global coordinates should be functions of (𝑥𝑖(𝜉1,𝜉2,𝜉3))ones. In order to these maps be one-to-one, it is necessary and sufficient that the Jacobians of transformations be nonzero||||𝐽=det𝜕𝑥𝑖𝜕𝜉𝑗||||0,𝐽1||||=det𝜕𝜉𝑖𝜕𝑥𝑗||||0.(7.6) The differential elements along coordinate axes are related by 𝑑𝑥𝑖=𝜕𝑥𝑖𝜕𝜉𝑗𝑑𝜉𝑗,𝑑𝐱=𝐽(𝜉)𝑑𝜉,𝑑𝜉𝑖=𝜕𝜉𝑖𝜕𝑥𝑗𝑑𝑥𝑗,𝑑𝜉=𝐽1(𝐱)𝑑𝐱.(7.7) The volume element in the 3 is transformed under the formula𝑑𝑉=𝑑𝑥1𝑑𝑥2𝑑𝑥3=𝐽(𝜉)𝑑𝜉1𝑑𝜉2𝑑𝜉3,(7.8) and the area element in the 2 is transformed under the formula𝑑𝐴=𝑑𝑥1𝑑𝑥2||||=det𝜕𝑥𝛼𝜕𝜉𝛽||||𝑑𝜉1𝑑𝜉2,𝛼,𝛽=1,2.(7.9) The differential of the surface located in the 3 is defined by the expression𝑛𝑑𝑆=21+𝑛22+𝑛231/2𝑑𝜉1𝑑𝜉2,(7.10) where𝑛1=𝜕𝑥1𝜕𝜉1𝜕𝑥3𝜕𝜉2𝜕𝑥2𝜕𝜉2𝜕𝑥3𝜕𝜉1,𝑛2=𝜕𝑥3𝜕𝜉1𝜕𝑥1𝜕𝜉2𝜕𝑥1𝜕𝜉1𝜕𝑥3𝜕𝜉2,𝑛3=𝜕𝑥1𝜕𝜉1𝜕𝑥2𝜕𝜉2𝜕𝑥2𝜕𝜉1𝜕𝑥1𝜕𝜉2.(7.11) The element of length of a contour in the 2 is defined by the expression𝑑𝑙=𝑑𝑥1𝑑𝜉12+𝑑𝑥2𝑑𝜉121/2𝑑𝜉1.(7.12)

It is important to point attention that it is quite enough to consider standard FE which can be transformed to the necessary form by suitable transformation of coordinates. The FE approximation has to be linear independent and compact in the corresponding functional space.

We consider here some examples of the FE approximation, which are frequently used in the BEM.

7.1. Piecewise Constant Approximation

The piecewise constant approximation is the simplest one. Interpolation functions in this case do not depend on the FE form and the dimension of the domain. They have the form𝜑𝑛𝑞(𝐱)=1𝐱𝑉𝑛,0𝐱𝑉𝑛.(7.13) This approximation is linear independent and compact in the functional spaces in which it is not necessary to approximate the derivative of the function.

7.2. Piecewise Linear Approximation

Interpolation functions in this case change linearly inside the FE and depend on its form and dimension of the domain. Therefore, interpolation functions are called shape functions.

In 1D case, the FEs are linear, and their shape functions are 𝜑𝑛𝑞=121+𝜉𝑞1𝜉1,𝜉𝑞𝑖=±1.(7.14)

In 2D case for the rectangular FE, shape functions are𝜑𝑛𝑞=141+𝜉𝑞1𝜉11+𝜉𝑞2𝜉2,𝜉𝑞𝑖=±1.(7.15)

For the triangular FE, it is convenient to introduce system of coordinates which connected with the Cartesian by relation 𝑥𝑖=𝜉𝑞𝑥𝑞𝑖,𝜉𝑞=1,𝑖=1,2,𝑞=1,2,3.(7.16) For the triangular FE, the interpolation functions are 𝜑𝑛𝑞=𝜉𝑞.(7.17) The local coordinates 𝜉𝑞 are also called area coordinates.

7.3. Piecewise Quadratic Approximation

Interpolation functions in this case change quadratically inside the FE and depend on its form and dimension of the domain.

In 1D case, the FEs are bilinear, and their shape functions are 𝜑𝑛𝑞=12𝜉𝑞1𝜉11+𝜉𝑞1𝜉1(𝑞=1,2,3)(7.18) at the end nodes and 𝜑𝑛𝑞=1𝜉21(7.19) at the midpoint node.

In 2D case for the rectangular FE, shape functions are𝜑𝑛𝑞=141+𝜉𝑞1𝜉11+𝜉𝑞2𝜉2𝜉𝑞1𝜉1𝜉𝑞2𝜉21(𝑞=1,,8)(7.20) at the angular nodes and 𝜑𝑛𝑞=121𝜉211+𝜉𝑞2𝜉2𝜉𝑞1𝜑=0,𝑛𝑞=121𝜉221+𝜉𝑞1𝜉1𝜉𝑞2=0(7.21) at the side nodes.

For the triangular FE, the interpolation functions are 𝜑𝑛𝑞=2𝜉𝑞𝜉1𝑞(𝑞=1,,6)(7.22) at the angular nodes and 𝜑𝑛𝑞=4𝜉1𝜉2.(7.23) at the side nodes.

In order to construct system of algebraic equation for the BEM (6.18), we have to represent global coordinates as function of local ones. The most convenient way to do that is to use the interpolation function 𝜑𝑛𝑞 defined in (7.13)–(7.23). In this case, the global coordinates have the form𝑥𝑖=𝑥𝑞𝑖𝜑𝑛𝑞(𝜉).(7.24) Here 𝑥𝑞𝑖 is the global coordinate of nodal points for the FE number 𝑛.

8. Regularization of the Divergent Integrals

As it was mentioned above in order to solve the BIE (5.9) numerically, we have to transform them to the finite dimensional BE equations (6.18). In order to do that transform, we have to calculate integrals (6.19). One of the main problems that occur in this situation is the presence of the divergent integrals. They can not be calculated in traditional way, for example, numerically using quadrature formulas. For example, the integrals with the kernels 𝑈𝑖𝑗(𝐱𝐲) are weakly singular (WS). They have to be considered as improper integrals. The integrals with the kernels 𝑊𝑖𝑗(𝐱,𝐲) and 𝐾𝑖𝑗(𝐱,𝐲) are singular. They have to be considered in the sense of Cauchy as principal values (PVs). The integrals with the kernels 𝐹𝑖𝑗(𝐱,𝐲) are hypersingular. They have to be considered in the sense of Hadamard as finite parts (FPs). Traditional approach to the divergent integral calculation may be found in [1, 2, 6, 8]. Following [1116], we will develop here and apply the method of the divergent integral calculation based on the theory of distribution. This approach consists in application of the second Green theorem and transformation of divergent integrals into regular ones. In [12] we have developed formulas for regularization of the divergent integrals 𝑟𝑚 in the form𝐹.𝑃.𝑎𝑎𝜑(𝑥)𝑟𝑚𝑑𝑥=𝑘1𝑖=0(1)𝑖+1𝑑𝑖𝑑𝑥𝑖𝑃𝑖𝑟𝑚𝑘𝑑𝑘1𝑖𝜑(𝑥)𝑑𝑥𝑘1𝑖||||𝑥=𝑎𝑥=𝑎+(1)𝑘𝑎𝑎𝑃𝑘𝑟𝑚𝑘𝑑𝑘𝜑(𝑥)𝑑𝑥𝑘,𝐹.𝑃.𝑉𝜑(𝐱)𝑟𝑚𝑑𝑉=𝑘1𝑖=0(1)𝑖+1𝜕𝑉Δ𝑘𝑖1𝜑(𝐱)𝜕𝑛𝑃𝑖𝑟𝑚2𝑖𝑃𝑖𝑟𝑚2𝑖𝜕𝑛Δ𝑘𝑖1𝜑(𝐱)𝑑𝑆+(1)𝑘𝑉1𝑟𝑚2𝑘Δ𝑘+1𝜑(𝐱)𝑑𝑉,(8.1) for 1D and 2D cases, respectively. We will demonstrate here how this approach works in the problems of elastostatics.

8.1.  1D Divergent Integrals

In 2D elastostatics after introducing local system of coordinates and simplification, all divergent integrals can be presented in the form 𝐽0=𝑏𝑎1𝜑(𝑥)ln𝑥𝑑𝑥,𝐽𝑘=𝑏𝑎𝜑(𝑥)𝑥𝑘𝑑𝑥,𝑘=1,2(8.2) Here 𝜑(𝑥) is the smooth function that depend on the shape of the BE and the interpolation polynomials.

We consider first the weakly singular integral 𝐽0. Because of logarithmical singularity, we can not use formula (8.1). Therefore, we start from the formula for integration by parts in the form𝑏𝑎𝑑𝑔(𝑥)||𝑑𝑥𝜑(𝑥)𝑑𝑥=𝜑(𝑥)𝑔(𝑥)𝑏𝑎𝑏𝑎𝑑𝜑(𝑥)𝑑𝑥𝑔(𝑥)𝑑𝑥.(8.3) Let in this formula be 𝑔(𝑥)=𝑥ln1/𝑥, 𝑑𝑔(𝑥)/𝑑𝑥=ln1/𝑥, then we obtain 𝐽0=𝑊.𝑆.𝑏𝑎1𝜑(𝑥)ln𝑥1𝑑𝑥=𝜑(𝑥)𝑥ln𝑥|||𝑏𝑎𝑏𝑎𝑑𝜑(𝑥)1𝑑𝑥𝑥ln𝑥𝑑𝑥.(8.4) Obviously the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation 𝜑(𝑥)=1, and we get 𝐽0=𝑊.𝑆.𝑏𝑎||||||||||||ln𝑥𝑦𝑑𝑥=(𝑏𝑦)ln𝑏𝑦(𝑎𝑦)ln𝑎𝑦(𝑎<𝑦<𝑏).(8.5)

For the singular integral 𝐽1 regularization, we will also use the formula for integration by parts (8.3). Let in this formula 𝑔(𝑥)=ln1/𝑥, (𝑑𝑔(𝑥))/𝑑𝑥=1/𝑥, then we obtain 𝐽1=𝑃.𝑉.𝑏𝑎𝜑(𝑥)𝑥𝑑𝑥=𝑑𝜑(𝑥)1𝑑𝑥𝑥ln𝑥1𝜑(𝑥)ln𝑥||||𝑏𝑎𝑏𝑎𝑑2𝜑(𝑥)𝑑𝑥21𝑥ln𝑥𝑑𝑥.(8.6) Here also the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation 𝜑(𝑥)=1, and we get 𝐽1=𝑃.𝑉.𝑏𝑎𝑑𝑥||||𝑥𝑦=ln𝑏𝑦||||,𝑎𝑦(𝑎<𝑦<𝑏).(8.7)

Finally for the hypersingular integral 𝐽2 regularization, we will use formula (8.1) and the above-obtained result for singular 𝐽1 regularization. Finally we get𝐽2=𝐹.𝑃.𝑏𝑎𝜑(𝑥)𝑥2𝑑𝑑𝑥=2𝜑(𝑥)𝑑𝑥21𝑥ln𝑥𝑑𝜑(𝑥)1𝑑𝑥ln𝑥𝜑(𝑥)𝑥||||𝑏𝑎𝑏𝑎𝑑3𝜑(𝑥)𝑑𝑥31𝑥ln𝑥𝑑𝑥.(8.8) Here also the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation 𝜑(𝑥)=1, and we get 𝐽2=𝐹.𝑃.𝑏𝑎𝑑𝑥(𝑥𝑦)21=+1(𝑏𝑦)(𝑎𝑦)(𝑎<𝑦<𝑏).(8.9) From this equation can be found Hadamard’s example of a function that is positive every were in the integration region, but its integral is a negative one 𝐹.𝑃.𝑎𝑎𝑑𝑦𝑦22=𝑎,𝑎>0.(8.10)

8.2.  2D Divergent Integrals

In 3D elastostatics after introducing local system of coordinates and simplification, all divergent integrals can be presented in the form 𝐽𝑘𝑙,𝑚=𝑆𝑛𝑥𝑙1𝑥𝑚2𝑟𝑘𝜑(𝑥)𝑑𝑆,𝑙,𝑚=0,1,2,𝑘=3,4,5.(8.11) Here 𝜑(𝑥) is the smooth function that depends on shape of the BE and interpolation polynomials.

8.2.1. Integrals with Kernels of the Type 𝑟𝑘,𝑘=1,2,3

From (8.1) for 𝑘=1, we get regularization for weakly singular integral 𝐽10,0=𝑊.𝑆.𝑉𝜑(𝐱)𝑟𝑑𝑉=𝜕𝑉𝑟𝜑(𝐱)𝑛2𝑟𝑟𝜕𝑛𝜑(𝐱)𝑑𝑆+𝑉𝑟Δ𝜑(𝐱)𝑑𝑉.(8.12) Here 𝑟𝑛=(𝑥𝛼𝑦𝛼)𝑛𝛼, and the summation convention applies to repeated indices 𝛼=1,2. The integral on the left is divergent and on the right are regular ones.

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate this integral analytically. Introducing polar, coordinates we will get 𝐽10,0=12𝜕𝑆𝑛𝑟𝑛𝑟𝑑𝑙=02𝜋𝑟𝑟𝑟𝑑𝜑=𝜋𝑟.(8.13) In order to regularize singular integral, we will use, relation 1/𝑟2=(1/2)Δ(ln𝑟)2. Then in the same way we get𝑉𝜑(𝐱)𝑟21𝑑𝑉=2𝜕𝑉𝜑(𝐱)2𝑟𝑛ln𝑟𝑟2(ln𝑟)2𝜕𝑛1𝜑(𝐱)𝑑𝑆+2𝑉(ln𝑟)2Δ𝜑(𝐱)𝑑𝑉.(8.14) The volume integral on the right is weakly singular. Taking into account the relation (ln𝑟)2=(𝑟2/6)Δ(ln𝑟)4, we obtain regularization for this weakly singular integral𝑉(ln𝑟)21Δ𝜑(𝐱)𝑑𝑉=6𝜕𝑉2Δ𝜑(𝐱)𝑟𝑛(ln𝑟)3𝑟2(ln𝑟)4𝜕𝑛+1Δ𝜑(𝐱)𝑑𝑆6𝑉𝑟2(ln𝑟)4Δ2𝜑(𝐱)𝑑𝑉.(8.15)

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate singular integral in (8.14) analytically. Introducing polar coordinates, we will get 𝐽20,0=𝜕𝑆𝑛𝑟𝑛ln𝑟𝑟2𝑑𝑙=02𝜋𝑟ln𝑟𝑟2𝑟𝑑𝜑=2𝜋ln𝑟.(8.16) Finally from (8.1) for 𝑘=3, we get regularization for hupersingular integral 𝑉𝜑(𝐱)𝑟3𝑑𝑉=𝜕𝑉𝑟Δ𝜑(𝐱)𝑛𝑟2𝑟𝜑(𝐱)𝑛𝑟31𝑟𝜕𝑛𝜑(𝐱)𝑟𝜕𝑛Δ𝜑(𝐱)𝑑𝑆+𝑉𝑟Δ2𝜑(𝐱)𝑑𝑉.(8.17)

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get 𝐽30,0=𝜕𝑆𝑛𝑟𝑛𝑟3𝑑𝑙=02𝜋𝑟𝑟3𝑟𝑑𝜑=2𝜋𝑟.(8.18)

Now using (8.1), any divergent integral with the kernels of the type 1/𝑟𝑘 for any positive integer 𝑘 can be calculated.

8.2.2. Integrals with Kernels of the Type 𝑥2𝛼/𝑟𝑘,𝑘=3,4,5

The weakly singular integral with kernel 𝑥2𝛼/𝑟3 is calculated taking into account the equation 𝑥2𝛼/𝑟3=(1/3)((2/𝑟)Δ(𝑥2𝛼/𝑟)). It is easy to show that𝐽32,0=𝑊.𝑆.𝑉𝑥𝜑(𝐱)21𝑟32𝑑𝑉=3𝑊.𝑆.𝑉𝜑(𝐱)𝑟1𝑑𝑉3𝑊.𝑆.𝑉𝑥𝜑(𝐱)Δ21𝑟𝑑𝑉.(8.19) The first integral here is already calculated in (8.12). The second one may be presented in the form𝑉𝑥𝜑(𝐱)Δ21𝑟𝑑𝑉=𝜕𝑉𝜑(𝐱)2𝑛1𝑥1𝑟𝑥12𝑟𝑛𝑟3𝑥21𝑟𝜕𝑛𝜑(𝐱)𝑑𝑆+𝑉𝑥21𝑟Δ𝜑(𝐱)𝑑𝑉.(8.20) Combining the last two equations, finally we will get 𝐽32,0=13𝜕𝑉𝜑(𝐱)2𝑛1𝑥1𝑟𝑥12𝑟𝑛𝑟3+𝑟𝑛𝑥2𝑟21𝑟𝜕+𝑟𝑛+1𝜑(𝐱)𝑑𝑆3𝑉𝑥21𝑟+𝑟Δ𝜑(𝐱)𝑑𝑉.(8.21)

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area we can calculate this integral analytically. Introducing polar coordinates, we will get 𝐽32,0=13𝜕𝑆𝑛𝑥21𝑟𝑛𝑟32𝑥1𝑛1𝑟+2𝑟𝑛𝑟=1𝑑𝑙302𝜋(𝑟cos𝜑)2𝑟𝑟3𝑟𝑑𝜑02𝜋2𝑟(cos𝜑)2𝑟𝑟𝑑𝜑+202𝜋𝑟𝑟𝑟𝑑𝜑=𝜋𝑟.(8.22)

The singular integral with kernel 𝑥2𝛼/𝑟4 is calculated taking into account that Δ𝑥2𝛼/𝑟4=(1/4)((2/𝑟2)Δ𝑥2𝛼/𝑟2). In this case,𝐽42,0=𝑃.𝑉.𝑉𝑥𝜑(𝐱)21𝑟21𝑑𝑉=4𝑉2𝜑(𝐱)𝑟2𝑥Δ21𝑟21𝑑𝑉=2𝑉𝜑(𝐱)𝑟21𝑑𝑉4𝑉𝑥𝜑(𝐱)Δ21𝑟2𝑑𝑉.(8.23) The first integral here is already calculated in (8.14). The second one may be presented in the form𝑉𝑥𝜑(𝐱)Δ12𝑟2𝑑𝑉=𝜕𝑉𝜑(𝐱)2𝑛1𝑥1𝑟22𝑥12𝑟𝑛𝑟4𝑥12𝑟2𝜕𝑛𝜑(𝐱)𝑑𝑆+𝑉𝑥12𝑟2Δ𝜑(𝐱)𝑑𝑉.(8.24) Combining the last two equations, finally we will get 𝐽42,0=12𝜕𝑉𝑥𝜑(𝐱)12𝑟𝑛𝑟4𝑛1𝑥1𝑟2+𝑟𝑛ln𝑟𝑟2(ln𝑟)22𝑥122𝑟2𝜕𝑛+1𝜑(𝐱)𝑑𝑆4𝑉(ln𝑟)2𝑥12𝑟2Δ𝜑(𝐱)𝑑𝑉.(8.25)

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get 𝐽42,0=12𝜕𝑆𝑛𝑥12𝑟𝑛𝑟4𝑥1𝑛1𝑟2+𝑟𝑛ln𝑟𝑟2𝑑𝑙=02𝜋(𝑟cos𝜑)2𝑟𝑟4𝑟𝑑𝜑02𝜋𝑟(cos𝜑)2𝑟2+𝑟𝑑𝜑02𝜋𝑟ln𝑟𝑟2𝑟𝑑𝜑=2𝜋ln𝑟.(8.26)

The hypersingular integral with kernel 𝑥2𝛼/𝑟5 is calculated taking into account that 𝑥2𝛼/𝑟5=(1/3)((2/𝑟3)Δ(𝑥2𝛼/𝑟3)). In this case, it is easy to show that𝐽52,0=𝐹.𝑃.𝑉𝑥𝜑(𝐱)21𝑟51𝑑𝑉=3𝑉2𝜑(𝐱)𝑟3𝑥Δ21𝑟32𝑑𝑉=3𝑉𝜑(𝐱)𝑟31𝑑𝑉3𝑉𝑥𝜑(𝐱)Δ21𝑟3𝑑𝑉.(8.27) The first integral here is already calculated in (8.17). The second one may be presented in the form𝑉𝑥𝜑(𝐱)Δ21𝑟3𝑑𝑉=𝜕𝑉𝜑(𝐱)2𝑛1𝑥1𝑟33𝑥12𝑟𝑛𝑟5𝑥21𝑟3𝜕𝑛𝜑(𝐱)𝑑𝑆+𝑉𝑥21𝑟3Δ𝜑(𝐱)𝑑𝑉.(8.28) Combining the last two equations, finally we will get 𝐽52,0=23𝜕𝑉𝑛𝜑(𝐱)1𝑥1𝑟3𝑟𝑛𝑟33𝑥12𝑟𝑛2𝑟51𝑟+𝑥21𝑟3𝜕𝑛+2𝜑(𝐱)𝑑𝑆3𝑊.𝑆.𝑉1𝑟𝑥21𝑟3Δ𝜑(𝐱)𝑑𝑉.(8.29) The volume integral here is weakly singular. Its regularization may be done using (8.12) and (8.21). For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get 𝐽52,0=𝜕𝑆𝑛2𝑟𝑛3𝑟3+2𝑥12𝑟𝑛3𝑟5𝑥1𝑛1𝑟32𝑑𝑙=302𝜋𝑟𝑟32𝑟𝑑𝜑+302𝜋(𝑟cos𝜑)2𝑟𝑟5𝑟𝑑𝜑02𝜋𝑟(cos𝜑)2𝑟3𝜋𝑟𝑑𝜑=𝑟.(8.30)

8.2.3. Integrals with Kernels of the Type (𝑥1𝑥2)/𝑟𝑘,𝑘=3,4,5

The weakly singular integral with kernel (𝑥1𝑥2)/𝑟3 is calculated using (8.1). Taking into account that (𝑥1𝑥2)/𝑟3=1/3Δ(𝑥1𝑥2)/𝑟, it is easy to show that𝐽31,1=𝑊.𝑆.𝑉𝑥𝜑(𝐱)1𝑥2𝑟31𝑑𝑉=3𝜕𝑉𝑥𝜑(𝐱)1𝑥2𝑟𝑛𝑟3𝑟𝑟+𝑥1𝑥2𝑟𝜕𝑛1𝜑(𝐱)𝑑𝑆3𝑉𝑥1𝑥2𝑟Δ𝜑(𝐱)𝑑𝑉.(8.31) Here 𝑟=𝑥1𝑛2+𝑥2𝑛1.

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get 𝐽31,1=13𝜕𝑆𝑛𝑥1𝑥2𝑟𝑛𝑟3𝑟𝑟𝑑𝑙=02𝜋𝑟3cos𝜑sin𝜑𝑟3𝑟𝑑𝜑02𝜋2𝑟cos𝜑sin𝜑𝑟𝑟𝑑𝜑=0.(8.32)

The singular integral with kernel 𝑥1𝑥2/𝑟4 is calculated using equation (8.1). Taking into account that 𝑥1𝑥2/𝑟4=(1/4)Δ𝑥1𝑥2/𝑟2, it is easy to show that𝐽41,1=𝑃.𝑉.𝑉𝑥𝜑(𝐱)1𝑥2𝑟41𝑑𝑉=4𝜕𝑉𝜑(𝐱)2𝑥1𝑥2𝑟𝑛𝑟4𝑟𝑟2𝑥1𝑥2𝑟2𝜕𝑛1𝜑(𝐱)𝑑𝑆4𝑉𝑥1𝑥2𝑟2Δ𝜑(𝐱)𝑑𝑉.(8.33)

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get 𝐽41,1=14𝜕𝑆𝑛2𝑥1𝑥2𝑟𝑛𝑟4𝑟𝑟21𝑑𝑙=402𝜋2𝑟3cos𝜑sin𝜑𝑟41𝑟𝑑𝜑402𝜋2𝑟cos𝜑sin𝜑𝑟2𝑟𝑑𝜑=0.(8.34)

The hypersingular integral with kernel 𝑥1𝑥2/𝑟5 is calculated using (8.1). Taking into account that 𝑥1𝑥2/𝑟5=(1/3)Δ𝑥1𝑥2/𝑟3, it is easy to show that𝐽51,1=𝐹.𝑃.𝑉𝑥𝜑(𝐱)1𝑥2𝑟5𝑑𝑉=𝜕𝑉𝑥𝜑(𝐱)1𝑥2𝑟𝑛𝑟5𝑟3𝑟3𝑥1𝑥2𝑟3𝜕𝑛1𝜑(𝐱)𝑑𝑆3𝑉𝑥1𝑥2𝑟3Δ𝜑(𝐱)𝑑𝑉.(8.35) The volume integral here is weakly singular. Its regularization may be done using (8.31).

For the linear BE and piecewise constant approximation 𝜑(𝑥)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get 𝐽51,1=𝜕𝑆𝑛𝑥1𝑥2𝑟𝑛𝑟5𝑟3𝑟3𝑑𝑙=02𝜋𝑟3cos𝜑sin𝜑𝑟5𝑟𝑑𝜑02𝜋2𝑟cos𝜑sin𝜑3𝑟3𝑟𝑑𝜑=0.(8.36)

9. Conclusions

The method of regularization of the weakly singular, singular and hypersingular integrals, based on the second Green theorem, is developed here. These divergent integrals are presented when the boundary integral equation methods are used to solve problems in elastostatics. The equations that permit easy calculation of the weakly singular, singular, and hypersingular integrals for different shape functions are presented here. The divergent integrals over the circular domain have been calculated analytically. The main equations related to the boundary integral equation and boundary elements methods formulation in 2D and 3D elastostatics have been discussed in details.