Boundary Element and Finite Element Methods

MTH-ME82 Boundary Element and Finite Element Methods
Course details
20 Lectures
3 hour exam. Choose 3 out of 4 questions.
Recommended Text
Pozrikidis, C. (1992) Boundary Integral & Singularity Methods for Linearised Viscous Flow,
Cambridge.
There is one copy of this book in the library.
Problem sheet
One exercise sheet contains assorted problems covering the whole course. One question
from the exercise sheet will appear verbatim on the summer exam paper.
Contents
1. Background material.
2. Boundary integral methods for potential flow: Laplaces equation, Poissons equation.
3. Boundary integral methods for Stokes flow
4. The Finite Element Method
MTH-ME82
Boundary Element and Finite Element Methods
1. Background
We will start by covering some background material, some of which you may already
be familiar with, to prepare for the main course material.
1.1 Index notation
In standard index notation, we represent the components of a vector using a subscript.
So
ui
for i = 1, . . . , 3 represents the three-dimensional vector u = (u1 , u2 , u3 ).
Definition: The Kronecker delta, ij , is defined by

ij =
0 if i 6= j
1 if i = j
So 12 = 0 and 22 = 1 and so on. If we write this as a matrix (assuming i = 1, 2 and

j = 1, 2) we have

1 0
,
0 1
we recognise it as the 2 2 identity matrix.
Definition: The alternating tensor, ijk , is defined as
ijk
0 if any of i, j, k are equal

=
1 if i, j, k are in cyclic order
1
otherwise.
So, for example,

123 = 231 = 1,
132 = 1,
133 = 0.
The alternating tensor provides a convenient way of expressing a vector (cross) product in
index notation.
Exercise: Show that ijk aj bk is equal to the ith component of a b by writing out the
components.
Einsteins summation convention: According to this convention we sum over an index whenever we see it repeated (the range of the sum will be obvious from the context). For example,
for a calculation in three dimensions, i = 1, 2, 3, and
aii
means
3
X
i=1
aii .
MTH-ME82
For a calculation in two dimensions, i = 1, 2, and

bii
means
2
X
bii .
i=1
For example, ii = 2 in two dimensions, and ii = 3 in three dimensions.

Index shifting: The Kronecker delta is useful for shifting the index of a vector. For example,
we can write
ui = ij uj .
The summation convention is being used here to sum over the repeated index j on the right
hand side. Note that the only non-zero contribution to the sum comes when i = j, and
thus we recover the term on the left hand side.
1.2 The delta function
The delta function is an example of a generalised function or distribution. Informally,
we express the delta function as

0 if x 6= 0
(x) =
if x = 0
We will not need the formal theory behind the delta function here1 . However, we will make
use of the following property.
Sifting property: When a delta function multiplies another function inside an integral, the
result is the value of that function at the point where the argument of the delta function
vanishes. So,
Z
(x)f (x) dx = f (0).
The delta function sifts out all values of f except that at x = 0. Also,
Z
(x a)f (x) dx = f (a).
The delta function sifts out all values of f except that at x = a.

1.3 The principal value of an integral
Consider the following divergent (or improper) integral,
Z
1
1
1
dx.
x
If you are interested, a good book is Applied Functional Analysis by D. Griffel.
MTH-ME82
There is evidently a problem with the integrand at x = 0. Proceeding carefully, we re-write

the integral as
Z 1
Z 1
h Z 1
i
1
1
dx = lim
dx +
dx .
0
1 x
x
1 x
Integrating, we have
lim
0
i
ln || ln | 1| + ln |1| ln || = 0.
However, we could equally have re-written the integral as

Z 1
Z 1
h Z a 1
i
1
1
dx = lim
dx +
dx ,
0
x
1 x
1
x
for any a > 0. Integrating, we find
h
i
lim ln |a| ln | 1| + ln |1| ln || = ln a,
0
and we can choose a to get any result we want!

Definition: The principal value of the integral
Z 1
1
dx
x
1
is defined to be
PV
1
h
1
dx = lim
0
x
1
dx +
x
Z

i
1
dx = 0.
x
The principal value of other improper integrals is defined in a similar way.
1.4 Applications of the Boundary Integral Method

Potential Flow
A potential flow is one which is both inviscid and irrotational. Consider a steady (timeindependent) flow with velocity field u. If the flow is inviscid, then u satisfies the Euler
equations
(u )u = p,
where p is the fluid pressure. The flow is irrotational if
u = 0.
In this case we can define a scalar potential, , such that
u = .
Conservation of mass requires that
u = 0.
4
MTH-ME82
Substituting u = , we find
2 = 0.
(1.1)
So a potential flow satisfies Laplaces equation.

For a potential flow, the boundary condition at a solid wall is
u n = 0,
where n is the normal to the wall. This is called the normal flow condition and it demands
that fluid not pass through the solid wall. Since u = , the condition may be written
n = 0
or
= 0,
n
(1.2)
on the solid wall.

A typical potential flow problem might require us to solve (1.1) subject to (1.2). If the
geometry of the flow domain is simple, like a circle say, then this may be straightforward.
Example: Solve the equation
2 = 0
in the circle 0 r 1 with = 1 on r = 1. In polar coordinates (assuming no
dependence), Laplaces equation becomes

1 d
d
r
= 0.
r dr
dr
Integrating, we find
= A log r + B.
We set A = 0 to remove the singularity at r = 0, and require B = 1 to satisfy the boundary
condition on the circle. So the solution is = 1.
Exercise: Think about how you would solve the same problem as above in a domain with
a more complex geometry. How would you satisfy the boundary condition = 1 on a
boundary with a complicated shape?
Electrostatics
A conducting metal sphere of given radius is charged to a potential = 1 and held above
flat ground at z = 0 where the potential is zero, = 0. The electric field outside the sphere
satisfies Laplaces equation,
2 = 0.
Compute the electric potential outside the sphere. This problem can be solved using the
Boundary Integral Method. Moreover, it can still be solved if the sphere is deformed to
some other shape.
MTH-ME82
Stokes flow
Stokes flow occurs at small values of the Reynolds number (we will discuss this more closely
later in the course). The Boundary Integral Method has been applied to a host of different
problems in Stokes flow. Examples include:
Oscillations of a gas bubble
Bursting of a bubble near to a wall or a free surface
Flow of a viscous fluid film over a shaped surface
Deformation of red blood cells passing through capillaries
Break-up of a viscous liquid jet
2. The boundary integral method (BIM) for potential flow

Over the last 30 or so years, a very popular method for solving Stokes flows in intricate
geometries has developed and matured. This nifty method calculates a Stokes flow solely by
reference to what happens at the flow boundaries. As we will see, we can write the velocity
at any point in the flow field purely in terms of the flow quantites on the boundaries. The
method is called the boundary integral method.
The historical development of the mathematical machinery sitting behind the boundary
integral method is nicely reviewed in the paper Cheng & Cheng (2005) Heritage and early
history of the boundary element method, Eng. Analysis Bound. Elem., 29, 268-302.
The main advantage of the boundary integral method for solving certain classes of PDEs is
as follows: Since reference is only made to boundary values, the dimension of the problem
is effectively reduced by one. So a three dimensional problem effectively becomes a two
dimensional problem. This brings about a tremendous saving in computing time.
Its is an extremely powerful method for this reason: It can cope with any geometry at
all. Since many practical applications involve very complex geometries, this constitutes a
major advantage over other methods such as finite differences which are very clumsy if the
geometry is not simple.
The boundary integral method can be applied to both potential flows and Stokes flows.
Since the implementation of the former is slightly easier, it is these which we shall discuss
first.
2.1 The BIM for potential flow in two dimensions
Aim: To solve Laplaces equation inside a domain D, with boundary C whose unit normal
n is defined to point inside D.
6
MTH-ME82
We take as our starting point Greens second identity,

2 2 = ( ).
(2.1)
We choose = G, where G satisfies the equation

2 G + (x x0 ) = 0,
(2.2)
where the second term is the delta function.

Definition: We call x0 the pole or the singular point or just the singularity.
nb: G is called a fundamental solution of Laplaces equation or Greens function2 (we will
discuss these in more detail later).
Let us define
r = |x x0 |,
so that r is the distance from x0 to x. When r 6= 0 (2.2) becomes
2 G = 0,
(2.3)
or (assuming G only depends on r)

1 d
r dr

dG
r
= 0.
dr
We notice that
G = log r,
(2.4)
for constant , is a solution. For reasons to be revealed, we choose = 1/2.

Definition: The fundamental solution to Laplaces equation,
G=
1
log r,
2
(2.5)
is called the free-space Greens function.

In the next step, we integrate (2.1) over the domain D. In fact, to avoid the problem with
the delta function at x = x0 , we delete from D a small disk, D , of radius centred at x0 .
Then we define D0 = D D and integrate (2.1) over D0 :
ZZ
ZZ
2
2
(G G) dS =
(G G) dS = 0,
D0
D0
since both 2 G = 0 and 2 = 0 within D0 . Using the divergence theorem we obtain

Z
(G G) n dl = 0,
(2.6)
C+C
2
Strictly speaking one should only refer to this as a Greens function when it satisfies a boundary value
problem, that is G adopts prescribed values on a given boundary
MTH-ME82
where C is the boundary to D and C is the boundary to D , and l measures arc length
along either C or C .
Consider now the integral around C ,
Z
I =
(G G) n dl
C
In the limit 0, the radius of the disk shrinks to zero. From (2.4)
G = log r.
Therefore
I =
C
log r
r
r

dl.
Since r = on C , letting 0, we find

Z
Z
(x0 )
(x0 )
I = (log ) (x0 )
dl
dl = 2( log ) (x0 )2
= 2(x0 ).
r

r

C
C
Rearranging (2.6) we have
Z
(G2 2 G) dS = 2(x0 ),
and so
(x0 ) =
1
2
(2 G G2 ) dS.
(2.7)
For tidiness, we choose = 1/2. Substituting for G we have

Z

Z
1
log r n dl
n log r dl .
(x0 ) =
2
C
C
(2.8)
This is called Greens third identity. In fact equation (2.7) applies for a general Greens
function G which satisfies (2.2) and any chosen boundary conditions.
Definition: The boundary integral equation (BIE) for potential flow is given by
Z
(x0 ) =
n G G n dl,
(2.9)
C
where the unit normal to the boundary C, n, points inside the solution domain D.
Note : Equation (2.8) is an example of an integral equation. Although appears as the
subject on the left hand side, the function and its derivatives also appear on the right
hand side inside the integrals.
Other Greens functions
8
MTH-ME82
It is worth emphasizing at this stage that a Greens function is simply a device to enable us
to solve problems. In the present scenario, we will use a Greens function to help us solve
Laplaces equation
2 = 0
subject to some boundary condition, for example
=0
on
where represents a boundary. The solution which we end up with does not depend on the
Greens function used to obtain it. In this sense, we can think of a Greens function simply
as a tool which is used to create an end product; but the end product itself is not changed
by the particular choice of tool used to put it together.
Above we derived the free space Greens function,
G=
1
log r,
2
where r = |x x0 | and x0 is the location of the singularity. In some cases, it is expedient

to select a Greens function which vanishes (or its derivative vanishes) on the boundary. A
simple example is a Greens function for a singularity x0 = (x0 , y0 ) (with y0 > 0) placed
above a plane wall at y = 0. The Greens function is required to vanish at the wall, that is
G = 0 when y = 0.
We can construct this Greens function by placing an image singularity of opposite
sign3 below the wall at y = y0 . We write down the free-space Greens function for each
singularity and add them together, thus
G=
1
1
log r +
log R,
2
2
where
r=
So
p
(x x0 )2 + (y y0 )2 ,
1
G = log
4
R=
p
(x x0 )2 + (y + y0 )2 .
(x x0 )2 + (y + y0 )2
(x x0 )2 + (y y0 )2
Now, setting y = 0 we obtain

1
G = log
4
(x x0 )2 + y02
(x x0 )2 + y02

= 0,
as required.
Alternatively, we may instead wish to construct a Greens function whose y-derivative
vanishes at the wall. After come experimenting, we notice that
G=
3
1
1
log r
log R,
2
2
By this we mean that we take G = (1/2) log r rather than G = (1/2) log r
MTH-ME82
fits the bill. To check, we differentiate with respect to y to obtain

Gy =
Now,
ry =
1 ry
1 Ry
.
2 r
2 R
y y0
,
r
So,
Gy =
Setting y = 0 we find
Gy =
ry =
y + y0
.
R
1 y y0
1 y + y0
.
2
2 r
2 R2
1 y0
1 y0
= 0,
2
2 r
2 R2
as required.
Periodic Greens functions
Sometimes it is useful to work with a periodic Greens function which satisfies a periodic
relation, for example
G(x, y; x0 , y0 ) = G(x + L, y; x0 , y0 ),
for some L > 0. One example is a periodic array of singularities located at the positions
(x0 nL, y0 )
for n = 0, 1, 2, . . .. After some working we find
n
o
1
G = log 2{cosh(k[y y0 ]) cos(k[x x0 ])} ,
4
where k = 2/L.
For a periodic array of singularities located above a wall at y = 0, with G required to vanish
at the wall, we have
n
o 1
n
o
1
G = log 2{cosh(k[yy0 ])cos(k[xx0 ])} + log 2{cosh(k[y+y0 ])cos(k[xx0 ])} .
4
4
Exercise: Check that the periodic Greens function above has the property G = 0 when
y = 0. How would you modify the formula so that Gy = 0 on y = 0?
The Boundary Integral Method
Our goal is to solve (2.9) for . To make this more explicit, we write out (2.9) again showing
explicitly all the dependencies:
Z
(x0 ) =
n(x) G(x, x0 ) G(x, x0 ) n(x) (x) dl(x).
C
Note that
10
MTH-ME82
x0 is the singular point in the interior of the domain D.

x is the so-called field point.
Integration is done with respect to arclength around C. Note that the local increment
of arc length, dl depends on x.
We notice that on the right hand side, the unknown function appears inside integrals
which are defined over the boundary C. So the right hand side is only concerned with
boundary values of . This suggests the following way to proceed:
(i) Take the point x0 to lie on the boundary C. Then the integral equation involves only
values of on the boundary.
(ii) Solve the integral equation for the boundary values of .
(iii) Once the boundary values of are known we can compute the right hand side of (2.8)
for any x0 inside D. Hence we can find everywhere in the domain D.
Points (i-iii) form the basic steps of the boundary integral method.
Step (i) requires us to take the point x0 to lie on C. This requires some careful thought.
Lets work with a general Greens function, G, and write out the boundary integral equation
(2.9)
Z
Z
(x) n G(x, x0 ) dl(x)
G(x, x0 ) n (x) dl(x).
(2.10)
(x0 ) =
C
nb: Weve written out explicitly the arguments on each function in (2.10) to keep clear the
functional dependencies.
On the right hand side of (2.10), there are two integrals,
Z
Z
G(x, x0 ) n (x) dl(x)
and
(x) n G(x, x0 ) dl(x).
C
We call the first integral the single layer potential (SLP).

We call the second integral the double layer potential (DLP).
The limit x0 C
It can be shown that the SLP is continuous as the point x0 approaches and then crosses C.
In contrast, as x0 crosses C, the double-layer potential undergoes a jump discontinuity. As
x0 approaches C from inside the domain D, we find (see the Problem Sheet)
Z
lim {
x0 C
PV
(x) [n G(x, x0 )] dl(x)} =

C
11
1
(x)[n G(x, x0 )] dl(x) + (x0 ), (2.11)
2
MTH-ME82
and when x0 approaches C from outside D, we find

Z PV
Z
1
(x)[n G(x, x0 )] dl(x) (x0 ). (2.12)
lim { (x) [n G(x, x0 )] dl(x)} =
x0 C
2
C
C
Note that on the right hand side of these equations x0 is located on C. (For multiplyconnected domains, inside is taken to mean the area into which the normal vector points.)
The superscript P V means that we are taking the principal value of the integral.
Example
Suppose that the domain D occupies the strip 1 x 1, y > 0, and label as L that part
of C lying on y = 0. Consider what happens to the DLP as x0 travels down the y axis
through L.
We have that x0 = (0, y0 ) as we are interested in what happens as y0 passes down
through zero from above. The unit normal on L pointing into D is the unit vector in the
positive y direction. So we have
G
n G =
.
y
Taking the free space Greens function we have
G
1
1 (y y0 )
=
(log r) =
.
y
2 y
2
r2
since r2 = x2 + (y y0 )2 and
(y y0 )
r
=
.
y
r2
So the DLP over L is
Z
Z
D=
(x) n G dx =
L
(x)
G
1
dx =

y y=0
2
(x) y0
dx.
x2 + y02
A difficulty will arise in the integral at x = 0 when y0 = 0. We rewrite as

Z 1
Z 1
(x) (0)
1
y0
1
D=
y0
dx +
(0)
2
2 dx,
2
2
2
2
x + y0
1
1 x + y0
Consider the first integral. For small x, (x) = (0) + x0 (0) + . Hence near to x = 0,
h (x) (0) i
h
x i
1
0
lim
lim
(0)
= 0 (0) .
2
2
2
2
+
+
x
x
+
y
x
+
y
y0 0
y0 0
0
0
Therefore in the limit y0 0+ the first integral is improper and behaves like the integral
Z 1
1
dx.
x
1
To make sense of it in this limit, we take the principal value of the original integral. Therefore
we write, for y0 0+ ,
Z 1
Z 1
1
(x) (0)
1
y0
D=
y0 P V
dx +
(0) P V
2
2 dx,
2
2
2
2
x + y0
1
1 x + y0
12
MTH-ME82
For the second integral,

Z
y0
dx = 2 tan1
2
x + y02
1
y0

.
We can see that the value of this integral is discontinuous at y0 = 0 since

Z 1
Z 1
y0
y0
dx = ,
lim
lim
2
2 dx = .
2
2
+
y0 0
y0 0
1 x + y0
1 x + y0
When y0 is on L (so that y0 = 0) the integrand vanishes so that
Z 1
y0
PV
2 dx = 0.
2
1 x + y0
Putting everything together, then, we have that in the limit y0 0+ ,
Z 1
Z 1
1
(x) (0)
y0
1
1
D=
y0 P V
dx +
(0) P V
dx +
(0),
2
2
2
2
2
2
2
x + y0
1
1 x + y0
and so
1
D=
PV
2
y0 (x)
1
dx + (0).
2
2
2
x + y0
Or, in the original form

Z
PV
D=
L
1
(x) n G dx + (0).
2
We can work along similar lines to show that as y0 approaches L from outside D, that is as
y0 0 ,
Z PV
1
D=
(x) n G dx (0).
2
L

To complete step (i) of the boundary integral method, we let the point x0 approach C from
the inside in (2.10) and then, using (2.11), we obtain
Z PV
Z
1
(x0 ) =
(x) [n G(x, x0 )] dl(x) + (x0 )
G(x, x0 ) n (x) dl(x),
2
C
C
where x0 lies on C. So,
Z PV
Z
1
(x0 ) =
(x) [n G(x, x0 )] dl(x)
G(x, x0 ) n (x) dl(x),
2
C
C
where the unit normal to the boundary C, n, points inside the solution domain D. We now
introduce the flux q, which is defined4 so that q n . Then we have
Z PV
Z
1
(x0 ) =
(x) [n G(x, x0 )] dl(x)
G(x, x0 ) q(x) dl(x),
(2.13)
2
C
C
4
Why do we call this the flux? Ficks law states that a substance flows down its concentration gradient
(like heat, for example). If represents the concentration of some substance, then is its flow rate
according to Ficks law. Then the flux out of a closed region with inward pointing normal n will be n .
13
MTH-ME82
where x0 lies on C.
If is prescribed on the boundary, then we have a Fredholm integral equation of the first
kind for the flux q,
Z
G(x, x0 ) q(x) dl(x) = F (x0 ),
C
where
Z
PV
F (x0 ) =
(x) [n G(x, x0 )] dl(x)

C
1
(x0 ).
2
It is called an integral equation of the first kind since the subject (in this case q) appears
inside the integral but not on the right hand side.
If instead q is prescribed on the boundary C, we have a Fredholm integral equation of the
second kind for the boundary distribution of ,
1
(x0 ) =
2
PV
(x) [n G(x, x0 )] dl(x) + (x0 ),

C
where
Z
(x0 ) =
G(x, x0 ) q(x) dl(x).

C
It is called an integral equation of the second kind because the subject (in this case )
appears both inside the integral and outside of it.
Example problem 2.1
The function satisifes Laplaces equation in the domain D inside the circle C of radius a,
2 = 0,
with the boundary condition = 1 on C. Formulate the solution for using the BIM.
First let us note that the exact solution, which is regular everywhere inside the circle,
is = 1. Now, we start with the BIE for a point x0 lying on C, namely (2.13),
1
(x0 ) =
2
PV
n G G q dl.
C
where q = n and noting that = 1 on C we have

1
=
2
PV
n G G q dl.
C
We will use polar coordinates (, ) with origin at the centre of the circle. The unit normal
points into D in the negative radial direction, and so

Z 2
Z 2
1
G
= PV
a d a
G q d.
2
=a
0
0
14
MTH-ME82
We will utilise the free space Greens function, G = (1/2) log r, where r = [(x x0 )2 +
(y y0 )2 ]. It will be useful to note that when both (x, y) and (x0 , y0 ) lie on C, using the
cosine rule we have
r2 = 2a2 (1 cos ),
(2.14)
where is the angle subtended between the points (x, y) and (x0 , y0 ) and the origin. Also,
when (x0 , y0 ) lies at a general point in space,
r2 = a2 + 2 2a cos ,
where again is the angle subtended between the points (x, y) and (x0 , y0 ) and the origin.
Hence
G
G r
1 a cos
1 a cos
.
=
=
=
r
2r
r
2
r2
So on = a,
G
a 1 cos
1
=
=
.
2
2
r
4a
Therefore the BIE becomes
Z 2
1
a
1
q log r d
= +
2
2 2 0
and so
Z 2
q log r d = 0.
(2.15)
0
This is satisfied if we take q = 0 on = a. According to the exact solution, = 1, we have

q( = a) = n
=
= 0,
=a
=a
and so the two solutions are consistent.

Cautionary warning: When we reformulate a partial differential equation, such as
Laplaces equation, as an integral equation, as we are doing here, we introduce the possibility of spurious solutions which satisfy the integral equation but not the original PDE.
To illustrate this, let us revisit the example above. We start by noting that when both x
and x0 are on C, using (2.14),
Z 2
Z
Z
1 2
1 2
log r d =
log[2a2 (1 cos )] d = log(2a2 ) +
log(1 cos ) d.
2 0
2 0
0
Exercise: Show that
log(1 cos ) d = 2 log 2.

0
So
Z
log r d = 2 log a.
0
15
MTH-ME82
For the example problem, we derived the constraint on the boundary flux (2.15),
Z
q log r d = 0.
0
Consider what happens if we take the circle to be of unit radius, so that a = 1. Then we
have that
Z 2
log r d = 0,
0
and so we may choose q = , where is any constant, to satisfy (2.15). But we know from
the exact solution that the only admissible value is = 0. This serves as a warning to be
on our guard when dealing with the integral formulation of a PDE.
2.2 The BIM for potential flow at a corner
The arguments in the previous section relied on the use of the divergence theorem,
which is applicable provided that the boundary C is piecewise smooth5 . However, there
is a difficulty with allowing the point x0 to approach a point on a curve, such as a sharp
corner, where the normal vector is not continuous. In this case, the argument needs some
modification.
In practice, we are often interested in computing solutions in domains which involve corners
or kinks, for example a rectangular domain. In this section, we address how to adapt the
boundary integral method for potential flow to such regions.
To proceed, we go back to the beginning of the argument: Greens second identity. This
states that
2 2 = ( ).
Integrating over the domain, D, we obtain
ZZ
Z
ZZ
2
2
dS =
( ) dS = ( n n ) dl
D
by the divergence theorem, which we note is still valid for a boundary with a corner.
As before, satisfies Laplaces equation, 2 = 0. For the sake of argument, we take to
be the free-space Greens function,
=G=
1
log r,
2
satisfying
2 G + (x x0 ) = 0.
5
A curve is smooth if its normal vector is a continuous function of position. Such a curve or surface is
also called Lyapunov. A curve is piecewise smooth if its normal vector is continuous at finitely many points.
One such example is a square (e.g., Kreysig 10.5).
16
MTH-ME82
So,
ZZ
2 G dS =
Z
(G n n G) dl.
C
Consider now a point x0 which lies outside of the domain D. It follows that 2 G = 0
leaving
Z
1
(G n n G) dl.
0=
2 C
Substituting for G,
1
0=
2
Z
(log r n n log r) dl.
C
Since r = |x x0 |,
log r =

x x y y x x
0
0
0
,
=
.
,
log r =
x y
r2
r2
r2
Thus,
1
0=
2
Z
(log r n n
C
(a)
x x0
) dl
r2
(2.16)
(b)
x
0
n
n
x0
Figure 1: A domain with a sharp corner of angle .

Now consider a domain including a sharp corner of local angle as shown in figure 1. A
close-up of the corner is shown to the right of the figure. To proceed with the formulation
of the method, we deform the contour locally to skirt round the inside of the corner through
a small arc of radius as shown in figure 1(b). We will call this circular arc . Note that
the point x0 lies outside of the deformed contour, so the most recent equation is applicable.
Over , the contribution from the right hand side is
Z
1
x x0
(log r n n
) dl.
2
r2
17
MTH-ME82
On the arc, r = , x x0 = cos and y y0 = sin , where is the angle shown, and
dl = d giving,
Z
x x0
1
(log n n
) d.

2
2
0
The unit normal to the arc
n=
x x0
,
r
leaving
1

2
Z
0
|x x0 |2
1
(log n
) d =

3
2
Also,
n =
giving
1

2
(log n
0
) d.

x x0
x x0 y y0
=
+
,
r
r
x
r
y

log cos
+ sin
d
x
y

(since x x0 = cos and y y0 = sin )

Z
1
=
log cos
+ sin
d.
2 0
x
y
Now we let the arc shrink to the corner point by taking the limit 0, i.e. we consider

Z
1
log cos
+ sin
d
lim
0 2 0
x
y
Note that lim0 log = 0. So,

Z
Z
1
x x0 y y0
1
lim
log
+
(x) d = (x0 )
d = (x0 ).
0 2 0
r
x
r
y
2
2
0
In summary,
1
2
So,
Z
(log r n n
1
2
x x0
) dl =
(x0 ).
2
r
2
x x0
) dl
r2
C
Z
Z
1
x x0
1
x x0
=
PV
(log r n n
) dl + lim
(log r n n
) dl,
2
0
2
r
2
r2
C
(log r n n
where the principal value P V refers to the fact that the small arc of radius has been
excluded from the integration.
Thus, (2.16) becomes
1
2
Z
(log r n n
C
18
x x0
) dl
r2
MTH-ME82

1
=
PV
2
Rearranging,
Z
(log r n n
C
(x0 ) =
PV
2
2
x x0
) dl
(x0 ) = 0.
2
r
2
Z
(log r n n
C
x x0
) dl
r2
The previous equation was derived working with the free-space Greens function. But it is
equally applicable to any Greens function. So we may write in general,
Z
Z PV
(x0 ) =
G(x, x0 ) n (x) dl(x) +
(x) n G(x, x0 ) dl(x).
2
C
C
Note that the P V is only required on the 2nd integral as the first is continuous as x0
approaches the contour.
nb: If = , then there is no corner, and we recover our original boundary integral equation
(2.13),
Z
Z PV
1
(x0 ) =
G(x, x0 ) n (x) dl(x) +
(x) n G(x, x0 ) dl(x).
2
C
C
Summary: Potential Flow

In this section we have:
starting from Greens second identity, derived the boundary integral equation (BIE)
for potential flow valid for a point x0 inside the domain of interest.
introduced the notion of a Greens function
using the notion of a principal value, obtained the form of the BIE valid for a point
x0 lying on the contour enclosing the domain of interest.
derived the form of the BIE valid for a contour with a sharp kink
3. The boundary element method (BEM) for potential flow

We have now established how to reformulate a partial differential equation as an integral
equation, and we have seen that the advantage in doing this is that the new formulation is
defined over a boundary of any shape. In contrast, traditional methods of solving PDEs,
for example separation of variables, usually require the boundary geometry to be relatively
simple.
Aim: The aim of this section is to see how to go about solving an integral equation for
potential flow numerically. The procedure is known as the Boundary Element Method
(BEM).
19
MTH-ME82
We demonstrate the method by way of a specific example. In short, the idea is to discretise
the boundary C of the solution domain D using a sequence of straight or curved boundary
elements over each of which the solution is approximated using constants or polynomial
approximations.
Example problem 3.1
Consider the inviscid, irrotational motion of fluid through a circle C of radius a. Inside the
circular domain, D, the velocity potential, , satisfies Laplaces equation 2 = 0. Fixing
polar coordinates (, ) with origin at the centre of the circle, fluid flows into and out of D
such that the normal component of velocity on C is given by

n =
= sin ,

=a
where n = e , where e is the unit vector pointing in the positive direction, points into
D. Thus fluid flows out of D through = a, 0 , and fluid flows into D through
= a, 2.
We reformulate the problem as an integral equation. Taking the BIE with the point,
x0 , placed on C, we make use of equation (2.13), namely
1
(x0 ) =
2
PV
Z
(x) [n G(x, x0 )] dl(x)
G(x, x0 ) q(x) dl(x).

C
We will use the free space Greens function,

G=
1
log r,
2
where r = |xx0 |. Using the information given in the question, we have q = n = sin ,
and so the BIE becomes

Z 2
Z 2
1
a
G
(a, ) = a P V
(a, )
d
log r sin d.
2
=a
2 0
0
Here is the angle subtended between x, x0 and the origin, and = is the location of
the point x0 on the circle C. Using the cosine rule we know that
r2 = a2 + 2 2a cos( ),
r2 = 2a2 [1 cos( )],
when takes a general value and = a respectively. Also

G
G r
1 a cos( )
1 a cos( )
=
=
=
.
r
2r
r
2
r2
and so

=
=a
a 1 cos( )
1
=
.
2
2
r
4a
So we have
20
MTH-ME82
Figure 2: Discretisation of a circle with straight elements: (left) N = 6 and (right) N = 12.
1
1
(a, ) =
2
4
Z
0
a
(a, ) d
2
log r sin d,
0
and so
Z 2
Z
1
a 2
(a, ) d (a, ) =
log r sin d.
2 0
0
This is a Fredholm integral equation of the second kind, because the unknown function
appears both inside and outside of the integral. For the integral on the right hand side,
Z 2
Z
1 2
sin log r d =
sin log{2a2 [1 cos( )]} d
2 0
0
Z 2
Z
1 2
= log(2a2 )
sin d +
sin log[1 cos( )] d
2 0
0
Z
1 2
=
sin log[1 cos( )] d.
2 0
So the BIE is given by
1
2
where
F () =
(a, ) d (a, ) = F (),

0
a
2
sin log[1 cos( )] d.

0
Now (see Problem Sheet),

F () = a sin .
Hence the BIE is
(a, )
1
2
(a, ) d = a sin .
(3.1)
For this problem, the exact solution is

= + sin ,
21
(3.2)
MTH-ME82
for any constant . Substituting this into the BIE, we confirm that it satisfies the equation.
Note in particular that we can add an arbitrary constant to in (3.1) without altering the
equation. We will return to this point below.
Assuming that we dont know the exact solution, we proceed by discretising the boundary with a sequence of boundary elements. We will use straight elements, and approximate
the circle with a regular polygon with N sides, with each vertex of the polygon lying on
C. Next we approximate the value of along the ith side/element with a constant i ,
for i = 1 . . . N . Then if i is the angle subtended by the ith element and the origin, we
approximate the BIE (3.1) by replacing it with
N Z
1 X i+1
(a, )
i d = a sin ,
2
i
i=1
where N +1 = 1 . Notice that we have approximated the integral with a sum over the N
elements. Since were assuming that the i are constant, we have
Z i+1
N
1 X
i
d = a sin ,
(a, )
2
i
i=1
and so
(a, )
N
1 X
i di = a sin ,
2
(3.3)
i=1
where
di = i+1 i .
Next we apply (3.3) at the mid-point of each element to obtain the N equations
N
1 X
j
i di = a sin j ,
2
(3.4)
i=1
for j = 1, . . . , N , where
1
j = (j + j+1 ),
2
and j = (a, j ). We may write the system of N equations in (3.4) as a matrix system
A x = b,
(3.5)
where
x = (1 , 2 , . . . , N )T ,
b = a (sin 1 , sin 2 , . . . , sin N )T ,
and A is a N N matrix which has the appropriate structure.

Problem: We noted above that the integral equation (3.1) is unchanged by adding a
constant to and so it has an infinite number of solutions. We should expect the linear
system (3.5) to have the same property, that is for the matrix A to be singular.
22
MTH-ME82
BEM solution with N=16 straight elements

BEM
Exact
1.5
0.5
(a, )
0.5
1.5
2.5
Figure 3: BEM solution to integral equation (3.3) with N = 16 elements and a = 2.0.
To get around the problem that detA = 0 and to obtain a solution to the problem, we will
fix the value of to be zero at the mid-point of the first element (this is tantamount to
fixing the constant in the exact solution 3.2), and solve the algebraic equations (3.4) for
j = 2, . . . , N . So now we have the adjusted linear system:
1 = 0,
x
= b,
A
(3.6)
where
x = (2 , . . . , N )T ,
b = a (sin 2 , . . . , sin N )T ,
is a (N 1) (N 1) matrix which has the appropriate structure (see Matlab

where A
code BEMexample1.m).
We solve the linear system in (3.6) for x
, using Gaussian elimination for example,
to determine the values of on the elements. This gives us an approximation to the
solution. As we increase N and take polygons with more and more sides, we more closely
approximate the circular boundary, and more closely approximate the true solution on the
boundary, = a sin . The result of an example calculation carried out using the Matlab
code BEMexample1.m with N = 16 elements for a = 2.0 is shown in figure 3.
3.1 The BEM for Poissons equation
Here we briefly discuss methods to solve Poissons equation, namely
2 = s(x),
23
(3.7)
MTH-ME82
using the boundary element method.

Starting from Greens second identity,
2 2 = ( ).
we proceed as in the formulation for Laplaces equation, choose G to satisfy
2 G + (x x0 ) = 0,
(3.8)
and derive
Z
ZZ
(x0 ) =
n(x) G(x, x0 ) G(x, x0 ) n(x) dl(x) +
s(x)G(x, x0 ) dA(x), (3.9)
C
which is the corresponding form of (2.9) for Poissons equation. The last term in (3.9)
threatens to compromise the method by introducing an integral over the entire domain D.
To get round this problem, we re-express this term as a line integral by making use of
Greens theorem in the plane, namely

ZZ
I
M
L
dA = (Ldx + M dy).
x
y
C
Example
We illustrate the last point for the simplest case when s is a constant and G is the free-space
Greens function, G = (1/2) log r. We select

1
x
1
y
M =
x
log r
,
L=
y log r
.
(3.10)
4
2
4
2
Noting that rx = x
/r, we have that

1
x
2 1
Mx =
log r + 2
,
4
r
2
Then
1
Mx Ly =
4
So we can now write

ZZ
1
Ly =
4
y2 1
log r + 2
r
2

x
2 + y2
1
1 = log r = G.
2 log r +
r2
2
1
G(x, x0 ) dA(x) =
4
D
Z
Ldx + M dy =
a dx,
C
where a = (L, M ). Alternatively,

Z
ZZ
Z
1
Ldx + M dy =
a t dl,
G(x, x0 ) dA(x) =
4 C
D
C
24
MTH-ME82
where t is the unit tangent to the curve C pointing in the direction of increasing arc length
l. We may then proceed using the BEM to solve the integral equation (3.9) now written in
the form
Z
Z
(x0 ) =
n(x) G(x, x0 ) G(x, x0 ) n(x) dl(x) + s
a t dl,
(3.11)
C
where a = (L, M ) and M , L are given by (3.10).

As an alternative to the preceding method we could instead break down into two parts
as follows,
= + w,
(3.12)
where 2 = 0 and 2 w = s(x). Then (3.12) solves (3.7) by construction. We would

then proceed by solving Laplaces equation for as above, remembering to modify the
boundary conditions to account for the decomposition (3.12). Depending on the form of
the inhomogeneity s(x) it may be a straightforward matter to obtain a solution for w in
closed form, in which case the solution to the problem is completed without difficutly.
Example problem 3.2
Solve the Poisson equation 2 = 1 on a circular disk D of radius a with boundary C
subject to the boundary condition = a2 /4 on C.
We re-formulate as an integral equation using (3.11). Working similarly to Example Problem
2.1, we arrive at the integral equation for the boundary flux q,
Z
Z
q G(x, x0 ) dl =
a t dl,
(3.13)
C
where, choosing to use the free-space Greens function,

G=
1
log r,
2
r = |x x0 |,
we have a = (L, M ) with L and M defined in (3.10). The exact solution to the problem is
=
r2
.
4
This gives q = r |r=a = a/2. Table 3.1 summarises results obtained from a BEM computation using N = 80 straight elements to discretise C:
Notice that the case a = 1 is excluded from the Table. This is because the integral equation
(3.13) does not have a unique solution for a unit disk (see Example Problem 2.1 for further
details).
4. The boundary integral method (BIM) for Stokes flow
25
MTH-ME82

a
2
3
4
q
0.9996
1.4992
1.9988
Exact
1
1.5
2.0
Table 3.1: Comparison of BEM-computed q (with N = 80 straight elements) and the

exact value for the Poisson problem of Example Problem 3.1.
The equations describing creeping fluid motion are called the Stokes equations:
0 = P + F + 2 u,
u = 0,
(4.1)
where is the fluid density, is the viscosity, and F is a body force acting per unit
mass (e.g., gravity). The first equation is the momentum equation for Stokes flow and is
an expression of Newtons second law of motion. The second equation is the continuity
equation expressing conservation of mass.
You will note that equation (4.1) is linear. This is a crucial property for the theory which
we shall develop in the following lectures.
4.1 Derivation of the governing equations
4.1.1 Internal stress in a fluid
What is the stress at any point in a fluid due to the surrounding fluid? Consider a fluid
and mentally isolate within it a parcel of fluid of volume V and surface S. The force acting
on an element of the surface dS with unit outward normal n is defined to be
(ij nj )dS
or
( n)dS.
So the total force acting on the surface S is

ZZ
ij nj dS.
S
Note that summation over j is implied in the above equation, according to the Einstein
summation convention discussed in section 1.1, so it really means
ZZ X
3
ij nj dS.
S j=1
We will be using the summation convention throughout.

4.1.2 Conservation of mass
26
MTH-ME82
Consider a closed volume V which is fixed in space inside a moving fluid. As the
fluid moves past, the rate of change of the mass of fluid within V must equal the net
amount of fluid coming into (or out of) it, assuming there are no sources of fluid within V .
Mathematically,
ZZ
ZZZ
d
u n dS = 0,
dV +
dt
S
V
where is the fluid density, S is the volume surface, and n is the normal to V pointing
outwards. Using the fact that V is fixed in space and applying the divergence theorem, this
becomes
ZZZ n
o
+ ( u) dV = 0.
t
V
Since our choice of volume V was arbitrary, it must be true that
+ (u) = 0.
t
(4.2)
u = 0.
(4.3)
If is constant, this reduces to
Equation (4.3) therefore expresses the principle of mass conservation. It is usually referred
to as the continuity equation.
4.1.3 Newtons second law
We now derive the equation of motion of a slow moving viscous fluid. Let F represent
some body force per unit mass acting on the fluid (e.g., gravity). We apply Newtons second
law of motion to a volume of fluid V (t), with surface area S(t) and unit outward normal n,
moving with the flow.
Following this line of argument, and ignoring any inertia present in the fluid, we eventually derive the Stokes equations of fluid motion.
u = 0,
0 = F P + 2 u.
(4.4)
These are the Stokes equations for an incompressible Newtonian fluid moving under the
conditions of Stokes flow.
Alternatively, it can be shown that these can be written as
u = 0,
0 = F ,
with the understanding that means

ij
.
xj
Here,
ij = P ij + 2eij ,
27
MTH-ME82
where
1
eij =
2
uj
ui
+
xj
xi

.
Note that in index notation, u is written as

ui
.
xi
4.2 The Lorentz reciprocal relation
The boundary integral method for potential flow discussed came from Greens second
identity (2.1). The boundary integral method for Stokes flow emereges from a similar, but
less familiar, identity which is valid for a velocity and pressure field which satisfy the Stokes
equations of motion (4.1). It is called the Lorentz reciprocal relation6 .
To derive this relation, consider two different, independent regular7 Stokes flows with ve0 . We have
locity fields ui and u0i , and respective stress tensors, ij and ij
u = 0,
=0
We form the product

u0i
u0 = 0,
and
0 = 0
ij
u0
=
(u0i ij ) ij i ,
xj
xj
xj
by the product rule of differentiation. Then

0
uj
ui
ui
0
0 ij
ui
=
(ui ij ) P ij +
+
,
xj
xj
xj
xi
xj
u0
(u0i ij ) + P i
xj
xi
uj
ui
+
xj
xi
u0i
.
xj
So
u0i
ij
=
(u0 ij )
xj
xj i
uj
ui
+
xj
xi
u0i
,
xj
by continuity. By swapping primed and unprimed variables we have

0

0
ij
u0j ui
ui
0
ui
=
(ui ij
)
+
.
xj
xj
xj
xi xj
(4.5)
(4.6)
Now (4.5)-(4.6) yields

ij
u0i
xj
0
ij
0
ui
=
(u0 ij ui ij
).
xj
xj i
6
This was derived by H. A. Lorentz in 1907 (his paper is reproduced in the Journal of Engineering
Mathematics 1996, volume 30).
7
By regular we mean there are no singular points in the domain of interest, that is neither of the velocity
fields or the stress fields blow up in the domain.
28
MTH-ME82
But, by the momentum equation (4.4) for each flow,

0
ij
= 0,
xj
ij
=0
xj
and we are left with
0
(u0 ij ui ij
) = 0.
xj i
(4.7)
This is the Lorentz reciprocal identity for Stokes flow. Note that, according to the Einstein
convention discussed in section 1.1, there is summation here over both i and j. In vector
form it is written as
(u0 u 0 ) = 0.
Note: To derive the above we assumed the two flows were for fluids with the same viscosity.
To generalize to flows with different viscosities, see the Problem Sheet.
4.3 Greens functions for Stokes flow
Before proceeding to the boundary integral formulation, we must first discuss the notion
of a Greens function for Stokes flow. As we saw for Laplaces equation, a Greens function
is the solution to a singularly forced differential equation (this is just a fancy way of saying
the equation has a delta function in it).
For Laplaces equation, we had
2 G + (x x0 ) = 0.
Including a singular forcing term in the Stokes equations we have
0 = P + b (x x0 ) + 2 u.
Here the singular forcing term is given by
b (x x0 ).
It represents a point force acting in the direction b at the point x0 . Here
b = bx i + by j + bz k
is a constant vector.
b
x0
29
(4.8)
MTH-ME82
We write the solution of (4.8) in the general form

ui (x) =
1
Gij (x, x0 )bj ,
8
(4.9)
where the Greens function Gij is yet to be determined. The factor of 1/8 is included
purely for convenience.
nb: The solution (4.9) gives the fluid velocity u at a general point in space x provoked by
the point force of strength b acting at the point x0 .
We write the solution for the pressure as
P =
1
pj bj
8
(4.10)
where pj is yet to be determined. The factor of 1/8 is purely for convenience.

In index notation, the momentum equation (4.4) is
0=
P
2 ui
+ bi (xi x0,i ) + 2 .
xi
xk
Substituting in the Greens function solution (4.9) and (4.10), we find

0=
1 pj
1 2 Gij
bj + bi (xi x0,i ) +
bj .
8 xi
8 x2k
This can be re-written as

0=
or
1 2 Gij
1 pj
bj + ij bj (xi x0,i ) +
bj ,
8 xi
8 x2k
h
1 pj
1 2 Gij i
,
0 = bj
+ ij (xi x0,i ) +
8 xi
8 x2k
Assuming bj is not the zero vector, we have leave

8(xi x0,i ) =
pj
2 Gij
+
.
xi
x2k
(4.11)
So we must ensure that the Greens function Gij and the pressure vector pj satisfy this
equation.
In due course we will consider specific examples of Greens functions. First, we discuss
generic properties common to all Greens functions.
4.4 Properties of Greens functions for Stokes flow
4.4.1 Conservation of mass
30
MTH-ME82
We have shown that a Stokes flow Greens function, Gij is a solution of equation (4.11),
together with a suitable choice for the pressure vector pj . The velocity field associated with
the Greens function was given above,
ui = Gij bj .
Taking the divergence, we have
ui
=
(Gij bj ) = 0,
xi
xi
since
ui
= 0,
xi
by continuity. If we integrate over a closed volume V with unit outward normal ni , then
ZZZ
ZZZ
ui
dV =
(Gij bj ) dV = 0.
V xi
V xi
Applying the divergence theorem8 , we obtain
ZZZ
ZZ
(Gij bj ) dV =
ni Gij bj dS = 0.
V xi
S
Since b is a constant vector, this shows that
ZZ
Gij (x, x0 )ni dS = 0.
(4.12)
This integral identity is satisfied by all Greens functions and effectively expresses conservation of mass of the flow.
4.4.2 Symmetry of the Greens function
A Greens function is symmetric in the following sense:
Gij (x, x0 ) = Gji (x0 , x).
i.e. swapping i and j and also x and x0 leaves the Greens function unchanged.
Proof: See Appendix C.
5. Examples of Greens functions for Stokes flow
Greens functions are the cornerstone of the boundary integral method. In this section,
we discuss some examples.
5.1 The free space Greens function or stokeslet
8
For a more formal argument, see Appendix D
31
MTH-ME82
Consider a point force in free space. We call the associated Greens function the free
space Greens function or stokeslet. The analysis depends on whether we are in two or three
dimensions.
5.1.1 The stokeslet in three-dimensional flow
A Stokes Greens function is a solution of the singularly forced equation
0 = P + b (x x0 ) + 2 u.
(5.1)
The solution is written in the form

ui =
1
Gij (x, x0 )bj ,
8
where the vector b is known and Gij is to be found.

Definition: The stokeslet , or free-space Greens function for Stokes flow (sometimes called
the Oseen-Burgers tensor), is a solution of (5.1) when there are no boundaries in the flow,
that is when the point force at x0 acts in an unbounded fluid.
It can be shown that (see Appendix E) the three-dimensional stokeslet has the form
Gij =
x
i x
j
ij
+ 3 .
r
r
(5.2)
where x
= x x0 and r = |x x0 | = |
x|. The stokeslet induces a pressure field in the fluid
given by
1
P =
pj bj ,
8
where
2
xj
pj = 3 .
r
Example: Suppose that x0 = (0, 0, 0) and that b = 8(1, 0, 0) so that the point force
points along the x axis. Then
uj = G1j .
So,
u1 =
1 x2
1
x2
+ 3 = 2
+
.
r
r
(x + y 2 + z 2 )1/2 (x2 + y 2 + z 2 )3/2
If we set x = z = 0, then
u1 |x=z=0 =
1
.
y
Thus the horizontal velocity component induced by the point force decays like 1/y. In
general the velocity induced by a point force decays like the inverse distance from the point
at which the force is acting.
5.1.2 The stokeslet in two-dimensional flow
32
MTH-ME82
In two dimensions the stokeslet takes the form

Gij = ij log r +
x
i x
j
,
r2
(5.3)
and the corresponding velocity field is

ui =
1
Gij (x, x0 )bj ,
4
where the 4 is included purely for convenience. The stokeslet induces a pressure field in
the fluid given by field is
1
pj (x, x0 )bj ,
P =
4
where
x
j
pj = 2 2 .
r
5.2 Greens function for flow above a plane wall
Sometimes it is convenient to work with a Greens function which vanishes over a boundary of the flow. In this section we compute the Greens function Gij such that
G(x, x0 ) = 0,
whenever the point x lies on the wall at x = 1.
The correct form of the Greens function is given by Blake (1971) to be
D
IM
SD
IM
Gij (x, x0 ) = Sij (x x0 ) Sij (x xIM
0 ) + 2Gij (x x0 ) 2Gij (x x0 )
(5.4)
where x0 is the singular point, and xIM

= (x0 2, y0 , z0 ) is its image reflected in the plane
0
wall, and
ij
xi xj
+
,
|x|
|x|3

ij
xi xj
D
Gij (x) =
3 5 ,
|x|3
|x|

j1 xi i1 xj
SD
D
Gij (x) = x1 Gij (x)
.
|x|3
Sij (x) =
The plus sign applies for j = 2, 3 and the minus sign for j = 1. Sij is the stokeslet discussed
above; GD is called a potential dipole and GSD is called a stokeslet doublet.
The construction of (5.4) is complicated and we will not go through it here. Instead we will
simply check that it has the desired property,
Gij (x, x0 ) = 0,
33
MTH-ME82
when x = (1, y, z), i.e. that the Greens function always vanishes on the wall at x = 1.
We will check by scrutinising the components of Gij . For simplicity, we consider the case
when x0 = (0, 0, 0). In this case,
xIM
= (2, 0, 0),
0
|x x0 | = |x xIM
0 |.
Henceforth we write |x x0 | = R. With our choice of x0 , we have simply

x x0 = (x, y, z) = x
and
x xIM
= (x + 2, y, z).
0
Since we are only interested in points lying on the wall at x = 1, we have
x xIM
= (1, y, z).
0
x x0 = (1, y, z),
Consider Gyy :
Syy (x x0 ) = Syy (x) =
1
y2
+ 3
R R
1
y2
+ 3 = Syy (x x0 )
R R

21 x2 21 x2
SD
IM
D
IM
Gyy (x x0 ) = 1 Gyy (x x0 ) +
,
R3
Syy (x xIM
0 )=
and so
IM
D
IM
GSD
yy (x x0 ) = Gyy (x x0 ).
The full Greens function gives

D
IM
SD
IM
Gyy (x x0 ) = Syy (x x0 ) Syy (x xIM
0 ) + 2Gyy (x x0 ) 2Gyy (x x0 )
IM
D
IM
= Syy (x x0 ) Syy (x x0 ) + 2GD
yy (x x0 ) 2Gyy (x x0 )
=0
as required.
Next, consider Gxy :
Sxy (x x0 ) = Sxy (x) =
1
y2
3
R R
since x x0 = (1, y, z).

Syy (x xIM
0 )=
34
y
1
+
R R3
MTH-ME82
since x xIM
= (1, y, z).
0
GSD
xy (x
xIM
0 )
=1
GD
ij (x
xIM
0 )

+
21 x1 11 x2
R3
so
IM
D
IM
GSD
xy (x x0 ) = Gxy (x x0 )

,
y
,
R3
So the full Greens function gives

D
IM
SD
IM
Gxy (x x0 ) = Sxy (x x0 ) Sxy (x xIM
0 ) + 2Gxy (x x0 ) 2Gxy (x x0 )

=
1
y2
3
R R
1
y
+ 3
R R
D
+ 2GD
xy 2Gxy + 2
y
.
R3
=0
as required.
By the symmetry property of the Greens function, Gyx will also be zero. See the Problem
Sheet for the remaining components.
Notes:
1. In two dimensions one may derive the expression for a stokeslet above a plane wall using
a complex variable formulation which makes use of so-called Goursat functions. A very nice,
and easy-to-follow description of this is given by Crowdy & Or (2010) in Physical Review
E, volume 81, manuscript number 036313.
2. A recently published paper in J. Fluid Mechanics gives a simpler solution for the stokeslet
above a plane wall than that given in this section. The simpler solution makes use of
Papkovich-Neuber potentials. The full citation is:
Gimbutas, Z., Greengard, L. & Veerapaneni, S. (2015) Simple and efficient representations
for the fundamental solutions of Stokes flow in a half-space, J. Fluid Mech., 776.
5.2 The Greens function stress tensor
In three-dimensions we expressed the velocity field due to a point force in the form,
ui =
1
Gij (x, x0 )bj ,
8
(5.5)
and we wrote the corresponding pressure field as

P =
1
pj (x, x0 )bj .
8
(5.6)
The stress in the fluid induced by the point force may be expressed in a similar way,
ik (x, x0 ) =
1
Tijk (x, x0 )bj .
8
35
(5.7)
MTH-ME82
For the stokeslet in three-dimensions, Tijk takes the form,

Tijk (x, x0 ) = ik pj (x, x0 ) +
Gkj
Gij
.
+
xk
xi
Substituting this into (5.7) and using (5.5), we can confirm by differentiation that

ui
uk
ik = P ik +
.
+
xk
xi
(5.8)
We note that the stress tensor, ik is symmetric (we can confirm this by swapping i and k
in 5.8). Since
1
ik (x, x0 ) =
Tijk (x, x0 )bj ,
8
and since ik = ki , it follows that
Tijk = Tkji .
5.2.1 Stress tensor identity in three dimensions

Using all of the information at our disposal, we can obtain the following identity
ZZ
8
Tijk (x, x0 )ni (x) dS(x) = 4 jk ,
(5.9)
D
0
where D is a closed surface, and the unit normal n points out of D. Note that
8 applies when x0 is inside D.
4 applies when x0 is on D.
0 applies when x0 is outside D.
Derivation of (5.9)
By definition (see section 4.1.1) the hydrodynamic force due to the surrounding fluid which
acts on a parcel of fluid of volume V , with bounding surface D and with unit outward
normal n, is
ZZ
ZZZ
ik
F =
ik ni dS =
dV,
D
V xi
by the divergence theorem. But the singularly-forced Stokes momentum equation states
that
ik
+ bk (x x0 ) = 0.
xi
and so,
ZZZ
F =
V
ik
dV = bk
xi
36
ZZZ
(x x0 ) dV.
V
MTH-ME82
Using (5.7),
ZZ
F =
1
ik ni dS =
bj
8
D
ZZ
ZZZ
Tijk ni dS = bk
D
(x x0 ) dV.
V
Writing bk = jk bj (see section 1.1), we have

ZZ
ZZZ
1
Tijk ni dS = bj jk
(x x0 ) dV.
bj
8
D
V
But since b is an arbitrary vector, it must be true that
ZZ
ZZZ
1
Tijk ni dS = jk
(x x0 ) dV,
8
D
V
and the result follows on using given properties of the delta function,
ZZZ
if x0 lies inside V
1
1
(x x0 ) dV =
if
x
0 lies exactly on S
2
V
0
if x0 lies outside V .
The middle property is explained in Appendix B.
5.2.2 Stress tensor identity in two dimensions
The two dimensional equivalent of the identity (5.9) is,
Z
4
Tijk (x, x0 )ni (x) dl(x) = 2 jk ,
C
0
(5.10)
where C is a closed contour, and the unit normal n points out of C. Note that
4 applies when x0 is inside C.
2 applies when x0 is on C.
0 applies when x0 is outside C.
5.3 The stress tensor for the stokeslet
The stress tensor corresponding to the three-dimensional stokeslet is
Tijk (x, x0 ) = 6
x
i x
j x
k
.
5
r
The stress tensor corresponding to the two-dimensional stokeslet is

Tijk (x, x0 ) = 4
5.4 Summary
In this section on Stokes flow, we have
37
x
i x
j x
k
.
4
r
(5.11)
MTH-ME82
Derived the governing equations

Computed the Greens function in free space
Discussed the Greens function with a plane wall boundary
Introduced the Greens function stress tensor Tijk (x, x0 )
6. The boundary integral equation for Stokes flow in three dimensions

We are now in a position to derive the boundary integral equation for Stokes flow. Our
starting point is the Lorentz reciprocal identity derived in section 4.2,
0
(u0 ik ui ik
) = 0.
xk i
(6.1)
Consider two different flows with velocity and stress fields (u, ) and (u0 , 0 ) respectively.
We choose u0 to be the velocity field induced by a point force of strength b at the point x0 ,
i.e.
1
1
0
(x) =
Gij bj ,
ik
Tijk bj .
u0i =
8
8
Substituting into (6.1), we obtain,
(Gij ik ui Tijk ) bj = 0.
xk
Dropping the arbitrary vector bj , we have
(Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )) = 0.

xk
(6.2)
We now subdivide the argument into two sections.

6.1 When x0 is outside V.
We integrate equation (6.2) over a closed volume V with boundary D, as shown in the
figure.
ZZ
[Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )] dV (x) = 0.

x
k
V
Applying the divergence theorem,
ZZ
[Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )] nk dS(x) = 0.
D
where the unit vector nj points out of the volume V as required by the divergence theorem.
38
MTH-ME82
D
V
D
n
n
For later convenience we redefine the unit normal to point into the volume V and write the
equation for convenience as
Z
Z
1
1
Gij (x, x0 )fi (x) dS(x) +
ui (x)Tijk (x, x0 )nk dS(x) = 0.
(6.3)
8 D
8 D
where fi = ik nk as usual.
nb: We were able to apply the divergence theorem because the integrand is analytic everywhere inside the domain of integration. This is a requirement of the theorem. The
integrand is analytic because x0 is outside V and so neither Gij nor Tijk become singular
inside V .
6.2 When x0 is inside V.
The divergence theorem requires that the function being integrated does not become
infinite anywhere inside the integration domain. If x0 is inside V , the Greens function
Gij (x, x0 ) diverges at x0 inside the integration domain. To get around this problem, we
integrate (6.2) over the volume enclosed by S excluding a a small sphere of radius 1,
with surface S and volume V , which is centred on the point x0 as shown in the figure. We
call this new region V 0 = V V .
D
V
S
x0
39
MTH-ME82
Performing the integration we have

ZZZ
(Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )) dV (x) = 0.

x
0
k
V
Notice that V 0 is a multiply-connected domain9 . Fortunately, the divergence theorem can
be easily extended to such domains (see the Problem Sheet). Applying the divergence
theorem,
ZZ
[Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )] nk dS(x) = 0
(6.4)
S +D
where the unit normal is chosen to point into the volume of integration. The idea now is
to take the limit 0.
Consider the integral over S ,
ZZ
[Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )] nk dS(x).
I
(6.5)
S
The point x lies on the surface S . As decreases, |x x0 | = decreases and so x moves

closer to x0 .
Examining the behaviour of Gij and Tijk as x x0 , it is found that these always behave
as a stokeslet Greens function and a stokeslet stress tensor as x x0 . Hence, as 0
and x x0 , we have
Gij
ij
x
i x
j
+ 3 ,

Tijk 6
x
i x
j x
k
.
5

(6.6)
The normal to the surface of S is n = x x0 , so the unit normal is

n=
x x0
(x x0 )
=
.
|x x0 |

, so,
With our usual notation, x x0 = x
n=
x
,

ni =
x
i

The integral (6.5) becomes, as 0,

Z Z
ij
x
i x
j
x
i x
j x
k
I =
+ 3 ik (x) ui (x) 6
nk dS(x).

5
S
Taking spherical polars (r, , ), so that dS = r2 sin d d, we have

Z Z
ij
x
i x
j
x
i x
j x
k
+ 3 ik (x) ui (x) 6
nk 2 sin d d.
I =
5

S
9
A simply connected domain is one for which we may shrink a closed loop lying inside the domain down
to a point without any point on the loop leaving the domain (e.g. a circular disk is simply-connected). If
this is not possible the domain is multiply-connected. For example an annulus is multiply-connected for a
circular loop inside the domain cannot be shrunk to a point without leaving the annulus.
40
MTH-ME82
So,

x
i x
j
x
i x
j x
k
ij +
ik (x) + 6 ui (x)
nk sin d d.

3
S
ZZ
I =
Consider the first term in I :

ZZ
S

x
i x
j
ij +
ik (x)nk sin d d

As 0, x x0 and ik (x) ik (x0 ), so we can write,

ZZ
ZZ
x
i x
j
x
i x
j
ij +
ij +
ik (x)nk sin d d ik (x0 )
nk sin d d

S
S
Also, as 0,
x
i x
j
2
0.

So the whole of the first term vanishes in the limit 0.
So this leaves

ZZ
ZZ
x
i x
j
x
i x
j x
k
nk sin d d =
6 ui (x)
x
k nk sin d d.
I =
6 ui (x)
3
3
S
S
But
1
1 2
x
k nk = x
k x
k = |
x| = .

Therefore the remaining term is

ZZ
x
i x
j
6 ui (x)
2

S

sin dd.
In the limit 0, ui (x) ui (x0 ), so it becomes

ZZ
ZZ
x
i x
j
6
sin
dd
=
u
(x
)
x
i x
j dS,
6 ui (x0 )
i 0
2
4
S
S
since dS = 2 sin d d.
Now (see the Problem Sheet)
ZZ
4
x
i x
j dS = ij 4 ,
3
S
and so the integral over the small sphere is

I = 8 uj (x0 ).
41
MTH-ME82
Summary: We have just shown that

ZZ
lim
[Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )] nk dS(x) = 8 uj (x0 ).
0
S
So (6.4) becomes
ZZ
[Gij (x, x0 )ik (x) ui (x)Tijk (x, x0 )] nk dS(x) = 8 uj (x0 ).
D
Rearranging, we have
ZZ
ZZ
1
1
uj (x0 ) =
Gij (x, x0 )ik (x)nk dS(x) +
ui (x)Tijk (x, x0 )nk dS(x).
8 D
8
D
Definition: We define f = n to be the surface stress or surface traction.
Substituting for the traction, we have the boundary integral equation (BIE) for Stokes
flow
ZZ
ZZ
1
1
uj (x0 ) =
fi (x) Gij (x, x0 ) dS(x) +
ui (x)Tijk (x, x0 )nk dS(x), (6.7)
8 D
8
D
where the unit normal n points into the domain enclosed by D.
Notes:
1. Equation (6.7) is quite remarkable. It states that the velocity at a general point x0 in
the flow field can be expressed purely in terms of boundary values of the velocity and stress.
2. The choice of Greens function Gij and Greens function stress tensor Tijk is arbitrary.
The only requirement is that they satisfy the singularly forced momentum equation (4.8).
3. In equation (6.7), the boundary surface D can take any smooth shape.
5. Equation (6.7) applies for a point x0 inside but not on the boundary S. We will say
more on this below.
Terminology: In (6.7)
ZZ
fi (x) Gij (x, x0 ) dS(x)
is the single-layer potential
and
ZZ
ui (x)Tijk (x, x0 )nk dS(x)
D
6.3 When x0 is on D.
42
is the double-layer potential
MTH-ME82
We turn to the case when the point x0 lies on the boundary of the domain D. The
discussion here is analogous to that for the integral formulation of Laplaces equation described in section 2.1. We start by discussing the behaviour of the two potentials in (6.7)
as the point x0 crosses D.
6.3.1 The single-layer potential
The single-layer potential is given by
ZZ
SL
Ij (x0 ) =
fi (x) Gij (x, x0 ) dS(x).
D
It can be shown that IjSL is continuous as x0 crosses D.

6.3.1 The double-layer potential
The double-layer potential is given by
ZZ
DL
Ij (x0 ) =
(6.8)
We ask the question: How does I DL behave as the point x0 crosses the boundary D?
Example:
To illustrate the behaviour of the the double-layer potential as x0 crosses the flow boundary
D, it is simplest to give an example in two dimensions. In this case the DLP is written as
Z
DL
Ij (x0 ) =
ui (x0 )Tijk (x, x0 )nk dl(x).
C
and the flow boundary is the curve C.

For simplicity, lets consider the case when C is the x axis, y = 0 and choose x0 = (0, y0 ).
Then
Z
DL
Ij (x0 , y0 ) =
ui (0, y0 )Tijk (0, y0 )nk dx.
In two dimensions the free space Greens function stress tensor has the form
Tijk (x, x0 ) = 4
So,
IjDL
= 4

ui (0, y0 )
x
i x
j x
k
.
4
r
x
i x
j x
k
4
r

nk dx.
y=0
The unit normal is n = (0, 1) and so ni = 2i . We assume that the flow is purely in the x
direction, so that u = (U, 0). Then we have

Z
Z
x
1 x
j x
k
x
1 x
j x
2
DL
Ij = 4
U (y0 )
2k dx = 4
U (y0 )
dx.
r4
r4
y=0
y=0
43
MTH-ME82
But
x
1 = x1 x0,1 = x,
so
IjDL (y0 )
x
2 = x2 x0,2 = y y0 ,

U (y0 )
=4
xx
j y0
r4

dx.
Consider now the limit as y0 approaches zero from above, namely y0 0+ . We rewrite the
last equation as

Z
Z
xx
j y0
xx
j y0
DL
dx + 4 U (0)
dx.
Ij (y0 ) = 4
[U (y0 ) U (0)]
r4
r4
For small y0 ,
U (y0 ) = U (0) + y0 U 0 (0) +
Therefore,
IjDL (y0 )
Z
4
U (0)
xx
j y02
r4
dx + 4 U (0)
xx
j y0
r4

dx,
for small y0 , where

r2 = (x x0 )2 + (y y0 )2 = x2 + y02 .
Now, setting j = 1,
I1DL (y0 )
4 U (0)
y02
x2
dx + 4 U (0) (y y0 )
(x2 + y02 )2
x2
dx.
(x2 + y02 )2
Carrying out the integrations, we find

Z
2y0
x2
=
(x2 + y02 )2

2y0
if y0 > 0
if y0 < 0
and so for y0 small and positive

I1DL (y0 ) 2 U 0 (0) y0 + 2 U (0),
but for y0 small and negative
I1DL (y0 ) 2 U 0 (0) y0 2 U (0).
Thus
lim I1DL (y0 ) = 2 U (0),
lim I1DL (y0 ) = 2 U (0),
y0 0
y0 0+
and so the double-later potential jumps in value as y0 crosses C.

Exercise: Finish the above argument by considering the case j = 2.
44
MTH-ME82

In general, returning to a three-dimensional setting, we simply state the result that as x0
crosses D we have
ZZ PV
DL
lim Ij (x0 ) = 4uj (x0 ) +
ui (x)Tijk (x, x0 )nk dS(x),
(6.9)
x0 D+
where D+ denotes the fluid side of D (i.e. the region into which n points; see figure below)
and
ZZ PV
DL
(6.10)
lim Ij (x0 ) = 4uj (x0 ) +
x0 D
where D denotes the side of D exterior to the fluid.

dm
dp
dp
F
dm
Remembering our intention to derive the form of the BIE when x0 is on D, we recall the
BIE (6.7)
uj (x0 ) =
1 SL
1 DL
I (x0 ) +
I (x0 ),
8 j
8 j
where x0 is inside D. Taking the limit x0 D+ and using (6.9), we obtain

ZZ PV
1 SL
1
uj (x0 ) =
I (x0 ) +
4uj (x0 ) +
ui (x)Tijk (x, x0 )nk dS(x) .
8 j
8
D
Rearranging, we get
1
1 SL
1 DLP V
uj (x0 ) =
I (x0 ) +
I
(x0 ).
2
8 j
8 j
Or, in full,
1
1
uj (x0 ) =
2
8
1
8
D
ZZ
45
ZZ
PV
ui (x)Tijk (x, x0 )nk dS(x).(6.11)

D
MTH-ME82
This is the form of the BIE when x0 is positioned precisely on the boundary S.
Notes:
1. Equation (6.11) is almost identical to the BIE, except for the half factor on the LHS.
2. Equation (6.11) could equally we have been derived by starting with (6.3) and then
taking x0 D .
Summary for three-dimensional flow:

When the point x0 is inside D, the boundary integral equation is
ZZ
ZZ
1
1

ui (x)Tijk (x, x0 )nk dS(x)
8 D
8
D
When the point x0 is outside D, the boundary integral equation is
ZZ
ZZ
1
1
0=
8 D
8
D
When the point x0 is on D, the boundary integral equation is
Z
Z PV
1
1
1
uj (x0 ) =
2
8 D
8 D
6.4 The BIEs for two-dimensional flow
The equivalent boundary integral equations for two-dimensional flow bounded by a
contour C are as follows.
When the point x0 is inside C, the boundary integral equation is
Z
Z
1
1
uj (x0 ) =
Gij (x, x0 )fi (x) dl(x) +
ui (x)Tijk (x, x0 )nk dl(x)
4 C
4 C
(6.12)
When the point x0 is outside C, the boundary integral equation is

Z
Z
1
1
0=
fi (x) Gij (x, x0 ) dl(x) +
ui (x)Tijk (x, x0 )nk dl(x).
4 C
4 C
When the point x0 is on C, the boundary integral equation is
Z
Z PV
1
1
1
uj (x0 ) =
2
4 C
4 C
7. The Boundary Element Method for Stokes flow

46
(6.13)
MTH-ME82
We have now derived the boundary integral equation (6.11) which is valid when x0
lies on the boundary D. The procedure for solving the integral equation is very similar to
that discussed for solutions of Laplaces equation in section 3. The boundary geometry is
discretised into a set of flat/straight or curved elements and the unknown velocity or traction
is approximated over each element with constants or with polynomial representations.
Example: The driven cavity problem
This is a well-known problem in Fluid Mechanics and is often used as a benchmark for
checking a new numerical method for computing fluid flow. Consider a two dimensional
square box, as shown in the figure, which contains viscous liquid. The lid of the box moves
at a prescribed speed U in the x direction. The fluid velocity on the remaining walls of the
box is zero. The length of each cavity wall is equal to one.
y
U
The moving lid drives a flow within the cavity. Our goal is to compute the traction on the
walls created by this flow. We will label the box lid L and the other three walls W .
We will use the two-dimensional BIE (6.13), namely,
1
1
uj (x0 ) =
2
4
1
4
C
PV

C
Taking the point x0 to lie on the cavity walls, we have

1
1
4
C
1
2
PV
ui (x)Tijk (x, x0 )nk dl(x)

C
U 1j
if x0 is on L
(7.1)
if x0 is on W ,
where the unit normal n points into the fluid inside the box. Note that the contour C
includes both L and W .
47
MTH-ME82
For the double layer integral we have

1
4
PV
1
ui (x)Tijk (x, x0 )nk dl(x) =
U 1i
4
PV
Tijk (x, x0 )nk dl(x),

L
and so
1
4
PV
1
ui (x)Tijk (x, x0 )nk dl(x) =
U
4
PV
T1jk (x, x0 )nk dl(x) = Dj (x0 ),
say.
So the BIE (7.1) becomes

1
4
Z
fi (x) Gij (x, x0 ) dl(x) =
C
Dj (x0 )
1
2
U 1j
if x0 is on L
(7.2)
Dj (x0 )
if x0 is on W ,
We will use the stokeslet or free space Greens function (5.3), namely
Gij = ij log r +
x
i x
j
,
r2
= x x0 and r = |x x0 |. The stress tensor is given by (5.11),

where x
Tijk (x, x0 ) = 4
x
i x
j x
k
.
4
r
We discretise each side of the contour C using straight elements of equal length
1
4
= ,
N/4
N
where N (a multiple of 4) is the total number of elements. We will label the elements
clockwise starting from the left end of the lid as element 1.
We approximate the traction f as a constant over each element, setting f k as the constant
on the kth element.
The discretised form of the BIE (7.2) is
1
4
N
X
fik (x) Iijk (x0 ) =
Dj (x0 )
k=1
1
2
if x0 E k for 1 k
U 1j
N
4
(7.3)
if x0 E k for
Dj (x0 )
where E k denotes the kth element, and where

Z
Z
Iijk (x0 ) =
Gij (x, x0 ) dl(x) =
Ek
Ek
Left vertical elements

48

ij log r +
N
4
+1k N
x
i x
j
r2

dl(x).
MTH-ME82
For a vertical element parallel to the y axis at x = 0, we have

Z yk+1
x
i x
j
k
ij log r + 2
Iij (x0 ) =
dy,
r
yk
x=0
where yk is the y coordinate of the lowermost point of the kth vertical element. Here
r2 = (x x0 )2 + (y y0 )2 = x20 + (y y0 )2
So,
Iijk (x0 )
yk+1
=
yk

x
i x
j
1
2
ij log r + 2
dy.
2
r
x=0
For example,
k
I11
(x0 )
yk+1
=
yk
x20
1
log[x20 + (y y0 )2 ] + 2
2
x0 + (y y0 )2
and
k
I12
(x0 )
yk+1
=
yk

dy,
x0 (y y0 )
dy.
x20 + (y y0 )2
These integrals can be done exactly. Most often, however, we end up with integrals which
we have to evaluate numerically.
Remaining elements
These may be written down in a similar way.
Boundary integral method
We apply the discretised BIE (7.3) at the mid-points of the elements. That is we take x0
to lie at the mid-point of element k for k = 1, . . . , N . In this way we produce 2N algebraic
equations for the 2N unknown traction values
(f11 , f21 ),
(f12 , f22 ),
(f1N , f2N ).
The left hand side of (7.3) is

N
1 X k
fi (x) Iijk (x0 ).
4
k=1
We express this in the matrix form
I11
1
I2
1
4
..
.
I1N
I21
I22
..
.
I2N
IN
1
IN
2
..
..
.
.
IN
N
49
f1
f2
..
fN
MTH-ME82
Here
Ikm (x0 ) =
k (x ) I k (x )
I11
0
12 0
when x0 lies on element m,
k (x )
I21
0
k (x )
I22
0
and
fk =
f1k
f2k
contains the two components of the constant traction on element k. So the complete linear
system to be solved is
0
1
I11 I21 IN
f
D
1
1
.
.
1
2
2
N
I2 I2 I2 f D2
1

1
1
(7.4)
0 ,
. = . 2 U
4
.
.
.
.
.
. .
.
.
.

.
.
.
. . .
..
N
I1N I2N IN
f
D
.
N
N
0
where
D1 (x0 )
Dm =
when x0 is on element m.
D2 (x0 )
We solve this linear system using Gaussian elimination.
Once we have computed the wall tractions, we can determine the fluid velocity at any point
in the flow using the BIE (6.12),
Z
Z
1
1
Gij (x, x0 )fi (x) dl(x) +
uj (x0 ) =
4 C
4 C
In the present case, using the identity (5.10) this reduces to
Z
Z
1
U
uj (x0 ) =
Gij (x, x0 )fi (x) dl(x)
T1j2 (x, x0 ) dl(x).
4 C
4 L
The Matlab code (see my webpage) cavity.m solves the driven cavity problem as described
above. (NB: The code uses a method of numerical integration called Gauss-Legendre
integration. Do not worry about this part.) Figure 4(a) shows the discretisation of the
cavity with N = 16 elements. Figure 4(b) shows the computed tractions using N = 64
elements. Note that the traction at a 90 corner in Stokes flow diverges, and this is why
we see the sharp spikes at the corners. This is therefore a property of the flowfield and not
50
MTH-ME82
(a)
(b)
200
fx
fy
1.5
150
100
50
0.5
0
50
0
100
150
0.5
0.5
0.5
0.5
1.5
1.5
2
s
2.5
3.5
Figure 4: (a) Cavity with N = 16 boundary elements (mid-points shown with red circles).
(b) The computed tractions with N = 64 boundary elements.
(a)
(b)
1.2
ux
uy
1.5
0.8
1
0.6
0.4
0.5
0.2
0.2
0.4
0.6
0.1
0.2
0.3
0.4
0.5
y
0.6
0.7
0.8
0.9
0.5
0.5
0.5
1.5
Figure 5: (a) Velocity profiles along x = 0.75L using N = 64 elements. (b) Computed
streamlines using N = 64 elements.
a failure of the method. In figure 5 velocity profiles of the horizontal and vertical velocity
components, ux and uy respectively, are plotted along the line x = 0.75L, and typical
streamlines are shown, both using N = 64 boundary elements.
8. Applications
The boundary integral method has two very powerful features
It reaches a solution using only the boundary data
51
MTH-ME82
It can handle complex flow geometries

In this section, we discuss some practical applications.
8.1 Shear flow over an obstacle
Simple shear flow over a wall at y = 0 is described by the solution
u = y.
The velocity profile is shown in the figure below.
Suppose we are interested in computing three-dimensional flow over a wall with an obstacle.
Practical examples of this type of flow include coating flows in manufacturing microelectronic circuitry or blood flow over a thrombus.
The figure illustrates the flow geometry. The wall lies at y = 0, and an obstacle is fixed on
the wall as shown.
y
P
Figure 6: Shear flow over an object attached to a plane wall

Let P represent the volume occupied by the obstacle.
Let C represent the wall area covered by the obstacle.
Let A represent the wall area left uncovered by the obstacle.
Aim: To derive a boundary integral equation for the velocity field at a general point x0 .
Method: We start by splitting the desired velocity field into two parts, writing
u = u + uD
where
u represents the background flow obtaining in the absence of the obstacle, i.e. the
undisturbed shear flow. So
u = y i.
uD represents the disturbance due to the obstacle. A long way from the obstacle this
disturbance flow decays to zero. If there is no obstacle, uD 0 everywhere.
52
MTH-ME82
Together, the background and disturbance flows combine to make up the total flow u.
Similarly the total traction is split up as follows:
f = f + f D,
where f and f D are the background and disturbance parts of the traction respectively.
From equation (6.7), we have the boundary integral equation for Stokes flow
ZZ
ZZ
1
1
uj (x0 ) =
8 S
8
S
(8.1)
where x0 is a point within the flow field and not on the boundary.
We apply the BIE (8.1) to the disturbance flow, uD :
ZZ
ZZ
1
1
D
uD
(x
)
=
f
(x)
G
(x,
x
)
dS(x)
+
uD
0
ij
0
j
i (x)Tijk (x, x0 )nk dS(x), (8.2)
8 S i
8
S
where S includes the whole of A, P and a surface extending from the wall to infinity and
back (indicated by the dotted line in the figure). Since the Greens function Gij and stress
tensor Tijk decay to zero at infinity the integrals over the surface extending to infinity equal
zero.
Now we can make some simplifications. The total flow satisfies no-slip, u = 0, on the wall A
and on the obstacle P . Therefore uD + u = 0 or uD = u on A and on P . Substituting
into (8.2),
ZZ
ZZ
1
1
D
D
uj (x0 ) =
f (x) Gij (x, x0 ) dS(x)
u
i (x)Tijk (x, x0 )nk dS(x).
8 A,P i
8
A,P
But, u
i = 0 on A so
ZZ
ZZ
1
1
D
D
uj (x0 ) =
u
i (x)Tijk (x, x0 )nk dS(x). (8.3)
8 A,P i
8
P
Note that the second integral is now taken only over the obstacle P . Equation (8.3) is valid
for a point x0 lying inside the flow.
For the next step, we recall the integral equation we derived for a point x0 lying outside
the integration domain, equation (6.3):
ZZ
ZZ
1
1

(8.4)
8 S
8
S
If we choose S to be the union of P and C and x0 to lie in the flow, then x0 is certainly
outside the domain of integration. Applying (8.4) for the background flow,
ZZ
ZZ
1
1

u
(8.5)
i (x)Tijk (x, x0 )nk dS(x) = 0,

8 P,C
8
P,C
53
MTH-ME82
n . But u = 0 on the wall C so

where fi = ik
k
i
ZZ
ZZ
1
1
u
i (x)Tijk (x, x0 )nk dS(x) = 0.

8 P,C
8
P
Rearranging,
ZZ
P
u
i (x)Tijk (x, x0 )nk dS(x) =
ZZ
Gij (x, x0 )fi (x) dS(x).
(8.6)
(8.7)
P,C
which is valid for a point x0 lying within the flow.

Since both (8.3) and (8.7) are valid for a point x0 within the flow, we can combine them.
Substituting (8.7) into (8.3) we get
ZZ
ZZ
1
1
D
D
uj (x0 ) =
u
i (x)Tijk (x, x0 )nk (x) dS(x)
8 A,P i
8
P
ZZ
ZZ
1
1
D
=
f Gij dS(x)
f Gij dS(x)
8 A,P i
8 P,C i
ZZ
ZZ
ZZ
1
1
1
D
=
f Gij dS(x)
fi Gij dS(x)
f Gij dS(x).
8 A i
8 P
8 C i
Note that in the last step the disturbance and background tractions in the integral over P
have combined to give the total traction.
Adding u
i to both sides, we obtain
ZZ
ZZ
1
1
D
uj (x0 ) = u
(x
)
f
(x)
G
(x,
x
)
dS(x)

0
ij
0
j
8 A i
8 P
ZZ
1
f (x)Gij (x, x0 ) dS(x).(8.8)

8 C i
This is a boundary integral equation for the velocity field. If we choose a Greens function
such that Gij (x, x0 ) = 0 when x is on the wall10 , it simplifies to
ZZ
1
uj (x0 ) = uj (x0 )
fi (x) Gij (x, x0 ) dS(x).
(8.9)
8 P
Thus it can sometimes be a good idea to use a particular Greens function in order to
simplify the formulation.
Notes:
1. Equation (8.9) is a boundary integral equation for the velocity involving the unknown
traction fi (x) and velocity ui (x). To obtain the solution, we must solve for these simultaneously.
10
The form of the required Greens function was given in section 5.2
54
MTH-ME82
8.2 Flow past a rotating object

Consider the slow flow of a viscous fluid around a three-dimensional rigid object rotating
about an axis passing through O perpendicular to the page as shown in the figure. Let x
be the position vector (relative to O) of a point P which is in or on the object.
n
x
O
solid
fluid
If the rigid body is rotating with angular velocity , where the vector points along the
axis of rotation (out of the page), the velocity of the point P is11
u(x) = x.
Making use of the alternating tensor defined in section 1.1, we can write this in index
notation as
ui (x) = ilm l xm .
(8.10)
To find the velocity of the fluid, we use the boundary integral formulation. The BIE equation
(6.7) states that
ZZ
ZZ
1
1
ui (x)Tijk (x, x0 )nk dS(x) (8.11)
uj (x0 ) =
8 S
8
S
when x0 lies in the fluid. Using (8.10), this becomes
ZZ
ZZ
l
1
uj (x0 ) =
fi (x) Gij (x, x0 ) dS(x) + ilm
xm Tijk (x, x0 )nk dS(x).
8 S
8
S
An identity states that
8
xm Tijk (x, x0 )nk (x) dS(x) = jlm x0,m 4
S
0
ZZ
ilm
(8.12)
when x0 is inside (8), right on (4) or outside S (0), and where the normal vector points
outside S as shown in the figure. Note that x0,m means the mth component of x0 . Hence,
since x0 is in the fluid outside S, the double layer potential vanishes leaving
ZZ
1
uj (x0 ) =
fi (x) Gij (x, x0 ) dS(x),
(8.13)
8 S
11
See any book on classical mechanics
55
MTH-ME82
for a point x0 lying inside the fluid.

In fact, (8.13) also applies when x0 lies on the surface S. To see this, we can either appeal
to section 6.3.1 and note that the single layer potential is a continuous function as x0
approaches S and hence (8.13) is equally valid on the surface, or we may simply apply the
BIE for a point x0 lying on S, namely equation (6.11),
1
1
uj (x0 ) =
2
8
ZZ
S
1
8
PV
ZZ

S
Recall that the P V on the second integral means that the integral is evaluated with x0 on
S. Using (8.10) and then identity (8.12) we find
ZZ
1
1
1
fi (x) Gij (x, x0 ) dS(x) jlm l x0,m .
uj (x0 ) =
2
8 S
2
But when x0 lies on S, uj (x0 ) = jlm l x0,m , and so we have
ZZ
1
uj (x0 ) =
8 S
again.
Note: Problem S14 on the Problem Sheet asks you to perform a similar analysis to the
above for an object moving at constant speed through a viscous liquid.
8.3 Flow past a liquid drop attached to a wall
Consider the flow of a viscous liquid past a drop of another liquid of the same viscosity
attached to a plane wall. We assume that the drop does not move along the wall, but can
be deformed by the flow. We also expect the passing flow to generate a fluid motion inside
the drop.
Label the outer fluid as fluid 1, and the fluid in the drop as fluid 2. Let n be the unit
normal vector pointing into the outer flow.
y
P
Let P represent the surface of the drop exposed to the fluid.

Let C represent the wall area covered by the drop.
Let A represent the wall area left uncovered by the drop.
56
MTH-ME82
We wish to derive a boundary integral equation for the velocity field at a general point x0
both in the outer flow and inside the drop.
To begin, we label the velocity and traction fields in fluid j = 1, 2 as
u(j) (x0 ),
f (j) (x0 ),
at a point x0 .
We split the velocity and traction fields in fluid 1 into a background and a disturbance part,
u(1) = u + u(1)D ,
f (1) = f + f (1)D .
Here,
u represents the background flow found in the absence of the drop, given by
u = y i.
u(1)D represents the disturbance outer flow in fluid 1 due to the drop. A long way
from the drop this disturbance flow decays to zero. If there is no drop, u(1)D 0
everywhere.
Now, apply the BIE (8.1) for a point inside fluid 1 to the disturbance flow, u(1)D :
ZZ
ZZ
1
1
(1)D
(1)D
(1)D
fi
(x) Gij dS(x) +
ui (x)Tijk nk dS(x),
uj (x0 ) =
8 S
8
S
(8.14)
where S includes the whole of A, P and a surface extending from the wall to infinity and
back (indicated by the dotted line in the figure). Since the Greens function Gij and stress
tensor Tijk decay to zero at infinity the integrals over the surface extending to infinity equal
zero.
The total flow in fluid 1 satisfies no-slip, u(1) = 0, on the wall A. Therefore u(1)D + u = 0
or u(1)D = u on A. Substituting into (8.14),
ZZ
ZZ
1
1
(1)D
(1)D
f
Gij dS(x)
u
uj (x0 ) =
i Tijk nk dS(x)
8 A,P i
8
A
ZZ
1
(1)D
+
ui Tijk nk dS(x).
8
P
But u
i = 0 on A so
(1)D
uj (x0 )
1
=
8
ZZ
A,P
(1)D
fi
1
Gij dS(x) +
8
ZZ
P
(1)D
ui
Tijk nk dS(x).
Adding u
i to both sides, we have
ZZ
ZZ
1
1
(1)
(1)D
(1)D
uj (x0 ) = uj (x0 )
f
Gij dS(x) +
ui Tijk nk dS(x)
8 A,P i
8
P
57
(8.15)
MTH-ME82
for a point x0 inside fluid 1.

For the next step, recall the integral equation we derived for a point x0 lying outside the
domain of integration, equation (6.3):
ZZ
ZZ
1
1
(8.16)
8 S
8
S
First, we take S to be the union of P and C and apply this equation for the background
flow u
i for a point x0 lying inside fluid 1 and therefore outside the domain of integration.
We obtain,
ZZ
ZZ
1
1
0=
fi Gij dS(x) +
u
i Tijk nk dS(x).
8 P,C
8
P,C
Since u
i = 0 on the wall C, this simplifies to
ZZ
ZZ
1
1
0=
f Gij dS(x) +
u
i Tijk nk dS(x).
8 P,C i
8
P
(8.17)
which is valid for a point x0 lying inside fluid 1.

Since both (8.15) and (8.17) apply for a point x0 inside fluid 1 we can combine them.
Adding (8.17) to (8.15) we find
ZZ
ZZ
1
1
(1)
(1)D
(1)D
uj (x0 ) = uj (x0 )
fi
Gij dS(x) +
ui Tijk nk dS(x)
8 A,P
8
P
ZZ
ZZ
1
1
fi Gij dS(x) +
u
i Tijk nk dS(x)
8 P,C
8
P
ZZ
ZZ
ZZ
1
1
1
(1)D
(1)
= uj (x0 )
fi
Gij dS(x)
fi Gij dS(x)
fi Gij dS(x)
8 A
8 C
8 P
ZZ
1
(1)
+
ui Tijk nk dS(x).
(8.18)
8
P
At this point, we decide to simplify matters by choosing a Greens function which vanishes
on the wall, so that Gij (x, x0 ) = 0 when x is on the wall. Then (8.18) reduces to
ZZ
ZZ
1
1
(1)
(1)
(1)
uj (x0 ) = u
(x
)
f
G
dS(x)
+
ui Tijk nk dS(x).
(8.19)
ij
0
j
8 P i
8
P
This equation is valid for a point x0 inside fluid 1.
Next we apply equation (8.16) with S as the union of P and C but this time for the flow
in the drop, fluid 2, and for a point x0 lying in fluid 1. We obtain,
ZZ
ZZ
1
1
(2)
(2)
0=
f Gij +
ui Tijk nk dS(x).
8 P,C i
8
P,C
58
MTH-ME82
(2)
But ui = 0 on the wall C by the no-slip condition. Also Gij has now been chosen to vanish
on the wall. So,
ZZ
ZZ
1
1
(2)
(2)
0=
fi Gij dS(x) +
ui Tijk nk dS(x)
(8.20)
8 P
8
P
for a point x0 lying inside fluid 1.
Since both (8.19) and (8.20) apply for a point x0 inside fluid 1 we can combine them.
Subtracting (8.20) from (8.19), we find
ZZ
ZZ
1
1
(1)
(1)
(1)
uj (x0 ) = uj (x0 )
f Gij dS(x) +
ui Tijk nk dS(x)
8 P i
8
P
ZZ
ZZ
1
1
(2)
(2)
+
f Gij dS(x)
ui Tijk nk dS(x)
(8.21)
8 P i
8
P
ZZ
ZZ h
i
1
1
(1)
(2)
= uj (x0 )
fi Gij dS(x) +
ui ui Tijk nk dS(x),
8 P
8
P
where
(1)
fi = fi
(2)
fi
is the jump in traction at the interface between the drop and the outer liquid.
Physically, we expect the velocity field to be continuous at the interface, so that
(1)
(2)
ui (x) = ui (x)
when x lies on the interface. Hence (8.21) reduces to
ZZ
1
(1)
fi (x) Gij (x, x0 ) dS(x),
uj (x0 ) = u
(x
)
0
j
8 P
and is valid for a point x0 inside fluid 1. This formula should be compared with (8.9) in
section 8.1 for flow over a solid object attached to a wall.
8.4 Final comment
The choice of boundary integral formulation is not necessarily unique and different
approaches may be taken. As a general rule, the aim is to end up with a Fredholm integral
equation of the second kind. The reason for this is that solving integral equations of the
first kind numerically usually results in unstable numerical schemes (e.g. Power & Miranda
SIAM J Applied Math V47(4), p.691). For a comparison of different boundary integral
formulations (single-layer, double-layer etc.) see Ingber & Mammoli (1999) A comparison
of integral formulations for the analysis of low Reynolds number flows, Eng. Analysis Bound.
Elem., 23, 307-315.
Appendix A: Alternative form of the delta function

59
MTH-ME82
In this appendix we note the alternative representation of the delta function:

(
x) = 2
1
.
4r
(8.22)
Pozrikidis (1992, page 22) uses this form as a starting point for deriving the three-dimensional
stokeslet. To see why (8.22) is true, note that the free-space Greens function G for Laplaces
equation in three-dimensions satsifies
2 G = (
x),
(8.23)
= x x0 . Away from x = x0 the delta function is zero, leaving

where x
2 G = 0.
Define spherical polar coordinates (r, , ) centred on x0 , so that
r = |x x0 |,
i.e. r is the radial distance from the field point to the pole. Then, working in spherical
polars and noting that we expect the solution to be radially symmetric we have

1
2 G
2
r
= 0,
G= 2
r r
r
away from x = x0 . Integrating, we see that one solution is
A
G= .
r
for constant A. To fix A we integrate (8.23) over a small sphere, volume V , radius , and
surface S, which is centred at x0 . So,
Z
Z
2
G dV =
(x x0 ) dV = 1
V
on using the sifting property of the delta function. Applying the divergence theorem12 in
fact, we find
Z
Z
G
(x x0 ) dV = 1
dS =
V
S r
Thus, since the surface area of the sphere is 42 , we find
4A = 1
and so A = 1/4.
12
Note that it is a little dicey using the divergence theorem over a domain containing a singularity! This
calculation should more properly be performed by first isolating the singularity inside a small ball. See
Stakgold Greens Functions and Boundary Value Problems, Chapter 4
60
MTH-ME82
Thus,
G=
1
.
4r
(8.24)
is a solution of 2 G = (
x). Accordingly, we can write

1
2
.
(
x) =
4r
Note: we can alternatively derive (8.24) using a Fourier transform method.
Appendix B: Properties of the delta function

We demonstrate the properties of the delta function given in section 5.3, namely
Z
if x0 lies inside V
1
1
(x x0 ) dV =
if
x
0 lies exactly on S
2
V
0
We will work in two dimensions, in which case we have
Z
if x0 lies inside V
1
1
(x x0 ) dS =
if x0 lies exactly on S
2
S
0
where S is a two-dimensional surface enclosed by a contour C with unit outward normal
n. The argument is readily extended to three dimensions. To proceed we re-write the
two-dimensional delta function in a manner analogous to that introduced13 in section 5.1.1,
=
Then, we integrate over S,
Z
Z
(x x0 ) dS =
S
1 2
log r.
2
1 2
1
log r dS =
2
2
Z
n log r dl,
C
by the divergence theorem.

Consider the case when the singular point, x0 , lies on the contour C. We indent C inwards
in the form of a circular arc, radius , so that x0 now lies outside S as shown. Now, since
x0 is outside S, we have
Z
Z
1
(x x0 ) dS = 0 =
n log r dl
2 C
S
13
1
Recall that in section 5.1.1, we wrote (
x) = 2 ( 4r
).
61
MTH-ME82
x0

Z PV
Z
1
1
n log r dl +
n log r dl ,
= lim
0 2 C
2
where the principal value P V means integration around C excluding the indented portion
. Integrating over the indented arc, , we find,
Z
Z
1
1
1
1
n log r dl, =
d = ,
2
2 0
2
since the outward normal n points in the inward radial direction normal to the arc. Thus,
Z PV
1
1
0= +
n log r dl
2 2 C
So,
Z PV
1
1
n log r dl = .
2 C
2
Hence, when x0 lies precisely on the contour C,
Z
Z PV
1
1
n log r dl = .
(x x0 ) dS =
2 C
2
S
The other results follow from the standard sifting property of the delta function.

Appendix C: Symmetry of the Greens function for Stokes flow

It was stated in section 4.4.1 that the Greens function is symmetric, that is
Gij (x1 , x2 ) = Gji (x2 , x1 ).
To show this, we apply a similar argument used to demonstrate the symmetry of the Greens
function for potential flow14 . We start with the Lorentz reciprocal relation (4.7) valid for two
0 ) which contain no singular points within the domain of interest,
flows (ui , ik ) and (u0i , ik
0
(u0 ik ui ik
) = 0.
(C.1)
xk i
Suppose that the domain in question is enclosed by a boundary S and that the Greens
function in question vanishes on the boundary, that is
Gij (x, x0 ) = 0
when x lies on S and where x0 is the location of the singularity. To apply the Lorentz
reciprocal identity, we consider the volume which consists of the inside of S but with two
small spheres of radius centred around two points x1 and x2 removed as shown. Call this
S S
volume V 0 and its bounding surfaces S 0 = S S1 S2 .
14
See, for example, Garabedian, Partial Differential Equations, p.244.
62
MTH-ME82
S1
S2
x1
n
x2
S
Note that the unit normals are taken to point out of V 0 . We take
ui (x) =
1
Gim (x, x1 )am ,
8
ik =
1
Tijk (x, x1 )aj ,
8
(C.2)
u0i (x) =
1
Gim (x, x2 )bm ,
8
0
ik
=
1
Tijk (x, x2 )bj ,
8
(C.3)
and
for arbitrary constant vectors a and b. Integrating (C.1) over the volume V 0 within S 0 ,
Z
0
)dV = 0,
(u0i ik ui ik
V 0 xk
or, using the divergence theorem15 ,
Z
0
nk )dS = 0.
(u0i ik nk ui ik
S0
Accordingly,
Z
0
nk )dS +
(u0i ik nk ui ik
0
(u0i ik nk ui ik
nk )dS = 0.
S1 ,S2
The first integral on the left hand side vanishes since both ui and u0i vanish on S. So,
Z
0
(u0i ik nk ui ik
nk ) dS = 0.
S1 ,S2
Consider now the integral
Z
I11 =
u0i ik nk dS.
S1
Using (C.2) and (C.3) this becomes

1
aj bm
64 2
15
Z
Gim (x, x2 ) Tijk (x, x1 )nk dS
S1
Note that this is permitted since the integrand is regular throughout S 0
63
MTH-ME82
Letting the radius of S1 tend to zero, this becomes

Z
x
i x
j x
k
1
nk dS.
a
b
G
(x
,
x
)
6
j
m
im
1
2
64 2
5
S1
Note that we have approximated Tijk by its Stokeslet form (compare section 6). From the
working in section 6, we have the result that
Z
x
i x
j x
k
6
Gim (x1 , x2 )
nk dS = 8 Gjm (x1 , x2 ).
5
S1
So,
I11 =
1
aj bm Gjm (x1 , x2 ).
8
Next consider
Z
0
nk dS
ui ik
I12 =
S1
1
=
am bj
16 2
Z
Gim (x, x1 )Tijk (x, x2 )nk dS.
S1
Letting 0,
1
am bj Tijk (x1 , x2 )
16 2
I12
Z
S1
im
x
i x
m
nk dS.
3
Note that we have approximated Gim by its Stokeslet form (compare section 6). From the
working in section 6, we have the result that
Z
im x
i x
m
+ 3 nk dS 0

S1
as 0. So I12 0.
Following similar reasoning we may determine the values of
Z
Z
0
0
nk dS.
I21 =
ui ik nk dS,
I22 =
ui ik
S2
S2
In fact, we find
I21 = 0,
I22 =
1
1
am bj Gjm (x2 , x1 ) =
aj bm Gmj (x2 , x1 ).
8
8
Note that in the last step we switched the roles of m and j. This is permissible since they
are both dummy indices. Finally, we find that
Z
0
(u0i ik nk ui ik
nk )dS = I11 + I12 I21 I22
S1 ,S2
1
1
aj bm Gjm (x1 , x2 )
aj bm Gmj (x2 , x1 ) = 0.
8
8
It therefore follows that

Gij (x1 , x2 ) = Gmj (x2 , x1 ),
64
MTH-ME82
which is the required symmetry property.
Appendix D: Conservation of mass for a point force in Stokes flow

In section 4.4.1, we argued using the divergence theorem that
Z
Gij (x, x0 )ni (x) dS = 0.
(D.1)
We applied the divergence theorem for an integral containing a singularity at the pole x0 .
Here we provide a more formal, classical argument.
The velocity field in response to a point force at x0 is
ui = Gij (x, x0 ) bj ,
(D.2)
where the Greens function satisfies the singularly forced Stokes equation (4.8).
We take the divergence of (D.2) and integrate over a volume V containing the singular point
x0 but excluding a small sphere, V , of radius centred at x0 . Call this punctured volume
V 0 = V V . Let S and S represent the surfaces of V and V respectively. We find that
Z
Z
ui
dV =
ui ni dS,
V 0 xi
S,S
on applying the divergence theorem (note that the integrand is regular over the punctured
region). Now,
Z
Z
Z
Z
x
j
ij
ij
x
i x
i x
j
+
dS =
+
dS
ui ni dS = bj
Gij ni dS =
4
2

S
S
S
S
on noting that ni = x
j / on S . Noting further that dS 2 and x
j 2 and taking the
limit as 0, we find that this integral vanishes and we are left with
Z
Z
ui
dV =
ui ni dS.
V 0 xi
S
Since u = 0 over V 0 we deduce that
Z
ui ni dS = 0
S
and so (4.12) holds.
Appendix E: Derivation of the stokeslet in 3D

The following is a rather long-winded argument showing how to derive the three-dimensional
stokeslet (or Oseen-Burgers tensor).
65
MTH-ME82
To proceed, we re-express the delta function in the more convenient form, (see Appendix
A)
1
(
x) = 2
.
4r
We take the singularly forced Stokes equation, namely,
0 = P + b (x x0 ) + 2 u,
and write it in the equivalent form
b 2 1
= P + 2 u,
4
r
(D.1)
How do we deal with the pressure term? The Stokes pressure is a harmonic function, i.e.
it satisfies Laplaces equation,
2 P = 0
(this is left as an exercise). If we write

1

1
1
1
P = b
, ,
= b
4
x y z
r
4
r
we can confirm that P is a harmonic function (see Problem Sheet 1).
For the next step, we introduce a scalar function H, and write the velocity field in the form

1

2
uj = bi
ij H,
x
i x
j
where ij is the Kronecker delta. Using a more compact vector notation, we write this as
u=
1
b ( I2 )H,
where I is the identity matrix,
1 0 0
0 1 0 ,
0 0 1
and

x
i x
j
is an example of a dyadic product. For any two vectors a and b, the dyadic product
ab = ai bj .
Note that ab is a matrix with ijth entry ai bj .
Summary: The singularly forced Stokes momentum equation to be satisfied is
b 2 1
= P + 2 u,
4
r
66
(D.2)
MTH-ME82
and thus far we have written

1
P = b
4
and
u=
Now, note that

1
,
r
(D.3)
1
b ( I2 )H.
2
2u
1
2 H
=
b
I
)
.
x
2
x
2
Consequently,
2 u =
1
b ( I2 )2 H.
(D.4)
Substituting (D.3) and (D.4) into (D.2), we have

1
1
1
2 1
b I
=
b
+ b ( I2 )2 H.
4
r
4
r
Rearranging,
1
0=
b
4

1
1
1
b I2
+ b ( I2 )2 H,
r
4
r
or
0 = b I2

2 H +
1
.
4r
We want this to be true for any choice of b, so

1
0 = I2 2 H +
.
4r
(D.5)
We can satisfy then (D.5) if we choose H such that

2 H =
1
.
4r
(D.6)
To find the solution, recall from above that we may express the delta function in the form
1
(
x) = 2
.
4r
Taking 2 of both sides of (D.6),
4 H = 2
1
= (
x).
4r
So now we know that H satisifies the forced biharmonic equation

4 H = (
x).
67
MTH-ME82
Recalling that in spherical polars (independent of and ),

f
1
r2
= 0,
2 f = 2
r r
r
we see that, away from x = x0 ,
1
H= 2
r r
4

1
2
2 H
r
r
= 0.
r
r2 r
r
Integrating,

1
2
2 H
r
r
= A,
r
r2 r
r
for constant A. Integrating again,

A
2 H
r
= + B,
r
r
1
r2 r
for constant B. Integrating once more,

r2
H
A
B
= r2 + r3 + C,
r
2
3
for constant C. Integrating for a fourth time,

A
B
C
H = r + r2 + D,
2
6
r
for constant D. As usual, we integrate the singularly forced equation over a sphere of radius
a centred at the singularity, obtaining
Z
Z
4
H dV =
(
x) dV = 1,
V
by the sifting property of the delta function. Now,

4 H = (2 H),
so
(2 H) dV = 1,
and on applying the divergence theorem,

Z
(2 H) n dS = 1.
S
But for the sphere, n = r, and so

Z
S
2
H dS = 1.
r
68
MTH-ME82
Now,
A
B
C
H = r + r2 + D,
2
6
r
as found above, and thus

2 H =
But
2
2 r = ,
r
1
r2 r

H
r2
.
r
2 r2 = 6,
and so
2 H =
1
r
= 0,
A
+ B.
r
Accordingly,
Z
S
2
H dS =
r
Z
S
A
dS =
r2
Z
S
A
A
dS = 2 4a2 = 4A = 1,
2
a
a
implying that A = 1/4, and the values of B, C and D are immaterial (since we have just
seen that we only need to fix A to be 1/4 to get a solution to the equation). Therefore,
we can set B = C = D = 0, leaving
r
H= .
8
Originally, we wrote
u=
So now we know that
1
b ( I2 )H.

r
1
b ( I2 )
.
8

r
1
r 1
= b ()
+ b (I2 )
8
u=
Since 2 r = 2/r, this gives

u=

1
1
r 1
b ()
+ bI
,
4r
i.e.,

2
1
u=
b (r) + I
,
8
r
In index notation,

2
1

uj =
bi
(r) + ij
,
8
x
i x
j
r
Noting that
we have
x
j
r
= ,
x
j
r

2
j
1
x
uj =
bi
+ ij
,
8
x
i r
r
69
MTH-ME82

2
rij x
j (
xi /r)
1
+ ij
,
bi
=
8
r2
r
using the quotient rule of differentiation and using the fact that
x
j
= ij .
x
i
Simplifying,

ij
x
i x
j
2ij
1
,
uj =
bi
+ 3 +
8
r
r
r
leaving,

x
i x
j
ij
1
.
uj =
bi
+
8
r3
r
We are now in a position to write the velocity field in the form of (4.9), i.e.
ui =
1
Gij (x, x0 )bj ,
8
where the Stokeslet Greens function

Gij =
x
i x
j
ij
+
.
r3
r
Note that the Stokeslet is sometimes called the Oseen-Burgers tensor.

We write the pressure associated with the Stokeslet as
P =
1
pj bj ,
8
(see equation 4.10). Above, we wrote

x
1
1

1
2
x
1
P = b
= b 3 =
b
4
r
4
r
8
r3

2
xj
1
=
bj
8
r3
in index notation. So the Stokeslet pressure is
pj = 2
x
j
.
r3
Appendix F: The discontinuity of the double-layer potential for Stokes flow

We wish to demonstrate that the double-layer potential has the following behaviour
ZZ PV
DL
lim Ij (x0 ) = 4uj (x0 ) +
(F.1)
x0 D+
70
MTH-ME82
and
lim IjDL (x0 )
x0 D
ZZ
PV
= 4uj (x0 ) +
(F.2)
We denote by D+ the side of D which faces in the fluid into which the normal vector points
(see the diagram below). We denote by D the other side of D, which is exterior to the
flow. Recall that in (6.8) the normal vector points into the fluid consistent with the figure.
dm
dp
dp
F
dm
To proceed, we re-write (6.8) in the equivalent form

ZZ
ZZ h
i
DL
ui (x) ui (x0 ) Tijk (x, x0 )nk dS(x) + ui (x0 )
Tijk (x, x0 )nk dS(x).
Ij (x0 ) =
D
Next we decompose I DL thus,

IjDL (x0 ) = Yj (x0 ) + Qj (x0 ),
where
Yj (x0 ) =
ZZ h
i
ui (x) ui (x0 ) Tijk (x, x0 )nk dS(x),
and
ZZ
Qj (x0 ) = ui (x0 )
Tijk (x, x0 )nk dS(x).

D
The reason for doing this is that the integral

ZZ h
i
Yj (x0 ) =
ui (x) ui (x0 ) Tijk (x, x0 )nk dS(x)
D
is continuous as x0 crosses D, and we can focus attention on Qj .
71
(F.3)
MTH-ME82
From section 5.2.1, we have the identity:
8
Tijk (x, x0 )ni (x) dS(x) = 4 jk ,
D
0
ZZ
where S is a closed surface, and the unit normal, n, points into16 S as shown in the figure
above. Note that
8 applies when x0 is located inside D.
4 applies when x0 is located on D.
0 applies when x0 is located outside D.
Recalling that Tijk = Tkji , the relation can also be written
8
Tijk (x, x0 )nk (x) dS(x) = 4 ji ,
D
0
ZZ
So, when x0 is inside D,

Qj (x0 ) = 8ij ui (x0 ) = 8uj (x0 ) = Q+
j ,
say.
When x0 is outside S
Qj (x0 ) = 0 = Q
j ,
say.
When x0 is precisely on S,
Qj (x0 ) = 4uj (x0 ) = QSj ,
say.
In summary,
Q+
j = 8uj (x0 ),
QSj = 4uj (x0 ),
Q
j = 0.
(F.4)
From (F.3),
IjDL (x0 ) = Yj (x0 ) + Qj (x0 ).
So as x0 D+ from inside the flow,
IjDL (x0 ) Yj+ + Q+
j ,
16
Note that for the identity quoted in section 5.3 the normal points out of S. Thus there is a minus sign
difference on the right hand side between the identities.
72
MTH-ME82
where Yj+ is the limiting value of Yj as x0 D+ . Also, as x0 S from outside of the

flow,
IjDL (x0 ) Yj + Q
j ,
where Yj is the limiting value of Yj as x0 D .
But since Yj is continuous as x0 crosses D, we have Yj+ = Yj = YjS , say, and so
lim IjDL (x0 ) = YjS + Q+
j ,
(F.5)
lim IjDL (x0 ) = YjS + Q

j .
(F.6)
x0 D+
x0 D
Define the principal value of the double-layer potential, IjDLP V as its value when x0 lies
on S, namely
Z PV
IjDLP V =
S
The superscript P V on the integral reminds us that x0 lies precisely on S.

Taking x0 to lie on S in (F.3) we have
IjDLP V (x0 ) = YjS + QSj
and so
YjS = IjDLP V (x0 ) QSj = IjDLP V (x0 ) 4uj (x0 )
(F.7)
using (F.4).
Therefore, putting (F.4), (F.5), (F.6), and (F.7) together, we have
lim IjDL (x0 )
x0 D+
YjS
Q+
j
PV
Z
= 4uj (x0 ) +
(F.8)
and
lim IjDL (x0 ) = YjS + Q
j = 4uj (x0 ) +
x0 D
PV

S
So across S, the double layer potential I DLP experiences a jump of size
Q+
j Qj = 8uj (x0 ).
73
(F.9)

Boundary Element and Finite Element Methods

Uploaded by

Copyright:

Available Formats

Boundary Element and Finite Element Methods

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Boundary Element and Finite Element Methods

Uploaded by

Copyright:

Available Formats

What are the main topics covered in the course?

What are the main topics covered in the course?

How is index notation used to represent vectors and tensors?

How is index notation used to represent vectors and tensors?

MTH-ME82 Boundary Element and Finite Element Methods