# Divergent Integrals in Elastostatics: General Considerations

Authors: V. V. Zozulya
Journal: ISRN Applied Mathematics
 (2011)

---

## 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.


---

## Body

## 1. Introduction

A huge amount of publications is devoted to the boundary integral equation methods (BIM) and its application in science and engineering [1–5]. 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, 6–8]. 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 [9–19], 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 V with smooth boundary ∂V. The region V is an open bounded subset of the three-dimensional Euclidean space ℜ3 with a C0,1 Lipschitzian regular boundary ∂V. The boundary contains two parts ∂Vu and ∂Vp such that ∂Vu∩∂Vp=∅ and ∂Vu∪∂Vp=∂V. On the part ∂Vu are prescribed displacement ui(x)of the body points, and on the part ∂Vp are prescribed traction pi(x), respectively. The body may be affected by volume forces bi(x). 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 εij(x). The strain tensor and the displacement vector are connected by Cauchy relations(2.1)εij=12(∂iuj+∂jui),
where ∂/∂xi is a derivative with respect to space coordinates xi. The components of the strain tensor must also satisfy the Saint-Venant’s relations(2.2)∂kl2εij-∂il2εkj=∂kj2εil-∂ij2εkl.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(2.3)∂jσij+bi=0i,∀x∈V.
Here and throughout the paper the summation convention applies to repeated indices.The tensor of deformationεij(x) and stress σij(x) are related by Hook’s law(2.4)σij=cijklεij.
Here cijkl are elastic modules. In the case of homogeneous anisotropic medium they are symmetric  (2.5)cijkl=cjikl=cklij.
and satisfy condition of ellipticity(2.6)cijklεijεkl≥α1εijεij,∀εij,∀α1>0.In the case of homogeneous isotropic medium, the elastic modules have the form(2.7)cijkl=λδijδkl+μ(δikδjl+δilδjk),
where λ and μ are Lame constants, μ>0 and λ>-μ, δij 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(2.8)Aijuj+bi=0,∀x∈V.The differential operatorAij for homogeneous anisotropic medium has the form(2.9)Aij=cikjl∂k∂l,
and for homogeneous isotropic medium has the form(2.10)Aij=μδij∂k∂k+(λ+μ)∂i∂j.If the problem is defined in an infinite region, then solution of (2.8) must satisfy additional conditions at the infinity in the form (2.11)uj(x)=O(r-1),σij(x)=O(r-2)forr⟶∞,
where r=x12+x22+x32 is the distance in the three-dimensional Euclidian space.If the body occupied a finite regionVwith the boundary ∂V, it is necessary to establish boundary conditions. We consider the mixed boundary conditions in the form (2.12)ui(x)=φi(x),∀x∈∂Vu,pi(x)=σij(x)nj(x)=Pij[uj(x)]=ψi(x),∀x∈∂Vp.The differential operatorPij:uj→pi is called stress operator.  It transforms the displacements into the tractions. For homogeneous anisotropic and isotropic medium, they have the forms (2.13)Pij=cikjlnk∂l,Pij=λni∂k+μ(δij∂n+nk∂i),
respectively. Here ni are components of the outward normal vector, ∂n=ni∂i is a derivative in direction of the vector n(x) normal to the surface ∂Vp.


## 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 displacementsu and u*(3.1)a(u,u*)=cijklεij(u)εkl(u*).
Obviously that (3.2)a(u,u*)=a(u*,u),a(u,u)=σij(u)εij(u)≥α1εij(u).Integrating the equality (3.1) over the volume V and applying the Gauss-Ostrogradskii formula, we will obtain(3.3)∫Va(u,u*)dV=∫Vσij(u)εkl(u*)dV=∫∂Vσijnjui*dS-∫Vui*∂jσijdV.
Taking into account that Aijuj=∂jσij and pi=σijnj=Pij[uj], we will find the first Betti’s theorem in the form(3.4)∫Vui*AijujdV=∫Va(u,u*)dV-∫∂Vui*Pij[uj]dS.
We will replace ui and ui* 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)∫V(ui*Aijuj-uiAijuj*)dV=∫∂V(uiPij[uj*]-ui*Pij[uj])dS.
Taking into the account the definition of differential operator Aij given in (2.8) and the definition of differential operator Pij given in (2.12), we obtain relation (3.6)∫V(biui*-bi*ui)dV=∫∂V(pi*ui-piui*)dS,
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 bi*(x)=δijδ(x-y)(3.7)AijUkj(x-y)+δki(x-y)=0,∀x,y∈R3.
Now considering that (3.8)ui*(x)=Uij(x-y),pi*(x)=Pij[uj*(x)]=Wij(x,y),from (3.6)  we obtain the integral representation for the displacements vector(3.9)ui(y)=∫∂V(pi(x)Uji(x-y)-uj(x)Wji(x,y))dS+∫Vbi(x)Uji(x-y)dV,
which is called Somigliana’s identity. The kernels Uji(x-y) and Wji(x,y) are called fundamental solutions for elastostatics.Applying to (3.9) the differential operator Pij, we will find integral representation for the traction in the form(3.10)pi(y)=∫∂V(pi(x)Kji(x,y)-uj(x)Fji(x,y))dS+∫Vbi(x)Kji(x,y)dV.
The kernels Kji(x,y) and Fji(x,y) may be obtained applying the differential operator Pij  to the kernels Uji(x-y) and Wji(x,y), 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 solutionsUij(x-y) for the differential operator Aij, 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(4.1)Uij(x-y)=116πμ(1-υ)r((3-4υ)δij+∂ir∂jr).In 2D elastostatics they have the form(4.2)Uij(x-y)=18πμ(1-υ)((3-4υ)δijln⁡1r+∂ir∂jr).
Here r=(x1-y1)2+(x2-y2)2+(x3-y3)2 and r=(x1-y1)2+(x2-y2)2 for 3D and 2D case, respectively, ∂ir=∂r/∂xi=-∂r/∂yi=(xi-yi)/r.The kernelsWij(x,y) from (3.9) may be obtained by applying to Uij(x-y) differential operator (4.3)Pik[⋅,(x)]=λni(x)∂k[⋅]+μ⌊δiknj(x)∂j[⋅]+nk(x)∂i[⋅]⌋,
as it is shown here(4.4)Wij(x,y)=λni(x)∂kUkj(x,y)+μnk(x)⌊∂kUij(x,y)+∂iUkj(x,y)⌋.
Then after some transformations and simplifications, the expression for the kernels Wij(x,y) will have the following form(4.5)Wij(x,y)=-14π(1-υ)rα(nk(x)∂kr((1-2υ)δij+β∂ir∂jr)+(1-2υ)(ni(x)∂jr-nj(x)∂ir)).The kernelsKij(x,y) and Fij(x,y) from (3.10) may be obtained by applying differential operator
(4.6)Pik[⋅,(y)]=λni(y)∂k+μ⌊δiknj(y)∂j+nk(y)∂i⌋
to Uij(x-y) and Wij(x,y), respectively(4.7)Kij(x,y)=λni(y)∂kUjk(x,y)+μnk(y)⌊∂kUji(x,y)+∂iUjk(x,y)⌋,Fij(x,y)=λni(y)∂kWjk(x,y)+μnk(y)⌊∂kWji(x,y)+∂iWjk(x,y)⌋.
Then after some transformations and simplifications the expression for the kernels Kij(x,y) will have the form(4.8)Kij(x,y)=14π(1-υ)rα(nk(y)∂kr((1-2υ)δij+β∂ir∂jr)+(1-2υ)(ni(y)∂jr-nj(y)∂ir)),
and for the kernels, Fij(x,y) will have the form(4.9)Fij(x,y)=μ2απ(1-υ)rβ(βnk(x)∂kr((1-2υ)ni(y)∂jr+υ(δijnk(y)∂kr+nj(y)∂ir)-γnk(y)∂kr∂ir∂jr)+βν(ni(x)nk(y)∂kr∂jr+nk(x)nk(y)∂ir∂jr)+(1-2υ)(βnj(x)nk(y)∂kr∂ir+nk(x)nk(y)δij+ni(x)nj(y))-(1-4υ)nj(x)ni(y)).
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 r→0. More detailed analysis of (4.1)–(4.9) gives us the following results.In the 3D case withx→y,(5.1)Uij(x-y)⟶r-1,Wij(x,y)⟶r-2,Kij(x,y)⟶r-2,Fij(x,y)⟶r-3.In the 2D case withx→y,(5.2)Uij(x-y)⟶ln⁡(r-1),Wij(x,y)⟶r-1,Kij(x,y)⟶r-1,Fij(x,y)⟶r-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 coordinatedx,y∈ℜn (where n=3 or n=2) and region V with smooth boundary ∂V of the class C0,1. The boundary integrals of the types
(5.3)∫∂VG(x,y)rαφ(x)dS,α>0,
where G(x,y) is a finite function in ℜn×∂V and φ(x) is a finite function in ∂V, are  weakly singular for  α<n-1, singular for α=n-1, and strongly  singular  of hypersingularfor α>n-1.Definition 5.2.
Let we consider two points with coordinatedx,y∈ℜn (where n=3 or n=2) and region V with smooth boundary ∂V of the class C0,1. The boundary integrals of the types
(5.4)∫∂VG(x,y)ln⁡(r)φ(y)dy,
where G(x,y) is a finite function in ℜn×∂V and φ(x) is a finite function in ∂V,  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 Uij(x-y) are weakly singular and must be considered as improper
(5.5)W.S.∫∂Vpi(x)Uij(x-y)dS=lim⁡ε→0∫∂V∖∂Vεpi(x)Uij(x-y)dS.
Here ∂Vε is a part of the boundary, projection of which on tangential plane is contained in the circle Cε(x) of the radio ε  with center in the point x.Definition 5.4.
Integrals in (3.9) and (3.10) with kernels Wij(x,y) and Kij(x,y) are singular and must be considered in sense of the Cauchy principal values as
(5.6)P.V.∫∂Vui(x)Wij(x,y)dS=lim⁡ε→0∫∂V∖∂V(r<ε)ui(x)Wij(x,y)dS,P.V.∫∂Vpi(x)Kij(x,y)dS=lim⁡ε→0∫∂V∖∂V(r<ε)pi(x)Kij(x,y)dS.
Here ∂V(r<ε) is a part of the boundary, projection of which on tangential plane is the circle Cε(x) of the radio ε with center in the point x.Definition 5.5.
Integrals in (3.10) with kernels Fij(x,y) are hypersingular and must be considered in sense of the Hadamard finite part as
(5.7)F.P.∫∂Vui(x)Fji(x-y)dS=lim⁡ε→0(∫∂V∖∂V(r<ε)ui(x)Wji(x-y)dS+2uj(x)fj(x)∂V(r<ε)).
Here functions fj(x) 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 Uij(x-y), are weakly singular, and, therefore, they are continuous everywhere in the ℜn and, therefore, may be continuously extended on the boundary ∂V. The potentials with the kernels Wij(x,y) and Kij(x,y) contain singular kernels, and they jump when crossing the boundary ∂V. The potential with the kernels Fij(x,y) contain hypersingular kernels. They continuously cross the boundary ∂V.
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
(5.8)(∫∂Vpi(x)Uji(x-y)dS)±=(∫∂Vpi(x)Uji(x-y)dS)0,(∫∂Vui(x)Wji(x,y)dS)±=∓12uj(y)+(∫∂Vui(x)Wji(x,y)dS)0,(∫∂Vpi(x)Kji(x,y)dS)±=±12pj(y)+(∫∂Vpi(x)Kji(x,y)dS)0,(∫∂Vui(x)Fji(x-y)dS)±=(∫∂Vui(x)Fji(x-y)dS)0.
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  ∂V 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. Tendingy in (3.9) and (3.10) to the boundary ∂V and taking into consideration boundary properties of the potentials (5.8), we obtain representation of the displacements and traction vectors on the boundary surface ∂V. On the smooth parts of the boundary, they have the following form
(5.9)±12ui(y)=∫∂V(pj(x)Uij(x-y)uj(x)-Wij(x,y))dS+∫Vpj(x)Uij(x-y)dV,∓12pi(y)=∫∂V(pj(x)Kij(x,y)uj(x)-Fij(x,y))dS+∫Vpj(x)Kij(x,y)dV.
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 spacesX and Y and functional equation in those spaces(6.1)A⋅u=f,u∈D(A)⊂X,f∈R(A)⊂Y.
Here A:X→Y is the linear operator mapping from Banach space X in Banach space Y, D(A) is a domain, and R(A) is a range of the operator A.The equation (6.1) is named the exact equation, and its solution is the exact solution. We denote L(X,Y) Banach space of the linear operators mapping from X in Y.Let inX and Y act sequences of projection operators Ph and Ph′ such that (6.2)Ph2=Ph,PhX=Xh,Xh⊂X,(Ph′)2=Ph′,Ph′Y=Yh,Yh⊂Y,
where  Xh and Yh are finite dimension subspaces of the Banach spaces  X and Y; h∈R1 is a parameter of discretization.Now we consider operatorAh∈L(Xh,Yh) mapping in finite dimensional subspaces Xh and Yh and an approximate equation(6.3)Ah⋅uh=fh,Ah=Ph′⋅A⋅Ph,uh=Ph⋅u,fh=Ph′⋅f.
Solution uh 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.(6.4)Now let us consider operatorAh∈L(Xh,Yh) mapping in finite-dimensional subspaces   Xh and Yh 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 operatorsPh and Ph′ converge to identity operators in X and Y, respectively. It means that(6.5)lim⁡h→0‖Ph⋅u-u‖X=0∀u∈X,lim⁡h→0‖Ph′⋅f-f‖Y=0∀f∈Y.Definition 6.1.
Let conditions (6.5) be satisfied and stating from some h=h0>0 for any f∈Y, (6.3) has unique solution uh. In this case if
(6.6)lim⁡h→0‖Ah⋅uh-A⋅u‖Y=0,
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{Ah} mapping from Xh into Yh there is a constant γ>0 such that stating from some h=h0>0(6.7)‖Ah⋅uh-A⋅u‖Y≥γ‖uh‖X∀uh∈Xh
then for sequence of the operators {Ah}, 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 operatorsPh and Ph′ converge to identity operators in the Banach spaces X and Y,  respectively, as it stated in (6.5);(2)
the sequence of approximate operators{Ah} converges to A on each exact solution;(3)
the condition of stability (6.7) is satisfied for the sequence of operators {Ah}.
Then the following consequences take place:(1)
the exact solution exists and it is unique;(2)
for all enough small  h   exists a unique solutionuh∈Xh of the approximate equation (6.3);(3)
the sequence of the approximate equation{uh} converges to the exact one and takes place in the estimation(6.8)‖uh-Ph⋅u‖Xh≤γ-1‖Ph′A⋅u-AhPhu‖Yh.Thus, using a projective method instead of the exact solution of (6.1) in functional space X, we can find the solution of the approximate equation (6.3) in finite dimensional space Xh. The functional spaces {X,Xh} and {Y,Yh} are related by means of the projection operators Ph∈L(X,Xh) and Ph′∈L(Y,Yh), respectively. It is also important to construct inverse operator Ph-1∈L(Xh,X) which maps the finite dimensional space Xh into the initial functional space X. Such operator refers to as the operator of interpolation. Because of Xh⊂X, the interpolation operator is not unique; moreover, for any two functional spaces Xh and X, 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 X(∂V) and H-1/2(∂V) that are trace of displacements and traction on the boundary of the region in the following way:
(6.9)Uij⋅pj=∫∂V(Uij(x,y)pj(x))dS:Y(∂V)⟶X(∂V),Wij⋅uj=∫∂V(Wij(x,y)uj(x))dS:X(∂V)⟶Y(∂V),Kij⋅pj=∫∂V(Kij(x,y)pj(x))dS:Y(∂V)⟶X(∂V),Fij⋅uj=∫∂V(Fij(x,y)uj(x))dS:X(∂V)⟶Y(∂V).We have to construct finite dimensional functional spaces that correspond toX(∂V) and Y(∂V) and the corresponding projection operators. To construct finite dimensional functional spaces, we shall apply approximation by finite functions and splitting ∂V into finite elements
(6.10)∂V=⋃n=1N∂Vn,∂Vn∩∂Vk=∅,ifn≠k.Because∂V is the boundary of the region, these elements are called boundary elements. On each boundary element, we shall choose Q nodes of interpolation. Local projection operators act from functional X(∂Vn) and Y(∂Vn) to the finite directional ones  Xq(∂Vn) and Yq∂Vn)(6.11)Pqu:X(∂Vn)⟶Xq(∂Vn)∀x∈∂Vn,Pqp:Y(∂Vn)⟶Yq(∂Vn)∀x∈∂Vn.
Global projection operators are defined as the sum of the local projection operators
(6.12)Pnqu=∑n=1NPqu,Pnqp=∑n=1NPqp.
They map X(∂V) and Y(∂V) to finite dimensional interpolations spaces
(6.13)Pnqu:X(∂V)⟶Xq(⋃n=1N∂Vn)∀x∈∂V,Pnqp:Y(∂V)⟶Yq(⋃n=1N∂Vn)∀x∈∂V.
The local projection operators Pnu and Pnp establish correspondence between vectors of displacements and traction and their value on the nodes of interpolation of the boundary elements ∂Vn(6.14)Pqu⋅ui(x)={uin(xq),q=1,…,Q}∀x∈∂Vn,Pqp⋅pi(x)={pin(xq),q=1,…,Q}∀x∈∂Vn.
Similarly for the global operators, we have
(6.15)Pnqu⋅ui(x)={uin(xq),q=1,…,Q;n=1,…,N}∀x∈∂V,Pnqp⋅pi(x)={pin(xq),q=1,…,Q;n=1,…,N}∀x∈∂V.
Let us construct local interpolation operators (Pqu)-1 and (Pqp)-1. For this purpose, we will introduce systems of shape functions φnq(x) and ψnq(x) in the finite dimensional functional spaces Xq(∂Vn) and Yq(∂Vn). Then the vectors of displacements and traction on the boundary element ∂Vn will be represented approximately in the form
(6.16)ui(x)≈∑q=1Quin(xq)φnq(x),x∈∂Vn,pi(x)≈∑q=1Qpin(xq)φnq(x),x∈∂Vn,
and on the whole crack surface ∂V in the form
(6.17)ui(x)≈∑n=1N∑q=1Quin(xq)φnq(x),x∈⋃n=1N∂Vn,pi(x)≈∑n=1N∑q=1Qpin(xq)φnq(x),x∈⋃n=1N∂Vn.Finite-dimensional analogies for the integral operators (6.9) are operators which map the finite-dimensional functional spaces Xq(⋃n=1N∂Vn) and Yq(⋃n=1N∂Vn), from one to another
(6.18)Uijnq=Pnqp⋅Uij⋅Pnqu:Yq(∂Vn)⟶Xq(∂Vn),Wijnq=Pnqp⋅Wij⋅Pnqu:Xq(∂Vn)⟶Yq(∂Vn),Kijnq=Pnqp⋅Kij⋅Pnqu:Yq(∂Vn)⟶Xq(∂Vn),Fijnq=Pnqp⋅Fij⋅Pnqu:Xq(∂Vn)⟶Yq(∂Vn).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
(6.19)12uim(yr)=∑n=1N∑q=1QUjin(yr,xq)pjn(xq)-Wjin(xr,xq)ujn(xq)+Ui(f,y,Vn),12pim(yr)=∑n=1N∑q=1QKjin(yr,xq)pjn(yq)-Fjin(yr,xq)ujn(xq)+Ki(f,y,Vn),
where
(6.20)Ujin(yr,xq)=∫∂VnUji(yr,x)ψnq(x)dS,Wjin(yr,xq)=∫∂VnWji(yr,x)φnq(x)dS,Kjin(yr,xq)=∫∂VnKji(yr,x)ψnq(x)dS,Fjin(yr,xq)=∫∂VnFji(yr,x)φnq(x)dS.
The volume potentials Ui(f,y,Vn) and Ki(f,y,Vn) depend on discretization of the V 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 areaV⊂Rn and a BE model of its boundary ∂V⊂Rn-1 (n=2,3). We fix in the area V finite number of points xq(g=1,…,Q); these points refer to as global node points V(q)={xq∈V:g=1,…,Q}. We shall divide the area V into finite number of subareas Vn(n=1,…,N) which are FEs. They have to satisfy the following conditions:  (7.1)Vn∪Vm=∅,m≠n,m,n=1,2,…,N,V=⋃n=1NVn.
On each FE we introduce a local coordinate system ξ. The nodal points xq∈Vn in the local system of coordinates we designate by ξq. They are coordinates of nodal points in the local coordinate system. Local and global coordinate are related in the following way:(7.2)xq=∑n=1NΛnξnq.
Functions Λn 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 V.Having constructed FE model of the areaV, we shall consider approximation of the function f(x) that belong to some functional space. The FE model of the area V is the domain of function which should be approximated. We denote function f(x) on the FE Vn by fn(x). Then(7.3)f(x)=∑n=1Nfn(x).
On each FE, the local functions fn(x) may be represented in the form (7.4)fn(x)≈∑q=1Qfn(xq)φnq(ξ),
where φnq(ξ) are interpolation polynomials or nodal functions of the FE with number n. In the nodal point with coordinates xq, 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 f(x) looks like(7.5)f(x)≈∑n=1N∑q=1Qfn(xq)φnq(ξ).
If the nodal point q 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ξi, in which the element has a simple appearance, and global xi where the FE represents more complex figure. Local coordinates ξi should be functions of global (ξi(x1,x2,x3)) ones, and on the contrary, global coordinates should be functions of (xi(ξ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(7.6)J=det⁡|∂xi∂ξj|≠0,J-1=det⁡|∂ξi∂xj|≠0.
The differential elements along coordinate axes are related by (7.7)dxi=(∂xi∂ξj)dξj,dx=J(ξ)dξ,dξi=(∂ξi∂xj)dxj,dξ=J-1(x)dx.
The volume element in the ℜ3 is transformed under the formula(7.8)dV=dx1dx2dx3=J(ξ)dξ1dξ2dξ3,
and the area element in the ℜ2  is transformed under the formula(7.9)dA=dx1dx2=det⁡|∂xα∂ξβ|dξ1dξ2,α,β=1,2.
The differential of the surface located in the ℜ3 is defined by the expression(7.10)dS=(n12+n22+n32)1/2dξ1dξ2,
where(7.11)n1=∂x1∂ξ1∂x3∂ξ2-∂x2∂ξ2∂x3∂ξ1,n2=∂x3∂ξ1∂x1∂ξ2-∂x1∂ξ1∂x3∂ξ2,n3=∂x1∂ξ1∂x2∂ξ2-∂x2∂ξ1∂x1∂ξ2.
The element of length of a contour in the ℜ2 is defined by the expression(7.12)dl=[(dx1dξ1)2+(dx2dξ1)2]1/2dξ1.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(7.13)φnq(x)={1∀x∈Vn,0∀x∉Vn.
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(7.14)φnq=12(1+ξ1qξ1),ξiq=±1.In 2D case for the rectangular FE, shape functions are(7.15)φnq=14(1+ξ1qξ1)(1+ξ2qξ2),ξiq=±1.For the triangular FE, it is convenient to introduce system of coordinates which connected with the Cartesian by relation(7.16)xi=ξqxiq,∑ξq=1,i=1,2,q=1,2,3.
For the triangular FE, the interpolation functions are (7.17)φnq=ξq.
The local coordinates ξq 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(7.18)φnq=12ξ1qξ1(1+ξ1qξ1)(q=1,2,3)
at the end nodes and (7.19)φnq=1-ξ12
at the midpoint node.In 2D case for the rectangular FE, shape functions are(7.20)φnq=14(1+ξ1qξ1)(1+ξ2qξ2)(ξ1qξ1-ξ2qξ2-1)(q=1,…,8)
at the angular nodes and (7.21)φnq=12(1-ξ12)(1+ξ2qξ2)ξ1q=0,φnq=12(1-ξ22)(1+ξ1qξ1)ξ2q=0
at the side nodes.For the triangular FE, the interpolation functions are(7.22)φnq=(2ξq-1)ξq(q=1,…,6)
at the angular nodes and (7.23)φnq=4ξ1ξ2⋯.
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 φnq defined in (7.13)–(7.23). In this case, the global coordinates have the form(7.24)xi=xiqφnq(ξ).
Here xiq is the global coordinate of nodal points for the FE number n.


## 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(7.13)φnq(x)={1∀x∈Vn,0∀x∉Vn.
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(7.14)φnq=12(1+ξ1qξ1),ξiq=±1.In 2D case for the rectangular FE, shape functions are(7.15)φnq=14(1+ξ1qξ1)(1+ξ2qξ2),ξiq=±1.For the triangular FE, it is convenient to introduce system of coordinates which connected with the Cartesian by relation(7.16)xi=ξqxiq,∑ξq=1,i=1,2,q=1,2,3.
For the triangular FE, the interpolation functions are (7.17)φnq=ξq.
The local coordinates ξq 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(7.18)φnq=12ξ1qξ1(1+ξ1qξ1)(q=1,2,3)
at the end nodes and (7.19)φnq=1-ξ12
at the midpoint node.In 2D case for the rectangular FE, shape functions are(7.20)φnq=14(1+ξ1qξ1)(1+ξ2qξ2)(ξ1qξ1-ξ2qξ2-1)(q=1,…,8)
at the angular nodes and (7.21)φnq=12(1-ξ12)(1+ξ2qξ2)ξ1q=0,φnq=12(1-ξ22)(1+ξ1qξ1)ξ2q=0
at the side nodes.For the triangular FE, the interpolation functions are(7.22)φnq=(2ξq-1)ξq(q=1,…,6)
at the angular nodes and (7.23)φnq=4ξ1ξ2⋯.
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 φnq defined in (7.13)–(7.23). In this case, the global coordinates have the form(7.24)xi=xiqφnq(ξ).
Here xiq is the global coordinate of nodal points for the FE number n.


## 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 Uij(x-y) are weakly singular (WS). They have to be considered as improper integrals. The integrals with the kernels Wij(x,y) and Kij(x,y) are singular. They have to be considered in the sense of Cauchy as principal values (PVs). The integrals with the kernels Fij(x,y) 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 [11–16], 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 r-m in the form(8.1)F.P.∫-aaφ(x)rmdx=∑i=0k-1(-1)i+1didxiPirm-kdk-1-iφ(x)dxk-1-i|x=-ax=a+(-1)k∫-aaPkrm-kdkφ(x)dxk,F.P.∫Vφ(x)rmdV=∑i=0k-1(-1)i+1∫∂V[Δk-i-1φ(x)∂nPirm-2i-Pirm-2i∂nΔk-i-1φ(x)]dS+(-1)k∫V1rm-2kΔk+1φ(x)dV,
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(8.2)J0=∫abφ(x)ln⁡1xdx,Jk=∫abφ(x)xkdx,k=1,2
Here φ(x) is the smooth function that depend on the shape of the BE and the interpolation polynomials.We consider first the weakly singular integralJ0. 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)∫abdg(x)dxφ(x)dx=φ(x)g(x)|ab-∫abdφ(x)dxg(x)dx.
Let in this formula be g(x)=xln1/x, dg(x)/dx=ln1/x, then we obtain (8.4)J0=W.S.∫abφ(x)ln⁡1xdx=φ(x)xln⁡1x|ab-∫abdφ(x)dxxln⁡1xdx.
Obviously the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation  φ(x)=1, and we get (8.5)J0=W.S.∫abln⁡|x-y|dx=(b-y)ln⁡|b-y|-(a-y)ln⁡|a-y|(a<y<b).For the singular integralJ1 regularization, we will also use the formula for integration by parts (8.3). Let in this formula g(x)=-ln1/x, (dg(x))/dx=1/x, then we obtain (8.6)J1=P.V.∫abφ(x)xdx=(dφ(x)dxxln⁡1x-φ(x)ln⁡1x)|ab-∫abd2φ(x)dx2xln⁡1xdx.
Here also the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation  φ(x)=1, and we get (8.7)J1=P.V.∫abdxx-y=ln⁡|b-ya-y|,(a<y<b).Finally for the hypersingular integralJ2 regularization, we will use formula (8.1) and the above-obtained result for singular J1 regularization. Finally we get(8.8)J2=F.P.∫abφ(x)x2dx=(d2φ(x)dx2xln⁡1x-dφ(x)dxln⁡1x-φ(x)x)|ab-∫abd3φ(x)dx3xln⁡1xdx.
Here also the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation  φ(x)=1, and we get (8.9)J2=F.P.∫abdx(x-y)2=-1(b-y)+1(a-y)(a<y<b).
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 (8.10)F.P.∫-aadyy2=-2a,a>0.


### 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(8.11)Jkl,m=∫Snx1lx2mrkφ(x)dS,l,m=0,1,2,k=3,4,5.
Here φ(x) is the smooth function that depends on shape of the BE and interpolation polynomials.


#### 8.2.1. Integrals with Kernels of the Typer-k,k=1,2,3

From (8.1) for k=1, we get regularization for weakly singular integral (8.12)J10,0=W.S.∫Vφ(x)rdV=∫∂V[φ(x)rn2r-r∂nφ(x)]dS+∫VrΔφ(x)dV.
Here rn=(xα-yα)nα, 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φ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar, coordinates we will get (8.13)J10,0=12∫∂Snrnrdl=∫02πrrrdφ=πr.
In order to regularize singular integral, we will use, relation 1/r2=(1/2)Δ(lnr)2. Then in the same way we get(8.14)∫Vφ(x)r2dV=12∫∂V(φ(x)2rnln⁡rr2-(ln⁡r)2∂nφ(x))dS+12∫V(ln⁡r)2Δφ(x)dV.
The volume integral on the right is weakly singular. Taking into account the relation (lnr)2=(r2/6)Δ(lnr)4, we obtain regularization for this weakly singular integral(8.15)∫V(ln⁡r)2Δφ(x)dV=16∫∂V(2Δφ(x)rn(ln⁡r)3-r2(ln⁡r)4∂nΔφ(x))dS+16∫Vr2(ln⁡r)4Δ2φ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate singular integral in (8.14) analytically. Introducing polar coordinates, we will get (8.16)J20,0=∫∂Snrnln⁡rr2dl=∫02πrln⁡rr2rdφ=2πln⁡r.
Finally from (8.1) for k=3, we get regularization for hupersingular integral (8.17)∫Vφ(x)r3dV=∫∂V[Δφ(x)rn2r-φ(x)rnr3-1r∂nφ(x)-r∂nΔφ(x)]dS+∫VrΔ2φ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.18)J30,0=-∫∂Snrnr3dl=-∫02πrr3rdφ=-2πr.Now using (8.1), any divergent integral with the kernels of the type 1/rk for any positive integer k can be calculated.


#### 8.2.2. Integrals with Kernels of the Typexα2/rk,k=3,4,5

The weakly singular integral with kernelxα2/r3 is calculated taking into account the equation xα2/r3=(1/3)((2/r)-Δ(xα2/r)). It is easy to show that(8.19)J32,0=W.S.∫Vφ(x)x12r3dV=23W.S.∫Vφ(x)rdV-13W.S.∫Vφ(x)Δx12rdV.
The first integral here is already calculated in (8.12). The second one may be presented in the form(8.20)∫Vφ(x)Δx12rdV=∫∂V[φ(x)(2n1x1r-x12rnr3)-x12r∂nφ(x)]dS+∫Vx12rΔφ(x)dV.
Combining the last two equations, finally we will get (8.21)J32,0=13∫∂V[φ(x)(2n1x1r-x12rnr3+rn2r)-(x12r+r)∂nφ(x)]dS+13∫V(x12r+r)Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area we can calculate this integral analytically. Introducing polar coordinates, we will get (8.22)J32,0=13∫∂Sn(x12rnr3-2x1n1r+2rnr)dl=13(∫02π(rcos⁡φ)2rr3rdφ-∫02π2r(cos⁡φ)2rrdφ+2∫02πrrrdφ)=πr.The singular integral with kernelxα2/r4 is calculated taking into account that Δxα2/r4=(1/4)((2/r2)-Δxα2/r2). In this case,(8.23)J42,0=P.V.∫Vφ(x)x12r2dV=14∫Vφ(x)(2r2-Δx12r2)dV=12∫Vφ(x)r2dV-14∫Vφ(x)Δx12r2dV.
The first integral here is already calculated in (8.14). The second one may be presented in the form(8.24)∫Vφ(x)Δx12r2dV=∫∂V[φ(x)(2n1x1r2-2x12rnr4)-x12r2∂nφ(x)]dS+∫Vx12r2Δφ(x)dV.
Combining the last two equations, finally we will get (8.25)J42,0=12∫∂V[φ(x)(x12rnr4-n1x1r2+rnln⁡rr2)-((ln⁡r)22-x122r2)∂nφ(x)]dS+14∫V((ln⁡r)2-x12r2)Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.26)J42,0=12∫∂Sn[x12rnr4-x1n1r2+rnln⁡rr2]dl=∫02π(rcos⁡φ)2rr4rdφ-∫02πr(cos⁡φ)2r2rdφ+∫02πrln⁡rr2rdφ=2πln⁡r.The hypersingular integral with kernelxα2/r5 is calculated taking into account that xα2/r5=(1/3)((2/r3)-Δ(xα2/r3)). In this case, it is easy to show that(8.27)J52,0=F.P.∫Vφ(x)x12r5dV=13∫Vφ(x)(2r3-Δx12r3)dV=23∫Vφ(x)r3dV-13∫Vφ(x)Δx12r3dV.
The first integral here is already calculated in (8.17). The second one may be presented in the form(8.28)∫Vφ(x)Δx12r3dV=∫∂V[φ(x)(2n1x1r3-3x12rnr5)-x12r3∂nφ(x)]dS+∫Vx12r3Δφ(x)dV.
Combining the last two equations, finally we will get (8.29)J52,0=23∫∂V(φ(x)(n1x1r3-rnr3-3x12rn2r5)-(1r+x12r3)∂nφ(x))dS+23W.S.∫V(1r-x12r3)Δφ(x)dV.
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  φ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.30)J52,0=∫∂Sn(2rn3r3+2x12rn3r5-x1n1r3)dl=23∫02πrr3rdφ+23∫02π(rcos⁡φ)2rr5rdφ-∫02πr(cos⁡φ)2r3rdφ=-πr.


#### 8.2.3. Integrals with Kernels of the Type(x1x2)/rk,k=3,4,5

The weakly singular integral with kernel(x1x2)/r3 is calculated using (8.1). Taking into account that (x1x2)/r3=-1/3Δ(x1x2)/r, it is easy to show that(8.31)J31,1=W.S.∫Vφ(x)x1x2r3dV=13∫∂V[φ(x)(x1x2rnr3-r*r)+x1x2r∂nφ(x)]dS-13∫Vx1x2rΔφ(x)dV.
Here r*=x1n2+x2n1.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.32)J31,1=13∫∂Sn[x1x2rnr3-r*r]dl=∫02πr3cos⁡⁡φsin⁡φr3rdφ-∫02π2rcos⁡⁡φsin⁡φrrdφ=0.The singular integral with kernelx1x2/r4 is calculated using equation (8.1). Taking into account that x1x2/r4=-(1/4)Δx1x2/r2, it is easy to show that(8.33)J41,1=P.V.∫Vφ(x)x1x2r4dV=14∫∂V[φ(x)(2x1x2rnr4-r*r2)-x1x2r2∂nφ(x)]dS-14∫Vx1x2r2Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.34)J41,1=14∫∂Sn[2x1x2rnr4-r*r2]dl=14∫02π2r3cos⁡⁡φsin⁡φr4rdφ-14∫02π2rcos⁡⁡φsin⁡φr2rdφ=0.The hypersingular integral with kernelx1x2/r5 is calculated using (8.1). Taking into account that x1x2/r5=-(1/3)Δx1x2/r3, it is easy to show that(8.35)J51,1=F.P.∫Vφ(x)x1x2r5dV=∫∂V[φ(x)(x1x2rnr5-r*3r3)-x1x2r3∂nφ(x)]dS-13∫Vx1x2r3Δφ(x)dV.
The volume integral here is weakly singular. Its regularization may be done using (8.31).For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.36)J51,1=∫∂Sn[x1x2rnr5-r*3r3]dl=∫02πr3cos⁡⁡φsin⁡φr5rdφ-∫02π2rcos⁡⁡φsin⁡φ3r3rdφ=0.


## 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(8.2)J0=∫abφ(x)ln⁡1xdx,Jk=∫abφ(x)xkdx,k=1,2
Here φ(x) is the smooth function that depend on the shape of the BE and the interpolation polynomials.We consider first the weakly singular integralJ0. 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)∫abdg(x)dxφ(x)dx=φ(x)g(x)|ab-∫abdφ(x)dxg(x)dx.
Let in this formula be g(x)=xln1/x, dg(x)/dx=ln1/x, then we obtain (8.4)J0=W.S.∫abφ(x)ln⁡1xdx=φ(x)xln⁡1x|ab-∫abdφ(x)dxxln⁡1xdx.
Obviously the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation  φ(x)=1, and we get (8.5)J0=W.S.∫abln⁡|x-y|dx=(b-y)ln⁡|b-y|-(a-y)ln⁡|a-y|(a<y<b).For the singular integralJ1 regularization, we will also use the formula for integration by parts (8.3). Let in this formula g(x)=-ln1/x, (dg(x))/dx=1/x, then we obtain (8.6)J1=P.V.∫abφ(x)xdx=(dφ(x)dxxln⁡1x-φ(x)ln⁡1x)|ab-∫abd2φ(x)dx2xln⁡1xdx.
Here also the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation  φ(x)=1, and we get (8.7)J1=P.V.∫abdxx-y=ln⁡|b-ya-y|,(a<y<b).Finally for the hypersingular integralJ2 regularization, we will use formula (8.1) and the above-obtained result for singular J1 regularization. Finally we get(8.8)J2=F.P.∫abφ(x)x2dx=(d2φ(x)dx2xln⁡1x-dφ(x)dxln⁡1x-φ(x)x)|ab-∫abd3φ(x)dx3xln⁡1xdx.
Here also the integral on the left is divergent and on the right is regular one. For the linear BE and piecewise constant approximation  φ(x)=1, and we get (8.9)J2=F.P.∫abdx(x-y)2=-1(b-y)+1(a-y)(a<y<b).
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 (8.10)F.P.∫-aadyy2=-2a,a>0.


## 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(8.11)Jkl,m=∫Snx1lx2mrkφ(x)dS,l,m=0,1,2,k=3,4,5.
Here φ(x) is the smooth function that depends on shape of the BE and interpolation polynomials.


### 8.2.1. Integrals with Kernels of the Typer-k,k=1,2,3

From (8.1) for k=1, we get regularization for weakly singular integral (8.12)J10,0=W.S.∫Vφ(x)rdV=∫∂V[φ(x)rn2r-r∂nφ(x)]dS+∫VrΔφ(x)dV.
Here rn=(xα-yα)nα, 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φ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar, coordinates we will get (8.13)J10,0=12∫∂Snrnrdl=∫02πrrrdφ=πr.
In order to regularize singular integral, we will use, relation 1/r2=(1/2)Δ(lnr)2. Then in the same way we get(8.14)∫Vφ(x)r2dV=12∫∂V(φ(x)2rnln⁡rr2-(ln⁡r)2∂nφ(x))dS+12∫V(ln⁡r)2Δφ(x)dV.
The volume integral on the right is weakly singular. Taking into account the relation (lnr)2=(r2/6)Δ(lnr)4, we obtain regularization for this weakly singular integral(8.15)∫V(ln⁡r)2Δφ(x)dV=16∫∂V(2Δφ(x)rn(ln⁡r)3-r2(ln⁡r)4∂nΔφ(x))dS+16∫Vr2(ln⁡r)4Δ2φ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate singular integral in (8.14) analytically. Introducing polar coordinates, we will get (8.16)J20,0=∫∂Snrnln⁡rr2dl=∫02πrln⁡rr2rdφ=2πln⁡r.
Finally from (8.1) for k=3, we get regularization for hupersingular integral (8.17)∫Vφ(x)r3dV=∫∂V[Δφ(x)rn2r-φ(x)rnr3-1r∂nφ(x)-r∂nΔφ(x)]dS+∫VrΔ2φ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.18)J30,0=-∫∂Snrnr3dl=-∫02πrr3rdφ=-2πr.Now using (8.1), any divergent integral with the kernels of the type 1/rk for any positive integer k can be calculated.


### 8.2.2. Integrals with Kernels of the Typexα2/rk,k=3,4,5

The weakly singular integral with kernelxα2/r3 is calculated taking into account the equation xα2/r3=(1/3)((2/r)-Δ(xα2/r)). It is easy to show that(8.19)J32,0=W.S.∫Vφ(x)x12r3dV=23W.S.∫Vφ(x)rdV-13W.S.∫Vφ(x)Δx12rdV.
The first integral here is already calculated in (8.12). The second one may be presented in the form(8.20)∫Vφ(x)Δx12rdV=∫∂V[φ(x)(2n1x1r-x12rnr3)-x12r∂nφ(x)]dS+∫Vx12rΔφ(x)dV.
Combining the last two equations, finally we will get (8.21)J32,0=13∫∂V[φ(x)(2n1x1r-x12rnr3+rn2r)-(x12r+r)∂nφ(x)]dS+13∫V(x12r+r)Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area we can calculate this integral analytically. Introducing polar coordinates, we will get (8.22)J32,0=13∫∂Sn(x12rnr3-2x1n1r+2rnr)dl=13(∫02π(rcos⁡φ)2rr3rdφ-∫02π2r(cos⁡φ)2rrdφ+2∫02πrrrdφ)=πr.The singular integral with kernelxα2/r4 is calculated taking into account that Δxα2/r4=(1/4)((2/r2)-Δxα2/r2). In this case,(8.23)J42,0=P.V.∫Vφ(x)x12r2dV=14∫Vφ(x)(2r2-Δx12r2)dV=12∫Vφ(x)r2dV-14∫Vφ(x)Δx12r2dV.
The first integral here is already calculated in (8.14). The second one may be presented in the form(8.24)∫Vφ(x)Δx12r2dV=∫∂V[φ(x)(2n1x1r2-2x12rnr4)-x12r2∂nφ(x)]dS+∫Vx12r2Δφ(x)dV.
Combining the last two equations, finally we will get (8.25)J42,0=12∫∂V[φ(x)(x12rnr4-n1x1r2+rnln⁡rr2)-((ln⁡r)22-x122r2)∂nφ(x)]dS+14∫V((ln⁡r)2-x12r2)Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.26)J42,0=12∫∂Sn[x12rnr4-x1n1r2+rnln⁡rr2]dl=∫02π(rcos⁡φ)2rr4rdφ-∫02πr(cos⁡φ)2r2rdφ+∫02πrln⁡rr2rdφ=2πln⁡r.The hypersingular integral with kernelxα2/r5 is calculated taking into account that xα2/r5=(1/3)((2/r3)-Δ(xα2/r3)). In this case, it is easy to show that(8.27)J52,0=F.P.∫Vφ(x)x12r5dV=13∫Vφ(x)(2r3-Δx12r3)dV=23∫Vφ(x)r3dV-13∫Vφ(x)Δx12r3dV.
The first integral here is already calculated in (8.17). The second one may be presented in the form(8.28)∫Vφ(x)Δx12r3dV=∫∂V[φ(x)(2n1x1r3-3x12rnr5)-x12r3∂nφ(x)]dS+∫Vx12r3Δφ(x)dV.
Combining the last two equations, finally we will get (8.29)J52,0=23∫∂V(φ(x)(n1x1r3-rnr3-3x12rn2r5)-(1r+x12r3)∂nφ(x))dS+23W.S.∫V(1r-x12r3)Δφ(x)dV.
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  φ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.30)J52,0=∫∂Sn(2rn3r3+2x12rn3r5-x1n1r3)dl=23∫02πrr3rdφ+23∫02π(rcos⁡φ)2rr5rdφ-∫02πr(cos⁡φ)2r3rdφ=-πr.


### 8.2.3. Integrals with Kernels of the Type(x1x2)/rk,k=3,4,5

The weakly singular integral with kernel(x1x2)/r3 is calculated using (8.1). Taking into account that (x1x2)/r3=-1/3Δ(x1x2)/r, it is easy to show that(8.31)J31,1=W.S.∫Vφ(x)x1x2r3dV=13∫∂V[φ(x)(x1x2rnr3-r*r)+x1x2r∂nφ(x)]dS-13∫Vx1x2rΔφ(x)dV.
Here r*=x1n2+x2n1.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.32)J31,1=13∫∂Sn[x1x2rnr3-r*r]dl=∫02πr3cos⁡⁡φsin⁡φr3rdφ-∫02π2rcos⁡⁡φsin⁡φrrdφ=0.The singular integral with kernelx1x2/r4 is calculated using equation (8.1). Taking into account that x1x2/r4=-(1/4)Δx1x2/r2, it is easy to show that(8.33)J41,1=P.V.∫Vφ(x)x1x2r4dV=14∫∂V[φ(x)(2x1x2rnr4-r*r2)-x1x2r2∂nφ(x)]dS-14∫Vx1x2r2Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.34)J41,1=14∫∂Sn[2x1x2rnr4-r*r2]dl=14∫02π2r3cos⁡⁡φsin⁡φr4rdφ-14∫02π2rcos⁡⁡φsin⁡φr2rdφ=0.The hypersingular integral with kernelx1x2/r5 is calculated using (8.1). Taking into account that x1x2/r5=-(1/3)Δx1x2/r3, it is easy to show that(8.35)J51,1=F.P.∫Vφ(x)x1x2r5dV=∫∂V[φ(x)(x1x2rnr5-r*3r3)-x1x2r3∂nφ(x)]dS-13∫Vx1x2r3Δφ(x)dV.
The volume integral here is weakly singular. Its regularization may be done using (8.31).For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.36)J51,1=∫∂Sn[x1x2rnr5-r*3r3]dl=∫02πr3cos⁡⁡φsin⁡φr5rdφ-∫02π2rcos⁡⁡φsin⁡φ3r3rdφ=0.


## 8.2.1. Integrals with Kernels of the Typer-k,k=1,2,3

From (8.1) for k=1, we get regularization for weakly singular integral (8.12)J10,0=W.S.∫Vφ(x)rdV=∫∂V[φ(x)rn2r-r∂nφ(x)]dS+∫VrΔφ(x)dV.
Here rn=(xα-yα)nα, 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φ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar, coordinates we will get (8.13)J10,0=12∫∂Snrnrdl=∫02πrrrdφ=πr.
In order to regularize singular integral, we will use, relation 1/r2=(1/2)Δ(lnr)2. Then in the same way we get(8.14)∫Vφ(x)r2dV=12∫∂V(φ(x)2rnln⁡rr2-(ln⁡r)2∂nφ(x))dS+12∫V(ln⁡r)2Δφ(x)dV.
The volume integral on the right is weakly singular. Taking into account the relation (lnr)2=(r2/6)Δ(lnr)4, we obtain regularization for this weakly singular integral(8.15)∫V(ln⁡r)2Δφ(x)dV=16∫∂V(2Δφ(x)rn(ln⁡r)3-r2(ln⁡r)4∂nΔφ(x))dS+16∫Vr2(ln⁡r)4Δ2φ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate singular integral in (8.14) analytically. Introducing polar coordinates, we will get (8.16)J20,0=∫∂Snrnln⁡rr2dl=∫02πrln⁡rr2rdφ=2πln⁡r.
Finally from (8.1) for k=3, we get regularization for hupersingular integral (8.17)∫Vφ(x)r3dV=∫∂V[Δφ(x)rn2r-φ(x)rnr3-1r∂nφ(x)-r∂nΔφ(x)]dS+∫VrΔ2φ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.18)J30,0=-∫∂Snrnr3dl=-∫02πrr3rdφ=-2πr.Now using (8.1), any divergent integral with the kernels of the type 1/rk for any positive integer k can be calculated.


## 8.2.2. Integrals with Kernels of the Typexα2/rk,k=3,4,5

The weakly singular integral with kernelxα2/r3 is calculated taking into account the equation xα2/r3=(1/3)((2/r)-Δ(xα2/r)). It is easy to show that(8.19)J32,0=W.S.∫Vφ(x)x12r3dV=23W.S.∫Vφ(x)rdV-13W.S.∫Vφ(x)Δx12rdV.
The first integral here is already calculated in (8.12). The second one may be presented in the form(8.20)∫Vφ(x)Δx12rdV=∫∂V[φ(x)(2n1x1r-x12rnr3)-x12r∂nφ(x)]dS+∫Vx12rΔφ(x)dV.
Combining the last two equations, finally we will get (8.21)J32,0=13∫∂V[φ(x)(2n1x1r-x12rnr3+rn2r)-(x12r+r)∂nφ(x)]dS+13∫V(x12r+r)Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area we can calculate this integral analytically. Introducing polar coordinates, we will get (8.22)J32,0=13∫∂Sn(x12rnr3-2x1n1r+2rnr)dl=13(∫02π(rcos⁡φ)2rr3rdφ-∫02π2r(cos⁡φ)2rrdφ+2∫02πrrrdφ)=πr.The singular integral with kernelxα2/r4 is calculated taking into account that Δxα2/r4=(1/4)((2/r2)-Δxα2/r2). In this case,(8.23)J42,0=P.V.∫Vφ(x)x12r2dV=14∫Vφ(x)(2r2-Δx12r2)dV=12∫Vφ(x)r2dV-14∫Vφ(x)Δx12r2dV.
The first integral here is already calculated in (8.14). The second one may be presented in the form(8.24)∫Vφ(x)Δx12r2dV=∫∂V[φ(x)(2n1x1r2-2x12rnr4)-x12r2∂nφ(x)]dS+∫Vx12r2Δφ(x)dV.
Combining the last two equations, finally we will get (8.25)J42,0=12∫∂V[φ(x)(x12rnr4-n1x1r2+rnln⁡rr2)-((ln⁡r)22-x122r2)∂nφ(x)]dS+14∫V((ln⁡r)2-x12r2)Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.26)J42,0=12∫∂Sn[x12rnr4-x1n1r2+rnln⁡rr2]dl=∫02π(rcos⁡φ)2rr4rdφ-∫02πr(cos⁡φ)2r2rdφ+∫02πrln⁡rr2rdφ=2πln⁡r.The hypersingular integral with kernelxα2/r5 is calculated taking into account that xα2/r5=(1/3)((2/r3)-Δ(xα2/r3)). In this case, it is easy to show that(8.27)J52,0=F.P.∫Vφ(x)x12r5dV=13∫Vφ(x)(2r3-Δx12r3)dV=23∫Vφ(x)r3dV-13∫Vφ(x)Δx12r3dV.
The first integral here is already calculated in (8.17). The second one may be presented in the form(8.28)∫Vφ(x)Δx12r3dV=∫∂V[φ(x)(2n1x1r3-3x12rnr5)-x12r3∂nφ(x)]dS+∫Vx12r3Δφ(x)dV.
Combining the last two equations, finally we will get (8.29)J52,0=23∫∂V(φ(x)(n1x1r3-rnr3-3x12rn2r5)-(1r+x12r3)∂nφ(x))dS+23W.S.∫V(1r-x12r3)Δφ(x)dV.
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  φ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.30)J52,0=∫∂Sn(2rn3r3+2x12rn3r5-x1n1r3)dl=23∫02πrr3rdφ+23∫02π(rcos⁡φ)2rr5rdφ-∫02πr(cos⁡φ)2r3rdφ=-πr.


## 8.2.3. Integrals with Kernels of the Type(x1x2)/rk,k=3,4,5

The weakly singular integral with kernel(x1x2)/r3 is calculated using (8.1). Taking into account that (x1x2)/r3=-1/3Δ(x1x2)/r, it is easy to show that(8.31)J31,1=W.S.∫Vφ(x)x1x2r3dV=13∫∂V[φ(x)(x1x2rnr3-r*r)+x1x2r∂nφ(x)]dS-13∫Vx1x2rΔφ(x)dV.
Here r*=x1n2+x2n1.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.32)J31,1=13∫∂Sn[x1x2rnr3-r*r]dl=∫02πr3cos⁡⁡φsin⁡φr3rdφ-∫02π2rcos⁡⁡φsin⁡φrrdφ=0.The singular integral with kernelx1x2/r4 is calculated using equation (8.1). Taking into account that x1x2/r4=-(1/4)Δx1x2/r2, it is easy to show that(8.33)J41,1=P.V.∫Vφ(x)x1x2r4dV=14∫∂V[φ(x)(2x1x2rnr4-r*r2)-x1x2r2∂nφ(x)]dS-14∫Vx1x2r2Δφ(x)dV.For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.34)J41,1=14∫∂Sn[2x1x2rnr4-r*r2]dl=14∫02π2r3cos⁡⁡φsin⁡φr4rdφ-14∫02π2rcos⁡⁡φsin⁡φr2rdφ=0.The hypersingular integral with kernelx1x2/r5 is calculated using (8.1). Taking into account that x1x2/r5=-(1/3)Δx1x2/r3, it is easy to show that(8.35)J51,1=F.P.∫Vφ(x)x1x2r5dV=∫∂V[φ(x)(x1x2rnr5-r*3r3)-x1x2r3∂nφ(x)]dS-13∫Vx1x2r3Δφ(x)dV.
The volume integral here is weakly singular. Its regularization may be done using (8.31).For the linear BE and piecewise constant approximationφ(x)=1 and circular area, we can calculate this integral analytically. Introducing polar coordinates, we will get (8.36)J51,1=∫∂Sn[x1x2rnr5-r*3r3]dl=∫02πr3cos⁡⁡φsin⁡φr5rdφ-∫02π2rcos⁡⁡φsin⁡φ3r3rdφ=0.


## 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.



---

*Source: 726402-2011-08-02.xml*