Download Boundary Element Method - Lecture Notes | MSC 321 and more Exams Earth Sciences in PDF only on Docsity!
Chapter 16
Boundary Element Method
We will proceed to learn about Object Oriented Programming and Design by look- ing at a specific example that revolves around the solution of the Laplace equation by the Boundary Element Method (BEM). The benefits of this approach is that it takes a 2D or 3D PDE and reduces it to an algorithm that involves 1D or 2D computations only, respectively. For simplicity we will focus on the 2D case only.
16.1 The equations and the boundary conditions
The Laplace equation takes the form
∇^2 φ = ∇ · (∇φ) = 0 (16.1)
subject to boundary conditions of the type
- Dirichlet where the unknown function is imposed on a portion of the bound- ary, i.e. φ(x) = a(s), s being a coordinate system along the boundary.
- Neumann where the normal derivative of the function is known ∇φ·n = q(s) where n is the outward unit normal.
- Robin where a relationship of the form αφ + β∇φ · n = r(s). A well-posed problem requires the specification of only one type of boundary conditions on any portion of the boundary, that is either the function value is specified, or the normal derivative but not both. The data imposed on the boundary form are part of the problem’s specifi- cation and can be quite arbitrary except for some minor constraints. For example, if Neumann conditions are applied over the entire boundary, then the solution can be determined only up to an additive constant, that is if φ is a solution of ∇^2 φ = 0 with ∇φ · n = q, then so is φ + C where C is an
155
156 CHAPTER 16. BOUNDARY ELEMENT METHOD
arbitrary constant. In that case the flux have to satisfy a compatability con- dition that can be derived by applying the divergence theorem to equation 16.1: ∫
Ω
∇^2 φ dA =
∫
∂Ω
∇φ · n ds =
∫
∂Ω
q ds = 0 (16.2)
16.2 From PDE to Integral Equation
Consider a function G, solution to a Poisson equation of the form:
∇^2 G = Q(x − xi) (16.3)
where Q(x − xi) is a function that depends on the space variable x and on the source point xi where an impulse forcing is applied. I will keep these terms vague for the time being, and will delay specifying them until later. Multiplying equation 16.3 by φ and 16.1 by G and substracting the resulting expressions we get:
φ∇^2 G − G∇^2 φ = φG(x − xi) (16.4)
Integrating the above equation over the area of the domain Ω we get: ∫
Ω
(φ∇^2 G − G∇^2 φ) dA =
∫
Ω
φG(x − xi) dA (16.5)
The chain rule of differentiation and the Gauss divergence theorem can be used to turn the left hand side into a boundary integral only. Indeed since for any differentiable functions a and b: ∇ · (a∇b) = ∇a · ∇b + a∇^2 b, we can then write:
φ∇^2 G − G∇^2 φ = ∇ · (φ∇G) − ∇φ · ∇G − ∇ · (G∇φ) + ∇φ · ∇G (16.6) = ∇ · (φ∇G) − ∇ · (G∇φ) = ∇ · (φ∇G − G∇φ) (16.7)
Furthermore the divergence form of the last equations permits the application of the Gauss divergence theorm to get: ∫
Ω
(φ∇^2 G − G∇^2 φ) dA =
∫
Ω
φQ(x − xi) dA (16.8) ∫
Ω
∇ · (φ∇G − G∇φ) dA =
∫
Ω
φQ(x − xi) dA (16.9) ∫ ∂Ω
n · (φ∇G − G∇φ) dS =
∫ Ω
φQ(x − xi) dA (16.10) ∫
∂Ω
( φ
∂G
∂n
− G
∂φ ∂n
) dS =
∫
Ω
φQ(x − xi) dA (16.11)
Further progress can be made to simplify equation 16.11 by imposing constraints on the choices of Q and G. Requiring that the function Q be the equivalent of the two-dimensional Dirac delta function, and G to be its associated free space
158 CHAPTER 16. BOUNDARY ELEMENT METHOD
One may wonder how to define a function which vanishes everywhere, spikes to infinity at a single point, and yet has unit area under the curve. A simple example is given by the function δa = max(0, (a − |x|)/a^2 ), and taking its limit as a → 0. In 2D the definition of the function differs a bit from its one-dimensional counterpart. In particular setting
Q(x, xi) =
δ(r) 2 πr
where r = |x − xi| is the distance of the point x from the source point xi. It is easy to verify that the area under Q is unity, indeed we have
∫ A
Q(x)dA =
∫ (^2) π 0
∫ (^) ∞ 0
δ(r) 2 πr
r dr dθ = 1 (16.18)
We are now ready to determine our free space Green function. Dropping the azimuthal dependency from the Laplacian in polar coordinates we have
∂^2 φ ∂r^2
r
∂G
∂r
δ(r) 2 πr
r
r
( r
∂G
∂r
)
δ(r) 2 πr
r
∂G
∂r
∫ (^) δ(r) 2 π
dr (16.21)
r
∂G
∂r
2 π
+ A (16.22)
∂G
∂r
2 πr
A
r
G =
ln |r| 2 π
where A and B are integration constants. The constant B is just an offset to the function, and we can arbitrarily set it to 0. The constant A can be determined by an application of the divergence theorem on a circle of radius R we then have
∫ A
Q dA =
∫ A
∇^2 GdA =
∫ ∂A
∂G
∂n
ds =
∫ (^2) π 0
∂G
∂r
Rdθ = 1 + 2πA (16.25)
Since the left hand side of the equation is equal to 1, the constant A = 0. A final note concerns the case where the source point is located on the boundary of the domain. The integral takes then a slightly different form. Figure 16.1 shows a sketch with the source point on the boundary. The area integral is now divided into a piece that excludes the point plus an infinitesimally small area near the source. The first integral is of couse zero while the second leads to
∫ (^) α
0
∫ (^) ǫ
0
φ(x)
δ(r) 2 πr
r dr dθ =
α 2 π
φ(xi) (16.26)
16.4. DISCRETIZATION 159
j
j+
r
i
s
n
Figure 16.1: BEM discretization of the geometry with straight line segments. The source point is located at xi while the field points cover segement sj ≤ s ≤ sj+1.
where α is the interior angle of the boundary. If the boundary is a smooth curve α = π. Our final form for equation 16.27 becomes ∫
∂Ω
( φ
∂G
∂n
− G
∂φ ∂n
) dS =
α 2 π
φ(xi), G =
ln |r| 2 π
or equivalently (^) ∫
∂Ω
( φ r
∂r ∂n
− ln |r|
∂φ ∂n
) dS = αφ(xi) (16.28)
where r is understood to be the disctance from the field point x to the source point xi.
16.4 Discretization
The discretization of the continuous problem proceeds by first dividing the bound- ary of the domain into line segments so that the geometry of the domain is well- represented. The choice of line type is rather arbitrary. In the following we will only consider the case where the line segments are straight lines as sketched in figure 16.1. Notice that for linear segments:
r^2 = s^2 + n^2 , and
∂r ∂n
n r
16.4. DISCRETIZATION 161
16.4.1 Evaluation of Integrals
All the integrals in 16.35-16.38 can be evaluated analytically; indeed it is easy to verify that:
A(sj+1, sj , n) =
n ∆s
∫ (^) sj+
sj
s n^2 + s^2
ds =
[ (^) n
2∆s
ln(s^2 + n^2 )
]sj+
sj
=
n 2∆s
ln
r j^2 + r j^2
B(sj+1, sj , n) =
n ∆s
∫ (^) sj+
sj
n^2 + s^2
ds =
∆s
[ tan−^1
s n
]sj+
sj
=
∆s
[ tan−^1
sj+ n
− tan−^1
sj n
] (16.40)
=
∆s
tan−^1
n∆s 1 + sj sj+
C(sj+1, sj , n) =
2∆s
∫ (^) sj+
sj
s ln(s^2 + n^2 ) ds
4∆s
[ (s^2 + n^2 ) ln(n^2 + s^2 ) − (s^2 + n^2 )
]sj+ sj =
4∆s
[( r^2 j+1 ln r^2 j+1 − r^2 j+
) −
( r j^2 ln r j^2 − r^2 j
)] (16.42)
D(sj+1, sj , n) =
2∆s
∫ (^) sj+
sj
ln(s^2 + n^2 ) ds
2∆s
[ s ln(n^2 + s^2 ) − 2 s + 2n tan−^1
s n
]sj+
sj
2∆s
{[ sj+1 ln r j^2 +1 − 2 sj+1 + 2n tan−^1
sj+ n
]
[ sj ln r^2 j − 2 sj + 2n tan−^1
sj n
]} (16.43)
=
2∆s
{[ sj+1 ln r j^2 +1 + 2n tan−^1
sj+ n
]
[ sj ln r^2 j + 2n tan−^1
sj n
]} − 1 (16.44)
2∆s
[ sj+1 ln r^2 j+1 − sj ln r j^2 + 2n tan−^1
n∆s 1 + sj sj+
] − 1
where we have used the fact that tan(x−y) = (^) 1+tantan^ x− xtan tan^ y y to simplify the expressions. Replacing these expressions into the definition of the matrix entries we have
Ki,j = sj+1B(sj+1, sj , n) − A(sj+1, sj , n) (16.45) Ki,j+1 = −sj B(sj+1, sj , n) + A(sj+1, sj , n) (16.46) Ri,j = sj+1D(sj+1, sj , n) − C(sj+1, sj , n) (16.47) Ri,j+1 = −sj D(sj+1, sj , n) + C(sj+1, sj , n) (16.48)
162 CHAPTER 16. BOUNDARY ELEMENT METHOD
Notice that special care must be taken when integrating along an edge that contains the source point, since in that case either rj or rj+1 is zero. To account for this case we notice first that when n = 0 we have
A(sj+1, sj , 0) = 0 (16.49) B(sj+1, sj , 0) = 0 (16.50)
C(sj+1, sj , 0) =
4∆s
[( s^2 j+1 ln s^2 j+1 − s^2 j+
) −
( s^2 j ln s^2 j − s^2 j
)] (16.51)
D(sj+1, sj , 0) =
2∆s
[ sj+1 ln s^2 j+1 − sj ln s^2 j
] − 1 (16.52)
In a second step it is now easy to take the limit when sj or sj+1 are zero, since x ln x → 0 when x → 0, and so we have:
C(0, sj , 0 , 0) =
sj (2 ln sj − 1) 4
C(sj+1, 0 , 0) =
sj+1(2 ln sj+1 − 1) 4
D(0, sj, 0) = ln sj − 1 (16.55) D(sj+1, 0 , 0) = ln sj+1 − 1 (16.56)
16.5 Algorithm
The development of a code for solving the Laplace equation on arbitrary domain using BEM requires several steps:
- Description of the geometry
- Description of the boundary conditions
- Building the matrices mostly by evaluating integrals
- Constructing of the linear system by assembly
- Solution of the linear system
- Evaluating the solution at interior points
- Plotting the solution and displaying it.
We will take up each of these tasks in the following subsections.
16.5.1 Geometry description
The geometry is essentially formed by a collection of nodes connected by edges. We can already envision a derived type called bemnode and a derived type called polygon made-up of a collection of bemnodes. The location of the nodes of each polygon must be read from files that contain the list of nodes and their attributes.
164 CHAPTER 16. BOUNDARY ELEMENT METHOD
16.5.2 Geometry module
It is clear that the development of the code will require a module that include some amount of analytical geometry manipulation. These include the calculation of the angle between 2 vectors, the inner product and the vector product of two vectors. The use of the following equations will be useful:
u · v = u 1 v 1 + u 2 v 2 (16.57) u × v = u 1 v 2 − u 2 v 1 (16.58) ‖u‖ = u^21 + u^22 (16.59) sin θ =
u × v ‖u‖ ‖v‖
cos θ =
u · v ‖u‖ ‖v‖
The normal distance of a point S to the edge spanned by the points A and B is given by
n =
AS~ × AB~
‖ AB~ ‖
= (xs − xa) sin β − (ya − ys) cos β (16.63)
where
cos β =
xb − xa ∆s
, sin β =
yb − ya ∆s
, ∆s =
√ (xb − xa)^2 + (yb − ya)^2 (16.64)
The point of intersection of the edge and the normal is given by
xh = (ys − ya) cos β sin β − (xs − xa) sin^2 β + xs (16.65) yh = (xs − xa) cos β sin β − (ys − ya) cos^2 β + ys (16.66)
The distance sa from A to H and sb from B to H are given by
sA =
SA~ · AB~
‖ AB~ ‖
= (xa − xs) cos β + (ya − ys) sin β (16.67)
sB =
SB~ · AB~
‖ AB~ ‖
= (xb − xs) cos β + (yb − ys) sin β (16.68)
