CHAPTER 3

INTRODUCTION TO THE BOUNDARY ELEMENT METHOD

3.1. INTRODUCTION

The boundary element method was developed at the University of Southampton by

combining the methodology of the finite element method with the boundary integral
method. The first international conference devoted to the boundary element method took
place in 1978 at Southampton [7]. Since that time, many books have been published ([8],
[9], [10], ...) and the numerous contributions to the annual conferences like BEM and
BE1ECH show the rapid development of the new method for all the engineering fields.
In this chapter, the weighed residual method will be used to develop the boundary element
method for the case of anisotropic Laplace problems. The weighed residual method is the
most general technique, because it can also be applied to develop the finite difference
method and the finite element method for instance.
The application of the boundary element method to Laplace problems is presented in some
well known references like the four books mentioned above. Because of this, only some
special points will be described in this chapter, especially those relating to the anisotropic
characteristics of the medium.
Finally, the application of the boundary element method to Laplace problems leads to
special numerical integration problems and problems connected with the subdivision of
the domain into sub-regions. These points will be investigated in chapters 4 and 6.

3.2. GOVERNING EQUATION

Let Q be a domain in which the following anisotropic Laplace equation applies, ie.

(3.1)

The values of the components of the tensor K depend on the characteristics of the domain
Kx and Ky. They are the principal values associated to the corresponding principal axis x
and y. In order to determine the value of the potential u on the domain Q, one has to
impose Dirichlet (along r 1) or Neumann (along r 2) conditions on the boundaries of the
domain:

E. K. Bruch, The Boundary Element Method for Groundwater Flow

© Springer-Verlag Berlin, Heidelberg 1991
 18

u = u along f, and q = q along f2 (3.2)

The bars on the terms in (3.2) indicate that those are the imposed values of the potential
(u) or the normal velocity (q) along the boundary.
y

x
Figure 3.1

In case of a Cauchy boundary condition ('YU + 'l'q + 't =0), one of the two unknowns will
be expressed in terms of the other and afterwards introduced in the same way as a
Dirichlet or a Neumann boundary condition into the linear system.

3.3. TIffi BASIC EQUATION OF TIffi BOUNDARY ELEMENT METHOD

3.3.1. The weighed residual method

When one applies a numerical method to solve the Laplace equation (3.1), the numerical
result will generally give an error EI' This error must be as small as possible in order to
obtain accurate results. Thus, the approximation of the potential field u in the Laplace
equation gives the following result:

K (~u) , :;to
= E (3.3)

In order to obtain an accurate solution, the error EI must be as small as possible. EI is also
called the residual and in order to minimize it, one can use a weighting function FI to
weight the residual on the domain. Thus, the weighed residual method applied to relation
(3.3) leads to :

(3.4)

The accuracy of the solution depends on the choice of the weighting function Fl'
The application of a numerical method to solve the Laplace equation does not allow to
satisfy it rigorously on the domain Q. Consequently, the solution will be only approached
 19

along the boundary where the boundary conditions will produce errors Ez and ~ :

(3.5)

Thus, one proceeds with the boundary conditions by using two weighting functions F2
and F3, which are generally different from Fl. So, one obtains the following
relationships for the two parts of the boundary r.

(3.6.a)

f ',F,dr, f
r2
=
r2
(3.6.b.)

As each of these three weighting residual relations (3.4) and (3.6) must be equal to zero,
one can group them together in a single relation, ie.

(3.7)

As each term of (3.7) must be separately equal to zero, one can choose their signs to
allow some simplifications.
Before proceeding with this developments, one needs to make a choice concerning the
weighting functions F l , F2 and F3. In order to obtain the boundary element method, one
can define,

(3.8)

The special choice concerning F2 will be justified showing that it produces a useful
relationship which leads to the cancellation of certain terms. One can also consider the
[Link] two terms of the right part of (3.7) to justify the special choice for F2.
After the introduction of (3.8) into (3.7), one obtains the following weighed residual
relation:

(3.9)
 20
-
dU dU
with: q = dii and q = dii

In order to continue the transformation of (3.9), one now requires the Green's theorem to
transfer the Laplace operator from the potential u to the weighting function w.

3.3.2. The Green's theorem

The Green's theorem is given by the following relationship:

I !l
W K (Au) dQ = I !l
uK (Aw) dQ + I~ I
r
w dr -
r
(3.10)

In the case of an anisotropic domain, the tensor K is given by equation (2.6) and it is
possible to show that the operator d/dn is as follows,

(3.11)

The coefficients nx and ny are the components of the unit outward normal vector to the
boundary.
In order to transfer the Laplace operator from the unknown potential u to the weighting
function w, one can use the Green's theorem (3.10) in equation (3.9).

3.3.3. The basic relation of the boundary element method

After the introduction of (3.10) into (3.9), one obtains after several steps the following
equation:

I!l
uK(~w)dQ=
I I
r
udildI"-
dw
r
au
Wdiidr (3.12)

with : u = u on r and u = u on r
2 1
(3.13)
q = q on r and q = q on r
1 2

The first term of (3.12) contains an integral in the domain Q, whereas its second term
contains two integrals on the boundary r of the domain. A judicious choice of the
weighting function w will allow to eliminate the integral on the domain Q.
 21

3.4. 1HE WEIGHTING FUNCTION

3.4.1. Choice of the weiihtini function

In order to eliminate the remaining integral in the domain in (3.12), one can choose as
weighting function w the fundamental solution of the Laplace equation with Sj, the Dirac
function or distribution. Thus, the weighting function and its derivate along the outward
normal vector can be written as,

(3.14)

~
2
-1 x y
w =21tjKx"Ky In Kx +"Ky (3.15)

(3.16)

Taking into account the properties of the Dirac distribution, the introduction of (3.14) into
(3.12) gives:

(3.17)

The two remaining integrals of (3.17) apply only on the boundary of the domain. This
relationship is the starting point to develop the boundary element method.

3.4.2. Localization of the collocation point

The relation (3.14) is the continuity equation corresponding to a unit source located at the
point i. This point is called the collocation point and its position can not be chosen
arbitrarily. Indeed, the weighting function (3.15) and its derivative (3.16) are infinite at
the collocation point i. Considering this property, one has three possibilities regarding the
localization of the collocation point:
- firstly, one can locate the collocation point i on the boundary r. In this case, the
integrands of (3.17) are infinite at i and one has to be careful when carrying out the
numerical integrations.
- secondly, one can locate the collocation point i outside the domain Q. In this caSe, the
integrands of (3.17) do not tend towards infinity. It is however difficult to define the
optimal position for i. The choice of the position of i will influence the accuracy of the
results [63], [60].
- finally, the collocation point can be located inside the domain, in which case the
integrals here have finite values. In this case, the potential at the collocation point
introduces an additional unknown into the system.
 22
In the present case the fIrst approach (which is the classical one) will be used. One will
see later on that the detailed study of the numerical integrations will allow to obtain very
accurate results. One also has to analyse the behaviour of the integrals of (3.17) in the
vicinity of the collocation point i. This point will be investigated in the following section.

3.5. ANALYSIS OF THE INTEGRALS

3.5.1. Introduction

In this section, the behaviour of the integrals of (3.17), in the vicinity of the collocation
point i will be studied. Indeed, x and y are equal to zero at i and consequently, the
weighting function w (3.15) and its derivative (3.16) are infInite there. Consequently,
one needs to verify if the two integrals of (3.17) are integrable and if the singularities at
this point could have any special contribution to equation (3.17).

3.5.2. The fIrst integral

Taking into account the relation (3.15) for the weighting function w, the second integral
in (3.17) is given by:

11 =
-1 I doau (3.18)
2njKxKy n
r
Y

Yj

Yi-r---------r.-r~~-----------------x~

Figure 3.2

The integrand of (3.18) tends towards infInite at the collocation point i, but as In r is
integrable on each interval [0 IL], with IL different from infinite, II is integrable. Thus,
the numerical integration of (3.18) do not pose problems, except in the vicinity of i where
one has to use a numerical integration appropriate to the behaviour of the integrand there.
In the case of this fIrst integral (3.18), the singularity does not give a specifIc
contribution. To show this, one can consider the collocation point as isolated from the
rest of the boundary with the help of a circular notch r e the radius e of which will tend
towards zero (fIgure 3.2).
 23
Thus, II gives the following relationship:

(3.19)

The first integral of (3.19) tends towards 11. To analyse the second one, one can change
from cartesian coordinates to polar ones. Thus, one obtains the following relationships:

x =£cosO y = £ sinO (3.20)

The introduction of (3.20) into the second integral of (3.19) gives:

lim (3.21)
£~O

This last relation shows that the singularity at the collocation point does not result in any
specific contribution to the integral 11.

3.5.3. The second integral

Taking into account the relation (3.16) for the derivative of the weighting function, the
second integral is given by :

(3.22)

The analysis of the integral (3.22) shows that its integrand is infinite at the collocation
point i, but does not tend towards infinite at this point. Consequently, the integrand is
continuous at the collocation point i. It is possible to show that the value of the derivative
of the weighting function at the collocation point is given by :

an x any
aw ny ar - nx ar
-1
lim (3.23)
x~O.y~O
do

Due to the special behaviour of the integrand of 12 in the vicinity of the collocation point i,
the singularity will give a non zero contribution to the integral. In order to determine this
contribution in the anisotropic case [13], [19], one can use the same technique as for the
first integral (see figure 3.2). In the present case, 12 gives:
 24

(3.24)

The first integral of (3.24) tends towards 12. If one introduces now (3.24) into the
relationship (3.17), one obtains the relation (3.25) because u tends towards Uj when E
goes to zero.

J J
au aw
CjU j = w- dr- u do dr (3.25)
an
r r

with : Cj =1+
lim
E-tO
I £
aw
- dr
an

FOIIDula (3.26) shows that the value of Cj depends on the following (see figure 3.2)
(3.26)

- the angle (8b - 8a );

- the characteristics Kx and Ky of the medium;
- the position of the angle a in the principal reference system (x, y) associated with the
characteristics of the medium.
The integration of (3.26) gives:

1-~[ g 8}f
cj =
2n
Arctg (
J[K:t
R; e
a
(3.27)

The primitive of the integral of the equation (3.26) given by (3.27) presents a
discontinuity at 8 =± n/2 which one needs to take into account when computing (3.27).
In the case of isotropic characteristics (K = Kx = Ky), (3.27) gives:

Cj = 1 - ro I 2n (3.28)

where ro = 8b - 8a. Notice that in the case of a smooth boundary (ro = n), equation (3.28)
gives'cj = 0,5.
As have been seen, the integrand of (3.22) does not tend towards infinite at the
collocation point i, but is continuous. Nevertheless, in order to integrate numerically
(3.22) in the vicinity of i, special integration methods will be required.
 25

3.6. DISCRETIZATION OF THE PROBLEM

To discretize in equation (3.25) the potential u and its nonnal derivative q on the boundary
r, one can use the classical isoparametric shape functions [8], [9], [10] for two, three and
four node elements. These elements correspond respectively to linear, quadratic and
cubic shape functions. The use of the shape functions entails a change of variables.
Thus, each element is transfonned from the general reference system (X,Y) to its own
non-dimensional reference system. This is done by the relation (3.29) with j varying
from 1 to Nn(e), the number of nodes of the element being e. As in the classical
isoparametric finite element, the different shape functions Nj are multiplied by the
corresponding node values 2j of the coordinates or the unknowns.

(3.29)

The figure 3.3 shows the change of variables in the case of a three node element

y
(2) (3)

(l)~ -1'
(l) (3)
1
.. ~

Figure 3.3

To discretize the whole boundary, Ne elements like those shown by the figure 3.3 will be
used. Thus, one obtains :
Ne
r = L re (3.30)
e=!

The introduction of the discretization of the boundary and the unknowns into the
boundary integral equation (3.25) gives the following discretized boundary element
equation:

Ne (Nn(e) e) Ne (Nn(e) e )
L L q.G..
J IJ
- L L u·J R ..IJ (3.31)
e=! j=! e=! j=!

W1·th : H IJ..e = i
re
dW
"'C"[Link]
an J e (3.32), and G~.D = ~(
re
w N. ctr
J e
(3.33)

Equation (3.31) contains one unknown at each node of the discretization. For this
reason, one locates successively the collocation point i at each node of the discretization
and obtains a linear system representing the boundary value problem.
 26

3.7. SPECIAL PROBLEMS

3.7.1. Determination of the coefficients oeij and Heij

In the discretized boundary element equation (3.31), the coefficients oeij and Hei} given
by the relations (3.32) and (3.33), need to be determined. These numencal integrations
will be considered in the next chapter. Indeed, their integration requires special
re
techniques when the collocation point i is located on the element under consideration.

3.7.2. Determination of the ci coefficient

As shown by the relation (3.26), the determination of the ci coefficient is more difficult
when the domain has anisotropic characteristics. In order to get round this problem, one
can use equation (3.31) to determine the ci coefficient. As shown by the relations (3.26),
(3.32) and (3.33), the ci coefficient, oeij and Heij depends only on the geometry of the
domain and its characteristics and not on the boundary conditions. Consequently, if one
considers that the potential u has a constant value on the whole domain, equation (3.31)
gives:
Ne (Nn(e) )
c·1 = - ~ ~ H~.
1J (3.34)
e=l j=l

3.7.3. The division of the domain into sub-regions

The method so far described in this chapter requires that the coefficients Kx and Ky are
constant on the whole domain. But in many practical cases, the domain needs to be
divided into several parts, each of them having different characteristics. In these cases,
one has to divide the domain into several sub-regions having different characteristics.
The sub-regions will be separated by interfaces along which one has to impose the
continuity of the potential u and the normal velocity q. The intersection points of several
interfaces or of an interface and a boundary lead to special numerical problems. The
division of the domain into sub-regions and these special problems will be considered in
chapter 6.

3.7.4. The junction point of two boundaries

At the junction point of two boundaries r A and r B (figure 3.4), there is uncertainty
regarding the continuity of the potential u and its derivative along the two outward normal
vectors corresponding to the two boundaries at this point.
In order to solve this problem, different methods have been suggested :
- to place a small curvature at the corner at the point P to obtain there one single
normal vector [39];
- the use of two physically separated nodes in the vicinity of P or the use of
discontinuous elements [64], [63];
- the use of an additional equation to describe the potential field in the vicinity of P
[2], [43], [55];
- the double node technique [70].
The double node technique consists in the localization of two nodes at the point P,
belonging respectively to the two boundaries r A and r B. This method avoids the
 27

ambiguity between the two boundary conditions at this point and is the one used here.
However, this technique can lead to problems at the junction point of several interfaces.
These problems will be analysed in chapter 6 which studies the division of the domain
into sub-regions .

...n p

Figure 3.4 Figure 3.5

3.7.5. Detennination of the unknowns at interior points

Mter the solution of the boundary value problem, the potential and the nonnal velocity at
every boundary point are known and equation (3.17) allows to determine the potential at
any interior point. If one locates the collocation point i at an interior point, the only
remaining unknown in (3.17) is the potential at this point. Finally, in order to obtain the
gradient at this point, one can differentiate the relation (3.17) [8].

3.8. CONCLUSIONS

In this chapter, the boundary element method for the case of a Laplace problem on an
anisotropic domain has been presented.
The chapter has shown that the numerical integrations poses some problems due to the
fact that the weighting function wand its nonnal derivative are infmite at the collocation
point i. The study of this problem is the subject of the next chapter.