FEM Lecture Notes by Peter Hunter, Andrew Pullian

FEM/BEM NOTES
Professor Peter Hunter

[email protected]
Associate Professor Andrew Pullan

[email protected]
Department of Engineering Science

The University of Auckland
New Zealand
April 7, 2003
c Copyright 1997-2003
Department of Engineering Science
The University of Auckland
Contents
1 Finite Element Basis Functions 1

1.1 Representing a One-Dimensional Field . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Linear Basis Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Basis Functions as Weighting Functions . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Quadratic Basis Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.5 Two- and Three-Dimensional Elements . . . . . . . . . . . . . . . . . . . . . . . 7
1.6 Higher Order Continuity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.7 Triangular Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.8 Curvilinear Coordinate Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.9 CMISS Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2 Steady-State Heat Conduction 23

2.1 One-Dimensional Steady-State Heat Conduction . . . . . . . . . . . . . . . . . . 23
2.1.1 Integral equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.1.2 Integration by parts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.1.3 Finite element approximation . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.4 Element integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.1.5 Assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.1.6 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1.7 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1.8 Fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 An x-Dependent Source Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3 The Galerkin Weight Function Revisited . . . . . . . . . . . . . . . . . . . . . . . 31
2.4 Two and Three-Dimensional Steady-State Heat Conduction . . . . . . . . . . . . . 32
2.5 Basis Functions - Element Discretisation . . . . . . . . . . . . . . . . . . . . . . . 34
2.6 Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.7 Assemble Global Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.8 Gaussian Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.9 CMISS Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3 The Boundary Element Method 43

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2 The Dirac-Delta Function and Fundamental Solutions . . . . . . . . . . . . . . . . 43
3.2.1 Dirac-Delta function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2.2 Fundamental solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
ii CONTENTS
3.3 The Two-Dimensional Boundary Element Method . . . . . . . . . . . . . . . . . . 48

3.4 Numerical Solution Procedures for the Boundary Integral Equation . . . . . . . . . 53
3.5 Numerical Evaluation of Coefficient Integrals . . . . . . . . . . . . . . . . . . . . 55
3.6 The Three-Dimensional Boundary Element Method . . . . . . . . . . . . . . . . . 57
3.7 A Comparison of the FE and BE Methods . . . . . . . . . . . . . . . . . . . . . . 58
3.8 More on Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.8.1 Logarithmic quadrature and other special schemes . . . . . . . . . . . . . 60
3.8.2 Special solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.9 The Boundary Element Method Applied to other Elliptic PDEs . . . . . . . . . . . 61
3.10 Solution of Matrix Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.11 Coupling the FE and BE techniques . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.12 Other BEM techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.12.1 Trefftz method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.12.2 Regular BEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.13 Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.14 Axisymmetric Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.15 Infinite Regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.16 Appendix: Common Fundamental Solutions . . . . . . . . . . . . . . . . . . . . . 72
3.16.1 Two-Dimensional equations . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.16.2 Three-Dimensional equations . . . . . . . . . . . . . . . . . . . . . . . . 72
3.16.3 Axisymmetric problems . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.17 CMISS Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4 Linear Elasticity 75
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.2 Truss Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3 Beam Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.4 Plane Stress Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.4.1 Notes on calculating nodal loads . . . . . . . . . . . . . . . . . . . . . . . 83
4.5 Three-Dimensional Elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.5.1 Weighted Residual Integral Equation . . . . . . . . . . . . . . . . . . . . 85
4.5.2 The Principle of Virtual Work . . . . . . . . . . . . . . . . . . . . . . . . 86
4.5.3 The Finite Element Approximation . . . . . . . . . . . . . . . . . . . . . 87
4.6 Linear Elasticity with Boundary Elements . . . . . . . . . . . . . . . . . . . . . . 89
4.7 Fundamental Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.8 Boundary Integral Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.9 Body Forces (and Domain Integrals in General) . . . . . . . . . . . . . . . . . . . 96
4.10 CMISS Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5 Transient Heat Conduction 99

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.2 Finite Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.2.1 Explicit Transient Finite Differences . . . . . . . . . . . . . . . . . . . . . 99
5.2.2 Von Neumann Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . 101
5.2.3 Higher Order Approximations . . . . . . . . . . . . . . . . . . . . . . . . 102
CONTENTS iii
5.3 The Transient Advection-Diffusion Equation . . . . . . . . . . . . . . . . . . . . 103

5.4 Mass lumping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.5 CMISS Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6 Modal Analysis 111

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.2 Free Vibration Modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.3 An Analytic Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.4 Proportional Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
6.5 CMISS Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
7 Domain Integrals in the BEM 117

7.1 Achieving a Boundary Integral Formulation . . . . . . . . . . . . . . . . . . . . . 117
7.2 Removing Domain Integrals due to Inhomogeneous Terms . . . . . . . . . . . . . 118
7.2.1 The Galerkin Vector technique . . . . . . . . . . . . . . . . . . . . . . . . 118
7.2.2 The Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
7.2.3 Complementary Function-Particular Integral method . . . . . . . . . . . . 120
7.3 Domain Integrals Involving the Dependent Variable . . . . . . . . . . . . . . . . . 120
7.3.1 The Perturbation Boundary Element Method . . . . . . . . . . . . . . . . 121
7.3.2 The Multiple Reciprocity Method . . . . . . . . . . . . . . . . . . . . . . 122
7.3.3 The Dual Reciprocity Boundary Element Method . . . . . . . . . . . . . . 124
8 The BEM for Parabolic PDES 135

8.1 Time-Stepping Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
8.1.1 Coupled Finite Difference - Boundary Element Method . . . . . . . . . . . 135
8.1.2 Direct Time-Integration Method . . . . . . . . . . . . . . . . . . . . . . . 137
8.2 Laplace Transform Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
8.3 The DR-BEM For Transient Problems . . . . . . . . . . . . . . . . . . . . . . . . 139
8.4 The MRM for Transient Problems . . . . . . . . . . . . . . . . . . . . . . . . . . 140
Bibliography 143
Index 147
Chapter 1
Finite Element Basis Functions
1.1 Representing a One-Dimensional Field

Consider the problem of finding a mathematical expression u (x) to represent a one-dimensional
field e.g., measurements of temperature u against distance x along a bar, as shown in Figure 1.1a.
u u
+ +
+ + + +
+ ++ + + ++ +
+ + + +
+ + + +
+ + + +
+ +
+ +
x x
(a) (b)
F IGURE 1.1: (a) Temperature distribution u (x) along a bar. The points are the measured
temperatures. (b) A least-squares polynomial fit to the data, showing the unacceptable oscillation
between data points.
One approach would be to use a polynomial expression u (x) = a + bx + cx 2 + dx3 + : : :

and to estimate the values of the parameters a, b, c and d from a least-squares fit to the data. As
the degree of the polynomial is increased the data points are fitted with increasing accuracy and
polynomials provide a very convenient form of expression because they can be differentiated and
integrated readily. For low degree polynomials this is a satisfactory approach, but if the polynomial
order is increased further to improve the accuracy of fit a problem arises: the polynomial can be
made to fit the data accurately, but it oscillates unacceptably between the data points, as shown in
Figure 1.1b.
To circumvent this, while retaining the advantages of low degree polynomials, we divide the
bar into three subregions and use low order polynomials over each subregion - called elements. For
later generality we also introduce a parameter s which is a measure of distance along the bar. u is
plotted as a function of this arclength in Figure 1.2a. Figure 1.2b shows three linear polynomials
in s fitted by least-squares separately to the data in each element.
2 F INITE E LEMENT BASIS F UNCTIONS
u u
+ +
+ +
+ +
+ + + + + + + + + +
+ + + + + +
+ + + + + +
+ +
+ +
s s
(a) (b)
F IGURE 1.2: (a) Temperature measurements replotted against arclength parameter s. (b) The s
domain is divided into three subdomains, elements, and linear polynomials are independently fitted
to the data in each subdomain.
1.2 Linear Basis Functions

A new problem has now arisen in Figure 1.2b: the piecewise linear polynomials are not continuous
in u across the boundaries between elements. One solution would be to constrain the parameters a,
b, c etc. to ensure continuity of u across the element boundaries, but a better solution is to replace
the parameters a and b in the first element with parameters u1 and u2 , which are the values of u at
the two ends of that element. We then define a linear variation between these two values by
u ( ) = (1 ; ) u1 + u2
where (0 1) is a normalized measure of distance along the curve.
We define
'1 ( ) = 1 ;
'2 ( ) =
such that
u ( ) = '1 ( ) u1 + '2 ( ) u2
and refer to these expressions as the basis functions associated with the nodal parameters u 1 and
u2. The basis functions '1 ( ) and '2 ( ) are straight lines varying between 0 and 1 as shown in
Figure 1.3.
It is convenient always to associate the nodal quantity u n with element node n and to map the
temperature U defined at global node onto local node n of element e by using a connectivity
matrix (n; e) i.e.,
un = U(n;e)
where (n; e) = global node number of local node n of element e. This has the advantage that the
1.2 L INEAR BASIS F UNCTIONS 3
'1 ( ) '2 ( )
1 11

0 1 0 1
F IGURE 1.3: Linear basis functions '1 ( ) = 1 ; and '2 ( ) = .
interpolation
u ( ) = '1 ( ) u1 + '2 ( ) u2
holds for any element provided that u 1 and u2 are correctly identified with their global counterparts,
as shown in Figure 1.4. Thus, in the first element
node 1 node 2 node 3 node 4

U1 U2 U3 U4
global nodes:
x
element
nodes: u1 u2 u1 u2 u1 u2
111111
000000 1111111
0000000 111111
000000
0 1 0 1 0 1
element 1 element 2 element 3
F IGURE 1.4: The relationship between global nodes and element nodes.
u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.1)
with u1 = U1 and u2 = U2 .
In the second element u is interpolated by
u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.2)
with u1 = U2 and u2 = U3 , since the parameter U2 is shared between the first and second elements
the temperature field u is implicitly continuous. Similarly, in the third element u is interpolated by
u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.3)
with u1 = U3 and u2 = U4 , with the parameter U3 being shared between the second and third
elements. Figure 1.6 shows the temperature field defined by the three interpolations (1.1)–(1.3).
node 1
+ node 3
+ node 2
+ +
+
node 4
+ + +
+ + +
+ + +
+
+
element 1 element 2 element 3
s
F IGURE 1.5: Temperature measurements fitted with nodal parameters and linear basis functions.
The fitted temperature field is now continuous across element boundaries.
1.3 Basis Functions as Weighting Functions

It is useful to think of the basis functions as weighting functions on the nodal parameters. Thus, in
element 1
at =0 u (0) = (1 ; 0) u1 + 0u2 = u1
which is the value of u at the left hand end of the element and has no dependence on u 2
1

1 = 1 ; 1 u + 1u = 3u + 1u

at = u
4 4 4 1 4 2 4 1 4 2
which depends on u1 and u2 , but is weighted more towards u1 than u2
1

1 = 1 ; 1 u + 1u = 1u + 1u

at = u
2 2 2 1 2 2 2 1 2 2
which depends equally on u1 and u2
3
3
at =
4
u 4 = 1 ; 43 u1 + 34 u2 = 14 u1 + 34 u2
1.3 BASIS F UNCTIONS AS W EIGHTING F UNCTIONS 5
which depends on u1 and u2 but is weighted more towards u2 than u1
at =1 u (1) = (1 ; 1) u1 + 1u2 = u2
which is the value of u at the right hand end of the region and has no dependence on u 1 .
Moreover, these weighting functions can be considered as global functions, as shown in Fig-
ure 1.6, where the weighting function wn associated with global node n is constructed from the
basis functions in the elements adjacent to that node.
w1
(a)
s
w2
(b)
s
w3
(c)
w4
(d)
s
F IGURE 1.6: (a) : : : (d) The weighting functions wn associated with the global nodes n = 1 : : : 4,
respectively. Notice the linear fall off in the elements adjacent to a node. Outside the immediately
adjacent elements, the weighting functions are defined to be zero.
For example, w2 weights the global parameter U2 and the influence of U2 falls off linearly in
the elements on either side of node 2.
We now have a continuous piecewise parametric description of the temperature field u ( ) but
in order to define u (x) we need to define the relationship between x and for each element. A
convenient way to do this is to define x as an interpolation of the nodal values of x.
For example, in element 1
x ( ) = '1 ( ) x1 + '2 ( ) x2 (1.4)
and similarly for the other two elements. The dependence of temperature on x, u (x), is therefore
defined by the parametric expressions

X
u ( ) = 'n ( ) un
n
X
x( ) = 'n ( ) xn
n
where summation is taken over all element nodes (in this case only 2) and the parameter (the
“element coordinate”) links temperature u to physical position x. x ( ) provides the mapping
between the mathematical space 0 1 and the physical space x1 x x2 , as illustrated in
Figure 1.7.
u1
u2
u

0 1
= 0:2 u (x) at = 0:2
u1

u2
x
x
0 x1 x2
x2
x1

0 1
= 0:2
F IGURE 1.7: Illustrating how x and u are related through the normalized element coordinate .
The values of x ( ) and u ( ) are obtained from a linear interpolation of the nodal variables and
then plotted as u (x). The points at = 0:2 are emphasized.
1.4 QUADRATIC BASIS F UNCTIONS 7
1.4 Quadratic Basis Functions

The essential property of the basis functions defined above is that the basis function associated
with a particular node takes the value of 1 when evaluated at that node and is zero at every other
node in the element (only one other in the case of linear basis functions). This ensures the linear
independence of the basis functions. It is also the key to establishing the form of the basis functions
for higher order interpolation. For example, a quadratic variation of u over an element requires
three nodal parameters u1 , u2 and u3
u ( ) = '1 ( ) u1 + '2 ( ) u2 + '3 ( ) u3 (1.5)
The quadratic basis functions are shown, with their mathematical expressions, in Figure 1.8. Notice
that since '1 ( ) must be zero at = 0:5 (node 2), '1 ( ) must have a factor ( ; 0:5) and since it
is also zero at = 1 (node 3), another factor is ( ; 1). Finally, since '1 ( ) is 1 at = 0 (node 1)
we have '1 ( ) = 2 ( ; 1) ( ; 0:5). Similarly for the other two basis functions.
'1 ( ) '2 ( )
1 1

0 0:5 1 0 0:5 1
(a) '1 ( ) = 2 ( ; 1) ( ; 0:5) (b) '2 ( ) = 4 (1 ; )
'3 ( )

0 0:5 1
(c) '3 ( ) = 2 ( ; 0:5)
F IGURE 1.8: One-dimensional quadratic basis functions.
1.5 Two- and Three-Dimensional Elements

Two-dimensional bilinear basis functions are constructed from the products of the above one-
dimensional linear functions as follows
Let
u (1; 2) = '1 (1; 2) u1 + '2 (1; 2) u2 + '3 (1; 2) u3 + '4 (1; 2) u4
where
'1 (1; 2) = (1 ; 1) (1 ; 2)
'2 (1; 2) = 1 (1 ; 2)
'3 (1; 2) = (1 ; 1) 2
(1.6)
'4 (1; 2) = 12

Note that '1 (1 ; 2 ) = '1 (1 )'1 (2 ) where '1 (1 ) and '1 (2 ) are the one-dimensional linear
basis functions. Similarly, ' 2 (1 ; 2 ) = '2 (1 )'1 (2 ) : : : etc.
These four bilinear basis functions are illustrated in Figure 1.9.
'1 '3
node 3
2 2
0
node 1 node 4
node 2 1 1
1
'2 '4
2
2
0 0
1 1 1 1
F IGURE 1.9: Two-dimensional bilinear basis functions.
Notice that 'n (1 ; 2 ) is 1 at node n and zero at the other three nodes. This ensures that the
temperature u (1 ; 2 ) receives a contribution from each nodal parameter un weighted by 'n (1 ; 2 )
and that when u (1 ; 2 ) is evaluated at node n it takes on the value un .
As before the geometry of the element is defined in terms of the node positions (x n ; yn ), n =
1.5 T WO - AND T HREE -D IMENSIONAL E LEMENTS 9
1; : : : ; 4 by
X
x= 'n (1; 2) xn
n
X
y= 'n (1; 2) yn
n
which provide the mapping between the mathematical space ( 1 ; 2 ) (where 0 1 ; 2 1) and
the physical space (x; y ).
Higher order 2D basis functions can be similarly constructed from products of the appropriate
1D basis functions. For example, a six-noded (see Figure 1.10) quadratic-linear element (quadratic
in 1 and linear in 2 ) would have
X
6
u= 'n (1; 2) un
n=1
where
'1 (1; 2) = 2 (1 ; 1) (1 ; 0:5) (1 ; 2) '2 (1; 2) = 41 (1 ; 1) (1 ; 2) (1.7)
'3 (1; 2) = 21 (1 ; 0:5) (1 ; 2) '4 (1; 2) = 2 (1 ; 1) (1 ; 0:5) 2 (1.8)
'5 (1; 2) = 41 (1 ; 1) 2 '6 (1; 2) = 21 (1 ; 0:5) 2 (1.9)
2
1
4 5 6
1 2 3
0 1
0 0:5 1
F IGURE 1.10: A 6-node quadratic-linear element (node numbers circled).
Three-dimensional basis functions are formed similarly, e.g., a trilinear element basis has eight
nodes (see Figure 1.11) with basis functions
'1 (1; 2; 3) = (1 ; 1) (1 ; 2) (1 ; 3) '2 (1; 2; 3) = 1 (1 ; 2) (1 ; 3) (1.10)
'3 (1; 2; 3) = (1 ; 1) 2 (1 ; 3) '4 (1; 2; 3) = 12 (1 ; 3) (1.11)
'5 (1; 2; 3) = (1 ; 1) (1 ; 2) 3 '6 (1; 2; 3) = 1 (1 ; 2) 3 (1.12)
'7 (1; 2; 3) = (1 ; 1) 23 '8 (1; 2; 3) = 123 (1.13)
3
7
5 8
6
2
3
1 4
2
1
F IGURE 1.11: An 8-node trilinear element.
1.6 Higher Order Continuity

All the basis functions mentioned so far are Lagrange 1 basis functions and provide continuity of u
across element boundaries but not higher order continuity. Sometimes it is desirable to use basis
functions which also preserve continuity of the derivative of u with respect to across element
du
boundaries. A convenient way to achieve this is by defining two additional nodal parameters
. The basis functions are chosen to ensure that
d n

du = du

du = du

= u01 and = u02
d =0 d 1 d =1 d 2
and since un is shared between adjacent elements derivative continuity is ensured. Since the num-
ber of element parameters is 4 the basis functions must be cubic in . To derive these cubic
Hermite2 basis functions let
u ( ) = a + b + c 2 + d 3;
du = b + 2c + 3d 2;
d
1
Joseph-Louis Lagrange (1736-1813).
2
Charles Hermite (1822-1901).
1.6 H IGHER O RDER C ONTINUITY 11
and impose the constraints
u (0) = a = u1
u (1) = a + b + c + d = u2
du (0) = b = u 0
d 1
du (1) = b + 2c + 3d = u0
d 2
These four equations in the four unknowns a, b, c and d are solved to give
a = u1
b = u01
c = 3u2 ; 3u1 ; 2u01 ; u02
d = u01 + u02 + 2u1 ; 2u2
Substituting a, b, c and d back into the original cubic then gives
u ( ) = u1 + u01 + (3u2 ; 3u1 ; 2u01 ; u02) 2 + (u01 + u02 + 2u1 ; 2u2) 3

or, rearranging,
u ( ) = 01 ( ) u1 + 11 ( ) u01 + 02 ( ) u2 + 12 ( ) u02 (1.14)
where the four cubic Hermite basis functions are drawn in Figure 1.12.
du
One further step is required to make cubic Hermite basis functions useful in practice. The
derivative defined at node n is dependent upon the element -coordinate in the two ad-
d n
jacent elements. It is much more useful to define a global node derivative
du where s is
ds n
arclength and then use
du du ds
d =
ds d (1.15)
n (n;e) n
ds
where
d n is an element scale factor which scales the arclength derivative of global node
to the -coordinate derivative of element node n. Thus duds is constrained to be continuous
across element boundaries rather than
du . A two- dimensional bicubic Hermite basis requires four
d
derivatives per node
@u ; @u
u; @ and
@2u
1 @ 2 @1 @2
01 ( ) = 1 ; 3 2 + 2 3
slope = 1
1
11 ( ) = ( ; 1)2
0
1
0
1
02 ( ) = 2(3 ; 2 )
1
12 ( ) = 2( ; 1)
0
1
0 slope = 1
1
F IGURE 1.12: Cubic Hermite basis functions.
The need for the second-order cross-derivative term can be explained as follows; If u is cubic in 1
and cubic in 2 , then
@u @u is cubic in and quadratic
is quadratic in 1 and cubic in 2 , and 1
@1 @2
in 2 . Now consider the side 1–3 Thecubic variation of u with 2 is specified by
@uinFigure 1.13. @u
the four nodal parameters u1 , , u3 and . But since
@u (the normal derivative) is
@2 1 @2 3 @1
also cubic in 2 along that side and is entirely independent of these four parameters,
@u four additional
@u
parameters are required to specify this cubic. Two of these are specified by and ,
@2u @2u @ 1 1 @ 1 3
and the remaining two by
@1 @2 1 and @1 @2 3.
1.6 H IGHER O RDER C ONTINUITY 13
2
@u node 3 node 4
@1 3
@u
@1 1
1 node 1 node 2
@u
F IGURE 1.13: Interpolation of nodal derivative along side 1–3.
@1
The bicubic interpolation of these nodal parameters is given by
u (1; 2) = 01 (1) 01 (2) u1 + 02 (1) 01 (2) u2

+ 01 (1) 02 (2) u3 + 02 (1) 02 (2) u4
@u @u
+ 1 (1) 1 (2)
1 0 + 2 (1) 1 (2)
1 0
@
@u 1 1
@@u1 2
+ 11 (1) 02 (2) @ + 12 (1) 02 (2) @
@u1 3 @u1 4
+ 01 (1) 11 (2) + 02 (1) 11 (2) (1.16)
@@u2 1 @@u2 2
+ 01 (1) 12 (2) @ + 02 (1) 12 (2) @
@22 u 3 2 @42 u
+ 11 (1) 11 (2) + 12 (1) 11 (2)
@ @
@12 u 2 1 @@12@u 2 2
+ 11 (1) 12 (2) @ @ + 12 (1) 12 (2) @ @
1 2 3 1 2 4
where
01 ( ) = 1 ; 3 2 + 2 3
11 ( ) = ( ; 1)2 (1.17)
02 ( ) = 2 (3 ; 2 )
12 ( ) = 2 ( ; 1)
are the one-dimensional cubic Hermite basis functions (see Figure 1.12).
As in the one-dimensional case above, to preserve derivative continuity in physical x-coordinate
space as well as in -coordinate space the global node derivatives need to be specified with respect
to physical arclength. There are now two arclengths to consider: s 1 , measuring arclength along the
1-coordinate, and s2, measuring arclength along the 2-coordinate. Thus
@u @u @s
@1 n = @s1 (n;e) @1 n
1
@u @u @s
@2 n = @s2 (n;e) @2 n
2
(1.18)
@2u @2u ds ds
@1 @2 n = 1

@s1 @s2 (n;e) d1 n d2 n
2
ds ds
1 2
where and
d1 n d2 n are element scale factors which scale the arclength derivatives of
global node to the -coordinate derivatives of element node n.
The bicubic Hermite basis is a powerful shape descriptor for curvilinear surfaces. Figure 1.14
shows a four element bicubic Hermite surface in 3D space where each node has the following
twelve parameters
@x ; @x ; @ 2 x ; y; @y ; @y ; @ 2 y ; z; @z ; @z and @ 2 z
x; @s
1 @s2 @s1 @s2 @s1 @s2 @s1 @s2 @s1 @s2 @s1 @s2
1.7 Triangular Elements

Triangular elements cannot use the 1 and 2 coordinates defined above for tensor product elements
(i.e., two- and three- dimensional elements whose basis functions are formed as the product of one-
dimensional basis functions). The natural coordinates for triangles are based on area ratios and are
called Area Coordinates . Consider the ratio of the area formed from the points 2, 3 and P (x; y )
in Figure 1.15 to the total area of the triangle
1 x y
L1 = Area < P 23 > = 1 1 x y = = (a + b x + c y) = (2)
Area < 123 > 2 1 x23 y23 1 1 1
1.7 T RIANGULAR E LEMENTS 15
z
12 parameters per node
2
y
1
x
F IGURE 1.14: A surface formed by four bicubic Hermite elements.
3 (x3 ; y3)
Area P 23
P(x,y )
2 (x2 ; y2)
1 (x1 ; y1)
L1 = 0
L1 = 13
L1 = 32
L1 = 1
F IGURE 1.15: Area coordinates for a triangular element.
1 x y
1 1
where = 21 1 x2 y2 is the area of the triangle with vertices 123, and a1 = x2 y3 ; x3 y2 ; b1 =
1 x3 y3
y2 ; y3; c1 = x3 ; x2 .
Notice that L1 is linear in x and y . Similarly, area coordinates for the other two triangles
containing P and two of the element vertices are
1 x y
L2 = Area < P 13 > = 1 1 x y = = (a + b x + c y) = (2)
Area < 123 > 2 1 x31 y31 2 2 2
1 x y
L3 = Area < P 12 > = 1 1 x y = = (a + b x + c y) = (2)
Area < 123 > 2 1 x12 y12 3 3 3
where a2 = x3 y1 ; x1 y3 ; b2 = y3 ; y1 ; c2 = x1 ; x3 and a3 = x1 y2 ; x2 y1 ; b3 = y1 ; y2 ; c3 = x2 ; x1 .
Notice that L1 + L2 + L3 = 1.
Area coordinate L1 varies linearly from L1 = 0 when P lies at node 2 or 3 to L1 = 1 when P
lies at node 1 and can therefore be used directly as the basis function for node 1 for a three node
triangle. Thus, interpolation over the triangle is given by
u (x; y) = '1 (x; y) u1 + '2 (x; y) u2 + '3 (x; y) u3

where '1 = L1 , '2 = L2 and '3 = L3 = 1 ; L1 ; L2 .
Six node quadratic triangular elements are constructed as shown in Figure 1.16.
'1 = L1 (2L1 ; 1)
'2 = L2 (2L2 ; 1)
'3 = L3 (2L3 ; 1)
4 6
'4 = 4L1 L2
'5 = 4L2 L3
'6 = 4L3 L1
2 3
5
F IGURE 1.16: Basis functions for a six node quadratic triangular element.
1.8 Curvilinear Coordinate Systems

It is sometimes convenient to model the geometry of the region (over which a finite element solu-
tion is sought) using an orthogonal curvilinear coordinate system. A 2D circular annulus, for ex-
1.8 C URVILINEAR C OORDINATE S YSTEMS 17
ample, can be modelled geometrically using one element with cylindrical polar (r; )-coordinates,
e.g., the annular plate in Figure 1.17a has two global nodes, the first with r = r 1 and the second
with r = r2 .
y 2
2
2 3
2 3
1 2 x
1 4 1 4
0 r 1
r1 r2
(a) (b) (c)
F IGURE 1.17: Defining a circular annulus with one cylindrical polar element. Notice that element
vertices 1 and 2 in (r; )-space or (1 ; 2 )-space, as shown in (b) and (c), respectively, map onto the
single global node 1 in (x; y )-space in (a). Similarly, element vertices 3 and 4 map onto global
node 2.
Global nodes 1 and 2, shown in (x; y )-space in Figure 1.17a, each map to two element vertices
in (r; )-space, as shown in Figure 1.17b, and in (1; 2)-space, as shown in Figure 1.17c. The
(r; ) coordinates at any (1; 2) point are given by a bilinear interpolation of the nodal coordinates
rn and n as
r = 'n (1; 2) rn
= 'n (1; 2) n
where the basis functions 'n (1 ; 2 ) are given by (1.6).
Three orthogonal curvilinear coordinate systems are defined here for use in later sections.
Cylindrical polar (r; ; z ) :
x = r cos
y = r sin (1.19)
z=z
Spherical polar (r; ; ) :
x = r cos cos
y = r sin cos (1.20)
z = r sin
Prolate spheroidal (; ; ) :
x = d cosh cos
y = d sinh sin cos (1.21)
z = d sinh sin sin
y

r
z
d
x
F IGURE 1.18: Prolate spheroidal coordinates.
The prolate spheroidal coordinates rae illustrated in Figure 1.18 and a single prolate spheroidal
element is shown in Figure 1.19. The coordinates (; ; ) are all trilinear in ( 1 ; 2 ; 3 ). Only four
global nodes are required provided the four global nodes map to eight element nodes as shown in
Figure 1.19.
1.8 C URVILINEAR C OORDINATE S YSTEMS 19
(a) (b)
24 y
z 2
1
1 3
3
x
(c) (d)
2
2 4 2 4
90o
2 1 4
1 3
2 4
0 3
1 3 1 3
2
1 3

F IGURE 1.19: A single prolate spheroidal element, shown (a) in (x; y; z )-coordinates, (c) in
(; ; )-coordinates and (d) in (1 ; 2 ; 3 )-coordinates, (b) shows the orientation of the
i -coordinates on the prolate spheroid.
1.9 CMISS Examples

1. To define a 2D bilinear finite element mesh run the CMISS example number 111. The nodes
should be positioned as shown in Figure 1.20. After defining elements the mesh should
appear like the one shown in Figure 1.21.
6
4 5
1 2
3
F IGURE 1.20: Node positions for example 111.
1 2
F IGURE 1.21: 2D bilinear finite element mesh for example 111.
2. To refine a mesh run the CMISS example 113. After the first refine the mesh should appear
like the one shown in Figure 1.22.
3. To define a quadratic-linear element run the cmiss example 115.
4. To define a 3D trilinear element run CMISS example 121.
5. To define a 2D cubic Hermite-linear finite element mesh run example 114.
6. To define a triangular element mesh run CMISS example 116 (see Figure 1.24).
7. To define a bilinear mesh in cylindrical polar coordinates run CMISS example 122.
1.9 CMISS E XAMPLES 21
1 3
7 9
2
1 3 2 4
8 5
4 10
6
F IGURE 1.22: First refined mesh for example 113
1 11 3
7 13 9
2
1 5 3 6 2 4
5
4 12 8 14 10
6
F IGURE 1.23: Second refined mesh for example 113
3
4
1 2
F IGURE 1.24: Defining a triangular mesh for example 116

Chapter 2
Steady-State Heat Conduction
2.1 One-Dimensional Steady-State Heat Conduction

Our first example of solving a partial differential equation by finite elements is the one-dimensional
steady-state heat equation. The equation arises from a simple heat balance over a region of con-
ducting material:
Rate of change of heat flux = heat source per unit volume
or
d (heat flux) + heat sink per unit volume = 0
dx
or

d ;k du + q (u; x) = 0
dx dx
where u is temperature, q (u; x) the heat sink and k the thermal conductivity (Watts=m= C).
Consider the case where q = u
d
du
; dx k dx + u = 0 0 < x < 1 (2.1)
subject to boundary conditions: u (0) = 0 and u (1) = 1.

This equation (with k = 1) has an exact solution
;
u (x) = e2 e; 1 ex ; e;x
(2.2)
with which we can compare the approximate finite element solutions.

To solve Equation (2.1) by the finite element method requires the following steps:
1. Write down the integral equation form of the heat equation.
2. Integrate by parts (in 1D) or use Green’s Theorem (in 2D or 3D) to reduce the order of
derivatives.
24 S TEADY-S TATE H EAT C ONDUCTION
3. Introduce the finite element approximation for the temperature field with nodal parameters
and element basis functions.
4. Integrate over the elements to calculate the element stiffness matrices and RHS vectors.
5. Assemble the global equations.
6. Apply the boundary conditions.
7. Solve the global equations.
8. Evaluate the fluxes.
2.1.1 Integral equation

Rather than solving Equation (2.1) directly, we form the weighted residual
Z
R!:dx = 0 (2.3)
where R is the residual

d k du + u

R = ; dx dx (2.4)
for an approximate solution u and ! is a weighting function to be chosen below. If u were an exact
solution over the whole domain, the residual R would be zero everywhere. But, given that in real
engineering problems this will not be the case, we try to obtain an approximate solution u for which
the residual or error (i.e., the amount by which the differential equation is not satisfied exactly at a
point) is distributed evenly over the domain. Substituting Equation (2.4) into Equation (2.3) gives
Z1 d du
; dx k dx ! + u! dx = 0 (2.5)
0
This formulation of the governing equation can be thought of as forcing the residual or error to
be zero in a spatially averaged sense. More precisely, ! is chosen such that the residual is kept
orthogonal to the space of functions used in the approximation of u (see step 3 below).
2.1.2 Integration by parts

A major advantage of the integral equation is that the order of the derivatives inside the integral can
be reduced from two to one by integrating by parts (or, equivalently for 2D problems, by applying
Green’s theorem - see later). Thus, substituting f
du into the integration by parts
= ! and g = ;k dx
2.1 O NE -D IMENSIONAL S TEADY-S TATE H EAT C ONDUCTION 25
formula
Z1 dg Z1 df
f dx dx = [f:g]0 ; g dx dx
1
0 0
gives
Z1 d du du 1 Z1 du d!
! dx ;k dx dx = ! ;k dx ; ;k dx dx dx
0
0 0
and Equation (2.5) becomes
Z1 du d! du 1
k + u! dx = k !
dx dx dx (2.6)
0
0
2.1.3 Finite element approximation

We divide the domain 0 < x < 1 into 3 equal length elements and replace the continuous field
variable u (x) within each element by the parametric finite element approximation
u ( ) = '1 ( ) u1 + '2 ( ) u2 = 'n ( ) un

x ( ) = '1 ( ) x1 + '2 ( ) x2 = 'n ( ) xn
(summation implied by repeated index) where ' 1 ( ) = 1 ; and '2 ( ) = are the linear basis
functions for both u and x.
We also choose ! = 'm (called the Galerkin1 assumption). This forces the residual R to be
orthogonal to the space of functions used to represent the dependent variable u, thereby ensuring
that the residual, or error, is monotonically reduced as the finite element mesh is refined (see later
for a more complete justification of this very important step) .
The domain integral in Equation (2.6) can now be replaced by the sum of integrals taken sepa-
rately over the three elements
Z1 Z31 Z32 Z1
dx = dx + dx + dx
0 0 1 2
3 3
1
Boris G. Galerkin (1871-1945). Galerkin was a Russian engineer who published his first technical paper on the
buckling of bars while imprisoned in 1906 by the Tzar in pre-revolutionary Russia. In many Russian texts the Galerkin
finite element method is known as the Bubnov-Galerkin method. He published a paper using this idea in 1915. The
method was also attributed to I.G. Bubnov in 1913.
and each element integral is then taken over -space
Zx2 Z1
dx = J d
x1 0
dx
where J = is the Jacobian of the transformation from x coordinates to coordinates.
d
2.1.4 Element integrals
The element integrals arising from the LHS of Equation (2.6) have the form
Z1 du d!
k + u! J d
dx dx (2.7)
0
where u = 'n un and ! = 'm . Since 'n and 'm are both functions of the derivatives with respect
to x need to be converted to derivatives with respect to . Thus Equation (2.7) becomes
Z1 d'm d d'n d
un k + 'm 'n J d
d dx d dx (2.8)
0
Notice that un has been taken outside the integral because it is not a function of . The term
d is
dx
1
evaluated by substituting the finite element approximation x ( ) = ' n :Xn . In this case x = or
3
d = 3 and the Jacobian is J = dx = 1 . The term multiplying the nodal parameters u is called
n
dx d 3
the element stiffness matrix, Emn
Z1 d'm d d'n d Z1 d'm d'n 1
Emn = k + 'm'n J d =
d dx d dx k 3 3 + 'm'n d d d 3
0 0
where the indices m and n are 1 or 2. To evaluate Emn we substitute the basis functions
'1 ( ) = 1 ; or
d'1 = ;1
d
'2 ( ) = or
d'2 = 1
d
Node 1 2 3 4
Node 1
X X 0 0 U1 X
Node 2
X X X 0 U2 X
=
Node 3
0 X X X U3 X
Node 4
0 0 X X U4 X
F IGURE 2.1: The rows of the global stiffness matrix are generated from the global weight
functions. The bar is shown at the top divided into three elements.
Thus,
Z d' 2 ! Z ;
1
E11 = 3
1
9k d 1
+ ('1) 2 1
1

d = 3 9k (;1)2 + (1 ; )2 d = 13 9k + 31
0 0
and, similarly,

E12 = E21 = 31 ;9k + 16
1
1
E22 = 3 9k + 3
1 ;9k + 1 1 ;;9k + 1
Emn = 13;;9k +31 31 ;9k + 16
3 6 3 3
Notice that the element stiffness matrix is symmetric. Notice also that the stiffness matrix, in this
particular case, is the same for all elements. For simplicity we put k = 1 in the following steps.
2.1.5 Assembly
The three element stiffness matrices (with k = 1) are assembled into one global stiffness matrix.
This process is illustrated in Figure 2.1 where rows 1; ::; 4 of the global stiffness matrix (shown here
multiplied by the vector of global unknowns) are generalised from the weight function associated
with nodes 1; ::; 4.
Note how each element stiffness matrix (the smaller square brackets in Figure 2.1) overlaps
with its neighbour because they share a common global node. The assembly process gives
2 28 ; 53 0 0
3 2U 3
66; 5318 289 + 289 ; 1853 0 77 66U1277
9 18
4 0 ; 5318 289 + 289 ; 5318 5 4U35

0 053 28 ; 18 9 U4
Notice that the first row (generating heat flux at node 1) has zeros multiplying U 3 and U4 since
nodes 3 and 4 have no direct connection through the basis functions to node 1. Finite element
matrices are always sparse matrices - containing many zeros - since the basis functions are local
to elements.
The RHS of Equation (2.6) is
du x=1 du du
k dx ! = k dx ! ; k dx !
(2.9)
x=0 x=1 x=0
To evaluate these expressions consider the weighting function ! corresponding to each global node
(see Fig.1.6). For node 1 !1 is obtained from the basis function ' 1 associated with the first node
of element 1 and therefore !1 jx=0 = 1. Also, since !1 is identically zero outside element 1,
!1 jx=1 = 0. Thus Equation (2.9) for node 1 reduces to
du x=1 du
k dx !1 = ; k dx = flux entering node 1.
x=0 x=0
Similarly,
du x=1
k dx !n =0 (nodes 2 and 3)
x=0
and
du x=1 du
k dx !4 = k = flux entering node 4.
dx
x=0 x=1
Note: k has been left in these expressions to emphasise that they are heat fluxes.
Putting these global equations together we get
2 du 3
2 28 ; 53 0 0
3 2U 3 6; k dx 7
66;95318 289 +18289 ; 1853 1 6
0 7 6U2 77 = 66
7 6 0
x=0 7
77
4 0 ; 53 28 + 28 ; 18 5 4U3 5 66 0 77
53 (2.10)
U4 4 k du 5
18 9 9
0 0 ; 5318 28
dx x=1
9
or
Ku = f
K u f
where is the global “stiffness” matrix, the vector of unknowns and the global “load” vector.
Note that if the governing differential equation had included a distributed source term that was
independent of u, this term would appear - via its weighted integral - on the RHS of Equation (2.10)
rather than on the LHS as here. Moreover, if the source term was a function of x, the contribution
from each element would be different - as shown in the next section.
2.1.6 Boundary conditions

The boundary conditions u (0) = 0 and u (1) = 1 are applied directly to the first and last nodal
values: i.e., U1 = 0 and U4 = 1. These so-called essential boundary conditions then replace the
first and last rows in the global Equation (2.10), where the flux terms on the RHS are at present
unknown
1st equation U1 =0
2nd equation ; 53
18 U1 + 9 U2 ; 18 U3
56 53 =0
3rd equation ; 1853 U2 + 569 U3 ; 1853 U4 =0
4 equation
th
U4 =1
Note that, if a flux boundary condition had been applied, rather than an essential boundary
condition, the known value of flux would enter the appropriate RHS term and the value of U at
that node would remain an unknown in the system of equations. An applied boundary flux of zero,
corresponding to an insulated boundary, is termed a natural boundary condition, since effectively
no additional constraint is applied to the global equation. At least one essential boundary condition
must be applied.
2.1.7 Solution
Solving these equations gives: U 2 = 0:2885 and U3 = 0:6098. From Equation (2.2) the exact
solutions at these points are 0:2889 and 0:6102, respectively. The finite element solution is shown
in Figure 2.2.
2.1.8 Fluxes
The fluxes at nodes 1 and 4 are evaluated by substituting the nodal solutions U 1 = 0, U2 = 0:2885,
U3 = 0:6098 and U4 = 1 into Equation (2.10)
du
flux entering node 1 = ; k = ;0:8496 (k = 1; exact solution 0:8509)
dx
du x=0
flux entering node 4 = k
dx x=1 = 1:3157 (k = 1; exact solution 1:3131)
These fluxes are shown in Figure 2.2 as heat entering node 4 and leaving node 1, consistent with
heat flow down the temperature gradient.
1 1:3157
0:6098
0:2885
0:8496 x
0 1
3
2
3 1
F IGURE 2.2: Finite element solution of one-dimensional heat equation.
2.2 An x-Dependent Source Term

Consider the addition of a source term dependent on x in Equation (2.1):

d k du + u ; x = 0 0 < x < 1
; dx dx
Equation (2.6) now becomes
Z1 du d! du 1 Z1
k + u! dx = k ! + x! dx
dx dx dx (2.11)
0
0 0
where the x-dependent source term appears on the RHS because it is not dependent on u. Replacing
the domain integral for this source term by the sum of three element integrals
Z1 Z13 Z23 Z1
x! dx = x! dx + x! dx + x! dx
0 0 1 2
3 3
and putting x in terms of gives (with

dx = 1 for all three elements)
d 3
Z1 Z1 Z1 (1 + ) ! d + 1
Z1 (2 + ) ! d
x! dx = 13 1
3 ! d + 3 3 3 3 (2.12)
0 0 0 0
2.3 T HE G ALERKIN W EIGHT F UNCTION R EVISITED 31
where ! is chosen to be the appropriate basis function within each element. For example, the first
1
Z1
term on the RHS of (2.12) corresponding to element 1 is 9 'm d , where '1 = 1 ; and
0
'2 = . Evaluating these expressions,
Z1 1 (1 ; ) d = 1
9 54
0
and
Z1 1 2 d = 1
9 27
0
1
Thus, the contribution to the element 1 RHS vector from the source term is 54
1 .
27
Similarly, for element 2,
Z1 1 (1 + ) (1 ; ) d = 2 and Z1 1 (1 + ) d = 5 gives
2
27
5
9 27 9 54 54
0 0
and for element 3,

Z1 1 (2 + ) (1 ; ) d = 7 and Z1 1 (2 + ) d = 5 gives
7
54
5
9 54 9 54 54
0 0
Assembling these into the global RHS vector, Equation (2.10) becomes
2 du 3
2 28 ; 53 0 0
3 2U 3 6; k dx 7 2 1 3
66;91853 56918 ; 5318 6
0 77 66U2 77 = 66
1
0
x=0 7
77 + 66 271 +54 272 77
4 0 ; 53 56 ; 1853 5 4U3 5 66 0 77 4 545 + 547 5
U4 4 k du 5
18 9
0 0 ; 5318 28 5
dx x=1
9 54
2.3 The Galerkin Weight Function Revisited

A key idea in the Galerkin finite element method is the choice of weighting functions which are
orthogonal to the equation residual (thought of here as the error or amount by which the equation
fails to be exactly zero). This idea is illustrated in Figure 2.3.
u
In Figure 2.3a an exact vector e (lying in 3D space) is approximated by a vector = 1 1 u u'
'
where 1 is a basis vector along the first coordinate axis (representing one degree of freedom
u
in the system). The difference between the exact vector e and the approximate vector is the u
(a) (b) '3 (c)
u3
u
ue u2 '2 R
R
u
u
u1
'1
u = u1'1 u = u1'1 + u2'2 u = u1'1 + u2'2 + u3'3
r '1 = 0 : : : + r '2 = 0 : : : + r '3 = 0
F IGURE 2.3: Showing how the Galerkin method maintains orthogonality between the residual
vector r and the set of basis vectors 'i as i is increased from (a) 1 to (b) 2 to (c) 3.
r u u
error or residual = e ; (shown by the broken line in Figure 2.3a). The Galerkin technique
minimises this residual by making it orthogonal to ' 1 and hence to the approximating vector . If u
a second degree of freedom (in the form of another coordinate axis in Figure 2.3b) is added, the
u ' '
approximating vector is = u1 1 + u2 2 and the residual is now also made orthogonal to ' 2
u
and hence to . Finally, in Figure 2.3c, a third degree of freedom (a third axis in Figure 2.3c) is
u
permitted in the approximation = u 1 '1 + u2 '2 + u3 '3 with the result that the residual (now
u u
also orthogonal to '3 ) is reduced to zero and = e . For a 3D vector space we only need three
u
axes or basis vectors to represent the true vector , but in the infinite dimensional vector space
associated with a spatially continuous field u (x) we need to impose the equivalent orthogonality
Z
condition R'dx = 0 for every basis function ' used in the approximate representation of
u (x). The key point is that in this analogy the residual is made orthogonal to the current set of basis
vectors - or, equivalently, in finite element analysis, to the set of basis functions used to represent
the dependent variable. This ensures that the error or residual is minimal (in a least-squares sense)
for the current number of degrees of freedom and that as the number of degrees of freedom is
increased (or the mesh refined) the error decreases monotonically.
2.4 Two and Three-Dimensional Steady-State Heat Conduction

Extending Equation (2.1) to two or three spatial dimensions introduces some additional complexity
which we examine here. Consider the three-dimensional steady-state heat equation with no source
terms:

@ k @u ; @ k @u ; @ k @u = 0

; @x x
@x @y y @y @z z @z
2.4 T WO AND T HREE -D IMENSIONAL S TEADY-S TATE H EAT C ONDUCTION 33
where kx ; ky and kz are the thermal diffusivities along the x, y and z axes respectively. If the
material is assumed to be isotropic, k x = ky = kz = k , and the above equation can be written as
;r (kru) = 0 (2.13)
and, if k is spatially constant (in the case of a homogeneous material), this reduces to Laplace’s
equation k r2 u = 0. Here we consider the solution of Equation (2.13) over the region
, subject
to boundary conditions on ; (see Figure 2.4).
Solution region boundary: ;
Solution region:
F IGURE 2.4: The region

and the boundary ;.
The weighted integral equation, corresponding to Equation (2.13), is

Z
;r (kru) ! d
= 0 (2.14)

The multi-dimensional equivalent of integration by parts is the Green-Gauss theorem:

Z Z @g
(f r rg + rf rg) d
= f @n d; (2.15)

;
(see p553 in Advanced Engineering Mathematics” by E. Kreysig, 7th edition, Wiley, 1993).
This is used (with f = ! , g = ;ku and assuming that k is constant) to reduce the derivative
order from two to one as follows:
Z Z Z @u
;r (kru) ! d
= kru r! d
; k @n ! d; (2.16)

;
Z Z du d! x2
cf. Integration by parts is ; dx k dx ! dx = k dx dx dx ; k du
d du
dx ! .
x x x 1
Using Equation (2.16) in Equation (2.14) gives the two-dimensional equivalent of Equation (2.6)
(but with no source term):

Z Z @u
kru r! d
= k @n ! d; (2.17)

;
subject to u being given on one part of the boundary and

@u being given on another part of the
@n
boundary.
The integrand on the LHS of (2.17) is evaluated using
@u @! = @u @i @! @j
ru r! = @x (2.18)
k @xk @i @xk @j @xk
where
@i
u = 'nun and ! = 'm, as before, and the geometric terms @x are found from the
k
inverse matrix
@ @x ;1
i
@xk= k @i
or, for a two-dimensional element,
2 @1 @1 3 2 @x @x 3;1 2 @y @x 3
64 @x @y 75 = 64 @1 @2 75 = 1 64 @2 ; @2 75
@2 @2 @y @y @x @y ; @x @y @y @x
; @
@x @y @1 @2 @1 @2 @2 @1 1 @1
2.5 Basis Functions - Element Discretisation

[I
Let
=
i , i.e., the solution region is the union of the individual elements. In each
i let
i=1
u = 'nun = '1u1 + '2u2 + : : : + 'N uN and map each
i to the 1; 2 plane. Figure 2.5 shows an
example of this mapping.
2.5 BASIS F UNCTIONS - E LEMENT D ISCRETISATION 35
2 2
y 1 7 1 8
8 9

3
4
7 8 9 0 4 5 0 5 6
0 1 1 0 1 1

4

3
5 6
4 2 2

2 1 4 1 5

1 5 6
3
1
2
2 0 1 2
1 0 2 3
1
1 0 1 0 1
x
F IGURE 2.5: Mapping each
to the 1 ; 2 plane in a 2 2 element plane.
For each element, the basis functions and their derivatives are:
'1 = (1 ; 1)(1 ; 2 ) @'1 = ;(1 ; )

2 (2.19)
@1
@'1 = ;(1 ; )
1 (2.20)
@2
(2.21)
'2 = 1(1 ; 2) @'2 = 1 ;

2 (2.22)
@1
@'1 = ;
1 (2.23)
@2
(2.24)
'3 = (1 ; 1)2 @'3 = ;

2 (2.25)
@1
@'3 = 1 ;
1 (2.26)
@2
(2.27)
'4 = 12 @'4 =

@1 2 (2.28)
@'4 =
@2 1 (2.29)
2.6 Integration
The equation is
Z Z @u
kru r! d
= k @n ! d; (2.30)

;
i.e.,
Z @u @! @u @! Z @u
k @x @x + @y @y d
= k @n ! d; (2.31)

;
u has already been approximated by 'n un and ! is a weight function but what should this be
chosen to be? For a Galerkin formulation choose ! = 'm i.e., weight function is one of the basis
functions used to approximate the dependent variable.
This gives
X Z @'n @'m @'n @'m Z @u

un k @x @x + @y @y d
= k @n 'm d; (2.32)
i
;
where the stiffness matrix is Emn where m = 1; : : : ; 4 and n = 1; : : : ; 4 and Fm is the (element)
load vector.
The names originated from earlier finite element applications and extension of spring systems,
i.e., F = kx where k is the stiffness of spring and F is the force/load.
This yields the system of equations E mn un = Fm . e.g., heat flow in a unit square (see Fig-
ure 2.6).
y (2)
x (1)
0 1
F IGURE 2.6: Considering heat flow in a unit square.
2.7 A SSEMBLE G LOBAL E QUATIONS 37
The first component E11 is calculated as
Z1 Z1
E11 = k (1 ; y)2 + (1 ; x)2 dxdy
0 0
= 23 k
and similarly for the other components of the matrix.
Note that if the element was not the unit square we would need to transform from (x; y ) to
(1; 2) coordinates. In this case we would have to include the Jacobian of the transformation and
also use the chain rule to calculate
@'i . e.g., @'n = @'n @1 + @'n @2 = @'n @i .
@xj @x @1 @x @2 @x @i @x
The system of Emn un = Fm becomes
2 2 ; 1 ; 1 ; 1 3 2u 3
1
6
k 64;
3
1 2
6 6
1
3
7 6 7
3 ; 3 ; 6 7 6u2 7 = RHS
1
; 6 ; 3 3 ; 6 5 4u35
6 (Right Hand Side) (2.33)
1 1 2 1
;3 ;6 ;6
1 1 1 2
3 u4
Note that the Galerkin formulation generates a symmetric stiffness matrix (this is true for self
adjoint operators which are the most common).
Given that boundary conditions can be applied and it is possible to solve for unknown nodal
temperatures or fluxes. However, typically there is more than one element and so the next step is
required.
2.7 Assemble Global Equations

Each element stiffness matrix must be assembled into a global stiffness matrix. For example,
consider 4 elements (each of unit size) and nine nodes. Each element has the same element stiffness
matrix as that given above. This is because each element is the same size, shape and interpolation.
2 2 ;1 ; 16 ; 31
3 2u 3
1
66; 6 3 + 3 ; 6 ; 13
3
1 2
6
2 1 ; 6 ; 61
1 ;31 7
7 6
6u277
66 ; 16 2 ; 31 ;61 77 66u377
66; 16 ; 31 3 23 + 23 ; 61 ; 61 ;6 ;3
1 1 77 66u477
66; 13 ; 16 ; 61 ; 32 ; 16 ; 16 2 + 23 + 32 + 23 ; 6 ; 6 ; 3 ; 6 ; 6 ; 3 77 66u577 = RHS
1 1 1 1 1 1
66 ; 31 ; 16
3
; 16 ; 61 3+3
2 2 ; 13 ; 16 77 66u677
66 ; 16 ; 31 2 ; 16 77 66u777
4 ; 13 ; 6 ; 61
1 ; 3 ; 6 3 + 3 ; 6 5 4u85
1
3
1 2 2 1
; 31 ; 16 ; 16 2
3 u9
(2.34)
y
7 8 9
global node numbering

3 4
element numbering
4 5 6
1 2
1 2 3 x
F IGURE 2.7: Assembling 4 unit sized elements into a global stiffness matrix.
This yields the system of equations

22 ; ; ;
3 2u13
77 666uu2777
1 1 1
66;316 4
6 6 3
; 61 ; ; ; 31
1 1
66 ; 16
3
2
3 3
; ; 61
1 77 66u377
66; 16 ; 13
3
4 ;
3
1 ; 16 ; 13 77 66u477
66; 13 ; 13
3 3
; 13 ; 31 38 ; 31 ; 13 ; 13 ; 13 77 66u5677 = RHS
66 ; 13 ; 16 ; 31 43 ; 31 ; 16 77 66u677
66 ;6 ;3
1 1 2 ; 16 77 66u 77
4 ; 6 5 64u7875
3
; 13 ; 13 ; 31 ; 16 43 1
; 13 ; 61 ; 16 2
3 u9
Note that the matrix is symmetric. It should also be clear that the matrix will be sparse if there is a
larger number of elements.
From this system of equations, boundary conditions can be applied and the equations solved.
To solve, firstly boundary conditions are applied to reduce the size of the system.
If at global node i, ui is known, we can remove the ith equation and replace it with the known
value of ui . This is because the RHS at node i is known but the RHS equation is uncoupled from
other equations so the equation can be removed. Therefore the size of the system is reduced. The
final system to solve is only as big as the number of unknown values of u.
As an example to illustrate this consider fixing the temperature (u) at the left and right sides of
the plate in Figure 2.7 and insulating the top (node 8) and the bottom (node 2). This means that
2.8 G AUSSIAN Q UADRATURE 39
there are only 3 unknown values of u at nodes (2,5 and 8), therefore there is a 3 3 matrix to solve.
The RHS is known at these three nodes (see below). We can then solve the 3 3 matrix and then
multiply out the original matrix to find the unknown RHS values.
The RHS is 0 at nodes 2 and 8 because it is insulated. To find out what the RHS is at node 5
Z @u
we need to examine the RHS expression
@n ! d; = 0 at node 5. This is zero as flux is always
;
0 at internal nodes. This can be explained in two ways.

1
2
n n
F IGURE 2.8: “Cancelling” of flux in internal nodes.
Correct way: ; does not pass through node 5 and each basis function that is not zero at 5 is zero
on ;
Other way:
@u is opposite in neighbouring elements so it cancels (see Figure 2.8).
@n
2.8 Gaussian Quadrature
The element integrals arising from two- or three-dimensional problems can seldom be evaluated an-
alytically. Numerical integration or quadrature is therefore required and the most efficient scheme
for integrating the expressions that arise in the finite element method is Gauss-Legendre quadra-
ture.
Consider first the problem of integrating f ( ) between the limits 0 and 1 by the sum of
weighted samples of f ( ) taken at points 1 ; 2 ; : : : ; I (see Figure 2.3):
Z1 X
I
f ( ) d = Wif (i) + E
0 i=1
Here Wi are the weights associated with sample points i - called Gauss points - and E is the
error in the approximation of the integral. We now choose the Gauss points and weights to exactly
integrate a polynomial of degree 2I ; 1 (since a general polynomial of degree 2I ; 1 has 2I
arbitrary coefficients and there are 2I unknown Gauss points and weights).
For example, with I = 2 we can exactly integrate a polynomial of degree 3:
f ( )
....

0 1 2 .... I 1
F IGURE 2.9: Gaussian quadrature. f ( ) is sampled at I Gauss points 1 ; 2 : : : I :
Z1
Let f ( ) d = W1 f (1) + W2 f (2)
0
and choose f ( ) = a + b + c 2 + d 3 . Then
Z1 Z1 Z1 Z1 Z1
f ( ) d = a d + b d + c d + d
2 3 d (2.35)
0 0 0 0 0
Since a, b, c and d are arbitrary coefficients, each integral on the RHS of 2.35 must be integrated
exactly. Thus,
Z1
d = 1 = W1 :1 + W2:1 (2.36)
0
Z1
d = 12 = W1 :1 + W2:2 (2.37)
0
Z 1
2 d = 13 = W1 :12 + W2:22 (2.38)
0
Z1
3 d = 14 = W1 :13 + W2:23 (2.39)
0
These four equations yield the solution for the two Gauss points and weights as follows:
2.8 G AUSSIAN Q UADRATURE 41
From symmetry and Equation (2.36),
W1 = W2 = 12 :
Then, from (2.37),
2 = 1 ; 1
and, substituting in (2.38),
12 + (1 ; 1)2 = 23
212 ; 21 + 31 = 0;
giving
1 = 21 p1 :
2 3
Equation (2.39) is satisfied identically. Thus, the two Gauss points are given by
1 = 12 ; p1 ;
2 3
1
2 = 2 + p1 ; (2.40)
2 3
W1 = W2 = 12
A similar calculation for a 5th degree polynomial using three Gauss points gives
r
1 = 21 ; 21 35 ; 5
W1 = 18
2 = 12 ; W2 = 94 (2.41)
r
3 = 2 + 2 35 ;
1 1 5
W3 = 18
2 For two- or three-dimensional Gaussian quadrature the Gauss point positions are simply the
values given above along each i -coordinate with the weights scaled to sum to 1 e.g., for 2x2 Gauss
1
quadrature the 4 weights are all . The number of Gauss points chosen for each i -direction is
4
governed by the complexity of the integrand in the element integral (2.8). In general two- and three-
dimensional problems the integral is not polynomial (owing to the
@i terms which come from the
@xj
@x
i
inverse of the matrix ) and no attempt is made to achieve exact integration. The quadrature
@j
error must be balanced against the discretization error. For example, if the two-dimensional basis
is cubic in the 1 -direction and linear in the 2 -direction, three Gauss points would be used in the
1-direction and two in the 2-direction.
2.9 CMISS Examples

1. To solve for the steady state temperature distribution inside a plate run CMISS example 311
2. To solve for the steady state temperature distribution inside an annulus run CMISS example
312
3. To investigate the convergence of the steady state temperature distribution with mesh refine-
ment run CMISS examples 3141, 3142, 3143 and 3144.
Chapter 3
The Boundary Element Method
3.1 Introduction
Having developed the basic ideas behind the finite element method, we now develop the basic ideas
of the boundary element method. There are several key differences between these two methods,
one of which involves the choice of weighting function (recall the Galerkin finite element method
used as a weighting function one of the basis functions used to approximate the solution variable).
Before launching into the boundary element method we must briefly develop some ideas that are
central to the weighting function used in the boundary element method.
3.2 The Dirac-Delta Function and Fundamental Solutions

Before one applies the boundary element method to a particular problem one must obtain a funda-
mental solution (which is similar to the idea of a particular solution in ordinary differential equa-
tions and is the weighting function). Fundamental solutions are tied to the Dirac 1 Delta function
and we deal with both here.
3.2.1 Dirac-Delta function

What we do here is very non-rigorous. To gain an intuitive feel for this unusual function, consider
the following sequence of force distributions applied to a large plate as shown in Figure 3.1
n jxj < 1
wn (x) = 2 n1
0 jxj > n
1
Paul A.M. Dirac (1902-1994) was awarded the Nobel Prize (with Erwin Schrodinger) in 1933 for his work in
quantum mechanics. Dirac introduced the idea of the “Dirac Delta” intuitively, as we will do here, around 1926-27.
It was rigorously defined as a so-called generalised function by Schwartz in 1950-51, and strictly speaking we should
talk about the “Dirac Delta Distribution”.
44 T HE B OUNDARY E LEMENT M ETHOD
Each has the property that

Z1
wn (x) dx = 1 (i.e., the total force applied is unity)
;1
but as n increases the area of force application decreases and the force/unit area increases.
2 w4
3
2 w3
1
w2
1
2 w1
1 1 1 1
4 3 2
F IGURE 3.1: Illustrations of unit force distributions wn .
As n gets larger we can easily see that the area of application of the force becomes smaller
and smaller, the magnitude of the force increases but the total force applied remains unity. If we
imagine letting n ! 1 we obtain an idealised “point” force of unit strength, given the symbol
(x), acting at x = 0. Thus, in a nonrigorous sense we have
(x) = nlim w (x)
!1 n
the Dirac Delta“function”.
This is not a function that we are used to dealing with because we have (x) = 0 if x 6= 0
and “ (0) = 1” i.e., the “function” is zero everywhere except at the origin, where it is infinite.
Z1 Z1
However, we have (x) dx = 1 since each wn (x) dx = 1.
;1 ;1
The Dirac delta “function” is not a function in the usual sense, and it is more correctly referred
to as the Dirac delta distribution. It also has the property that for any continuous function h (x)
Z1
(x) h (x) dx = h (0) (3.1)
;1
3.2 T HE D IRAC -D ELTA F UNCTION AND F UNDAMENTAL S OLUTIONS 45
A rough proof of this is as follows

Z1 Z1
(x) h (x) dx = nlim
!1
wn (x) h (x) dx by definition of (x)
;1 ;1
n Zn
1
= nlim
!1 2
h (x) dx by definition of wn (x)
; n1
n h ( ) 2 1 1

= nlim by the Mean Value Theorem, where 2 ; ;
!1 2 n 1 1 n n
since 2 ; ;
= h (0) n n and as n ! 1; ! 0
The above result (Equation (3.1)) is often used as the defining property of the Dirac delta in
more rigorous derivations. One does not usually talk about the values of the Dirac delta at a
particular point, but rather its integral behaviour. Some properties of the Dirac delta are listed
below
Z1
( ; x) h (x) dx = h ( ) (3.2)
;1
(Note: ( ; x) is the Dirac delta distribution centred at x = instead of x = 0)
( ; x) = H 0 ( ; t) (3.3)
(
0 if <t
where H ( ; t) = (i.e., the Dirac Delta function is the slope of the Heaviside 2
1 if >t
step function.)
( ; x; ; y) = ( ; x) ( ; y) (3.4)
(i.e., the two dimensional Dirac delta is just a product of two one-dimensional Dirac deltas.)
3.2.2 Fundamental solutions

We develop here the fundamental solution (also called the freespace Green’s3 function) for Laplace’s
Equation in two variables. The fundamental solution of a particular equation is the weighting func-
tion that is used in the boundary element formulation of that equation. It is therefore important to
be able to find the fundamental solution for a particular equation. Most of the common equations
2
Oliver Heaviside (1850-1925) was a British physicist, who pioneered the mathematical study of electrical circuits
and helped develop vector analysis.
3
George Green (1793-1841) was a self-educated miller’s son. Most widely known for his integral theorem (the
Green-Gauss theorem).
have well-known fundamental solutions (see Appendix 3.16). We briefly illustrate here how to find
a simple fundamental solution.
Consider solving the Laplace Equation
@ 2 u + @ 2 u = 0 in some domain
2 <2.
@x2 @y2
The fundamental solution for this equation (analogous to a particular solution in ODE work) is
a solution of
@ 2 ! + @ 2 ! + ( ; x; ; y) = 0 (3.5)
@x2 @y2
in <2 (i.e., we solve the above without reference to the original domain
or original boundary
conditions). The method is to try and find solution to r 2 ! = 0 in <2 which contains a singularity
at the point (; ). This is not as difficult as it sounds. We expect the solution to be symmetric
about the point (; ) since ( ; x; ; y ) is symmetric about this point. So we adopt a local
polar coordinate system about the singular point (; ).
Let
q
r = ( ; x)2 + ( ; y)2
Then, from Section 1.8 we have

@ r @! + 1 @ 2 !

r ! = 1r @r
2
@r r2 @2 (3.6)
> 0; ( ; x; ; y) = 0 and owing to symmetry, @@!2 is zero. Thus Equation (3.6) becomes
2
For r

1 @ r @! = 0
r @r @r
This can be solved by straight (one-dimensional) integration. The solution is
! = A log r + B (3.7)
Note that this function is singular at r = 0 as required.

To find A and B we make use of the integral property of the Delta function. From Equa-
tion (3.5) we must have
Z Z
r ! dD = ; dD = ;1
2 (3.8)
D D
where D is any domain containing r = 0.
We choose a simple domain to allow us to evaluate the above integrals. If D is a small disk of
3.2 T HE D IRAC -D ELTA F UNCTION AND F UNDAMENTAL S OLUTIONS 47
"
D
(; )
F IGURE 3.2: Domain used to evaluate fundamental solution coefficients.
radius " > 0 centred at r = 0 (Figure 3.2) then from the Green-Gauss theorem
Z Z @!
r ! dD = @n dS
2 @D is the surface of the disk D
D @D
Z @!
= @r dS since D is a disk centred at r = 0 so n and r are in the same direction
@D
= A" 2" from Equation (3.7), and the fact that D is a disc of radius "
= 2A
Therefore, from Equation (3.8)
A = ; 21 :
So we have
! = ; 21 log r + B
B remains arbitrary but usually put equal to zero, so that the fundamental solution for the two-
dimensional Laplace Equation is

! = ; 21 log r = 21 log 1r (3.9)
q
where r = ( ; x)2 + ( ; y )2 (singular at the point (; )).
The fundamental solution for the three-dimensional Laplace Equation can be found by a similar
technique. The result is
1
! = 4r
where r is now a distance measured in three-dimensions.
3.3 The Two-Dimensional Boundary Element Method

We are now at a point where we can develop the boundary element method for the solution of
r2 u = 0 in a two-dimensional domain
. The basic steps are in fact quite similar to those used for
the finite element method (refer Section 2.1). We firstly must form an integral equation from the
Laplace Equation by using a weighted integral equation and then use the Green-Gauss theorem.
From Section 2.4 we have seen that
Z Z @u Z
0= r u:! d
= @n ! d; ; ru:r! d
2 (3.10)

@

This was the starting point for the finite element method. To derive the starting equation for
the boundary element method we use the Green-Gauss theorem again on the second integral. This
gives
Z @u Z
0 = @n ! d; ; ru:r! d
Z @u Z
@! Z (3.11)
= @n ! d; ; u @n d; + ur ! d
@
@

For the Galerkin FEM we chose ! , the weighting function, to be ' m , one of the basis functions
used to approximate u. For the boundary element method we choose ! to be the fundamental
solution of Laplace’s Equation derived in the previous section i.e.,
! = ; 21 log r
q
where r = ( ; x)2 + ( ; y )2 (singular at the point (; ) 2
).
Then from Equation (3.11), using the property of the Dirac delta
Z Z
ur ! d
= ;
2 u ( ; x; ; y) d
= ;u (; ) (; ) 2
(3.12)

i.e., the domain integral has been replaced by a point value.

3.3 T HE T WO -D IMENSIONAL B OUNDARY E LEMENT M ETHOD 49
Thus Equation (3.11) becomes

Z @! Z @u
u (; ) + u @n d; = @n ! d; (; ) 2
(3.13)
@
@
This equation contains only boundary integrals (and no domain integrals as in Finite Elements)
and is referred to as a boundary integral equation. It relates the value of u at some point inside
the solution domain to integral expressions involving u and
@u
@n over the boundary of the solution
domain. Rather than having an expression relating the value of u at some point inside the domain
to boundary integrals, a more useful expression would be one relating the value of u at some point
on the boundary to boundary integrals. We derive such an expression below.
The previous equation (Equation (3.13)) holds if (; ) 2
(i.e., the singularity of Dirac Delta
function is inside the domain). If (; ) is outside
then
Z Z
ur ! d
= ;
2 u ( ; x; ; y) d
= 0

since the integrand of the second integral is zero at every point except (; ) and this point is
outside the region of integration. The case which needs special consideration is when the singular
point (; ) is on the boundary of the domain
. This case also happens to be the most important
for numerical work as we shall see. The integral expression we will ultimately obtain is simply
Equation (3.13) with u (; ) replaced by 12 u (; ). We can see this in a non-rigorous way as
follows. When (; ) was inside the domain, we integrated around the entire singularity of the
Dirac Delta to get u (; ) in Equation (3.13). When (; ) is on the boundary we only have half of
the singularity contained inside the domain, so we integrate around one-half of the singularity to
get
1 u (; ). Rigorous details of where this coefficient 1 comes from are given below.
2 2 Z
Let P denote the point (; ) 2
. In order to be able to evaluate ur2! d
in this case we

enlarge
to include a disk of radius " about P (Figure 3.3). We call this enlarged region
0 and
let ;0 = ;;" [ ;" .
Now, since P is inside the enlarged region
0 , Equation (3.13) holds for this enlarged domain
i.e.,
Z Z @u
u (P ) + u @!
@n d; = @n ! d; (3.14)
;;" [;" ;;" [;"
We must now investigate this equation as lim"#0 . There are 4 integrals to consider, and we look at
each of these in turn.
;"

0
"
;;"
P

F IGURE 3.3: Illustration of enlarged domain when singular point is on the boundary.
Firstly consider
Z @! Z
@ ; 1 log r d;
u @n d; = u @n 2 by definition of !
;" ;"
Z
@ ; 1 log r d; @ @ on ;
= u @r 2
since
@n @r "
;"
1 Z u
= ; 2 r d;
1 1 Z;"
= ; 2 " u d; since r = " on ;"

;"
1 1
! ; 2 " u (P ) "
by the mean value theorem for a surface with a unique tangent at P .
Thus
Z @! 1 u (P ) u (P )
lim
"#0
u @n d; = lim
"#0
; 2 " " = ; 2 (3.15)
;"
By a similar process we obtain

Z @u 1 @u
lim
"#0
! @n d; = lim
"#0
; 2 @n (P ) " log " = 0 (3.16)
;"
since lim" log "#0 as lim"#0 .

3.3 T HE T WO -D IMENSIONAL B OUNDARY E LEMENT M ETHOD 51
It only remains to consider the integrand over ; ;" . For “nice” integrals (which includes the
integrals we are dealing with here) we have
0Z 1 Z
lim
"#0
@ (nice integrand) d;A = (nice integrand) d;
;;" ;
since ;;" ! ; as lim"#0 .

Note: If the integrand is too badly behaved we cannot always replace ;;" by ; in the limit and
one must deal with Cauchy Principal Values. (refer Section 4.8)
Thus we have
0Z 1 Z
lim @ @u ! d;A = @u ! d; (3.17)
"!0 @n @n
0Z;;"
1 Z;
lim @ @! u d;A = @! u d; (3.18)

"!0 @n @n
;;" ;
Combining Equations (3.14)–(3.18) we get

Z @! 1 Z @u
u (P ) + u @n d; = 2 u (P ) + @n ! d;
; ;
or
1 u (P ) + Z u @! d; = Z @u ! d;
2 @n @n
; ;
where P = (; ) 2 @
(i.e., singular point is on the boundary of the region).
Note: The above is true if the point P is at a smooth point (i.e., a point with a unique tangent) on
the boundary of
. If P happens to lie at some nonsmooth point e.g. a corner, then the coefficient
1 is replaced by where is the internal angle at P (Figure 3.4).
2 2
P

F IGURE 3.4: Illustration of internal angle .

Thus we get the boundary integral equation.

Z @! Z @u
c (P ) u (P ) + u @n d; = @n ! d; (3.19)
; ;
where
! = ; 21 log r
q
r = ( ; x)2 + ( ; y)2
8
>
< 1 if P 2
c (P ) = > internal2 angle if P 2 ; and ; smooth at P

1
: if P 2 ; and ; not smooth at P
2
For three-dimensional problems, the boundary integral equation expression above is the same,
with
1
! = 4r
q
r = ( ; x)2 + ( ; y)2 + ( ; z)2
8
>
< 1 if P 2
1
c (P ) = > inner solid
2 if P 2 ; and ; smooth at P
: angle
if P 2 ; and ; not smooth at P
4
Equation (3.19) involves only the surface distributions of
@u
u and @n and the value of u at a
point P . Once the surface distributions of u and
@u
@n are known, the value of u at any point P
inside
can be found since all surface integrals in Equation (3.19) are then known. The procedure
is thus to use Equation (3.19) to find the surface distributions of u and
@u
@n and then (if required)
use Equation (3.19) to find the solution at any point P 2
. Thus we solve for the boundary data
first, and find the volume data as a separate step.
Since Equation (3.19) only involves surface integrals, as opposed to volume integrals in a finite
element formulation, the overall size of the problem has been reduced by one dimension (from
volumes to surfaces). This can result in huge savings for problems with large volume to surface
ratios (i.e., problems with large domains). Also the effort required to produce a volume mesh of a
complex three-dimensional object is far greater than that required to produce a mesh of the surface.
Thus the boundary element method offers some distinct advantages over the finite element method
in certain situations. It also has some disadvantages when compared to the finite element method
and these will be discussed in Section 3.6. We now turn our attention to solving the boundary
integral equation given in Equation (3.19).
3.4 N UMERICAL S OLUTION P ROCEDURES FOR THE B OUNDARY I NTEGRAL E QUATION 53
3.4 Numerical Solution Procedures for the Boundary Integral

Equation
The first step is to discretise the surface ; into some set of elements (hence the name boundary
elements).
[N
;= ;j (3.20)
j =1
(a) (b)
F IGURE 3.5: Schematic illustration of a boundary element mesh (a) and a finite element mesh (b).
Then Equation (3.19) becomes

N Z
X N Z
X
c (P ) u (P ) + @!
u @n d; = @u ! d; (3.21)
j =1 ;j @n j =1 ;j
Over each element ;j we introduce standard (finite element) basis functions

X X
uj = ' uj and qj @uj
@n = ' qj (3.22)

where uj ; qj are values of u and q on element ;j and uj ; qj are values of u and q at node on
element ;j .
These basis functions for u and q can be any of the standard one-dimensional finite element
basis functions (although we are dealing with a two-dimensional problem, we only have to inter-
polate the functions over a one-dimensional element). In general the basis functions used for u and
q do not have to be the same (typically they are) and these basis functions can even be different to
the basis functions used for the geometry, but are generally taken to be the same (this is termed an
isoparametric formulation).
This gives
X
N X Z @! XN X Z
c (P ) u (P ) + uj ' @n d; = qj '! d; (3.23)
j =1 ;j j =1 ;j
This equation holds for any point P on the surface ;. We now generate one equation per node by
putting the point P to be at each node in turn. If P is at node i, say, then we have
X
N X Z @! i XN X Z
ciui + uj ' @n d; = qj '!i d; (3.24)
j =1 ;j j =1 ;j
where !i is the fundamental solution with the singularity at node i (recall ! is ;

1
2 log r , where
r is the distance from the singularity point). We can write Equation (3.24) in a more abbreviated
form as
X
N X X
N X
ciui + ujaij = qjbij (3.25)
j =1 j =1
where
Z Z
aij = ' @! i
@n d; and bij = '!i d; (3.26)
;j ;j
Equation (3.25) is for node i and if we have L nodes, then we can generate L equations.
We can assemble these equations into the matrix system
Au = Bq (3.27)
(compare to the global finite element equations Ku = f ) where the vectors u and q are the vectors
of nodal values of u and q . Note that the ij component of the A matrix in general is not a ij and
th
similarly for B .
At each node, we must specify either a value of u or q (or some combination of these) to have a
well-defined problem. We therefore have L equations (the number of nodes) and have L unknowns
to find. We need to rearrange the above system of equations to get
Cx = f (3.28)
x
where is the vector of unknowns. This can be solved using standard linear equation solvers,
although specialist solvers are required if the problem is large (refer [todo : Section ???]).
A B C
The matrices and (and hence ) are fully populated and not symmetric (compare to the
finite element formulation where the global stiffness matrix K
is sparse and symmetric). The
size of the A Band matrices are dependent on the number of surface nodes, while the matrix
K is dependent on the number of finite element nodes (which include nodes in the domain). As
3.5 N UMERICAL E VALUATION OF C OEFFICIENT I NTEGRALS 55
mentioned earlier, it depends on the surface to volume ratio as to which method will generate the
smallest and quickest solution.
The use of the fundamental solution as a weight function ensures that the and A
matrices B
A
are generally well conditioned (see Section 3.5 for more on this). In fact the matrix is diagonally
dominant (at least for Laplace’s equation). The matrix C
is therefore also well conditioned and
Equation (3.28) can be solved reasonably easily.
x u q
The vector contains the unknown values of and on the boundary. Once this has been
u q
found, all boundary values of and are known. If a solution is then required at a point inside the
domain, then we can use Equation (3.25) with the singular point P located at the required solution
point i.e.,
X
N X X
N X
u (P ) = qjbPj ; ujaPj (3.29)
j =1 j =1
The right hand side of Equation (3.29) contains no unknowns and only involves evaluating the
surface integrals using the fundamental solution with the singular point located at P .
3.5 Numerical Evaluation of Coefficient Integrals

We consider in detail here how one evaluates the aij and bij integrals for two-dimensional problems.
These integrals typically must be evaluated numerically, and require far more work and effort than
the analogous finite element integrals.
Recall that
Z Z
aij = ' @! i
@n d; and bij = '!i d;
;j ;j
where
!i = ; 21 log ri
ri = distance measured from node i
In terms of a local coordinate we have
Z1
bij = ' ( ) !i ( ) jJ ( )j d (3.30)
0
Z Z 1
aij = ' ( ) @n jJ ( )j d = ' ( ) @!
@!i ( ) i
( ) dri jJ ( )j d
@ri dn (3.31)
;j 0
The Jacobian J ( ) can be found by

s 2 2
J ( ) = d; = ds = dx + dy
d d d d (3.32)
where s represents the arclength and dx d and

dy
d can be found by straight differentiation of the
interpolation expression for x ( ) and y ( ).
The fundamental solution is
!i = ; 21 log (ri ( ))

q
ri ( ) = (x ( ) ; xi )2 + (y ( ) ; yi)2
where (xi ; yi ) are the coordinates of node i.
To find
dri we note that
dn
dri = rr n^
i (3.33)
dn
n
where ^ is a unit outward normal vector. To find a unit normal vector, we simply rotate the tangent

vector (given by (x0 ( ) ; y 0 ( )) ) by in the appropriate direction and then normalise.
2
Thus every expression in the integrands of the a ij and bij integrals can be found at any value of
, and the integrals can therefore be evaluated numerically using some suitable quadrature schemes.
If node i is well removed from element ;j then standard Gaussian quadrature can be used to
evaluate these integrals. However, if node i is in ;j (or close to it) we see that ri approaches 0
and the fundamental solution ! i tends to 1. The integral still exists, but the integrand becomes
singular. In such cases special care must be taken - either by using special quadrature schemes,
large numbers of Gauss points or other special treatment.
The integrals for which node i lies in element ; j are in general the largest in magnitude and
lead to the diagonally dominant matrix equation. It is therefore important to ensure that these
integrals are calculated as accurately as possible since these terms will have most influence on the
solution. This is one of the disadvantages of the BEM - the fact that singular integrands must be
accurately integrated.
A relatively straightforward way to evaluate all the integrals is simply to use Gaussian quadra-
ture with varying number of quadrature points, depending on how close or far the singular point is
from the current element. This is not very elegant or efficient, but has the benefit that it is relatively
easy to implement. For the case when node i is contained in the current element one can use special
quadrature schemes which are designed to integrate log-type functions. These are to be preferred
when one is dealing with Laplace’s equation. However, these special log-type schemes cannot be
so readily used on other types of fundamental solution so for a general purpose implementation,
Gaussian quadrature is still the norm. It is possible to incorporate adaptive integration schemes
that keep adding more quadrature points until some error estimate is small enough, or also to sub-
divide the current element into two or more smaller elements and evaluate the integral over each
3.6 T HE T HREE -D IMENSIONAL B OUNDARY E LEMENT M ETHOD 57
;j node i ;j
ri
ri
node i
(a) (b)
F IGURE 3.6: Illustration of the decrease in ri as node i approaches element ;j .
subelement. It is also possible to evaluate the “worst” integrals by using simple solutions to the
governing equation, and this technique is the norm for elasticity problems (Section 4.8). Details
on each of these methods is given in Section 3.8. It should be noted that research still continues in
an attempt to find more efficient ways of evaluating the boundary element integrals.
3.6 The Three-Dimensional Boundary Element Method

The three-dimensional boundary element method is very similar to the two-dimensional bound-
ary element method discussed above. As noted above, the three-dimensional boundary integral
equation is the same as the two-dimensional equation (3.19), with ! and c (P ) being defined as
in Section 3.3. The numerical solution procedure also parallels that given in Section 3.4, and the
expressions given for aij and bij apply equally well to the three-dimensional case. The only real
difference between the two procedures is how to numerically evaluate the terms in each integrand
of these coefficient integrals.
As in Section 3.5 we illustrate how to evaluate each of the terms in the integrand of a ij and bij .
The relevant expressions are

Z
aij = ' @! i
@n d;
;j
Z Z1
1
= ' (1; 2) @!
@r
i
( 1 ; 2 )
dri jJ ( ; )j d d
dn 1 2 1 2 (3.34)
i
Z0 0
bij = ' !i d;
;j
Z Z1
1
= ' (1; 2) !i (1; 2) jJ (1 ; 2)j d1 d2 (3.35)
0 0
The fundamental solution is
!i (1; 2) = 4r (1 ; )

qi 1 2
where ri (1 ; 2 ) = (x (1 ; 2 ) ; xi )2 + (y (1 ; 2 ) ; yi )2 + (z (1 ; 2 ) ; zi )2
where (xi ; yi ; zi ) are the coordinates of node i. As before we use

dri = rr ^ to find dri . n
i
dn dn
@x @y @z n @x @y @z
The unit outward normal ^ is found by normalising the cross product of the two tangent vectors
t1 = @1 ; @1 ; @1 and t2 = @ ; @ ; @ (it relies on the user of any BEM code to
2 2 2
ensure that the elements have been defined with a consistent set of element coordinates 1 and 2 ).
The Jacobian J (1 ; 2 ) is given by kt1 t2 k (where t1 and t2 are the two tangent vectors).
Note that this is different for the determinant in a two-dimensional finite element code - in that
case we are dealing with a two-dimensional surface in two-dimensional space, whereas here we
have a (possibly curved) two-dimensional surface in three-dimensional space.
The integrals are evaluated numerically using some suitable quadrature schemes (see Sec-
tion 3.8) (typically a Gauss-type scheme in both the 1 and 2 directions).
3.7 A Comparison of the FE and BE Methods

We comment here on some of the major differences between the two methods. Depending on the
application some of these differences can either be considered as advantageous or disadvantageous
to a particular scheme.
1. FEM: An entire domain mesh is required.

BEM: A mesh of the boundary only is required.
Comment: Because of the reduction in size of the mesh, one often hears of people saying
that the problem size has been reduced by one dimension. This is one of the major pluses of
3.7 A C OMPARISON OF THE FE AND BE M ETHODS 59
the BEM - construction of meshes for complicated objects, particularly in 3D, is a very time
consuming exercise.
2. FEM: Entire domain solution is calculated as part of the solution.

BEM: Solution on the boundary is calculated first, and then the solution at domain points (if
required) are found as a separate step.
Comment: There are many problems where the details of interest occur on the boundary, or
are localised to a particular part of the domain, and hence an entire domain solution is not
required.
3. FEM: Reactions on the boundary typically less accurate than the dependent variables.
u q
BEM: Both and of the same accuracy.
4. FEM: Differential Equation is being approximated.

BEM: Only boundary conditions are being approximated.
Comment: The use of the Green-Gauss theorem and a fundamental solution in the formu-
lation means that the BEM involves no approximations of the differential Equation in the
domain - only in its approximations of the boundary conditions.
5. FEM: Sparse symmetric matrix generated.

BEM: Fully populated nonsymmetric matrices generated.
Comment: The matrices are generally of different sizes due to the differences in size of
the domain mesh compared to the surface mesh. There are problems where either method
can give rise to the smaller system and quickest solution - it depends partly on the volume
to surface ratio. For problems involving infinite or semi-infinite domains, BEM is to be
favoured.
6. FEM: Element integrals easy to evaluate.

BEM: Integrals are more difficult to evaluate, and some contain integrands that become
singular.
Comment: BEM integrals are far harder to evaluate. Also the integrals that are the most
difficult (those containing singular integrands) have a significant effect on the accuracy of
the solution, so these integrals need to be evaluated accurately.
7. FEM: Widely applicable. Handles nonlinear problems well.

BEM: Cannot even handle all linear problems.
Comment: A fundamental solution must be found (or at least an approximate one) before the
BEM can be applied. There are many linear problems (e.g., virtually any nonhomogeneous
equation) for which fundamental solutions are not known. There are certain areas in which
the BEM is clearly superior, but it can be rather restrictive in its applicability.
8. FEM: Relatively easy to implement.

BEM: Much more difficult to implement.
Comment: The need to evaluate integrals involving singular integrands makes the BEM at
least an order of magnitude more difficult to implement than a corresponding finite element
procedure.
3.8 More on Numerical Integration

The BEM involves integrals whose integrands in generally become singular when the source point
is contained in the element of integration. If one uses constant or linear interpolation for the
geometry and dependent variable, then it is possible to obtain analytic expressions to most (if not
all) of the integrals that will appear in the BEM (at least for two-dimensional problems). The
expressions can become quite lengthy to write down and evaluate, but benefit from the fact that
they will be exact. However, when one begins to use general curved elements and/or solve three-
dimensional problems then the integrals will not be available as analytic expressions. The basic
tool for evaluation of these integrals is quadrature. As discussed in Section 2.8 a one-dimensional
integral is approximated by a sum in which the integrand is evaluated at certain discrete points or
abscissa
Z1 X
N
f ( ) d f (i) wi
0 i=1
where wi are the weights and i are the abscissa.

The weights and abscissa for the Gaussian quadrature scheme of order N are chosen so that the
above expression will exactly integrate any polynomial of degree 2N ; 1 or less. For the numerical
evaluation of two or three-dimensional integrals, a Gaussian scheme can be used of each variable
of integration if the region of integration is rectangular. This is generally not the optimal choice
for the weights and abscissae but it allows easy extension to higher order integration.
3.8.1 Logarithmic quadrature and other special schemes

Low order Gaussian schemes are generally sufficient for all FEM integrals, but that is not the
case for BEM. For a two-dimensional BEM solution of Laplace’s equation, integrals of the form
Z1
log ( ) f ( ) d will be required. It is relatively common to use logarithmic schemes for this.
0
These are obtained by approximating the integral as
Z1 X
N
log ( ) f ( ) d f (i) wi
0 i=1
i.e., the log function has been factored out.

In the same way as Gaussian quadrature schemes were developed in Section 2.8, log quadrature
schemes can be developed which will exactly integrate polynomial functions f ( ). Tables of these
are given below
It is possible to develop similar quadrature schemes for use in the BEM solution of other PDEs,
which use different fundamental solutions to the log function. The problem with this approach is
the lack of generality - each new equation to be used requires its own special quadrature scheme.
3.9 T HE B OUNDARY E LEMENT M ETHOD A PPLIED TO OTHER E LLIPTIC PDE S 61
Abscissas = ri Weight Factors = wi

n i ;wi n i ;wi n i ;wi
2 0.112009 0.718539 3 0.063891 0.513405 4 0.041448 0.383464
0.602277 0.281461 0.368997 0.391980 0.245275 0.386875
0.766880 0.094615 0.556165 0.190435
0.848982 0.039225
TABLE 3.1: Abscissas and weight factors for Gaussian integration for integrands with a
logarithmic singularity.
3.8.2 Special solutions

Another approach, particularly useful if Cauchy principal values are to be found (see Section 4.8) is
to use special solutions of the governing equation to find one or more of the more difficult integrals.
For example u = x is a solution to Laplaces’ equation (assuming the boundary conditions
are set correctly). Thus if one sets both u and q in Equation (3.27) at every node according to
the solution u = x, one can then use this to solve for some entry in either the A B or matrix
(typically the diagonal entry since this is the most important and difficult to find). Further solutions
to Laplaces equation (e.g., u = x2 ; y 2 ) can be used to find the other matrix entries (or just used
to check the accuracy of the matrices).
3.9 The Boundary Element Method Applied to other Elliptic

PDEs
Helmholtz, modified Helmholtz (CMISS example) Poisson Equation (domain integral and MRM,
DRM, Monte-carlo integration.
3.10 Solution of Matrix Equations

The standard BEM approach results in a system of equations of the form Cx f
= (refer (3.28)).
C
As mentioned above the matrix is generally well conditioned, fully populated and nonsymmet-
ric. For small problems, direct solution methods, based on LU factorisations, can be used. As the
problem size increases, the time taken for the matrix solution begins to dominate the matrix assem-
bly stage. This usually occurs when there is between 500 and 1000 degrees of freedom, although it
is very dependent on the implementation of the BE method. The current technique of favour in the
BE community for solution of large BEM matrix equations is a preconditioned Conjugate Gradient
solver. Preconditioners are generally problem dependent - what works well for one problem may
not be so good for another problem. The conjugate gradient technique is generally regarded as a
solution technique for (sparse) symmetric matrix equations.

B

f
;B
;f
;I
F IGURE 3.7: Coupled finite element/boundary element solution domain.
3.11 Coupling the FE and BE techniques

There are undoubtably situations which favour FEM over BEM and vice versa. Often one problem
can give rise to a model favouring one method in one region and the other method in another
region, e.g., in a detailed analysis of stresses around a foundation one needs FEM close to the
foundation to handle nonlinearities, but to handle the semi-infinite domain (well removed from the
foundation), BEM is better. There has been a lot of research on coupling FE and BE procedures -
we will only talk about the basic ideas and use Laplace’s Equation to illustrate this. There are at
least two possible methods.
1. Treat the BEM region as a finite element and combine with FEM
2. Treat the FEM region as an equivalent boundary element and combine with BEM
Note that these are essentially equivalent - the use of one or the other depends on the problem,
in the sense of which part is more dominant FEM or BEM)
Consider the region shown in Figure 3.7, where

f = FEM region

B = BEM region
;f = FEM boundary
;B = BEM boundary
;I = interface boundary
The BEM matrices for
B can be written as
Au = Bq (3.36)
u
where is a vector of the nodal values of u and
@u
q is a vector of the nodal values of @n
The FEM matrices for
F can be written as
Ku = f (3.37)
3.11 C OUPLING THE FE AND BE TECHNIQUES 63
where K f
is the stiffness matrix and is the load vector.
To apply method 1 (i.e., treating BEM as an equivalent FEM region) we get (from Equa-
tion (3.36))
B;1Au = q (3.38)
f
If we recall what the elements of in Equation (3.37) contained, then we can convert in q
q
Equation (3.38) to an equivalent load vector by weighting the nodal values of by the appropriate
basis functions, producing a matrix M f Mq
i.e., B =
Therefore Equation (3.38) becomes
M ;B;1A u = Mq = f B
i.e.,
KBu = f B
where K B = MB ;1 A
an equivalent stiffness matrix obtain from BEM.

Therefore we can assemble this together with original FEM matrix to produce an FEM-type
system for the entire region
B .
Notes:
1. K B is in general not symmetric and not sparse. This means that different matrix equation
solvers must be used for solving the new combined FEM-type system (most solvers in FEM
codes assume sparse and symmetric). Attempts have been made to “symmetricise” the K B
matrix - of doubtful quality. (e.g., replace K B by
1 ;K ; K T - often yields inaccurate
2 B B
results).
2. On ;I nodal values of u and q are unknown. One must make use of the following
uIB = uIF (u is continuous)
@ uIB = ; @ uIF (q is continuous, but ; = ;; )
@ nB @ nF B F
To apply method 2 (i.e., to treat the FEM region as an equivalent BEM region) we firstly note
that, as before, = f Mq. Applying this to (3.37) yields =Ku Mq an equivalent BEM system.
This can be assembled into the existing BEM system (using compatability conditions) and use
existing BEM matrix solvers.
Notes:
1. This approach does not require any matrix inversion and is hence easier (cheaper) to imple-
ment
2. Existing BEM solvers will not assume symmetric or sparse matrices therefore no new matrix
solvers to be implemented
3.12 Other BEM techniques

What we have mentioned to date is the so-called singular (direct) BEM. Given a BIE there are
other ways of solving the Equation although these are not so widely used.
3.12.1 Trefftz method

Trefftz was the first person to perform a BEM calculation (in 1917 - calculated the value (numer-
ical) of the contraction coefficient of a round jet issuing from an infinite tank - a nonlinear free
surface problem). This method basically uses a “complete” set of solutions instead of a Funda-
mental Solution. e.g., Consider Laplaces Equation in a (bounded) domain
Z @u Z @!
weighted residuals ) ! @n d; = u @n d; if r2 ! =0
@
@
The procedure is to express u as a series of (complete) functions satisfying Laplace’s equation

with coefficients which need to be numerically determined through utilisation of the boundary
conditions.
Notes:
1. Doesn’t introduce singular functions so integrals are easy to evaluate
2. Must find a (complete) set of functions (If you just use usual approximations for u matrix
system is not diagonally dominant so not so good)
3. Method is not so popular - Green’s functions more widely available that complete systems
3.12.2 Regular BEM

Consider the BIE for Laplace’s equation
Z @! Z @u
c (P ) u (P ) + u @n d; = @n ! d;
@
@
with ! = ; 21 log r

The usual procedure is to put point P at each solution variable node - creating an equation for each
node. This leads to singular integrands.
Another possibility is to put point P outside of the domain
- this yields
Z @!p Z @u
u @n d; = @n d;
p
@
@
3.13 S YMMETRY 65
Following discretisation as before gives
X
N X Z @! p XN X Z
uj ' @n d; = qj '!p d;
j =1 ;j j =1 ;j
- an equation involving u and q at each surface node.

By placing the point P (the singular point) at other distinct points outside
one can generate
as many equations as there are unknowns (or more if required).
Notes:
1. This method does not involve singular integrands, so that integrals are inexpensive to calcu-
late.
2. There is considerable choice for the location of the point P . Often the set of Equations
generated are ill-conditioned unless P chosen carefully. In practise P is chosen along the unit
outward normal of the surface at each solution variable node. The distance along each node
is often found by experimentation - various research papers suggesting “ideal” distances
(Patterson & Shiekh).
3. This method is not very popular.
4. The idea of placing the singularity point P away from the solution variable node is often of
use in other situations e.g., Exterior Acoustic Problems. For an acoustic problem (governed
by Helmholtz Equation r 2 u + k 2 u = 0) in an unbounded region the system of Equations pro-
duced by the usual (singular) BEM approach is singular for certain “fictitious” frequencies
(i.e., certain values of k ). To overcome this further equations are generated (by placing the
singular point P at various locations outside
). The system of equations are then overde-
termined and are solved in a least squares sense.
3.13 Symmetry
Consider the problem given in Figure 3.8 (the domain is outside the circle). Both the boundary
conditions and the governing Equation exhibit symmetry about the vertical axis. i.e., putting x to
;x makes no difference to the problem formulation. Thus the solution H (x; z) has the property
that H (x; z ) = H (;x; z ) 8x. This behaviour can be found in many problems and we can make
use of this as follows. The Boundary Element Equation is (with N = 2M (i.e., N is even) constant
elements)
1u + X
N Z @!i XN Z
2 i j=1 uj @n d; = j=1 qj !i d; i = 1; : : : ; N (3.39)
;j ;j
z
H = e;sz
r2H = s2H
x
F IGURE 3.8: A problem exhibiting symmetry.
We have N Equations and N unknowns (allowing for the boundary conditions). From symmetry
we know that (refer to Figure 3.9).
ui = un+1;i i = 1; : : : ; M (3.40)
So we can write
8 9 8 9
1u + X u
M >
< Z @!i d; + Z >
=
@!i d; = X q
M >
<Z Z >
=
2 i j=1 j > ! d; + ! d;
:;j @n ;N+1;j @n > ; j=1 j >
:;j i ;N+1;j i > ;
(3.41)
for nodes i = 1; : : : ; M . (The Equations for nodes i = M +1; : : : ; N are the same as the Equations
for nodes i = 1; : : : ; M ). The above M Equations have only M unknowns.
If we define
Z @!i Z @!i
aij = @n d; + @n d; (3.42)
;j ;N +1;j
Z Z
bij = !i d; + !i d; (3.43)
;j ;N +1;j
then we can write Equation (3.41) as
1u + X
M XM
2 i j=1 aij uj = j=1 bij qj i = 1; : : : ; M (3.44)
and solve as before. (This procedure has halved the number of unknowns.)
3.14 A XISYMMETRIC P ROBLEMS 67
;N +1;i ;i
;N ;1
F IGURE 3.9: Illustration of a symmetric mesh.
Note: Since i = 1; : : : ; M this means that the integrals over the elements ; M +1 to ;N will never
contain a singularity arising from the fundamental solution, except possibly on the axis of symme-
try if linear or higher order elements are used.
An alternative approach to the method above arises from the implied no flux across the z axis.
This approach ignores the negative x axis and considers the half plane problem shown.
However now the surface to be discretised extends to infinity in the positive and negative z
directions and the resulting systems of equations produced is much larger.
Further examples of how symmetry can be used (e.g., radial symmetry) are given in the next
section.
3.14 Axisymmetric Problems

If a three-dimensional problem exhibits radial or axial symmetry (i.e., u (r; 1 ; z ) = u (r; 2 ; z )) it
is possible to reduce the two-dimensional integrals appearing in the standard boundary Equation
to one-dimensional line integrals and thus substantially reduce the amount of computer time that
would otherwise be required to solve the fully three-dimensional problem. The first step in such a
procedure is to write the standard boundary integral equation in terms of cylindrical polars (r; ; z )
i.e.,
Z 0Z2 @!p 1 Z 0Z2 1
c (P ) u (P ) + u @ @n dq A rq d; = q @ !p dq A rq d; (3.45)
; 0 ; 0
where (rp ; p ; zp ) and (rq ; q ; zq ) are the polar coordinates of P and Q respectively, and ; is the
intersection of ; and = 0 semi-plane (Refer Figure 3.10). (n.b. Q is a point on the surface being
integrated over.)
z ;
F IGURE 3.10: Illustration of surface ; for an axisymmetric problem.
Q
r
r2
P
r1
rp
rq r
F IGURE 3.11: The distance from the source point (P ) to the point of interest (Q) in terms of
cylindrical polar coordinates.
For three-dimensional problems governed by Laplace’s equation

1
!p = 4r
where rp is the distance from P to Q. From Figure 3.11
r12 = rp2 + rq2 ; 2rprq cos (p ; q )

q2 2
r = r1 + r2
2
q
r = rp2 + rq2 ; 2rprq cos (p ; q ) + (zp ; zq )2
q
= a ; b cos (p ; q ) + (zp ; zq )2
(3.46)
3.15 I NFINITE R EGIONS 69
We define
Z 2
!p = 4 !p dq Kp (m)
1 where m =
2b
a+b (3.47)
0
a+b
and K (m) is the complete elliptic integral of the first kind.
!p is called the axisymmetric fundamental solution and is the Green’s function for a ring source
as opposed to a point source. i.e., ! p is a solution of
r2! + (r ; rp) = 0 (3.48)
instead of
r2 ! + p = 0 (3.49)
where p is the dirac delta centered at the point P and (r ; rP ) is the dirac delta centered on the
ring r = rp .
Unlike the two- and three-dimensional cases, the axisymmetric fundamental solution cannot be
written as simply a function of the distance between two points P and Q, but it also depends upon
the distance of these points to the axis of revolution.
We also define
Z @!p
2
qp = 4 @n dq @!
1
@n
p
(3.50)
0
For Laplace’s equation Equation (3.50) becomes
1
1 r2 ; r2 + (zp ; zq ) z ; z

p p q p q
q =
p a;b E (m) ; K (m) nr (Q) + a ; b E (m) nz (Q)
a + b 2rq
(3.51)
where E (m) is the complete elliptic integral of the second kind.

Using Equation (3.47) and Equation (3.50) we can write Equation (3.45) as
Z @!p Z
c (P ) u (P ) + u @n d; = q!p d; (3.52)
; ;
and the solution procedure for this Equation follows the same lines as the solution procedure given
previously for the two-dimensional version of boundary element method.
3.15 Infinite Regions

The boundary integral equations we have been using have been derived assuming the domain
is bounded (although this was never stated). However all concepts presented thus far are also
;
R

R
x0
;
F IGURE 3.12: Derivation of infinite domain boundary integral equations.
valid for infinite regular (i.e., nice) regions provided the solution and its normal derivative behave
appropriately as ; ! 1.
Consider the problem of solving r 2 u = 0 outside some surface ;.
; is the centre of a circle (or sphere in three dimensions) of radius centred at some point x 0 on
; and surrounding ; (see Figure 3.12). The boundary integral equations for the bounded domain
R can be written as
Z @!P Z @!P Z Z
c (P ) u (P ) + u @n d; + u @n d; = q!P d; + q!P d; (3.53)
; ; ; ;
If we let the radius R ! 1 Equation (3.53) will only be valid for the points on ; if
Z @!P
lim
R!1
u @n ; q!P d; = 0 (3.54)
;
If this is satisfied, the boundary integral Equation for

will be as expected i.e.,
Z @!P Z
c (P ) u (P ) + u @n d; = q!P d; (3.55)
; ;
3.15 I NFINITE R EGIONS 71
For three-dimensional problems with !

1
= 4r
d; = jJ j dd
;
jJ j = O R2
;
! = O R;1
where
@! = O ;R;2
@n
where jJ j is the Jacobian and O () represents the asymptotic behaviour of the function as R !
1. In this case Equation (3.53) will be satisfied if u behaves at most as O (R ;1) so that q =
O (R;2). These are the regularity conditions at infinity and these ensure that each term in the
integral Equation (3.53) behaves at most as O (R ;1 ) (i.e., each term will ! 0 as R ! 1)
For two-dimensional problems with ! = O (log (R)) we require u to behave as log (R) so that
q = O (R;1 ). For almost all well posed infinite domain problems the solution behaves appropri-
ately at infinity.
3.16 Appendix: Common Fundamental Solutions

3.16.1 Two-Dimensional equations
Here r
p
= (x21 + x22 ).
@ 2 u + @ 2 u + = 0
Laplace Equation
@x1 2 @x2 2 0
Solution u = 21 log 1r
@ 2 u + @ 2 u + 2u + = 0
Helmholtz Equation 0
@x1 2 @x2 2
Solution u = 41i H0(2) (r)
where H is the Hankel funtion.
@2u @2u @2 u
Wave Equation c @x 2 + @x 2 ; @t2 + 0 (t) = 0
2
1 2
where c is the wave speed.
Solution u = ; 2cH((cct2t2;;r)r2)
Diffusion Equation
@ 2 u + @ 2 u ; 1 @u = 0
@x1 2 @x2 2 k @t
where k is the diffusivity.
r2
u =; 1
3 exp ; 4kt
Solution
(4kt) 2
@jk
Navier’s Equation
@xj + l = 0 for a point load in direction l.
Solution pi = pjiej
pji = ; 8 (1 ;1 2 ) r2
@r
@n [(1 ; 2 ) ij + 3r;ir;j ] + (1 ; 2 ) (nj r;i ; nir;j ) ej
for a traction in direction k where is Poisson’s ratio.
3.16.2 Three-Dimensional equations

Here r
p
= (x21 + x22 + x23 ).
@ 2 u + @ 2 u + @ 2 u + = 0
Laplace Equation
@x1 2 @x2 2 @x3 2 0
Solution u = 4r1
@ 2 u + @ 2 u + @ 2 u + 2 u + = 0
Helmholtz Equation 0
@x1 2 @x2 2 @x3 2
Solution u = 4r1 exp (;ir)
@2u @2u @2 u @2u
Wave Equation c @x 2 + @x 2 + @x 2 ; @t2 + t = 0
2
1 2 3
where cis the wave
r
speed.
t; c
Solution u = 4r
@jk
Navier’s Equation
@xj + l = 0 for a isotropic homogenenous Kelvin
solution for a point load in direction l.
Solution uk = ulk el 3 ; 4
1
ulk = 16G (1 ; ) @r @r
r lk + @x1 @x2
for a displacement in direction k where is Poisson’s
ratio and G is the shear modulus.
3.16.3 Axisymmetric problems

Laplace For u see Equation (3.47) and for q see Equation (3.51)
3.17 CMISS Examples

1. 2D steady-state heat conduction inside an annulus To determine the steady-state heat con-
duction inside an annulus run the CMISS example 324.
2. 3D steady-state heat conduction inside a sphere. To determine the steady-state heat conduc-
tion inside a sphere run the CMISS example 328.
3. CMISS comparison of 2-D FEM and BEM calculations To determine the CMISS comparison
of 2-D FEM and BEM calculations run examples 324 and 312.
4. CMISS biopotential problems C4 and C5.

Chapter 4
Linear Elasticity
4.1 Introduction
To analyse the stress in various elastic bodies we calculate the strain energy of the body in terms of
nodal displacements and then minimize the strain energy with respect to these parameters - a tech-
nique known as the Rayleigh-Ritz. In fact, as we will show later, this leads to the same algebraic
equations as would be obtained by the Galerkin method (now equivalent to virtual work) but the
physical assumptions made (in neglecting certain strain energy terms) are exposed more clearly in
the Rayleigh-Ritz method. We will first consider one-dimensional truss and beam elements, then
two-dimensional plane stress and plane strain elements, and finally three-dimensional elasticity.
In all cases the steps are:
1. Evaluate the components of strain in terms of nodal displacements,
2. Evaluate the components of stress from strain using the elastic material constants,
3. Evaluate the strain energy for each element by integrating the products of stress and strain
components over the element volume,
4. Evaluate the potential energy from the sum of total strain energy for all elements together
with the work done by applied boundary forces,
5. Apply the boundary conditions, e.g., by fixing nodal displacements,
6. Minimize the potential energy with respect to the unconstrained nodal displacements,
7. Solve the resulting system of equations for the unconstrained nodal displacements,
8. Evaluate the stresses and strains using the nodal displacements and element basis functions,
9. Evaluate the boundary reaction forces (or moments) at the nodes where displacement is
constrained.
76 L INEAR E LASTICITY
4.2 Truss Elements

Consider the one-dimensional truss of undeformed length L in Figure 3.1 with end points (0; 0)
and (x; y ) and making an angle with the x-axis. Under the action of forces in the x- and y -
directions the right hand end of the truss displaces by u in the x-direction and v in the y -direction,
relative to the left hand end.
(X + u; Y + v)
(X; Y ) u v
L
l

x
F IGURE 4.1: A truss of initial length L is stretched to a new length l. Displacements of the right
hand end relative to the left hand end are u and v in the x- and y - directions, respectively.
The new length is l with axial strain

q
l (X + u)2 + (Y + v)2
e= L ;1= p 2 2 ;1
X + Y
pL2 + 2 (Xu + Y v) + u2 + v2
= ;1
r L 2 2
= 1 + 2 cos : Lu + sin : Lv + u L+2 v ; 1
using
X = cos and Y = sin , where is defined to be positive in the anticlockwise direction.
L L p 1
Neglecting second order terms in the binomial expansion (1 + ") = 1 + " + O ("2 ), the strain
2
for small displacements u and v is
e u v
= cos : L + sin : L (4.1)
4.2 T RUSS E LEMENTS 77
The strain energy associated with this uniaxial stretch is
1 Z 1 Z L
1 Z L
SE = e dV = A e dx = EAe 2 dx = 1 ALEe2 (4.2)
2 2 2 2
0 0
where = Ee is the stress in the truss (of cross-sectional area A), linearly related to the strain e
via Young’s modulus E . We now substitute for e from Equation (4.1) into Equation (4.2) and put
u = u2 ; u1 and v = v2 ; v1 , where (u1; v1) and (u2; v2) are the nodal displacements of the two
ends of the truss
1
u2 ; u1 + sin : v2 ; v1
2
SE = ALE cos : (4.3)
2 L L
The potential energy is the combined strain energy from all trusses in the structure minus the
work done on the structure by external forces. The Rayleigh-Ritz approach is to minimize this
potential energy with respect to the nodal displacements once all displacement boundary conditions
have been applied.
For example, consider the system of three trusses shown in Figure 4.2. A force of 100 kN
is applied in the x-direction at node 1. Node 2 is a sliding joint and has zero displacement in the
y-direction only. Node 3 is a pivot and therefore has zero displacement in both x- and y - directions.
The problem is to find all nodal displacements and the stress in the three trusses.
node 1
100 kN
1
30 2
30
node 3 node 2
3
F IGURE 4.2: A system of three trusses.
The strain in truss 1 (joining nodes 1 and 3) is

p
u1 cos 30 + v1 sin 30 = 3 u1 + 1 v1
L L 2 L 2L
(u1 ; u2) cos 90 + v1 sin 90 = v1

L L L

p
u2 cos (;30) = 3 u2
L 2 L
Since a force of 100 kN acts at node 1 in the x-direction, the potential energy is
2 p !2 p !2 3
X1 2 ; 100u = 1 AE 4 3 u + 1 v + 3 u + (v )25 ; 100u
PE =
2
ALEe 1
2 L 2 1 2 1 2 2 1 1
trusses
[Note that if the force was applied in the negative x-direction, the final term would be +100u 1 ]
Minimizing the potential energy with respect to the three unknowns u 1 , v1 and u2 gives
p !p
@ PE = AE 3u + 1v 3 ; 100 = 0
@u1 L 2 1 2 1 2 (4.4)
" p ! #
@ PE = AE 3 u 1 1 +v =0
1 + v1
2 2 1
(4.5)
@v1 L 2
p !p
@ PE = AE 3u 3 =0
@u2 L 2 2 2 (4.6)
If we choose A= 5 10;3 m2, E = 10 GPa and L = 1 m (e.g., 100 mm 50 mm timber

truss) then
AE = 5 10;3 m2 107 kPa=m = 5 104 kN m;1.
L
Equation (4.6) gives
u2 = 0
Equation (4.4) gives
p
3u1 + 3v1 = 4 102= 5 104
;
Equation (4.5) gives for two dimensions
p
v1 = ; 53 u1
p these last two equations gives u1 = 3:34 mm and v1 = ;1:15 mm. Thus the strain in truss
Solving
1 is ( 23 3:34 ; 12 1:15) 10;3 = 0:232%, in truss 2 is ;0:115% and in truss 3 is zero.
The tension in truss 1 is A = AEe = 5 10 ;3 m2 107 kPa 0:232 10;2 = 116 kN (tensile),
in truss 2 is ;57:5 kN (compressive) and in truss 3 is zero. The nodal reaction forces are shown in
Figure 4.3.
4.3 B EAM E LEMENTS 79
100 kN
100 kN
57:7 kN
57:7 kN
F IGURE 4.3: Reaction forces for the truss system of Figure 4.2.
4.3 Beam Elements

Simple beam theory ignores all but axial strain e x and stress x = Eex (E = Young’s modulus)
along the beam (assumed here to be in the x-direction). The axial strain is given by e x =
z ,
R
where z is the lateral distance from the neutral axis in the plane of the bending and R is the radius
Z
of curvature in that plane. The bending moment is given by M = xz dA , where A is the beam
crossectional area. Thus
x = Eex = E Rz (4.7)
Z Z
M= E
xz dA = R z2 dA = EI
R (4.8)
Z
where I = z2 dA is the second moment of area of the beam cross-section. Thus, E =M
R I and
Equation (4.7) becomes
x = Mz
I (4.9)
The slope of the beam is

dw = and the rate of change of slope is the curvature
dx
K = dxd = d2w = 1 (4.10)
dx2 R
Thus the bending moment is
M = EI ddxw2 = EIw00
2
(4.11)
and a force balance gives the shear force
V = ; dM
dx = ; d (EIw00)
dx (4.12)
and the normal force (per unit length of beam)
p = dV
dx = ; d2 (EIw00 )
dx2 (4.13)
This last equation is the equilibrium equation for the beam, balancing the loading forces p with the
axial stresses associated with beam flexure

d2 EI d2w = p

; dx 2 dx2 (4.14)
The elastic strain energy stored in a bent beam is the sum of flexural strain energy and shear
strain energy, but this latter is ignored in the simple beam theory considered here. Thus, the
(flexural) strain energy is
1 ZL Z 1 Z Z L
SE =
2 xex dA dx = 2 E e2x dA dx
x=0 A x=0 A
1Z Z z 2 ZL
L
E 1 00 2
=2 R dA dx = 2 EI (w ) dx
x=0 A x=0
where x is taken along the beam and A is the cross-sectional area of the beam.
The external work associated with forces p acting normal to the beam and moving through a
LZ
transverse displacement w is pw dx. The potential energy is therefore
0
1 Z LZ L
2 EI (w ) dx ; pw dx:
PE = 00 2
(4.15)
0 0
The finite element approximation for the transverse displacement w must be able to represent
the second derivative w 00 . A linear basis function has a zero second derivative and therefore cannot
represent the flexural strain. The natural choice of basis function for beam deflection is in fact cubic
Hermite because the inter-element slope continuity of this basis ensures transmission of bending
moment as well as shear force across element boundaries.
The boundary conditions associated with the 4 th order equilibrium Equation (4.14) or the equa-
4.4 P LANE S TRESS E LEMENTS 81
tions arising from minimum potential energy Equation (4.15) (which contain the square of 2 nd
derivative terms) are more complex than the simple temperature or flux boundary conditions for
the (second order) heat equation. Three possible combinations of boundary condition with their
associated reactions are
Boundary conditions Reactions
(i) Simply supported zero displacement w = 0 shear force V

zero moment M = EIw 00 =0 slope (= w 0 )
(ii) Cantilever zero displacement w = 0 shear force V
zero slope = w 0 = 0 moment M
(iii) Free end = ; d (EIw00 ) = 0
zero shear force V displacement w
dx
zero moment M = EI 00 = 0 slope
4.4 Plane Stress Elements u

For two-dimensional problems, we define the displacement vector u = v , strain vector e =
2e 3 2 3
4 exy 5 and stress vector = 4 xy 5. The stress-strain relation for two-dimensional plane stress:
exy xy
x = 1 ;E 2 (ex + ey )
y = 1 ;E 2 (ey + ex ) (4.16)
E (e )
xy = 1 + xy
can be written in matrix form
= Ee
21 0
3
where E = 1; E 4 1 0 5. The strain components are given in terms of displacement
2
0 0 1;
gradients by
ex = @u
@x
ey = @v
@y
1 @u @v
(4.17)
exy = 2 @y + @x
The strain energy is
1 Z 1 Z
SE =
2 T e dV = 2 (exx + ey y + exy xy ) dV
V V
= 1Z eT 1 Z E
Ee dV = 2 1 ; 2 e2x + e2y + 2exey + (1 ; ) e2xy dV
2
V V
The potential energy is
1 Z Z
PE = SE ; external work =2 e Ee dV ; uT l dA
T (4.18)
V A
l
where represents the external loads (forces) acting on the elastic body.
Following the steps outlined in Section 4.1 we approximate the displacement field u with a
finite element basis u = 'n un , v = 'n vn and calculate the strains
ex = @u = @'n un
@x @x
@v = @'n v
ey = @y
@y n @'
(4.19)
1 @u @v 1 n @'n
exy = 2 @y + @x = 2 @y un + @x vn
or
2 @'n 3
2 e 3 66 @x 0 77
e = 4 ey 5 = 666 0 @' n 7 un
x
@y 77 vn = Bu (4.20)
exy 4 1 @'n 1 @'n 5
2 @y 2 @x
From Equation (4.18) the potential energy is therefore
1 Z Z
PE = (Bu) E (Bu) dV ;
T uT l dA
2
V A
Z Z
= 12 uT BT EB dV:u ; uT l dA
V A
Z
= 12 uT Ku ; uT l dA
A
4.4 P LANE S TRESS E LEMENTS 83
Z
where K = BT EB dV is the element stiffness matrix.
V
We next minimize the potential energy with respect to the nodal parameters u n and vn giving
Ku = f (4.21)
Z
where f = l dA is a vector of nodal forces.
A
4.4.1 Notes on calculating nodal loads

If a known stress acts normal to a given surface (e.g., a surface pressure), it may be applied by
calculating equivalent nodal forces. For example, consider a uniform load p kN m ;1 applied to the
edge of the plane stress element in Figure 4.4a.
f
The nodal load vector in Equation (4.21) has components
Z Z1
fn = p'n dx = pL 'n d (4.22)
x 0
where is the normalized element coordinate along the side of length L loaded by the constant
stress p kN m;1 . If the element side has a linear basis, Equation (4.22) gives
Z1 Z1
f1 = pL '1 d = pL (1 ; ) d = 21 pL
0 0
Z 1 Z1
f2 = pL '2 d = pL d = 12 pL
0 0
as shown in Figure 4.4b. If the element side has a quadratic basis, Equation (4.22) gives
Z1 Z1 1
f1 = pL '1 d = pL 2 2 ; (1 ; ) d = 61 pL
0 0
Z 1 Z1
f2 = pL '2 d = pL 4 (1 ; ) d = 23 pL
0 0
Z 1 Z1
f3 = pL '3 d = pL 2 ; 12 d = 16 pL
0 0
as shown in Figure 4.4c. A node common to two elements will receive contributions from both
elements, as shown in Figure 4.4d.
p kN m;1 1 pL 1 pL 1 ;L
p
2p ;L
2 6 6 2
L 2 pL
3
1p L;
2
2p ;L
3 3 2
1 pL 1 pL
;
3
1p L 2
2 6 6 2
(a) (b) (c) (d)
F IGURE 4.4: A uniform boundary stress applied to the element side in (a) is equivalent to nodal
loads of 12 pL and 12 pL for the linear basis used in (b) and to 16 pL, 23 pL and 16 pL for the quadratic
basis used in (c). Two adjacent quadratic elements both contribute to a common node in (d), where
the element length is now L2 .
4.5 Three-Dimensional Elasticity

Consider a surface ; enclosing a volume
of material of mass density . Conservation of linear
momentum over the domain
results in the governing equilibrium equations
ij;j + bi = 0 i; j = 1; 2; 3 (4.23)
where ij are the components of the stress tensor (ij is the component of the traction or stress
vector in the ith direction which is acting on the face of a rectangle whose normal is in the j th
b g
direction), and bi is the body force/unit volume (e.g., = ). Note that the notation ij;j =
@ij
@xj
has been introduced to represent the derivative.
Recall that the components of the linear (or small) strain tensor are
eij = 12 (ui;j + uj;i) i; j = 1; 2; 3 (4.24)
u u
where is the displacement vector (i.e., is the difference between the final and initial positions
of a material point in question). Note: we are assuming here that the displacement gradients are
small compared to unity, which is appropriate for many materials in solid mechanics. However, for
soft materials, such as rubber or biological tissue, then we need to use the exact finite strain tensor.
The object of solving an elasticity problem is to find the distributions of stress and displacement
in an elastic body, subject to a known set of body forces and prescribed stresses or displacements
at the boundaries. In the general three-dimensional case, this means finding 6 stress components
ij (= ji which arises from the conservation of angular momentum) and 3 displacements u i each
as a function of position in the body. Currently we have 15 unknowns (6 stresses, 6 strains and 3
displacements), but only 9 equations (3 equilibrium equations and 6 strain-displacement relations).
To progress, we require an equation of state, i.e., a stress-strain relation or constitutive law. For
a linear elastic material we may propose that the components of stress ij depend linearly on eij .
4.5 T HREE -D IMENSIONAL E LASTICITY 85
i.e.,
ij = cijklekl
where cijkl are the 81 components of a 4th order tensor, although symmetry of the strain and stress
tensors reduces the number of independent components to 21.
If the material is assumed to be isotropic (i.e., the material response is independent of orienta-
tion of the material element), then we end up with the generalized Hooke’s Law.
ij = ekk ij + 2eij (4.25)
or inversely
eij = 21 ij ; 2 (3+ 2) kk ij

where , are Lamés constants.
Note: , are related to Young’s modules E and Poisson’s ratio by
E = (3++2)
= 2 (+ )
Providing that the displacements are continuous functions of position, then Equation (4.23),
Equation (4.24) and Equation (4.25) are sufficient to determine the 15 unknown quantities. This
can often work with some smaller grouping or simplification of these equations, e.g., if all bound-
ary conditions are expressed in terms of displacements, substituting Equation (4.24) into Equa-
tion (4.25) then into Equation (4.23) yields Navier’s equation of motion.
ui;kk + ( + ) uk;ki + bi = 0 i; k = 1; 2; 3
These 3 equations can be solved for the unknown displacements. Then Equation (4.24) can be used
to determine the strains and Equation (4.25) to calculate the stresses.
4.5.1 Weighted Residual Integral Equation

Using weighted residuals as before we can write
Z
(ij;j + bi ) ui d
= 0 (4.26)

u
where = (ui ) is a (vector) weighting field. The ui are usually interpreted as a consistent set of
virtual displacements (hence we use the notation u instead of w ).
By the chain-rule
(ij ui );j = ij;j ui + ij ui;j

Therefore, the first term in the integrand of Equation (4.26) can be re-written
Z Z Z
ij;j u d
=
i (ij i ;j d
;
u) ij ui;j d
Z
= r (ij u)
i d
; ij ui;j d
Z
Z
= ij ui nj d; ; ij ui;j d

(4.27)
@

where the domain integral involving “r = @x@ ” has been transformed into a surface integral
j
using the divergence theorem
Z Z Z Z
r g d
= g n d; or gj;j d
= gj nj d;

@

@
n i
where = nj j is the outward normal vector to the surface ;.
Thus, combining Equation (4.26) and Equation (4.27) we have
Z Z Z
ij u
i;j d
= b u d
+
i i ij nj ui d;
@
Z
Z
= b u d
+
i i tiui d; (4.28)

@
t
where ti are the components of the internal stress vector ( ) and are related to the components of
the stress tensor (ij ) by Cauchy’s formula
t = ij nj ii (4.29)
To arrive at this point, we have used weighted residuals to tie in with Chapter 2, however
Equation (4.28) is more usually derived using the principle of virtual work (below). Note that the
weighted integral Equation (4.28) is independent of the constitutive law of the material.
4.5.2 The Principle of Virtual Work

The governing equations for elastostatics can also be derived from a physically appealing argument.
s
Let be the external traction vector (i.e., force per unit surface area). For equilibrium, the work
s i
done by the external surface forces = si i , in moving through a virtual displacement = ui i u i
4.5 T HREE -D IMENSIONAL E LASTICITY 87
t i
is equal to the work done by the stress vector = t i i in moving through a compatible set of virtual
u
displacements . In mathematical terms, the principle of virtual work can be written
Z Z Z
s u d; =
i i t u d; =
i i ij nj ui d; (4.30)
@
@
@
using Cauchy’s formula (Equation (4.29)).

The Green-Gauss theorem (Equation (2.15)) is now used to replace the right hand surface
integral in Equation (4.30) by a volume integral, giving
Z Z ;
s u d; =
i i ij;j ui + ij ui;j d
(4.31)
@

Substituting the equilibrium relation (Equation (4.23)) into the first integrand on the right hand
side, yields the virtual work equation
Z Z Z
ij u
i;j d
= b u d
+
i i siui d; (4.32)

@
where the internal work done due to the stress field is equated to the work due to internal body
forces and external surface forces. Note that Equation (4.32) is equivalent to Equation (4.28) via
Equation (4.30). In practice, Equation (4.32) is in a more useful form than Equation (4.28), because
the right hand side integrals can be expressed in terms of the known body forces and the applied
boundary conditions (surface traction forces or stresses).
4.5.3 The Finite Element Approximation

S

Let
= e and interpolate the virtual displacements u i from their nodal values. i.e.,
ui = 'm(_ umi)

so ui;j = @'m_ m
@xj (ui ) (4.33)
= 'm;k @k (_ umi)

@xj
where (um
i )
= (Ui4(m;e)), 4(m; e) is the global node number of local node m on element e, and
the shorthand 'm;k =
@'m has been introduced.
@ k
Substituting this into Equation (4.32) gives
0 1 0Z 1
X @Z X Z
@k d
A U 4(m;e) = @ b ' d
+ s ' d;A U 4(m;e)
ij 'm;k @x i i m i m i
e j e

e
e @
e
and since the virtual displacements are arbitrary we get

0 1
XZ @k d
= X @Z Z
ij 'm;k @x bi 'm d
+ si'm d;A (4.34)
e
e j e
e @
e
The next step is to express the stress components ij in terms of the virtual displacements
and their finite element approximation by substituting Equation (4.33) into Equation (4.24) (the
strain-displacement relation) and in turn into Equation (4.25) (the generalized Hooke’s law).
We first introduce the finite element approximation for the displacement field u j = 'n unj which
gives
;
eij = 12 @x@ ('nuni) + @x@ 'nunj = 12 @' n @l n @'n @l n
@l @xj ui + @l @xi uj (4.35)
j i
and
ekk = uk;k = @' n @l n

@ @x uk l k
Thus

ij = ij @'n @l n
@l @xk u k + 2 1 @'n @l un + 1 @'n @l un
2 @l @xj i 2 @l @xi j
which, due to symmetry of the stress tensor, simplifies to
ij = ij @' n @l n

u k + 2 @'n @l un
@l @xi j
@l @x@ k
@l un
= i(j) 'n;l @xl + 2'n;l @x j (4.36)
j i
where the summation index k has been replaced with j , but the parenthesis in i(j ) implies that
there is no sum with respect to that particular index.
Substituting this expression into Equation (4.34) and simplifying, we get for each element
Z @l ' @k + 2' @l ' @k d
= f

unj 'n;l @x m;k
@xi n;l
@xi m;k @xj im (4.37)
j

e
where fim denotes the right hand side terms in Equation (4.34). (Note that there has been some
careful manipulation of summation indices with the substitution of Equation (4.36) to arrive at
Equation (4.37).)
So for each element
Eimjnunj = fim
4.6 L INEAR E LASTICITY WITH B OUNDARY E LEMENTS 89
where
ZZZ
1 @ @ @

l k l @k
Eimjn = + 2
@xj @xi 'n;l 'm;k J (1; 2; 3)d1d2d3
@xi @xj
0
ZZZ
1 ZZ1 (4.38)
fim = bi 'mJ (1 ; 2; 3)d1d2d3 + si'm J2D (1; 2)d1d2

0 0
where the Jacobians J (1 ; 2 ; 3 ) and J2D (1 ; 2 ) have been used to transform volume and sur-
face integrals so that they can be can be calculated using -coordinates. (Note: without loss of
generality, the above definition of fim assumes that (1 ; 2 ) are defined to lie in the surface ;.)
So in summary, the finite element approximation leads to element stiffness matrix components
that can be calculated from the known material parameters, the chosen interpolation functions, and
the geometry of the material (note that the element stiffness components are independent of the
unknown displacement parameters). Element stiffness components are then assembled into the
global stiffness matrix in the usual manner (as described previously). Note that this is implicitly a
Galerkin formulation, since the unknown displacement fields are interpolated using the same basis
functions as those used to weight the integral equations.
4.6 Linear Elasticity with Boundary Elements

Equation (4.28) is the starting point for the general finite element formulation (Section 4.5). In
the above derivation, we have essentially used the Green-Gauss theorem once to move from Equa-
tion (4.26) to Equation (4.28) (as was done for the derivation of the FEM equation for Laplace’s
equation). To continue, we firstly note that
ij eij = 12 ij ui;j + 12 ij uj;i

= 1 ij ui;j + 1 jiuj;i
2 2
= ij ui;j + 1 ij ui;j
1
2 2
= ij ui;j
where eij are the virtual strains corresponding to the virtual displacements.
Using the constitutive law for linearly elastic materials (Equation (4.25)) we have
Z Z
ij u
i;j d
= ij eij d
Z Z
= ekk e
ij ij d
+ 2 eij eij d
Z
Z
= ekk ekk d
+ 2 eij eij d
= eij ij d
due to symmetry.
Thus from the virtual work statement, Equation (4.28) and the above symmetry we have
Z Z Z Z
b u d
+
i i t u d; =
i i b u
i i d
+ ti ui d; (4.39)

@

@
This is known as Betti’s second reciprical work theorem or the Maxwell-Betti reciprocity relation-
ship between two different elastic problems (the starred and unstarred variables) established on the
same domain.
Note that bi = ;ij;j
(i.e., + b = 0). Therefore Equation (4.39) can be written as
ij;j i
Z ; Z Z Z
ij;j ui d
+ bi ui d
= ti ui d; ; tiui d;
(4.40)

@
@
(ij ; eij ; ti represents the equilibrium state corresponding to the virtual displacements u i ).
Note: What we have essentially done is use integration of parts to get Equation (4.28), then use
it again to get Equation (4.39) above (after noting the reciprocity between ij and eij ).
Since the body forces, bi , are known functions, the second domain integral on the left hand
side of Equation (4.40) does not introduce any unknowns into the problem (more about this later).
The first domain integral contains unknown displacements in
and it is this integral we wish to
remove.
We choose the virtual displacements such that
+e =0
ij;j i (4.41)
(or equivalently ;bi + ei = 0), where ei is the ith component of a unit vector in the ith direction
x
and ei = ei ( ; P ). We can interpret this as the body force components which correspond to a
positive unit point load applied at a point P 2
in each of the three orthogonal directions.
Therefore
4.7 F UNDAMENTAL S OLUTIONS 91
Z Z
ij;j ui d
= ;
(x ; P ) eiui d
= ;ui(P )ei

i.e., the volume integral is replaced with a point value (as for Laplace’s equation).
Therefore, Equation (4.40) becomes
Z Z Z
ui (P ) ei = tj u d; ;
j t u
j j d; + bj uj d
P 2
(4.42)
@
@

If each point load is taken to be independent then u j and tj can be written as
uj = uij (P; x) ei (4.43)

tj = tij (P; x) ei (4.44)
where uij (P; x) and tij (P; x) represent the displacements and tractions in the j th direction at x
corresponding to a unit point force acting in the i th direction (ei ) applied at P . Substituting these
into Equation (4.42) (and equating components in each e i direction) yields
Z Z
ui (P ) = u
ij (P; x) tj (x) d; (x) ; tij (P; x) uj (x) d; (x)
@
@
Z
+ uij (P; x) bj (x) d
(x) (4.45)

where P 2
(see later for P 2 @
).
This is known as Somigliana’s 1 identity for displacement.
4.7 Fundamental Solutions

Recall from Equation (4.41) that ij satisfied
+ (x ; P ) e = 0
ij;j i (4.46)
or equivalently
bi = ei (x ; P )
Navier’s equation for the displacements ui is
G ui;kk + 1 ;G2 uk;ki + bi = 0

1
Somigliana was an Italian Mathematician who published this result around 1894-1902.
where G = shear Modulus.

Thus ui satisfy
G ui;kk + 1 ;G2 uk;ki + (x ; P ) ei = 0 (4.47)
The solutions to the above equation in either two or three dimensions are known as Kelvin 2 ’s
fundamental solutions and are given by
uij (P; x) = 16 (1 1; ) Gr f(3 ; 4 ) ij + r;ir;j g (4.48)
for three-dimensions and for two-dimensional plane strain problems,
uij (P; x) = 8 (1;;1 ) G f(3 ; 4 ) ij log r ; r;ir;j g (4.49)
and
; 1
@r

t ij (P; x) = 4 (1 ; ) r ((1 ; 2 ) ij + r;ir;j ) @n ; (1 ; 2 ) (r;i nj ; r;j ni ) (4.50)
where = 1; 2; = 2; 3 for two-dimensional plane strain and three-dimensional problems respec-

tively.
x x
Here r r (P; ), the distance between load point (P ) and field point ( ), ri = xi ( ) ; xi (P ) x
and r;i =
@r = ri .
@xi ( ) r x
x
In addition the strains at an point due to a unit point load applied at P in the i th direction are
given by
ejki (P; x) = 8 (1;;1 ) Gr f(1 ; 2 ) (r;k ij + r;j ik ) ; r;ijk + r;ir;j r;k g
and the stresses are given by
(P; x) = ;1
ijk 4 (1 ; ) r f(1 ; 2 ) (r;k ij + r;j ki ; r;ijk ) + r;ir;j r;k g
where and are defined above.
The plane strain expressions are valid for plane stress if is replaced by = 1 + (This is a
mathematical equivalence of plane stress and plane strain - there are obviously physical differences.
What the mathematical equivalence allows us to do is to use one program to solve both types of
problems - all we have to do is modify the values of the elastic constants).
Note that in three dimensions
1 1
uij =O r t
ij = O r2
2
Lord Kelvin (1824-1907) Scottish physicist who made great contributions to the science of thermodynamics
4.8 B OUNDARY I NTEGRAL E QUATION 93
and for two dimensions

1
u ij = O (log r) t
ij =O
r :
Somigliana’s identity (Equation (4.45)) is a continuous representation of displacements at any
point P 2
. Consequently, one can find the stress at any P 2
firstly by combining derivatives
of (4.45) to produce the strains and then substituting into Hooke’s law. Details can be found in
Brebbia, Telles & Wrobel (1984b) pp 190–191, 255–258.
This yields
Z Z
ij (P ) = u
ijk (P; x) tk (x) d; (x) ; tijk (P; x) uk (x) d; (x)
Z; ;
+ uijk (P; x) bk (x) d

(x)

Note: One can find internal stress via numerical differentiation as in FE/FD but these are not
as accurate as the above expressions.
Expressions for the new tensors uijk and tijk are on page 191 in (Brebbia et al. 1984b).
4.8 Boundary Integral Equation

Just as we did for Laplace’s equation we need to consider the limiting case of Equation (4.45) as
P is moved to @
. (i.e., we need to find the equivalent of c (P ) (in section 3) - called here cij (P ).)
We use the same procedure as for Laplace’s equation but here things are not so easy.
If P 2 @
we enlarge
to
0 as shown.
;"

0
"
;;"
P

F IGURE 4.5: Illustration of enlarged domain when singular point is on the boundary.
Then Equation (4.45) can be written as

Z Z
ui (P ) = uij (P; x) tj (x) d; (x) ; tij (P; x) uj (x) d; (x)
;;" +;" ;;" +;"
Z
+ uij (P; x) bj (x) d
(x) (4.51)

0
We need to look at each integral in turn as "# 0 (i.e., " ! 0 from above). The only integral that
presents a problem is the second integral. This can be written as
Z Z
t
ij (P; x) uj (x) d; (x) = tij (P; x) uj (x) d; (x)
;;" +;" ;"
Z
+ tij (P; x) uj (x) d; (x) (4.52)
;;"
The first integral on the right hand side can be written as

Z Z
tij (P; x) uj (x) d; (x) = tij (P; x) [uj (x) ; uj (P )] d;(x)
;"
|
;"
{z }
0 uj (x)
Z
by continuity of
+ uj (P ) tij (P; x) d; (x) (4.53)

;"
Let
Z
cij (P ) = ij + lim
"#0
tij (P; x) d; (x) (4.54)
;"
Z
As "# 0, ;;" ! ; and we write the second integral of Equation (4.52) as tij (P; x) uj (x) d; (x)
;
where we interpret this in the Cauchy Principal Value3 sense.
3
What is a Cauchy Principle Value?
Consider f (x) =
1 on ; = [;2; ;") [ ("; 2]
x
;"
4.8 B OUNDARY I NTEGRAL E QUATION 95
Thus as "# 0 we get the boundary integral equation

Z
cij (P ) uj (P ) + tij (P; x) uj (x) d; (x)
;
Z Z
= u
ij (P; x) tj (x) d; (x) + uij (P; x) bj (x) d
(4.55)
;
Z Z
(or, in brief (no body force), cij uj + t u
ij j d; = uij tj d;) where the integral on the left hand
; ;
side is interpreted in the Cauchy Principal sense. In practical applications c ij and the principal value
integral can be found indirectly from using Equation (4.55) to represent rigid-body movements.
The numerical implementation of Equation (4.55) is similar to the numerical implementation
of an elliptic equation (e.g., Laplace’s Equation). However, whereas with Laplace’s equation the
unknowns were u and
@u
@n (scalar quantities) here the unknowns are vector quantities. Thus it is
more convenient to work with matrices instead of indicial notation.
i.e., use
2u 3 2t 3
1 1
u = 4u25 ; t = 4t25
u3 t3
2u u u 3 2t t t 3
u = 4u21 u22 u23 5 ; t = 4t21 t22 t23 5
11 12 13 11 12 13
u31 u32 u33 t31 t32 t33

Then
Z Z;" 1 Z2 1 dx = ln jxjj;" + ln jxjj2
f (x) dx = dx + ;2 "
x x
;;" ;2 "
= ln " ; ln 2 + ln 2 ; ln " = 0 8" > 0

Z Z
) "lim
!0 f (x) dx =0 This is the Cauchy Principle Value of f (x) dx
;;" ;
But if we replace ;;" by

!0;;" = [;2; 2] = ; then
lim
"
Z 0 1 0 1
1 dx = Z 1 dx = @ lim Z" Z2
2 1
1 dxA + @ lim 1 dxA(by definition of improper integration)
x x "1 !0 x "2 !0 x
; ;2 ;2 ;"2
which does NOT exist. i.e., the integral does not exist in the proper sense, but it does in the Cauchy Principal Value
sense. However, if an integral exists in the proper sense, then it exists in the Cauchy Principal Value sense and the two
values are the same.
Then (in absence of a body force) we can write Equation (4.55) as

Z Z
cu + tu d; = ut d; (4.56)
; ;
We can discretise the boundary as before and put P , the singular point, at each node (each
node has 6 unknowns - 3 displacements and 3 tractions - we get 3 equations per node). The overall
matrix equation
Au = Bt (4.57)
2 3 2 3
u1 t
66u2 77 66t12 77
where u = 6 .. 7 and t = 6 .. 7 where n is the number nodes.
4.5 4.5
un tn
The diagonal elements of the A matrix in Equation (4.57) (for three-dimensions, a 3 x 3 matrix)
contains principal value components. If we have a rigid-body displacement of a finite body in any
one direction then we get
Ail = 0
i
( l = vector defining a rigid body displacement in direction l)
X
) aii = ; aij (no sum on i)
i6=j
i.e., the diagonal entries of A (the cij ’s) do not need to be determined explicitly. There is a similar
result for an infinite body.
4.9 Body Forces (and Domain Integrals in General)

The body force gives rise to a domain integral although it does not give rise to any further unknowns
in the system of equations. (This is because the body force is known - the fundamental solution
was chosen so that it removed all unknowns appearing in domain integrals).
Thus Equation (4.55) is still classed as a Boundary Integral Equation. Integrals over the domain
containing known functions (eg body force integral) appear in many situations e.g., the Poisson
equation r2 u = f yields a domain integral involving f .
The question is how do we evaluate domain integrals such as those appearing in the boundary
integral forumalation of such equations? Since the functions are known a coarse domain mesh
may work.(n.b. Since the integral also contains the fundamental solution and
may not be a
“nice” region it is unlikely that it can be evaluated analytically). However, a domain mesh nullifies
one of the advantages of BEM - that of having to prepare only a boundary mesh.
In some cases domain integrals must be used but there are techniques developing to avoid many
of them. In some standard situations a domain integral can be transformed to a boundary integral.
4.9 B ODY F ORCES ( AND D OMAIN I NTEGRALS IN G ENERAL ) 97
e.g., a body force arising from a constant gravitational load, or a centrifugal load due to rotation
about a fixed axis or the effect of a steady state thermal load can all be transformed to a boundary
integral.
Firstly, let Gij (the Galerkin tension) be related to uij by
uij = Gij;kk ; 2 (1 1; ) Gik;kj

8 1
>
< rij (3D)
) Gij = > 8G 1 r2 log 1
: 8G r ij (2D)
Then
Z Z 1

Bi = u b
ij j d
= Gij;kk ; 2 (1 ; ) G ik;kj bj d
Under a constant gravitational load g = (gj )

bj = gj
Z 1

) Bi = gj G ik;j ; 2 (1 ; )
Gik;kj d
Z
( )
= gj Gij;k ; 2 (1 ; 1g) G nk d;
;
ik;j
which is a boundary integral.

Unless the domain integrand is “nice” the above simple application of Green’s theorem won’t
work in general. There has been a considerable amount of research on domain integrals in BEM
which has produced techniques for overcoming some domain methods. The two integrals of note
are the DRM, dual reciprocity method, developed around 1982 and the MRM, multiple reciprocity
method, developed around 1988.
4.10 CMISS Examples

1. To solve a truss system run CMISS example 411 This solves the simple three truss system
shown in Figure 4.2.
2. To solve stresses in a bicycle frame modelled with truss elements run CMISS example 412.
Chapter 5
Transient Heat Conduction
5.1 Introduction
In the previous discussion of steady state boundary value problems the principal advantage of the
finite element method over the finite difference method has been the greater ease with which com-
plex boundary shapes can be modelled. In time-dependent problems the solution proceeds from
an initial solution at t = 0 and it is almost always convenient to calculate each new solution at a
constant time (t > 0) throughout the entire spatial domain
. There is, therefore, no need to use
the greater flexibility (and cost) of finite elements to subdivide the time domain: finite difference
approximations of the time derivatives are usually preferred. Finite difference techniques are intro-
duced in Section 5.2 to solve the transient one dimensional heat equation. A combination of finite
elements for the spatial domain and finite differences for the time domain is used in Section 5.3 to
solve the transient advection-diffusion equation - a slight generalization of the heat equation.
5.2 Finite Differences

5.2.1 Explicit Transient Finite Differences
Consider the transient one-dimensional heat equation
@u = D @ 2 u ; (0 < x < L; t > 0) (5.1)

@t @x2
where D is the conductivity and u = u (x; t) is the temperature, subject to the boundary conditions
u (0; t) = u0 and u (L; t) = u1 and the initial conditions u (x; 0) = 0. A finite difference approxi-
mation of this equation is obtained by defining a grid with spacing x in the x-domain and t in
the time domain, as shown in Figure 5.1.
Grid points are labelled by the indices i = 0; 1; : : : ; I (for the x-direction) and n = 0; 1; : : : ;N
(for the t-direction). The temperature at the grid point (i; n) is therefore labelled as
u (x; t) = u (ix; nt) = uni: (5.2)
Finite difference equations are derived by writing Taylor Series expansions for uni+1 ; uni;1 uni +1
100 T RANSIENT H EAT C ONDUCTION
about the grid point (i; n)

@u n
n n ;
+ 1 x2 : @ u2 + 1 x3 : @ u3 + O x4
2 3
uni+1 = uni + x: (5.3)
@x
@u
i 2
n 1
@x i 6
@ 2u n
1
@x i
@ 3 u n ;
ui;1 = ui ; x: @x + 2 x : @x2 ; 6 x : @x3 + O x4
n n 2 3 (5.4)
@u ni ; i i
uni +1 = uni + t: @t + O t2 (5.5)
i
where O (x4 ) and O (t2 ) represent all the remaining terms in the Taylor Series expansions.
Adding Equations (5.3) and (5.4) gives
@ 2 u n ;
uni+1 + uni;1 = 2uni + x2 : @x2 + O x4
i
or
@2 u n un ; 2un + un
= i+1 i i;1 + O ;x2 ; (5.6)
@x2 i x2
which is a “central difference” approximation of the second order spatial derivative.
Rearranging Equation (5.5) gives a “difference” approximation of the first order time derivative
@u n un+1 ; un
= i i + O (t) : (5.7)
@t i t
Substituting Equation (5.6) and Equation (5.7) into the transient heat equation Equation (5.1)
gives the finite difference approximation
uni +1 ; uni + O (t) = D uni+1 ; 2uni + uni;1 + O ;x2

t x2
which is rearranged to give an expression for uni +1 in terms of the values of u at the nth time step
; ;
uni +1 = uni + D xt2 uni+1 ; 2uni + uni;1 + O t2 ; x2 :
(5.8)
Given the initial values of uni at n = 0 (i.e., t = 0), the values of uni +1 for the next time step
are found from Equation (5.8) with i = 1; 2; : : : ; I . Applying Equation (5.8) iteratively for time
steps n = 1; 2; : : : etc. yields the time dependent temperatures at the grid points (see Figure 5.1).
This is an explicit finite difference formula because the value of uni depends only on the values of
uni (i = 1; 2; : : : ; I ) at the previous time step and not on the neighbouring terms u ni+1
+1 and un+1 at
i;1
the latest time step. The accuracy of the solution depends on the chosen values of x and t and
in fact the stability of the scheme depends on these satisfying the Courant condition:
D xt2 12 : (5.9)
5.2 F INITE D IFFERENCES 101
n+1 x
n x x x
:
1
00 1 i;1 i i+1 I x
2 ... ...
F IGURE 5.1: A finite difference grid for the solution of the transient 1D heat equation. The
O
equation is centred at grid point (i; n) shown by the . The lightly shaded region shows where the
solution is known at time step n. With central differences in x and a forward difference in t an
explicit finite difference formula gives the solution at time step n + 1 explicitly in terms of the
solution at the three points below it at step n, as indicated by the dark shading.
5.2.2 Von Neumann Stability Analysis

The concept behind the Von Neumann analysis is that all Fourier components decay as time ad-
vances or as they are processed by an iterative solver. Considering Equation (5.8), we can rearrange
this to be of the form,
uni +1 = uni+1 + (1 ; 2)uni + uni;1 (5.10)
where = D
t . By subsituting the general Fourier component un = Anei( kjLx ), we obtain,
x2 j k
kj x h k(j+1)x kj x k(j ;1)x i
Ank +1ei( L ) = Ank ei( L ) + (1 ; 2) ei( L ) + ei( L ) (5.11)
kj x
If divide Equation (5.11) by, Ank ei( L ) we obtain (no sum on k ),
Ank +1 = (1 ; 2) + ei( kjLx ) + e;i( kjLx )

Ank kj x ;ix
(x) = e + e
ix
= 1 ; 2 + 2cos using cos (5.12)
kj x L 2
= 1 ; 4sin 2 using cos(2x) = 1 ; 2sin2 (x)
2L
Equation (5.12) predicts the growth of any component (specified by k or j ) admitted by the
system. If all components are to decay,

An+1
k n 1 for stability (no sum on k ) (5.13)
Ak
As the sin2 term in Equation (5.12) is always between 0 and 1, we effective have the stablity
criteria that
1 ; 4 1 and 1 ; 4 ;1 (5.14)
The first inequality is trivially satisfied, since 0 for positive values of t and D, and the
second condition will always hold if
= D t2 1 (5.15)
x 2
Thus, to ensure stability, the time step should be chosen such that
t x
2
The Courant condition (5.16)
2D
5.2.3 Higher Order Approximations
An improvement in accuracy and stability can be obtained by using a higher order approximation
@u
;
for the time derivative. For example, if a central difference approximation is used for
centering the equation at (ix; n + 12 t) rather than (ix; nt) we get
@t by
@u n+ 12 un+1 ; un ;
= i i + O t2 (5.17)
@t i t
in place of Equation (5.7) and Equation (5.1) is approximated with the “Crank-Nicolson”formula
( n+1 n )
uni +1 ; uni = D 1 @ 2 u + 1 @ u2
2
(5.18)
t 2 @x2 i 2 @x i
in which the spatial second derivative term is weighted by 12 at the old time step n and by 12 at the
new time step n + 1. Notice that the finite difference time derivative has not changed - only the
time position at which it is centred. The price paid for the better accuracy (for a given t) and
unconditional stability (i.e., stable for any t) is that Equation (5.18) is an implicit scheme - the
equations for the new time step are now coupled in that u ni +1 depends on the neighbouring terms
uni+1
+1 and un+1 . Thus each new time step requires the solution of a system of coupled equations.
i;1
A generalization of (5.18) is
( n+1 @ 2 u n )
uni +1 ; uni = D @ 2 u + (1 ; ) @x2 (5.19)
t @x2 i i
5.3 T HE T RANSIENT A DVECTION -D IFFUSION E QUATION 103
n+1 x x x
n x x x
:
1
00 1 2 ... i;1 i i+1 ... I x

F IGURE 5.2: An implicit finite difference scheme based on central differences in t, as well as x,
1
x
which tie together the 6 points shown by . The equation is centred at the point (i; n + ) shown
2
O
by the . The lightly shaded region shows where the solution is known at time step n. The dark
shading shows the region of the coupled equations.
in which the spatial second derivative of Equation (5.1) has been weighted by at the new time step
and by (1 ; ) at the old time step. The original explicit forward difference scheme Equation (5.8)
is recovered when = 0 and the implicit central difference (Crank-Nicolson) scheme (5.19) when
= 21 . An implicit backward difference scheme is obtained when = 1.
In the following section the transient heat equation is approximated for numerical analysis
by using finite differences in time and finite elements in space. We also generalize the partial
differential equation to include an advection term and a source term.
5.3 The Transient Advection-Diffusion Equation

Consider a linear parabolic equation
@u + v ru = Dr2u + f (5.20)
@t
where u is a scalar variable (e.g., the advection-diffusion equation, where u is concentration or
v v
temperature; ru then represents advective transport by a velocity field ; D is the diffusivity
and f is source term. The ratio of advective to diffusive transport is characterised by the Peclet
v
number V L=D where V = k k2 and L is a characteristic length).
Applying the Galerkin weighted residual method to Equation (5.20) with weight gives !
Z @u
@t + v ru ; Dr u ; f ! d
= 0
2
or
Z @u Z Z @u
@t + v ru ! + Dru r! d
= f! d
+ D @n ! d; (5.21)

@
@
where
@n is the normal derivative to the boundary @
.
Putting u = 'n un and ! = 'm and summing the element contributions to the global equations,
Equation (5.21) can be represented by a system of first order ordinary differential equations,
M ddtu + Ku = Ku1 (5.22)
where M is the global mass matrix, K the global stiffness matrix and u a vector of global nodal
unknowns with steady state values (t ! 1) u1 . The element contributions to M and K are
given by
Z1
Mmne = 'm'nJ d (5.23)
0
and
Z1 @'m @'n @i @j Z1
Kmne = D @ @ @x @x J d + vj 'm @'n @i
@i @xj J d (5.24)
i j k k
0 0
If the time domain is now discretized (t = nt; n = 0; 1; 2; : : :) Equation (5.22) can be re-
placed by
M u ;t u + K un+1 + (1 ; ) un = Ku1

n+1 n
01 (5.25)
1
where is a weighting factor discussed in Section 5.2. Note that for = the method is known
2
as the Crank-Nicolson-Galerkin method and errors arising from the time domain discretization are
O (t2 ). Rearranging Equation (5.25) as
[M + tK ] un+1 = [M ; (1 ; ) tK ] un + tKu1 (5.26)
gives a set of linear algebraic equations to solve at the new time step (n + 1) t from the known
u
solution n at the previous time step nt.
u
The stability of the above scheme can be examined by expanding (assumed to be smoothly
s
continuous in time) in terms of the eigenvectors i (with associated eigenvalues i ) of the matrix
X
A= M ;1K . Writing the initial conditions u(0) = aisi and steady state solution u1 =
i
5.3 T HE T RANSIENT A DVECTION -D IFFUSION E QUATION 105
X
bi si , the set of ordinary differential equations Equation (5.22) has solution
i
X s
u= bi + (ai ; bi) e;it i (5.27)
i
u
The time-difference scheme Equation (5.26) on the other hand, with now replaced by a set
u
of discrete values n at each time step nt, can be written as the recursion formula
[I + tA] un+1 = [I ; (1 ; ) tA] un + tAu1 (5.28)
with solution
X 1 ; t (1 ; ) n
u= bi + (ai ; bi ) 1 + ti
i
si (5.29)
i
(You can verify that Equation (5.27) and Equation (5.29) are indeed the solutions of Equation (5.22)
and Equation (5.25), respectively, by substituting and using As
i = i i .) s
Comparing Equation (5.27) and Equation (5.29) shows that replacing the ordinary differential
equations (5.22) by the finite difference approximation Equation (5.25) is equivalent to replacing
the exponential e;i t in Equation (5.27) by the approximation
1 ; t (1 ; ) n
i
e;it 1 + ti (5.30)
or, with t = nt,
e;it 1 ;1+t (1 ; ) i = 1 ; ti

t 1 + t
(5.31)
i i
The stability of the numerical time integration scheme can now be investigated by examining
the behaviour of this approximation to the exponential. For stability we require
;1 1 ; 1 +tt
i
1 (5.32)
i
since this term appears in Equation (5.29) raised to the power n. The right hand inequality in
Equation (5.32) is trivially satisfied, since t; i and are all positive, and the left hand inequality
gives
ti 2 or ti (1 ; 2) 2 (5.33)
1 + ti
A consequence of Equation (5.33) is that the scheme is unconditionally stable if 12 1.
For < 12 the stability criterion is
ti < 1 ;2 2 (5.34)

If the exponential approximation given by Equation (5.31) is negative for any i the solution
will contain components which change sign with each time step n. This oscillatory noise can be
avoided by choosing
t < (1 ; 1) ; (5.35)

max
A
where max is the largest eigenvalue in the matrix , but in practice this imposes a limit which
is too severe for t and a small amount of oscillatory noise, associated with the high frequency
vibration modes of the system, is tolerated. Alternatively the oscillatory noise can be filtered out
by averaging.
These theoretical results are explored numerically with a Crank-Nicolson-Galerkin scheme
( = 21 ) in Figure 5.3, where the one-dimensional diffusion equation
@u = D @ 2 u on 0x1
@t @x2
(5.36)
subject to initial conditions u (x; 0) = 0
and boundary conditions u (0; t) = 0; u (1; t) = 1
is solved for various time increments (t) and element lengths (x) for both linear and cubic
Hermite elements.
Decreasing x from 0:25 to 0:1 with linear elements produces more oscillation because the
system has more degrees of freedom and leads to greater oscillation. At a sufficiently small t the
oscillations are negligible (bottom right, Figure 5.3). With this value of t (0:01 s) the numerical
results agree well with the exact solution (top, Figure 5.3) given by
X1
u (x; t) = x + (;n1) e;n22 t sin (nx)
2 n
(5.37)
n=1
5.4 Mass lumping

A technique known as mass lumping is sometimes used in which the mass matrix Mis replaced
by a diagonal matrix having diagonal terms equal to the row sums. For example, consider the mass
5.4 M ASS LUMPING 107
u(x; t) Exact solution

x x Linear CNG with x = 0:1; t = 0:01
1:0
t = 0:05 x
t = 0:02 x
0:5 x
t = 0:1 x
x
x
t=1 x
x
x x
x
x x
x
x
xx x
x x x
0 1
u u
(b) x = 0:75 (c) x = 0:75
x x x x x x
x x x x x x x
x x x
x linear elements x cubic elements
= 0 25
x : = 0 25
x :
=01t :
t =01
t :
t
0 1 0 1
x
u (d) x = 0:9 u (e) x = 0:9
x
x
x
x x x x x x x
x x
x
x
x
x
linear elements x linear elements
x =01x :
x =01
x :
=01t :
t = 0 01
t :
t
0 1 0 1
F IGURE 5.3: Analytical and numerical solutions of the transient 1D heat equation showing the
effects of element size x and time step size t. The top graph shows the exact and approximate
solutions as functions of x at various times. The lower graphs show the solution through time at the
specified x positions and with various choices of x and t as indicated.
matrix ((5.23)) for a bilinear element (see Figure 1.9 and (1.6)).
ZZ 3
3

M11 = (1 ; 1)2 (1 ; 2 )2 12 = ; (1 ;31) 1 ; (1 ;32) 1 = 13 13 = 19
ZZ 1 1 1
0 0
M22 = 12 (1 ; 2)2 12 = 3 3 = 9 and similarly M33 and M44 .
ZZ 1 1 1 1
M12 = 1 (1 ; 1) (1 ; 1) 12 = ; 2 ; 3 3 = 18
2
ZZ 1 and similarly M and M .
M13 = (1 ; 1)2 2 (1 ; 2)12 = 18 34 24
ZZ 1 and similarly M .
M14 = 1 (1 ; 1) 2 (1 ; 2)12 = 36 23
21 1 1 1 3 21 0 0 0
3
6
therefore M = 6
9
1
18
1
18
1
36
1 77 ;! 6604 4 0 07
1 7
4118 9 36 18
5 40 0 4 05
mass lumping
1 1 1 1
18 136 9 18
1
36
1
18
1
18
1
9 0 0 0 41
The element mass is effectively lumped at the element vertices. Such a scheme has computa-
tional advantages when = 0 in Equation (5.26) because each component of the vector n+1 is u
obtained directly without the need to solve a set of coupled equations. This explicit time integration
scheme, however, is only conditionally stable (see (5.34)) and suffers from phase lag errors - see
below. For evenly spaced elements the finite element scheme with mass lumping is equivalent to
finite differences with central spatial differences.
In Figure 5.4, the finite element and finite differences (lumped f.e. mass matrix) solutions of the
one-dimensional advection-diffusion equation (5.20) with V = 5 m s ;1 , D = 0:1 m2 s;1 , f = 0
are compared for the propogation and dispersion of an initial unit mass pulse at x = 0. The length
of the solution domain is sufficient to avoid reflected end effects.
The exact solution is a Gaussian distribution whose variance increases with time:
(x ; V t)2
u (x; t) = p M e 4Dt
;
(5.38)
4t
The finite element solution, using the Crank-Nicolson-Galerkin technique, shows excellent
amplitude and phase characteristics when compared with the exact solution. The finite difference,
or lumped mass, solution also using centered time differences, reproduces the amplitude of the
pulse very well but shows a slight phase lag.
5.5 CMISS Examples

1. To solve for the transient heat flow in a plate run CMISS example 331
2. To investigate the stability of time integration schemes run CMISS examples 3321 and 3322.
u (x; t)
V = 5 m s;1 D = 0:1 m2 s;1

t = 0:01 s
Exact solution
x x Finite element solution
t = 0:05 s
o o Finite difference solution
x
x x t = 0:2 s
t = 0:4 s
x
x x
x x ox
x x o o
o x oxx
x x x o x o x
x
x ox x o x
F IGURE 5.4: Advection-diffusion of a unit mass pulse. The finite element solutions (at t=0:01 s,
0:05 s, 0:2 s and 0:4 s) and finite difference solutions (at t=0:4 s only) are compared with the exact
solution. x= 0.1, t = 0:001 s for 0< t <0:01 s and t = 0.01 for t 0:01 s.
Chapter 6
Modal Analysis
6.1 Introduction
The system of ordinary differential equations which results from the application of the Galerkin
finite element (or other) discretization of the spatial domain to linear parabolic or hyperbolic equa-
tions can either be integrated directly - as in the last section for parabolic equations - or analysed
by mode superposition. That is, the time-dependent solution is expressed as the superposition of
the natural (or resonant) modes of the system. To find these modes requires the solution of an
eigenvalue problem.
6.2 Free Vibration Modes

Consider an extension of Equation (5.22) which includes second order time derivatives (e.g., nodal
point accelerations)
M u (t) + C u_ (t) + Ku (t) = f (t) (6.1)
M ; C and K are the mass, damping and stiffness matrices, respectively, f (t) is the external load
vector and u (t) is the vector of n nodal unknowns. In direct time integration methods u (t) and
u_ (t) are replaced by finite differences and the resulting system of algebraic equations is solved at
successive time steps. For a small number of steps this is the most economical method of solution
but, if a solution is required over a long time period, or for a large number of different load vectors
f (t), a suitable transformation
u (t) = Px (t) (6.2)
applied to Equation (6.1) can result in the matrices of the transformed system
P T MP x (t) + P T CP x_ (t) + P T KPx (t) = P T f (6.3)
having a much smaller bandwidth than in the original system and hence being more economical
P
to solve. In fact, if damping is neglected, can be chosen to diagonalize M
and K
and thereby
uncouple the equations entirely. This transformation (which is still applicable when damping is
112 M ODAL A NALYSIS
included but does not then result in an uncoupled system unless further simplications are made) is
found by solving the free vibration problem
M u (t) + Ku (t) = 0 (6.4)
Proof: Consider a solution to Equation (6.4) of the form
u (t) = s sin ! (t ; t0) ; (6.5)
s
where ! and t0 are constants and is a vector of order n. Substituting Equation (6.5) into Equa-
tion (6.4) gives the generalized eigenproblem
Ks = !2Ms (6.6)
having n eigensolutions (! 12 s1 ) ; (!22 s2 ) ; : : : ; (!n2 sn ). If K is a symmetric matrix (as is the case

when the original partial differential operator is self-adjoint) the eigenvectors are orthogonal and
can be “normalized” such that
1 i=j
si Msj = 0 i 6= j
T (6.7)
(the eigenvectors are said to be M -orthonormalised). Combining the n eigenvectors into a matrix
S s s s
= [ 1; 2 ; : : : ; n] - the modal matrix - rewriting Equation (6.7) as
ST MS = I (6.8)
where I is the identity matrix, (6.6) becomes

KS = MS (6.9)
where
22 3
! 0
66 1 !22 77
= 64 ..
.
75 (6.10)
0 !n2
or
ST KS = ST MS = I = (6.11)
Thus the modal matrix - whose columns are the M -orthonormalised eigenvectors of (i.e., K
P
satisfying Equation (6.6)) - can be used as the transformation matrix in Equation (6.2) required
to reduce the original system of equations (6.1) to the canonical form
x (t) + ST CS x_ (t) + x (t) = S T f (t) (6.12)

6.3 A N A NALYTIC E XAMPLE 113
With damping neglected equation Equation (6.12) becomes a system of uncoupled equations
xi (t) + !i2xi (t) = ri (t) i = 1; 2; : : : ; n (6.13)
x
where xi is the ith component of and ri is the ith component of the vector S T f . The solution of
this system is given by the Duhamel integral
1 Z t
xi (t) = ! ri ( ) sin !i (t ; ) d + i sin !it + i cos !id (6.14)
i
0
where the constants i and i are determined from the initial conditions
xi (0)jt=0 = siT M u (0)jt=0

x_ i (0)jt=0 = xi (0)jt=0 = siT M u (0)jt=0 siT M u (0)jt=0
(6.15)
6.3 An Analytic Example

M u + Ku = f where
As an example, consider the equilibrium equations
2 0 6 ;2 0
M = 0 1 ; K = ;2 4 and f = 10
To find the solution by modal analysis we first solve the generalised eigenproblem Ks = ! 2Ms
i.e.,
6 ; 2!2 ;2
;2 2 s=0
4;!
K M
has a solution when det [ ; ! 2 ] = 0 or ! 4 ; 7! 2 + 10 = 0. This characteristic polynomial
s
s
has solutions ! 2 = 2; 5 with corresponding eigenvectors T1 = a 1 1 ; T1 = b ;1 2 . To find
the magnitude of the eigenvectors we use Equation (6.7), i.e.,
2 01 1 2 0;1
a 1 1 0 1 1 =1)a= p
2 b ;1 2 0 1 2 = 1 ) b = p1
2
3 3

2 0 ;1 = 0).
(Notice that the orthogonality condition is satisfied: ab 1 1
1 01 1 2 1 2
The M -orthognormalised eigenvectors are now sT1 = p p and sT2 = ; p p ,
21 3 3 3 6 6
p ;p 7 1
giving the modal matrix S = 4 13
6 6
p 2 5 which, when used as the transformation matrix,
p
3 6
reduces the stiffness matrix to

2 1 13 21 13
p p p ; p
ST KS = 64 13 23 75 ;62 ;42 64 13 2 6 75 = 20 05 =
;p p p p
6 6 3 6
and the mass matrix to
2 1 13 21 13
p p p ; p
S T MS = 64 13 23 75 20 01 64 13 2 6 75 = 10 01 = I
;p p p p
6 6 3 6
Thus the natural modes of the system are given by
213 2 13
p ;p
u1 (t) = 64 13 75 sin p2 (t ; t0) and u2 (t) = 64 2 6 75 sin p5 (t ; t00) :
p p
3 6
The solution of the non-homogeneous system, subject to given initial conditions, is found by solv-
ing the uncoupled equations
2 1 13 2 10 3
2 0 6 p p 7 0
6 p7
x (t) + 0 5 x (t) = 4 1 2 5 10 = 4 203 5
3 3
;p p p
6 6 6
r
by means of the Duhamel integral (6.14) (in this case with constant) and then, from Equation (6.2)
P S s s
with = [( 1 ; 2 ; : : : ; n ] s
X
n
u (t) = Sx (t) = sixi (t) (6.16)
i=1
Notice that the solution is expressed in Equation (6.16) as the superposition of the natural
modes (eigenvectors) of the homogeneous equations. If the forcing function (load vector) is close
to one of these modes the corresponding coefficient xi will be large and will dominate the response
- if it coincides then resonance will occur. Very often it is unnecessary to evaluate all n eigenvectors
of the system; the higher frequency modes can be ignored and the solution adequately represented
by superposition of the p eigenvectors associated with the p lowest eigenvalues, where p < n.
6.4 Proportional Damping

When element damping terms are included in the original dynamic equations (6.1) the transforma-
S
tion to a lower bandwidth system is still based on the model matrix but Equation (6.12) is then
not a system of uncoupled equations. One simplification often made in order to retain the diago-
nal nature of Equation (6.12) is to approximate the overall energy dissipation of the finite element
system with proportional damping
siT Csj = 2!iiij ; (6.17)
where i is a modal damping parameter and ij is the Kronecker delta. Equation (6.12) now reduces
to n equations of the form
xi (t) + 2!iix_ i (t) + !12x1 (t) = ri (t) (6.18)
with solution (the Duhamel integral)
1 Z t
xi (t) = ! ri ( ) :ei!i (t; ) : sin !i (t ; ) dt + e;i!it fi sin !it + i cos !itg (6.19)
i
0
p
where !i = !i 1 ; i2 . i and i are calculated from the initial conditions Equation (6.15).
Once the components xi (t) have been found from Equation (6.19) (or alternative time integration
methods applied to (6.18)), the solution u (t) is expressed as a superposition of the mode shapes
si by Equation (6.16).
6.5 CMISS Examples
1. To analyse a plane stress modal analysis run CMISS example 451
2. To analyse a clamped beam modal analysis run CMISS example 452
3. To analyse a steel-framed building modal analysis run CMISS example 453

Chapter 7
Domain Integrals in the BEM
7.1 Achieving a Boundary Integral Formulation

The principal advantage of the BEM over other numerical methods is the ability to reduce the
problem dimension by one. This property is advantagous as it reduces the size of the solution
system leading to improved computational efficiency. This reduction of dimension also eases the
burden on the engineer as it is only necessary to construct a boundary mesh to implement the BEM.
To achieve this reduction of dimension it is necessary to formulate the governing equation as a
boundary integral equation. To achieve a boundary integral formulation it is necessary to find an
appropriate reciprocity relationship for the problem and to determine an appropriate fundamental
solution. If either of these requirements cannot be satisfied then a boundary integral formulation
cannot be achieved. The most common difficulty in applying the BEM is in determining an appro-
priate fundamental solution.
A linear differential equation can be expressed in operator form as Lu = where L is a linear
operator, is an inhomogeneous source term and u is the dependent variable. The fundamental
solution for this equation is a solution of
L! (x; ) + () = 0 (7.1)
where * indicates the adjoint of the operator L and is the Dirac delta function. No specific
boundary conditions are prescribed but in some cases regularity conditions at infinity need to be
satisfied. The fundamental solution is a Green’s function which is not required to satisfy any
boundary conditions and is therefore also commonly termed the free-space Green’s function.
The mathematical theory required to determine the fundamental solution of a constant coef-
ficient PDE is well-developed and has been used successfully to determine the fundamental so-
lutions for a wide range of constant coefficient equations (Brebbia & Walker 1980) (Clements &
Rizzo 1978) (Ortner 1987). Fundamental solutions are known and have been published for many
of the most important equations in engineering such as Laplace’s equation, the diffusion equation
and the wave equation (Brebbia, Telles & Wrobel 1984a). However, by no means can it be guaran-
teed that the fundamental solution to a specific differential equation is known. In particular, PDEs
with variable coefficients do not, in general, have known fundamental solutions. If the fundamental
solution to an operator cannot be found then domain integrals cannot be completely removed from
the integral formulation. Domain integrals will also arise for inhomogeneous equations.
118 D OMAIN I NTEGRALS IN THE BEM
Wu (1985) argued that the BEM has several advantages over other numerical methods which
justify its use for many practical problems - even in cases where domain integration is required.
He argued that for problems such as flow problems a wide range of phenomena are described by
the same governing equations. What distinguishes these phenomena is the boundary conditions of
the problem. For this reason accurate description of the boundary conditions is vital for solution
accuracy. The BEM generates a formulation involving both the dependent variable u and the flux
q. This allows flux boundary conditions to be applied directly which cannot be achieved in either
the finite element or finite difference methods.
Another advantage of the BEM over other numerical methods is that it allows an explicit ex-
pression for the solution at an internal point. This allows a problem to be subdivided into a number
of zones for which the BEM can be applied individually. This zoning approach is suited to prob-
lems with significantly different length scales or different properties in different areas.
Domain integration can be simply and accurately performed in the BEM. However, the pres-
ence of domain integrals in the BEM formulation negates one of the principal advantages of the
BEM in that the problem dimension is no longer reduced by one. Several methods have been de-
veloped which allow domain integrals to be expressed as equivalent boundary integrals. In this
section these methods will be discussed.
7.2 Removing Domain Integrals due to Inhomogeneous Terms

Inhomogeneous PDEs occur for a large number of physical problems. An inhomogeneous term
may arise due to a number of factors including a source term, a body force term, or due to ini-
tial conditions in time-dependent problems. An inhomogeneous linear PDE can be expressed in
operator form as Lu = where is a known function of position or a non-zero constant. If the
fundamental solution is known for the operator L, the resulting BEM formulation will be
Z
Hu ; Gq = ; ! d
(7.2)

The domain integral in this formulation does not involve any unknowns so domain integration can
be used directly to solve this equation. This requires discretising the domain into internal cells
in much the same way as for the finite element method. As the domain integral does not involve
any unknown values accurate results can generally be achieved using a fairly coarse mesh. This
method is simple and has been shown to produce accurate results (Brebbia et al. 1984a). This
approach, however, requires a domain discretisation and a numerical domain integration procedure
which reduces the attraction of the BEM over domain-based numerical methods.
7.2.1 The Galerkin Vector technique

For some particular forms of the inhomogeneous function the domain integral can be transformed
directly into boundary integrals.
Consider the Poisson equation r2 u = . Applying the BEM gives an equation of the form of
7.2 R EMOVING D OMAIN I NTEGRALS DUE TO I NHOMOGENEOUS T ERMS 119
Equation (7.2). Using Green’s second identity

Z ; Z @v @
r v ; vr d
=
2 2 @n ; v @n d; (7.3)

;
domain integration can be avoided for certain forms of . If a v can be found which satisfies
r2v = !, where ! is the fundamental solution of Laplace’s equation, then for the specific case of
being harmonic (r2 = 0) Green’s second identity can be reduced to
Z Z @v @
! d
= @n ; v @n d; (7.4)

;
Therefore if a Galerkin vector can be found and is harmonic the domain integral in Equation (7.2)
can be expressed as equivalent boundary integrals.
Fairweather, Rizzo, Shippy & Wu (1979) determined the Galerkin vector for the two-dimensional
Poisson equation and Monaco & Rangogni (1982) determined the Galerkin vector for the three-
dimensional Poisson equation. Danson (1981) showed how this method can be applied successfully
for a number of physical problems involving linear isotropic problems with body forces. He con-
sidered the practical cases where the body force term arose due to either a constant gravitational
load, rotation about a fixed axis or steady-state thermal loading. In each of these cases the domain
integral can be expressed as equivalent boundary integrals.
This Galerkin vector approach provides a simple method of expressing domain integrals as
equivalent boundary integrals. Unfortunately, it only applies to specific forms of the inhomoge-
neous term (i.e., is required to be harmonic).
7.2.2 The Monte Carlo method

Domain discretisation could be avoided by using a Monte Carlo technique. This technique approx-
imates a domain integral as a sum of the integrand at a number of random points. Specifically, in
two dimensions, a domain integral I is approximated as
A XN
I N f (xi ; yi) (7.5)
i=1
where f (xi ; yi ) is the value of the integrand at random point (x i ; yi ), N is the number of random
points used and A is the area of the region over which the integration is performed. This approxi-
mation allows a domain integral to be approximated by a summation over a set of random points
so domain integration can be performed without requiring a domain mesh. This method has the
secondary advantage of allowing the integration to be performed over a simple geometry enclosing
the problem domain - if a random point is not in the problem domain its contribution is ignored.
The method was proposed by Gipson (1987). Gipson has successfully applied this method to a
number of Poisson-type problems. Unfortunately this method often proves to be computationally
expensive as a large number of integration points are needed for accurate domain integration.
Gipson argues however that, as this method removes the burden of preparing a domain mesh,
the extra computational expense is justified.
7.2.3 Complementary Function-Particular Integral method

A more general approach can be developed using particular solutions. Consider the linear problem
Lu = . u can be considered as the sum of the complementary function uc, which is a solution of
the homogeneous equation Luc = 0, and a particular solution up which satisfies Lup = but is
not required to satisfy the boundary conditions of the problem. Applying BEM to the governing
equation using the expansion u = u c + up gives
Hu ; Gq = Hup ; Gqp (7.6)
If a particular solution up can be found, all values on the right-hand-side of Equation (7.6) are
known - reducing the problem to
Hu ; Gq = d (7.7)
d
where is a vector of known values. This linear system can be solved by applying boundary
conditions.
This approach can be applied in a situation where an analytic expression for a particular solu-
tion can be found. Unfortunately particular solutions are generally only known for simple operators
and for simple forms of . Alternatively an approximate particular solution could be calculated nu-
merically. Zheng, Coleman & Phan-Thien (1991) proposed a method where a particular solution
is determined by approximating the inhomogeneous source term using a global interpolation func-
tion. This approach is a special case of a more general method known as the dual reciprocity
boundary element method.
7.3 Domain Integrals Involving the Dependent Variable

Consider the linear homogeneous PDE Lu = 0. For many operators the fundamental solution to
the operator L may be unobtainable or may be in an unusable form. This is especially likely if L
involves variable coefficients for which case it has been shown that it is particularly difficult to find
a fundamental solution. Instead, a BEM formulation can be derived based on a related operator L ^
with known fundamental solution. A BEM formulation for Lu = 0 based on the operator L will ^
be of the form
Z
Hu ; Gq = ; ^ L ; L u! d
(7.8)

where ! is the fundamental solution corresponding to the operator L ^ . This integral equation is
similar to Equation (7.2). However in this case the domain integral term involves the dependent
variable u. This problem could be solved using domain integration where the internal nodes are
treated as formal problem unknowns.
7.3 D OMAIN I NTEGRALS I NVOLVING THE D EPENDENT VARIABLE 121
7.3.1 The Perturbation Boundary Element Method

Rangogni (1986) proposed solving variable coefficient PDEs by coupling the boundary element
method with a perturbation method. He considered the two-dimensional generalised Laplace equa-
tion
r ( (x; y) rV (x; y)) = 0 (7.9)
1;
Using the substitution V (x; y ) = 2 u (x; y ) Equation (7.9) can be recast as a heterogeneous
Helmholtz equation
r2 u + f (x; y) u = 0 (7.10)
where f is a known function of position.

Rangogni treated this equation as a perturbation about Laplace’s equation. He considered the
class of equations
r2 u + "f (x; y) u = 0 where 0 " 1 (7.11)
for which he sought a solution of the form

X
1
u = u0 + "u1 + " u2 + : : : =
2 uj "j (7.12)
j =0
Substituting Equation (7.12) into Equation (7.11) and grouping powers of " gives
; ;
r2 u0 + " r2u1 + fu0 + "2 r2 u2 + fu1 + : : : = 0
(7.13)
A solution will only exist for all values of " if the terms at each power of " equal zero. This allows
Equation (7.13) to be treated as an infinite series of distinct problems which can be solved using
the boundary element method. u0 can be found by solving r2 u0 = 0 which Rangogni assumes
will satisfy the boundary conditions of the original problem. Each successive u j can then be found
by solving a Poisson equation with homogeneous boundary conditions as u j ;1 has been previously
determined. Rangogni used a domain discretisation to solve these Poisson problems.
Equation (7.10) is a particular member of this family of equations for which " = 1. The
X
1
solution to Equation (7.10) is therefore given by uj . Rangogni reported that in practice this
j =0
series converged rapidly and in his numerical examples he achieved accurate results using only u 0
and u1 .
Rangogni (1991) extended this coupled perturbation - boundary element method to the general
second-order variable coefficient PDE
@u + g (x; y) @u = h (x; y)
r2 u + f (x; y) @x (7.14)
@y
He considered the family of equations

@u + g (x; y) @u = h (x; y) (0 " 1)
r u + " f (x; y) @x
2
@y (7.15)
Applying the perturbation method to this family of equations allows Equation (7.15) to be ex-
pressed as an infinite series of distinct Poisson equations which can be solved using the boundary
element method. Again Rangogni used an domain mesh to solve these Poisson equations. Ran-
gogni found that in practice convergence was rapid and accurate results were produced.
Gipson, Reible & Savant (1987) considered a class of hyperbolic and elliptic problems which
can be transformed into an inhomogeneous Helmholtz equation. They used the perturbation method
to recast this as an infinite sequence of Poisson equations. They avoided domain discretisation by
using a Monte Carlo integration technique (Gipson 1987) to evaluate the required domain integrals.
Lafe & Cheng (1987) used the perturbation method to solve steady-state groundwater flow
problems in heterogeneous aquifers. They showed the method produced accurate results for sim-
ply varying hydraulic conductivities with convergence after two or three terms. Lafe & Cheng
investigated the convergence of the perturbation method. They found that for rapidly varying hy-
draulic conductivity convergence is not guaranteed. From this investigation they concluded that
accurate results can be obtained so long as the hydraulic conductivity does not vary by more than
one order of magnitude within the solution domain. If the hydraulic conductivity variation is more
significant they recommend using the perturbation method in conjunction with a subregion tech-
nique so that the variation of conductivity within each subregion satisfies their requirements. This
process could become computationally expensive, particularly if convergence is not rapid, as the
solution of multiple subproblems will be required within each subregion.
7.3.2 The Multiple Reciprocity Method

The multiple reciprocity method (MRM) was initially proposed by Nowak (1987) for the solution
of transient heat conduction problems. Since then the method has been successfully applied to a
wide range of problems. The MRM can be viewed as a generalisation of the Galerkin vector ap-
proach. Instead of using one higher-order fundamental solution, the Galerkin vector, to convert the
remaining domain integrals to equivalent boundary integrals a series of higher-order fundamental
solutions is used.
Consider the Poisson equation
r2u = b0 (7.16)
x
where b0 = b0 ( ) is a known function of position. Applying BEM to this equation, using the
fundamental solution to the Laplace operator, gives
Z @!0 Z Z @u
c () u () + u @n d; + b0 !0 d
= !0 @n d; (7.17)
;
;

where !0 is the known fundamental solution to Laplace’s equation applied at point . To avoid
domain discretisation the domain integral in Equation (7.17) needs to be expressed as equivalent
boundary integrals. Using MRM this is achieved by defining a higher-order fundamental solution
!1 such that
r2!1 = !0 (7.18)
Using this higher-order fundamental solution the domain integral in Equation (7.17) can be written
as
Z Z
b0 !0 d
= b0 r2 !1 d
(7.19)

or
Z Z @!1 @u Z
b0 !0 d
= u @n ; !1 @n d; + !1r2 b0 d
(7.20)

;
This formulation has generated a new domain integral. b 0 is a known function so we can introduce
a new function b1 which can be determined analytically from the relationship
b1 = r2 b0 (7.21)
giving
Z Z
!1 r b0 d
=
2 !1b1 d
(7.22)

This process can be repeated by introducing a new higher-order fundamental solution ! 2 such that
r2!2 = !1 (7.23)
and continuing until convergence is reached.

This procedure is based on the recurrence relationships
bj+1 = r2 bj for j = 0; 1; 2; : : : (7.24)

r2 !j+1 = !j for j = 0; 1; 2; : : : (7.25)
Using these recurrence relationships gives the boundary integral formulation

Z @!0 @u X1 Z
@! @b

c () u () + u @n ; !0 @n d; + j +1 j
bj @n ; !j+1 @n d; = 0 (7.26)
; j =0 ;
which is an exact formulation if the infinite series converges. Errors are only introduced at the
stage of boundary discretisation.
Introducing interpolattion functions and discretising the boundary gives the matrix system
X
1
H0u ; G0q = (HJ +1pj ; GJ +1rj ) (7.27)
j =0
H G
where J +1 and J +1 are influence coefficient matrices corresponding to the higher-order fun-
p r
damental solutions and j and j contain the nodal values of bj and its normal derivative.
The MRM can be applied based on operators other than the Laplace operator. This approach
relies on knowledge of the higher-order fundamental solutions necessary for application of the
method. These solutions have been determined and successfully used for the Laplace operator
in both two and three dimensions but the extension of the method to other equation types needs
further research. Itagaki & Brebbia (1993) have determined the higher order fundamental solutions
for the two-dimensional modified Helmholtz equation.
The MRM can be extended to other equations by allowing the forcing function b 0 to be a general
x
function such that b0 = b0 ( ; u; t). The MRM will be restricted to cases where the recurrence
relationships - Equations (7.24) and (7.25) - can be employed. Brebbia & Nowak (1989) have
applied the MRM to the Helmholtz equation r 2 u + 2 u = 0 where b0 = ;2 u and the recurrence
relationship defined by Equation (7.24) becomes simply
uj+1 = r2 uj = ;2j u (7.28)
In this case the boundary integral formulation will be

1 Z
X @! @u

c () u () + j
u @n ; !j @n d; = 0
2j (7.29)
j =0 ;
7.3.3 The Dual Reciprocity Boundary Element Method

Equation Derivation
The dual reciprocity boundary element method (DR-BEM) was developed to avoid the need for
domain integration in cases where the fundamental solution of the governing differential equation
is unknown or is impractical to apply. Instead the DR-BEM is applied using an appropriate related
operator with known fundamental solution. The most common choice is the Laplace operator
(Partridge, Brebbia & Wrobel 1992) and in this chapter the DR-BEM will be illustrated for this
choice.
Consider a second-order PDE which can be expressed in the form
r2u = b (7.30)
x
The forcing function b can be completely general. If b = b ( ) then b is a known function of posi-
tion and the differential equation described is simply the Poisson equation. For potential problems
x x
b = b ( ; u) and for transient problems b = b ( ; u; t). Applying the BEM to Equation (7.30) will
give
Z
Hu ; Gq = ; b! d
(7.31)

where ! is the known fundamental solution to Laplace’s equation. The aim of the DR-BEM is to
express the domain integral due to the forcing function b as equivalent boundary integrals.
The DR-BEM uses the idea of approximating b using interpolation functions. A global approx-
imation to b of the form
X
M
b= j fj (7.32)
j =1
is proposed. j are unknown coefficients and fj are approximating functions used in the interpo-
lation and are generally chosen to be functions of the source point and the field point of the fun-
damental solution. The approximating functions f j are applied at M different collocation points
- called poles - generally most, but not all, of which are located on the boundary of the problem
domain.
As discussed in the previous chapter the solution to a linear PDE Lu = can be constructed as
the sum of a complimentary function u c (which satisfies the homogeneous equation Lu c = 0) and a
particular solution up to the equation Lup = . Instead of using a single particular solution, which
may be difficult to determine, the DR-BEM employs a series of particular solutions u ^ j which are
related to the approximating functions f j as shown in Equation (7.33).
r2 u^j = fj j = 1; : : : ; M (7.33)
By substituting Equations (7.32) and (7.33) into Equation (7.30) the forcing function b is approxi-
mated by a weighted summation of particular solutions to the Poisson equation.
X
M
r2u = j r2u^j (7.34)
j =1
The DR-BEM essentially constructs an approximate particular solution to the governing PDE as a
summation of localised particular solutions.
With the governing equation rewritten in the form of Equation (7.34) the standard boundary
element approach can be applied. Equation (7.34) is multiplied by a weighting function ! and
integrated over the domain. Green’s theorem is applied twice and the fundamental solution of the
Laplacian is used to remove the remaining domain integrals. The name dual reciprocity BEM is
derived from the application of reciprocity relationships to both sides of Equation (7.34). After
applying these steps Equation (7.35) is obtained, where the fundamental solution pole is applied at

point .
Z @! @u
c () u () + u @n ; ! @n d;
;
0 Z 1
X
M
= j @c () uj () + u^j @!
@n ; ! @hatuj d;A
@n (7.35)
j =1 ;
In implementing a numerical solution of this equation similar steps are taken as for the standard
BEM. The boundary is discretised into elements and interpolation functions are introduced to ap-
proximate the dependent variable within each element.
The form of each u^j is known from Equation (7.33) once the approximating functions f j have
been defined. It is not necessary to use interpolation functions to approximate each u^ j . However
by using the same interpolation functions to approximate u and u ^ j the numerical implementation
will generate the same matrices H G
and on both sides of Equation (7.35). The error generated
by approximating each u ^j in this manner has been found to be small and can be justified by the
improved computational efficiency of the method (Partridge et al. 1992).
The application of this method results in the system
X
N +I
Hu ; Gq = j (H u^j ; Gq^j ) (7.36)
j =1
where the M poles were chosen to be the N boundary nodes plus I internal points so that M =
N + I . Although it is not generally necessary to include poles at internal points it has been found
that in general improved accuracy is achieved by doing so (Nowak & Partridge 1992). It has
been shown that for many problems (Partridge et al. 1992) (Huang & Cruse 1993) using boundary
points only in this procedure is insufficient to define the problem. In general using internal points
is likely to improve the solution accuracy as it increases the number of degrees of freedom. No
theory has been developed of how many internal collocation points should be used for optimal
accuracy, or where these points should be positioned within the problem domain. Using internal
poles in this interpolation does not require domain discretisation - it is only necessary to specify
the coordinates of the internal collocation points. The internal points can be chosen to be locations
where the solution is of interest.
The ^ j and ^j vectors can be treated as columns of the matrices ^ and ^ respectively. This
u q U Q
allows Equation (7.36) to be rewritten as

Hu ; Gq = H U^ ; GQ^ (7.37)

where is a vector containing the nodal values of . To solve this system it is necessary to evaluate

. is defined by Equation (7.32) which, for the nodal values, can be expressed in matrix form as
b F F
= . If the matrix is nonsingular this expression can be rearranged to give Equation (7.38)
which provides an explicit expression for .
= F ;1b (7.38)
in Equation (7.37) gives

Including this explicit expression for

Hu ; Gq = H U^ ; GQ^ F ;1b (7.39)
The approach taken to solve this equation will depend on the form of b.
The Approximating Function f

The accuracy of the DR-BEM hinges on the accuracy of the global approximation to the forcing
function b (defined by Equation (7.32)). Therefore the choice of the approximating functions f j is
a key consideration when implementing the DR-BEM. The only requirement so far prescribed on
F
the form of the approximating functions f j is that the matrix generated should be nonsingular
and that the related particular solutions u
^ j can be determined and can be expressed in a practical
closed form. Some work has been conducted into investigating what form of f j should be used in
a given situation to provide the highest accuracy and computational efficiency.
Usually a form of fj is defined and this can be used, applying Equation (7.33), to specify u ^
and q^. The fundamental solution of Laplace’s equation is ! ( x; ) = ; 21 ln r in two-dimensional
x
space and ! ( ; ) =
1
4r in three-dimensional space - where r is the Euclidean distance between
x
the field point and the source point of the fundamental solution. Due to the dependence of this
fundamental solution only on r the approximating function is generally chosen to be some radial
function i.e., fj = fj (r ). Several other options for fj have been tried (Partridge et al. 1992) but it
has been found that in general the most accurate results were generated using some radial function.
For both two and three-dimensional problems Wrobel, Brebbia & Nardini (1986) recommended
choosing fj from the series
fj = 1 + rj + rj2 + : : : + rjm (7.40)
where rj is the distance between the field point (node j ) and the DR-BEM collocation point (node
i). They showed that accurate results can be achieved using some combination of terms from this
series. Generally the same approximating function f j is used at all the collocation points so in this
thesis, for simplicity, the form of approximating functions f j will be referred to by a single f .
Choosing f to be a function of only one variable simplifies the process of determining u ^ and q^.
For two-dimensional problems, if f = f (r ) then the relationship
r2 u^ = f (r) (7.41)
can be reduced to the ordinary differential equation
d2u^ + 1 du^ = f (7.42)

dr2 r dr
Using f defined by Equation (7.40) the corresponding forms of u^ and q^, for two-dimensional
problems, can be shown to be
u^ = r4 + r9 + : : : + r 2
2 3 m+2
(7.43)
@x @y (m + 2)
1 r r m
q^ = rx @n + ry @n 2 + 3 + : : : + m + 2 (7.44)
where rx = xj ; xi and ry = yj ; yi .
Any combination of terms from Equation (7.40) can be used for specifying f . It has been found
that in general including higher-order terms leads to little improvement in accuracy (Partridge
et al. 1992). The most commonly used form is f = 1 + r as this approximation will generally give
accurate results with greater computational efficiency than other choices.
Equation (7.40) was recommended as a basis for the approximating function f due to the
particular form of the fundamental solution of Laplace’s equation and its dependence on r only. If
a different operator is used as the basis of the DR-BEM then it is likely a different form of f will
be more appropriate. The choice of f in this case will be discussed in Section 7.3.3.
The performance of the DR-BEM hinges on the choice of the approximating function f . The
theory of how to determine the best approximating function is therefore a vital component of
the DR-BEM. Unfortunately the approximating function has generally been chosen and used in a
rather ad-hoc manner. Recently some more formal analysis of the use of approximating functions
has been undertaken.
Golberg & Chen (1994) argued that a formal analysis of the approximating function f can be
undertaken using the theory of radial basis functions. Radial basis functions are a generalisation
of cubic splines in multi-dimensions. Cubic splines are known to be optimal for one-dimensional
interpolation. Therefore, rather than being an arbitrary choice, it seems that choosing f to be
a radial function is a logical extension for two or three-dimensional problems. Golberg & Chen
showed that, for the Poisson equation, choosing f to be a radial basis function ensures convergence
of the DR-BEM.
They also demonstrated that f = 1 + r is a specific member of the group of radial basis
functions. The theory of using radial basis functions for multi-dimensional approximation is fairly
advanced. It has been shown that f = r is optimal for three-dimensional problems which justifies
the use of f = 1 + r when applying the DR-BEM to three-dimensional problems - the constant
F
is included to ensure a non-zero diagonal for . However for two-dimensional problems it has
been shown that optimal approximation is attained using the thin plate spline f = r 2 logr . This
observation lead Golberg & Chen to suggest that choosing f to be a thin plate spline may improve
the accuracy of the DR-BEM in two dimensions. Recently Golberg (1995) has published a review
of the DR-BEM concentrating on developments since 1990 concerning the numerical evaluation
of particular solutions.
Inhomogeneous Equations
If the forcing function b is a function of position only then the differential equation under consid-
F
eration is simply Poisson’s equation. In this case it is not necessary to invert the matrix as can
simply be calculated from =b F using Gaussian elimination. Equation (7.39) can be rewritten
as
^ ^
Hu ; Gq = d where d = H U ; GQ (7.45)
By applying boundary conditions Equation (7.45) can be reduced to a linear system Ax

=
which can be solved to give the unknown nodal values of u and q .
Zheng et al. (1991) and Coleman, Tullock & Phan-Thien (1991) have proposed a method which
uses a global shape function to construct an approximate particular solution. As discussed by
Polyzos, Dassios & Beskos (1994) this method is essentially equivalent to the DR-BEM. However,
Zheng et al. and Coleman et al. suggested several alternative ways of determining the unknown
coefficients j for inhomogeneous equations. Zheng et al. (1991) used a least-squares method
where they minimised the sum of squares
X
M X
N !
S= b (rm ) ; j fj (rm) (7.46)
m=1 j =1
using singular value decomposition. For large systems they found the computational efficiency
could be improved by employing the conjugate gradient method. Coleman et al. (1991) success-
fully solved inhomogeneous potential and elasticity problems which are governed by operators
other than the Laplacian.
Elliptic Problems

If b is a function of the dependent variable then will also be a function of the dependent variable.
Consider, for example, the linear second-order differential equation
r2 u + u = 0 (7.47)
F ;u
In this case b = ;u so = ;1 . Applying the DR-BEM to Equation (7.47), based on the
fundamental solution to Laplace’s equation, gives

Hu ; Gq = ; H U^ ; GQ^ F ;1u (7.48)
which can be rearranged to give

^ ^ ;1
(H + S ) u = Gq where S = H U ; GQ F (7.49)
Again, by applying boundary conditions Equation (7.49) can be reduced to a linear system AX =
which can be solved to determine the unknown nodal values.
F
Due to the presence of the fully-populated ;1 matrix in Equation (7.49) it is not possible to
solve the boundary problem and internal problem separately. Instead the solution can be treated as
a coupled problem and the solutions at boundary and internal nodes are generated simultaneously.
Derivative Terms The DR-BEM can also be applied for elliptic problems where b involves
derivatives of the dependent variable (Partridge et al. 1992). Consider, for example, the differ-
ential equation
r2u + @u
@x = 0 (7.50)
In this case applying DR-BEM, using the Laplace fundamental solution, gives

Hu ; Gq = ; H U^ ; GQ^ F ;1 @u
@x (7.51)
To solve this problem it is necessary to relate the nodal values of u to the nodal values of
@u
@x . This
u
is achieved by using interpolation functions to approximate in a similar manner as was used to
approximate b in Equation (7.32). A global approximation function of the form
X
M
u= j (x; y) j (7.52)
j =1
can be used to approximate u where j are the chosen interpolation functions and j are the un-
known coefficients. In system form this can be expressed as
u = (7.53)
F
Although it is not necessary, equating to improves the computational efficiency of the method
as only one matrix inversion procedure is required. Differentiating Equation (7.53) gives
@u = @
@x @X (7.54)
F
Choosing = and inverting Equation (7.53) to give an explicit expression for allows Equa-
tion (7.54) to be rewritten as
@u = @F F ;1u
@x @X (7.55)
Equation (7.39) can now be rewritten as

@F ;1
(H + R) u = Gq where R = H U^ ; GQ^ F ;1 @X F (7.56)
By applying boundary conditions Equation (7.56) can be reduced to a linear system which can be
solved to give the unknown nodal values.
As mentioned earlier, the approximating function f is generally chosen to be f = 1 + r . This
can lead to numerical problems if derivative terms are included in the forcing function b. As shown
in Equation (7.55) derivative terms require derivatives of f to be evaluated. For example, evaluating
@F matrix requires calculation of @f . Using the approximating function f = 1 + r gives

the
@X @x
@f = @f @r = @f rx (7.57)
@x @r @x @r r
This derivative function can become singular, so - as shown by Zhang (1993) - significant numerical
error may result. This will especially be the case in problems where collocation points are located
close together.
Zhang (1993) suggested two possibilities for avoiding this problem. The first suggestion in-
volved using a mapping procedure to map the governing equation to an equation without convec-
tive terms. This method was shown to produce accurate results but is somewhat cumbersome and
can only be applied to linear problems. A simpler approach is to choose an approximating func-
tion which does not lead to singularities for convective terms. Zhang recommended use of either
f = 1 + r3 or f = 1 + r2 + r3. These approximating functions produce accurate results and can be
simply applied for both linear and nonlinear problems. Zhang recommended the adoption of these
approximating functions for all use of the DR-BEM.
The same idea of using Equation (7.53) to allow nodal values of u to be associated to its
derivatives can be applied to extend the DR-BEM to cases involving higher-order derivatives or
cross derivatives of the dependent variable. Appropriate approximating functions need to be chosen
to avoid the problem of singularities.
Variable Coefficients The DR-BEM can be readily extended to equations with variable coeffi-
cients. Consider the variable coefficient Helmholtz equation
r2u + (x) u = 0 (7.58)
where is a function of position - = (x; y ) in two dimensions. If the DR-BEM is applied using
the known fundamental solution to the Laplace operator then the forcing function is b = ;u.
Applying the DR-BEM gives

Hu ; Gq = H U^ ; GQ^ F ;1b (7.59)
b
where is a vector of the nodal values of the forcing function b. The relationship b = ;u can be
b
written in matrix form as = ; Ku where K
is a diagonal matrix containing the nodal values of
(x; y) i.e.,
2 3
(x ; y ) 0 0
66 01 1 (x2; y2) 0 77
K = 64 ... ..
.
..
.
..
.
75 (7.60)
0 0 (xM ; yM )
where M is the number of collocation points used in applying the DR-BEM.
b Equation (7.59) can be rearranged to give

Using this matrix expression for

(H + SK ) u = Gq whereS = H U^ ; GQ^ F ;1 (7.61)
which is a boundary-only expression for the variable coefficient Helmholtz equation. This method
is general and can easily be extended to accommodate variable coefficient derivative terms and a
sum of variable coefficient terms.
Formulating the DR-BEM for a General Elliptic Problem In this section it has been shown
how the DR-BEM can be applied for elliptic problems with varying forms of b. The DR-BEM can
be applied in cases where b involves a sum of terms due to the basic property
Z Z Z
(b1 + b2 ) d
= b1 d
+ b2 d
(7.62)

Consider a two-dimensional equation of the form

@u + m (x; y) @u + n (x; y)
r2 u (x; y) = k (x; y) u + l (x; y) @x (7.63)
@y
Applying the DR-BEM to this equation gives a matrix system of the form
(H ; R) u = Gq + Sn (7.64)
where

S = H U^ ; GQ^ F ;1 (7.65)
@F @F
R = S K + L @X + M @Y F ;1 (7.66)
K , L and M are diagonal matrices where the diagonals contain the nodal values of k, l and m
respectively. n is a vector containing the nodal values of n.
The DR-BEM Using Other Operators

The DR-BEM has been presented in this chapter based on the Laplace operator. However the DR-
BEM can be be applied using essentially any operator of appropriate order with known fundamental
solution. If an appropriate operator can be found the complexity of the forcing function b can be
reduced. This should improve the accuracy of the method. The problem with applying the DR-
BEM based on another operator is in choosing the approximating function f . A choice of f which
produces accurate results is required but it is also necessary to choose an f for which a particular
solution u^ can be determined.
Zhu (1993) has determined the particular solutions necessary for applying the DR-BEM based
on the two-dimensional Helmholtz operator.
r2 u + 2 u = b (x; y; u; t) (7.67)
Radial functions have generally been used when applying the DR-BEM. Along the lines of Wrobel
et al. (1986), Zhu chose an approximating function of the form f = r m where m is a positive
integer. Determining the particular solution u
^ requires solving the ordinary differential equation
d2u^ + 1 du^ + 2u = rm (7.68)
dr2 r dr
which can be achieved using a variation of coefficients method.
Partridge et al. (1992) applied the DR-BEM to the transient convection diffusion equation
Dr2u ; vx @u
@x ; v y
@u ; ku = @u
@y @t (7.69)
where the material parameters D , vx , vy and k are all assumed to be homogeneous. They applied
the DR-BEM based on the steady-state convection-diffusion operator
Dr2u ; vx @u
@x ; v y
@u ; ku = 0
@y (7.70)
which has a known fundamental solution.

This analysis requires the determination of a particular solution u
^ which satisfies
Dr2u^ ; vx @@xu^ ; vy @@yu^ ; ku^ = f (7.71)
Instead of defining a form of the approximating function f and solving for u

^ Partridge et al. chose
to define u
^ and use Equation (7.71) to determine the corresponding approximating function. Al-
though somewhat ad-hoc this approach was found to produce accurate results.
Chapter 8
The BEM for Parabolic PDES
8.1 Time-Stepping Methods

Several approaches have been proposed for applying the BEM to parabolic problems. These meth-
ods can be broadly classified into two main approaches. Either some form of time-stepping proce-
dure is used to advance the solution in time, or a semi-analytic technique is used which can directly
calculate a solution at a specified time. In this section time-stepping procedures will be considered.
Time-stepping approaches discretise the time domain in some manner and use some form of
time marching scheme to advance the solution from one discrete time to the next. The two most
commonly used time-stepping methods are the coupled finite difference - BEM and the direct time
integration method. These two methods will be outlined in this section for the diffusion equation
r2 u (x; t) = 1 @u (@tx; t) (8.1)
where the diffusivity is a material parameter which can be a constant or a function of position.
8.1.1 Coupled Finite Difference - Boundary Element Method

This approach discretises the time-domain in a finite difference form. Consider the variation be-
tween a time tm and a time tm+1 = tm + t. The most common approach (Brebbia et al. 1984b) is
to assume that, for sufficiently small t, the time derivative can be approximated using a first-order
fully implicit finite difference scheme
@u (x; tm+1 ) = u (x; tm+1 ) ; u (x; tm ) (8.2)

@t t
which allows the diffusion equation in this time-range to be approximated as
; ;
r2u x; tm+1 ; 1 t u x; tm+1 + 1 t u (x; tm) = 0 (8.3)
Using this finite difference approximation the original parabolic equation has been reduced to an
elliptic equation. Using the weighted residuals method an integral equation can be generated from
136 T HE BEM FOR PARABOLIC PDES
Equation (8.3).
Z Z 1 Z um ! d
c ( ) um+1 + um+1 @! d; =
@n qm+1 !d; + t (8.4)
;

x
where um+1 = u ( ; tm+1 ) and um = u (x; tm). The fundamental solution ! is a solution of the
modified Helmholtz equation
r2 ! (x; ) ; 1 t ! (x; ) + () = 0 (8.5)
applied at some source point . The fundamental solution of the modified Helmholtz equation
is known in both two and three dimensions. If an internal solution is required at a specific time
this can be determined explicitly from Equation (8.4) where the fundamental solution is applied at

internal point and c ( ) = 1.
Unfortunately Equation (8.4) contains a domain integral. This integral is generally evaluated by
using a domain mesh (Brebbia et al. 1984b). The domain integral does not include any problem un-
knowns so a fairly coarse domain mesh will generally suffice. Applying the BEM to Equation (8.4)
gives
Hum+1 ; Gqm+1 = Bum (8.6)
B
where is a matrix containing the influence coefficients due to the domain integral. Using Equa-
tion (8.6) the solution can be advanced in time. U
0 is known from the initial conditions so a
solution can be calculated at t = t0 + t. A solution at internal nodes can then be calculated. The
time-stepping procedure can be repeated using the internal solution at t = t 0 +t as pseudo-initial
conditions for the next time-step.
If a constant time-step is used the matrices , HG and B
can be calculated once and stored.
The boundary conditions can be applied to form a solution system of the form Ax
m+1 = where
x
m+1 is the vector of unknown nodal values at time t m+1 and is a vector constructed from
known nodal values from the previous time-step. For a problem with time-independent boundary
x
conditions at each time-step it is only necessary to update and solve the system for m+1 . If a
problem has time-dependent boundary conditions the solution system needs to be reformed at each
time-step.
This coupled finite difference - boundary element method (FD-BEM) was first proposed by
Brebbia & Walker (1980) for the diffusion equation. It was implemented and investigated by
Curran, Cross & Lewis (1980). They found that this method will only produce accurate results
if Equation (8.2) accurately approximates the time derivative. This will generally require small
time-steps to be adopted. Curran et al. investigated the use of a higher-order approximation to the
time-derivative. They found that this improved the accuracy of the method. Unfortunately it lead
to a deterioration in convergence behaviour.
Tanaka, Matsimoto & Yang (1994) proposed a generalised version of this time-stepping scheme.
They approximated the time variation within an interval as
;
u (x; t) = u x; tm+1 + (1 ; ) u (x; tm ) (8.7)
8.1 T IME -S TEPPING M ETHODS 137
where , termed the time-scheme parameter, is a constant in the range 0 < 1. Substituting
this approximation and a first-order finite difference approximation of the time derivative into the
diffusion equation gives
r2 um+1 + (1 ; ) r2um = u ;t u

m+1 m
(8.8)
If = 1 this approximation of the diffusion equation is equivalent to the standard FD-BEM dis-
cussed earlier. An integral equation can be derived from Equation (8.8). Tanaka et al. implemented
this method and found it gave accurate results for a range of diffusion problems. They tested the
accuracy for a Crank-Nicolson scheme ( = 21 ), a Galerkin scheme ( = 23 ) and a fully implicit
scheme ( = 1). They found that the best results were achieved using a Crank-Nicolson scheme.
8.1.2 Direct Time-Integration Method

Instead of converting the original parabolic equation to an elliptic equation the problem can be
treated directly in the time domain by directly integrating over both time and space. The weighted
residual statement using this approach is
ZtF Z
(x; t) ! (; x; t ; t) d
dt = 0
r u (x; t) ; 1 @u @t
2
F (8.9)
t0
Integrating in time once and in space twice gives
ZtF Z @! ZtF Z Z
c ( ) u (; tF ) + u @n d; dt = q! d; dt + u (x; t0) ! d
(8.10)
t0 ; t0 ;
where the fundamental solution ! satisfies

x; tF ; t) + () (tF ) = 0
r2! (; x; tF ; t) + @! (;@t (8.11)
This time dependent fundamental solution is known in two and three dimensions. Physically this
x
fundamental solution represents the effect at a field point at time t of a unit point source applied
at a point at time tF . If an internal solution is required at a specific time this can be determined
from Equation (8.10) with c ( ) = 1.
The variation of u and q with time is unknown so it is still necessary to step in time. However,
as the time dependence is included in the fundamental solution, accurate results can be achieved
using larger time-steps than with the FD-BEM. Two different time-stepping schemes can be used.
Similarly to the FD-BEM, each time-step can be treated as a new problem so that an internal
solution is constructed at the end of each time-step to be used as pseudo-initial conditions for the
next time-step. Alternatively the time integration process can be restarted at t 0 with increasing
numbers of intermediate steps used. These two time-stepping approaches are discussed in detail in
Brebbia et al. (1984b).
The first method requires a new domain integral to be calculated after each time-step due to
the updated pseudo-initial conditions. The second time-stepping procedure involves only a domain
integral at t0 so, ideally, a domain integral only needs to be calculated once. This, however, will still
require the user to create a domain mesh. As mentioned by Brebbia et al. (1984b), in many practical
cases the domain integral can be avoided. If the initial conditions are u 0 = 0 throughout the body
the domain integral equals zero. If the initial conditions satisfy Laplace’s equation r 2 u0 = 0 then
a Galerkin vector can be found and the domain integral can be expressed as equivalent boundary
integrals. This includes many practical cases such as constant initial temperature or an initial linear
temperature profile.
Unfortunately, in practice it is not always feasible to restart the integration process at t 0 . At
each time-step new H G
and matrices are required so if many time-steps are required the storage
capacity of the computer is likely to be exceeded. This requires the procedure to be restarted
at some time where an internal solution is constructed and used as pseudo-initial conditions to
repeat the process. Therefore, in practice, both time-stepping methods are likely to require domain
integration.
8.2 Laplace Transform Method

An alternative approach which avoids time-stepping is to solve the problem in a transform domain
which removes the time dependence of the problem. The parabolic PDE is thus converted to an
elliptic problem for which the boundary element method has been shown to generally produce
accurate results. Once the solution to the elliptic problem is determined in the transform space
a solution in the original space can be attained using an inverse transform procedure. The most
appropriate transform approach for parabolic problems is the Laplace transform.
Consider the diffusion equation
(x; t)
r2u (x; t) = 1 @u @t (8.12)
with appropriate boundary and initial conditions. The Laplace transform of u ( x; t) will be sym-
x
bolised as U ( ; ) and is defined as
Z1
U (x; ) = e;t u (x; t) dt (8.13)
0
Applying Laplace transforms to Equation (8.12) gives
r2 U (x; ) = 1 (U (x; ) ; u0 (x)) (8.14)
x
with transformed boundary conditions. u 0 ( ) is the initial conditions of u. Equation (8.14) is an
elliptic PDE which can be readily solved using the boundary element method. Once the solution
is determined in Laplace transform space this solution can be inverted to give a solution in the
time-domain. This inversion procedure requires solutions to be generated for several values of the
transform parameter .
This method was first proposed by Rizzo & Shippy (1970) and has since been successfully
8.3 T HE DR-BEM F OR T RANSIENT P ROBLEMS 139
used by other practitioners (Moridis & Reddell 1991) (Zhu, Satravaha & Lu 1994). Liggett & Liu
(1979) compared the Laplace transform method with the time-dependent Green’s function method.
They noted that the direct method is simpler to apply. However, due to its greater efficiency, they
recommended the Laplace transform method for solving diffusion problems.
One limitation of the Laplace transform method is that Equation (8.14) is inhomogeneous so
that applying the standard BEM will generate a domain integral involving the initial conditions.
Traditionally this domain integral has been calculated by using a domain discretisation (Brebbia
et al. 1984b). However, recently Zhu et al. (1994) proposed using the DR-BEM to convert this
domain integral term to equivalent boundary integrals. They chose to apply the DR-BEM based
on the known fundamental solution to the Laplace operator. Considering Equation (8.14) this
means that the DR-BEM will be used to convert the right-hand-side to equivalent domain integrals.
Therefore the required DR-BEM approximation is
1 (U (x; ) ; u (x)) = NX +L

fj j
0 (8.15)
j =1
The DR-BEM can now be applied to Equation (8.14), giving a matrix system of the form

H ; S u ; Gq = ; 1 Su0 (8.16)
which can be reduced to a square system by applying boundary conditions. Once the solution
is determined for this elliptic equation in the transform space a solution at a given time can be
constructed using an inversion process.
This Laplace transform dual reciprocity method (LT-DRM) can easily be extended to equations
of the form
(x; t) + b (x; u)
r2 u (x; t) = 1 @u @t (8.17)
in which case a matrix expression of the form

H ; R ; S u ; Gq = ; 1 Su0 (8.18)
is generated. Zhu and his colleagues have successfully extended the LT-DRM to solve diffusion
problems with nonlinear source terms.
8.3 The DR-BEM For Transient Problems

The DR-BEM can also be applied to parabolic problems. Consider, for example, the diffusion
equation
r2 u = 1 @u
@t (8.19)
where the thermal diffusivity, , is a constant. In this case the global approximation of b implies a
separation of variables such that
@u = X
M
@t j=1 fj (x) j (t) (8.20)
Using Equation (8.20), Equation (7.39) becomes

Hu ; Gq = 1 H U^ ; GQ^ F ;1 @u
@t (8.21)
or

C @u
@t + Hu = Gq where C = ; 1 H u^ ; GQ^ F ;1 (8.22)
Equation (8.22) can be solved using a standard direct time-integration method.

Partridge & Brebbia (1990) recommended using a first-order finite difference approximation to
the time derivative
@u = um+1 ; um (8.23)
@t t
and linear approximations to u and q within a time-step.
u = (1 ; u) um + uum+1 (8.24)

q = (1 ; q ) qm + q qm+1 (8.25)
where u and q are weighting parameters with values in the range (0; 1] and the time-step is
between times tm and tm+1 = tm + t. Substituting these approximations into Equation (8.22) an
expression at tm+1 can be derived in terms of values at tm .
1 hc i
t C + uH um+1 ; q Gqm+1 = u) H u + (1 ; q ) Gq
; (1 ; m m (8.26)
t
The values of u0 and q 0 are known from the initial conditions so a time-stepping procedure can be
CH
used. If a constant time-step is used the matrices , G
and only need to be constructed once.
Using this two-level time-integration scheme the most common choice of time-scheme parameters
is u = 0:5; q = 1:0.
8.4 The MRM for Transient Problems

The MRM can be applied to the diffusion equation r 2 u = 1 @u
@t using the fundamental solution
of Laplace’s equation. In this case the forcing function becomes b 0 =
1 @u and the recurrence
@t
8.4 T HE MRM FOR T RANSIENT P ROBLEMS 141
relationship defined by Equation (7.24) becomes
uj+1 = r2 uj = @@tuj
j
(8.27)
The higher-order fundamental solutions are known for Laplace’s equation. In this case the MRM
formulation becomes
1 Z
X ju 1 Z
X
c ( ) u ( ) + @! @ d ; = ! @ j q d; (8.28)
j =0 ; @n @tj @tj
j =0 ;
The standard BEM numerical procedure can be applied to this boundary integral equation. This
gives the matrix system
H0u + H1u_ + H2u + : : : = G0q + G1q_ + G2q + : : : (8.29)
where the matrices H1 ; G1 etc are the influence coefficient matrices relating to the higher-order
fundamental solutions. This equation can be solved using a time-integration procedure.
The most common approach is to solve this system numerically by discretising the time domain
and using a time-stepping procedure. This requires some interpolation between the two time-levels
marked by m and m + 1. This most common approach is to use a linear approximation to u and q
in this time-range
u = (1 ; ) um + um+1 (8.30)
q = (1 ; ) qm + qm+1 (8.31)
where has a value in the range 0 to 1. Differentiating these linear approximations gives
u_ = u t;m u
m+1 m
(8.32)
q m+1 ; q m
q_ = tm (8.33)
and all the other derivatives vanish.

This allows Equation (8.29) to be simplified to
HLum+1 ; GLqm+1 = ;HRum + GRqm (8.34)
where HL = 1tm H1 + H0, HR = ; 1tm H1 + (1 ; ) H0, GL = 1tm G1 + G0, GR =

; 1tm G1 + (1 ; ) G0 . This approach is termed a first-order approach as it removes all but the
first derivatives. A second order approach can be formulated by using quadratic interpolation of u
and q within the time-range.
Using Equation (8.34) the solution can be advanced in time. If a constant time-step is used the
H H G G
matrices L ; R ; L and R only need to be constructed once outside the time-stepping loop.
If the boundary conditions are not time-dependent the boundary conditions only need to be applied
once.
Bibliography
Brebbia, C. A. & Nowak, A. (1989), A new approach for transforming domain integrals to the
boundary, in R. Gruber, J. Periaux & R. P. Shaw, eds, ’Proceedings of the Fifth International
Symposium on Numerical Methods in Engineering’, Computational Mechanics Publications,
Springer-Verlag, pp. 73–84.
Brebbia, C. A., Telles, J. C. F. & Wrobel, L. C. (1984a), Boundary Element Techniques, Springer-
Verlag.
Brebbia, C. A., Telles, J. C. F. & Wrobel, L. C. (1984b), Boundary element techniques: Theory
and applications in engineering, Springer-Verlag, New York.
Brebbia, C. A. & Walker, S. (1980), Boundary Element Techniques In Engineering, Newnes–

Butterworths.
Clements, D. L. & Rizzo, F. (1978), ’A method for the numerical solution of boundary value
problems governed by second-order elliptic systems’, Journal of the Institute of Mathematics
Applications 22, 197–202.
Coleman, C. J., Tullock, D. L. & Phan-Thien, N. (1991), ’An effective boundary element method
for inhomogeneous partial differential equations’, Journal of Applied Mathematics and
Physics (ZAMP) 42, 730–745.
Curran, D. A. S., Cross, M. & Lewis, B. A. (1980), ’Solution of parabolic differential equations by
the boundary element method using discretisation in time’, Applied Mathematical Modelling
4, 398–400.
Danson, D. J. (1981), A boundary element formulation of problems in linear isotropic elasticity

with body forces, in C. A. Brebbia, ed., ’Boundary Element Methods’, Computational Me-
chanics Publications and Springer-Verlag, pp. 105–122.
Fairweather, G., Rizzo, F. J., Shippy, D. J. & Wu, Y. S. (1979), ’On the numerical solution of
two-dimensional problems by an improved boundary integral equation method’, J. Comput.
Phys. 31, 96–112.
Gipson, G. S. (1987), Boundary Element Fundamentals - Basic Concepts And Recent Develop-
ments In The Poisson Equation, Computational Mechanics Publications.
Gipson, G. S., Reible, D. D. & Savant, S. A. (1987), Boundary elements and perturbation theory
for certain classes of hyperbolic and parabolic problems, in C. A. Brebbia, W. L. Wendland &
144 BIBLIOGRAPHY
G. Kuhn, eds, ’Boundary Elements IX’, Computational Mechanics Publications and Springer-
Verlag, pp. 115–127.
Golberg, M. A. (1995), ’The numerical evaluation of particular solutions in the BEM - a review’,
Boundary Element Communications .
Golberg, M. A. & Chen, C. S. (1994), ’The theory of radial basis functions applied to the BEM for
inhomogeneous partial differential equations’, BE Communications 5(2), 57–61.
Huang, Q. & Cruse, T. A. (1993), ’Some remarks on particular solution used in BEM formulation’,
Computational Mechanics 13, 68–73.
Itagaki, M. & Brebbia, C. A. (1993), ’Generation of higher order fundamental solutions to the two-
dimensional modified Helmholtz equation’, Engineering Analysis With Boundary Elements
11(1), 87–90.
Lafe, O. E. & Cheng, A. H.-D. (1987), ’A perturbation boundary element code for steady state
groundwater flow in heterogeneous aquifers’, Water Resources Research 23(6), 1079–1084.
Liggett, J. A. & Liu, P. L.-F. (1979), ’Unsteady flow in confined aquifers: A comparison of two
boundary integral methods’, Water Resources Research 15(4), 861–866.
Monaco, A. D. & Rangogni, R. (1982), Boundary element method: Processing of the source term
of the Poisson equation by means of boundary integrals only, in K. P. Holz, U. Meissner,
W. Zielke, C. A. Brebbia, G. Pinder & W. Gray, eds, ’Finite Elements In Water Resources
IV’, Springer-Verlag, pp. 19.29–19.36.
Moridis, G. L. & Reddell, D. L. (1991), The Laplace transform boundary element (LTBE) method
for the solution of diffusion-type problems, in ’Boundary Elements XIII’, Computational
Mechanics Publications and Springer-Verlag, pp. 83–97.
Nowak, A. J. (1987), Solution of transient heat conduction problems using boundary-only formu-
lation, in C. A. Brebbia, W. L. Wendland & G. Kuhn, eds, ’Boundary Elements IX Vol 3’,
Computational Mechanics Publications and Springer-Verlag, pp. 265–276.
Nowak, A. J. & Partridge, P. W. (1992), ’Comparison of the dual reciprocity and the multiple
reciprocity methods’, Engineering Analysis With Boundary Elements 10, 155–160.
Ortner, N. (1987), Construction of fundamental solutions, in C. A. Brebbia, ed., ’Topics In Bound-
ary Elements Research’, Springer-Verlag, Berlin and New York.
Partridge, P. W. & Brebbia, C. A. (1990), The BEM dual reciprocity method for diffusion prob-
lems, in ’Computational Methods In Water Resources VIII’, Computational Mechanics Pub-
lications and Springer-Verlag, pp. 397–403.
Partridge, P. W., Brebbia, C. A. & Wrobel, L. C. (1992), The Dual Reciprocity Boundary Element
Method, Computational Mechanics Publications and Elsevier Applied Science.
Polyzos, D., Dassios, G. & Beskos, D. E. (1994), ’On the equivalence of dual reciprocity and
particular integral approaches in the BEM’, BE Communications 5(6), 285–288.
BIBLIOGRAPHY 145
Rangogni, R. (1986), ’Numerical solution of the generalized Laplace equation by coupling the
boundary element method and the perturbation method’, Applied Mathematical Modelling
10, 266–270.
Rangogni, R. (1991), Solution of variable coefficients PDEs by means of BEM and perturba-
tion technique, in ’Boundary Elements XIII’, Computational Mechanics Publications and
Springer-Verlag.
Rizzo, F. J. & Shippy, D. J. (1970), ’A method of solution for certain problems of heat conduction’,
AIAA Journal 8(11), 2004–2009.
Tanaka, M., Matsimoto, T. & Yang, Q. F. (1994), ’Time-stepping boundary element method applied
to 2-D transient heat conduction problems’, Applied Mathematical Modelling 18, 569–576.
Wrobel, L. C., Brebbia, C. A. & Nardini, D. (1986), The dual reciprocity boundary element for-
mulation for transient heat conduction, in ’Finite Elements In Water Resources VI’, Compu-
tational Mechanics Publications and Springer-Verlag.
Wu, J. C. (1985), Boundary element methods and inhomogeneous parabolic equations, in C. A.

Brebbia & B. J. Noye, eds, ’BETECH 85’, Springer-Verlag, pp. 19–30.
Zhang, Y. (1993), On the dual reciprocity boundary element method, Master’s thesis, The Univer-
sity of Wollongong.
Zheng, R., Coleman, C. J. & Phan-Thien, N. (1991), ’A boundary element approach for non-
homogeneous potential problems’, Computational Mechanics 7, 279–288.
Zhu, S. (1993), ’Particular solutions associated with the Helmholtz operator used in DRBEM’, BE
Abstracts 4(6), 231–233.
Zhu, S., Satravaha, P. & Lu, X. (1994), ’Solving linear diffusion equations with the dual reciprocity
method in Laplace space’, Engineering Analysis With Boundary Elements 13, 1–10.
Index
Advection-diffusion equation, 103 Helmholtz, 72

Area Coordinates, 14 Kelvin, 91
Laplace, 46
Basis functions Navier, 72
Hermite wave equation, 72
cubic, 12
Basis functions, 2–16, 53 Galerkin formulation, 25, 31
Hermite, 10–14 Galerkin vector, 119
bicubic, 12, 14 Gaussian quadrature, 39–42, 56
cubic, 10 Global stiffness matrix, 27
Lagrange, 10
bilinear, 7–9 integration by parts, 24
linear, 2–4 Isoparametric formulation, 53
quadratic, 7 Laplace transform method, 138
Beam elements, 79
Boundary conditions Mass lumping, 106
application of, 29, 54 Multiple reciprocity method, 122–124
diffusion equation, 140
Coupled finite difference - boundary element Helmholtz equation, 124
method, 135 Poisson equation, 122
Curvilinear coordinate systems, 17–18
Cylindrical polar, 17 Perturbation BEM, 121
Prolate spheroidal, 18 Plane stress elements, 81–83
Spherical polar, 18 Potential energy, 82
Dirac-Delta function, 43–45 Rayleigh-Ritz method, 75, 77
Direct time-integration method, 137 Regular BEM, 64
Dual reciprocity BEM
approximating function, 127 Strain energy, 82
derivative terms, 129
Trefftz method, 64
elliptic problems, 129
Triangular elements, 14–16
transient problems, 139
Truss elements, 76
variable coefficients, 131
Dual reciprocity BEM, 124, 133 Weighted residual, 24
Weighting function, 24
element stiffness matrix, 26
Fundamental solution, 43, 45–48, 72, 117

diffusion equation, 72

FEM Lecture Notes by Peter Hunter, Andrew Pullian

Uploaded by

Copyright:

Available Formats

FEM Lecture Notes by Peter Hunter, Andrew Pullian

Uploaded by

Document Information

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

FEM Lecture Notes by Peter Hunter, Andrew Pullian

Uploaded by

Copyright:

Available Formats

FEM/BEM NOTES

Professor Peter Hunter

Associate Professor Andrew Pullan

Department of Engineering Science

1 Finite Element Basis Functions 1

2 Steady-State Heat Conduction 23

3 The Boundary Element Method 43

3.3 The Two-Dimensional Boundary Element Method . . . . . . . . . . . . . . . . . . 48

5 Transient Heat Conduction 99

5.3 The Transient Advection-Diffusion Equation . . . . . . . . . . . . . . . . . . . . 103

6 Modal Analysis 111

7 Domain Integrals in the BEM 117

8 The BEM for Parabolic PDES 135

Finite Element Basis Functions

1.1 Representing a One-Dimensional Field

One approach would be to use a polynomial expression u (x) = a + bx + cx 2 + dx3 + : : :

1.2 Linear Basis Functions

node 1 node 2 node 3 node 4

u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.1)

u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.2)

u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.3)

1.3 Basis Functions as Weighting Functions

which depends on u1 and u2 but is weighted more towards u2 than u1

x ( ) = '1 ( ) x1 + '2 ( ) x2 (1.4)

defined by the parametric expressions

1.4 Quadratic Basis Functions

u ( ) = '1 ( ) u1 + '2 ( ) u2 + '3 ( ) u3 (1.5)

F IGURE 1.8: One-dimensional quadratic basis functions.

1.5 Two- and Three-Dimensional Elements

'4 (1; 2) = 12

1.6 Higher Order Continuity

and impose the constraints

u ( ) = u1 + u01 + (3u2 ; 3u1 ; 2u01 ; u02)  2 + (u01 + u02 + 2u1 ; 2u2)  3

u ( ) = 01 ( ) u1 + 11 ( ) u01 + 02 ( ) u2 + 12 ( ) u02 (1.14)

The bicubic interpolation of these nodal parameters is given by

u (1; 2) = 01 (1) 01 (2) u1 + 02 (1) 01 (2) u2

1.7 Triangular Elements

u (x; y) = '1 (x; y) u1 + '2 (x; y) u2 + '3 (x; y) u3

1.8 Curvilinear Coordinate Systems

Prolate spheroidal (; ;  ) :

1.9 CMISS Examples

F IGURE 1.20: Node positions for example 111.

F IGURE 1.21: 2D bilinear finite element mesh for example 111.

3. To define a quadratic-linear element run the cmiss example 115.

4. To define a 3D trilinear element run CMISS example 121.

5. To define a 2D cubic Hermite-linear finite element mesh run example 114.

F IGURE 1.22: First refined mesh for example 113

F IGURE 1.23: Second refined mesh for example 113

F IGURE 1.24: Defining a triangular mesh for example 116

Steady-State Heat Conduction

2.1 One-Dimensional Steady-State Heat Conduction

subject to boundary conditions: u (0) = 0 and u (1) = 1.

with which we can compare the approximate finite element solutions.

5. Assemble the global equations.

6. Apply the boundary conditions.

7. Solve the global equations.

8. Evaluate the fluxes.

2.1.1 Integral equation

where R is the residual

2.1.2 Integration by parts

u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.1)

u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.2)

u ( ) = '1 ( ) u1 + '2 ( ) u2 (1.3)

x ( ) = '1 ( ) x1 + '2 ( ) x2 (1.4)

u ( ) = '1 ( ) u1 + '2 ( ) u2 + '3 ( ) u3 (1.5)

'4 (1; 2) = 12

u ( ) = u1 + u01 + (3u2 ; 3u1 ; 2u01 ; u02) 2 + (u01 + u02 + 2u1 ; 2u2) 3

u ( ) = 01 ( ) u1 + 11 ( ) u01 + 02 ( ) u2 + 12 ( ) u02 (1.14)

u (1; 2) = 01 (1) 01 (2) u1 + 02 (1) 01 (2) u2

Prolate spheroidal (; ; ) :

u ( ) = '1 ( ) u1 + '2 ( ) u2 = 'n ( ) un

and each element integral is then taken over -space

and putting x in terms of gives (with

'1 = (1 ; 1)(1 ; 2 ) @'1 = ;(1 ; )

'2 = 1(1 ; 2) @'2 = 1 ;

'3 = (1 ; 1)2 @'3 = ;

'4 = 12 @'4 =

X Z @'n @'m @'n @'m Z @u

and choose f ( ) = a + b + c 2 + d 3 . Then

= ; 2 " u d; since r = " on ;"

The Jacobian J ( ) can be found by

where s represents the arclength and dx d and

!i = ; 21 log (ri ( ))

!i (1; 2) = 4r (1 ; )